ecearth3/sources/tm5mp/proj/cbm4/emission_data.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_DATA
!
! !DESCRIPTION: Provides general methods and fields needed for the emission
!               routines. 
!
!    emission_vdist_by_sector : distribute vertically emissions according to source sector
!    msg_emis                 : print formated emission total
!    do_add_2d                : add 2D emission to tracer mass array
!    do_add_3d                : add 3D emission to tracer mass array
!    do_add_3d_cycle          : add 3D emission to tracer mass array with a daily cycle
!    scale_cycle              : evaluate a daily cycle
!    distribute_l44           : scale input field with NH (SH) fires distributed over JAS (JFM) 
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_DATA
  !
  ! !USES:
  !
  use GO,           only : gol, goPr, goErr
  use Dims,         only : nregions, lm
  use global_types, only : emis_data, d2_data, d23_data, d3_data

#ifdef with_m7 
  use mo_aero_m7,   only : nmod
  !use mo_aero,      only : zom2oc
  use global_types, only : modal_emissions
#endif

#ifdef with_budgets
  use budget_global,only : budg_dat, nzon_vg
#endif

  IMPLICIT NONE
  !
  ! !PUBLIC DATA MEMBERS & TYPES:
  !
  public :: emis_data, d2_data, d23_data
  !
  type(emis_data), target    :: plandr(nregions)    ! Land fraction
  type(emis_data), target    :: emis2D(nregions)    ! 2D emiss prior vertical redistribution
  type(emis_data), target    :: emis_temp(nregions) ! temp array to store emission

#ifdef with_m7
  type(modal_emissions), dimension(nregions,nmod), target :: emis_mass
  type(modal_emissions), dimension(nregions,nmod), target :: emis_number

  ! The value of 1.4 for the POM to OC mass ratio is an outdated estimate,
  ! see e.g. Turpin and Lim (Aerosol Sci. Technol., 2001) and
  ! Aiken et al. (Environ. Sci. Technol., 2008).
  ! Turpin and Lim estimate a ratio of 1.6 +- 0.2 for urban aerosol,
  ! and 2.1 +- 0.2 for aged (nonurban) aerosol;
  ! They also note that aerosols heavily impacted by woodsmoke can
  ! have an even higher ratio (2.2 to 2.6).
  ! According to Reid et al. (ACP, 2005),
  ! the POM to OC ratio in fresh biomass burning smoke is very uncertain,
  ! somewhere in between 1.4 and ~2.
  ! Aiken et al. measure ambient aerosol values 
  ! between 1.25 for O/C = 0 to 2.44 for O/C = 1.0.
  ! For OA from biomass burning, they measure 1.56-1.70,
  ! lower than the estimates from Turpin and Lim (~2.0).
  ! They find the highest ratios for 
  ! aged and freshly formed SOA (~2.4 and ~1.9, respectively)
  ! and lowest values for primary OA from urban combustion.
  ! Based on these studies, we apply different values
  ! for emissions from biomass burning versus other emissions,
  ! as is also done in some other models
  ! (see Table 1 in Tsigaridis et al., ACP, 2014),
  ! For SOA (seem emission_pom.F90)
  ! we apply a relatively high value valid for aged SOA.
  ! This compensates for the lack of SOA formation from isoprene,
  ! and improves the agreement with aerosol optical depth (AOD)
  ! derived from satellite observations (MODIS).
  !
  ! We could go even further and apply different values for the 
  ! water soluble and insoluble fractions (Turpin and Lim).
  ! 
  ! It should be acknowledged that the representation of OA 
  ! with a single tracer is very simplistic.
  ! In particular, increase in OA mass due to ageing
  ! is not properly accounted for.
  ! 
  !real, parameter  ::  oc2pom = zom2oc    !factor for conversion of OC mass to POM
  real, parameter  ::  oc2pom_ff = 1.6
  real, parameter  ::  oc2pom_vg = 1.6
  real, parameter  ::  oc2pom_soa = 2.4
#endif

  ! AR5 - EDGAR4 - RETRO - GFED3 - MACC - LPJ - HYMN - MEGAN
  logical            :: LAR5, LAR5BMB, LEDGAR4, LRETROF, LGFED3, LMACC, LLPJ, LHYMN, LMACCITY, LMEGAN
  character(len=256) :: emis_input_dir
  character(len=256) :: emis_input_dir_ar5
  character(len=256) :: emis_input_dir_megan
  character(len=256) :: emis_input_dir_ed4
  character(len=256) :: emis_input_dir_mac
  integer            :: emis_input_year
  character(len=256) :: emis_input_dir_gfed
  character(len=256) :: emis_input_dir_retro
  character(len=256) :: emis_input_dir_aerocom
  character(len=256) :: emis_input_dir_dms
  character(len=256) :: emis_input_dir_rn222
  character(len=256) :: emis_input_dir_dust, emis_input_dust
#ifdef with_ch4_emis
  character(len=256) :: emis_input_dir_natch4
#endif
  logical            :: emis_ch4_single
  logical            :: emis_ch4_fix3d
  real               :: emis_ch4_fixed_ppb
  character(len=256) :: emis_zch4_fname
  
  ! biomass burning levels:
  integer, parameter ::  bb_nlev = 5
  integer, parameter ::  bb_lm   = lm(1)    ! to be reduced to strip stratosphere!
  integer, parameter ::  bl_lm   = lm(1)    ! to be reduced to the BL (around 8)
  real,    parameter ::  bb_hlow (0:bb_nlev) = (/ 0.,   100.,  500., 1000., 2000., 3000./)
  real,    parameter ::  bb_hhigh(0:bb_nlev) = (/ 100., 500., 1000., 2000., 3000., 6000./)

  ! biomass burning levels in tropics (-20 - 20)::
  integer, parameter ::  bb_nlev_trop = 3
  real,    parameter ::  bb_hlow_trop (0:bb_nlev_trop) = (/ 0.,   100.,  500., 1000./)
  real,    parameter ::  bb_hhigh_trop(0:bb_nlev_trop) = (/ 100., 500., 1000., 2000./)

  ! -----------------------------------------------------------------------------------
  ! Biomassburning time splitting factors (now globally, instead of by constituent)
  logical :: emis_bb_trop_cycle

  type bb_cycle_data
     real, pointer ::  scalef(:)
  end type bb_cycle_data
  type(bb_cycle_data), dimension(nregions), target :: bb_cycle

  ! -----------------------------------------------------------------------------------
  ! Budgets
#ifdef with_budgets
  real, dimension(nregions)                    :: sum_emission

  type budemi_data
     real,dimension(:,:,:,:),pointer           :: emi    ! [i, j, nbud_vg, ntracet]
  end type budemi_data
  type(budemi_data),dimension(nregions),target :: budemi_dat
#endif

  logical :: use_tiedkte ! read from rc file, used to scale LiNOx
  
  ! ----------------------------------------------------
  ! Vertical splitting : used now for all sectors, for biomassburning and other categories
  ! (tables in Dentener, 2006, ACP)
  ! 
  ! ATTENTION: changes here have impact on the routine emission_vdist_by_sector (contained)
  ! Note: vdist_* variable could be hold private
  !
  integer, parameter                      ::  vd_class_name_len = 12
  ! biomassburning
  integer, parameter                      ::  vdist_bb_nlev  = 6
  real,dimension(vdist_bb_nlev),parameter ::  vdist_bb_hlow  = (/ 0.,   100.,  500., 1000., 2000., 3000./)
  real,dimension(vdist_bb_nlev),parameter ::  vdist_bb_hhigh = (/ 100., 500., 1000., 2000., 3000., 6000./)
  ! other sources (related to surface)
  integer, parameter                      ::  vdist_nn_nlev  = 6
  real,dimension(vdist_nn_nlev),parameter ::  vdist_nn_hlow  = (/ 0.,  30. , 100., 300.,  600., 1000. /)
  real,dimension(vdist_nn_nlev),parameter ::  vdist_nn_hhigh = (/ 30., 100., 300., 600., 1000., 2000. /)
  !
  ! !PRIVATE TYPES:
  !
  character(len=*), parameter, private  ::  mname = 'emission_data'
  !
  ! !REVISION HISTORY: 
  !    1 Oct 2010 - Achim Strunk - adapted for AR5
  !   28 Jun 2012 - Ph. Le Sager - made everything public by default
  !                              - adapted for lon-lat MPI domain decomposition
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_VDIST_BY_SECTOR
  !
  ! !DESCRIPTION: Return vertically distributed emissions depending on given
  !               sector and constituent.
  !               Distribution is done as 'compromise' between 
  !               1)  Dentener et al., 2006, ACP, slightly modified for
  !                   emissions supposed to enter the surface layer. 
  !               2)  De Meij et al., ACP, 2006
  !               3)  EMEP model (docu, model code, publications)
  !               4)  Bieser et al., GMDD, 2010 
  !               5)  Small updates to biomassburning following 
  !                   Huijnen et al., GMDD, 2010
  !\\
  !\\
  ! !INTERFACE:
  !  
  SUBROUTINE EMISSION_VDIST_BY_SECTOR( splitname, constituent, region, emis2D, emis3D, status )
    !
    ! !USES:
    !
    use toolbox,         only : distribute_emis2D
    !
    ! !INPUT PARAMETERS:
    !
    character(len=vd_class_name_len), intent(in)          :: splitname
    character(len=*),                 intent(in)          :: constituent
    integer,                          intent(in)          :: region
    type(emis_data),                  intent(in)          :: emis2D
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(d3_data),                    intent(out)         :: emis3D
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,                          intent(out)         :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - v0
    !
    ! !REMARKS: Splitting depending on constituent is still to code. 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer                         :: ilev
    character(len=2)                :: splittype
    integer, parameter              :: maxvdist_nlev = max(vdist_bb_nlev,vdist_nn_nlev)
    real, dimension(maxvdist_nlev)  :: sfglobal, sfintrop, sfintemp

    ! --- begin --------------------------------------
    ! control implicit settings here in this routine
    if( vdist_bb_nlev /= 6  .or. vdist_nn_nlev /= 6 ) then 
       write(gol, '("ERROR: wrong number of layers: bb_nlev /= 6  or  nn_lev /= 6")'); call goErr
       status=1; return
    end if

    ! initialise
    sfglobal = 0.0 
    sfintrop = 0.0 
    sfintemp = 0.0 
    emis3d%d3= 0.0
    
    ! -------------------------------------------------
    ! This here is selection by emission source sectors 
    ! -------------------------------------------------
    select case( trim(splitname) )

    case( 'combenergy' ) 
       ! Combustion from energy production and transformation, power-plants
       ! no zonal differences, no dependency on species
       ! --
       ! assumed to be related to EMEP SNAP category [1]
       ! AEROCOM equivalent: power-plants
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 0.0, 0.1, 0.7, 0.2, 0.0, 0.0 /) 

    case( 'combrescom' )
       ! Residential and commercial combustion
       ! no zonal differences, no dependency on species
       ! --
       ! assumed to be related to EMEP SNAP category [2]
       ! AEROCOM equivalent: domestic/industry
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 0.4, 0.4, 0.2, 0.0, 0.0, 0.0 /) 

    case( 'industry' )
       ! Industrial processes and combustion
       ! no zonal differences, no dependency on species
       ! --
       ! assumed to be an aggregate of EMEP SNAP categories [3,4,5]
       ! AEROCOM equivalent: industry
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 0.1, 0.2, 0.6, 0.1, 0.0, 0.0 /)

    case( 'waste' )
       ! Waste treatment and disposal
       ! no zonal differences, no dependency on species
       ! --
       ! assumed to be related to EMEP SNAP category [9]
       ! AEROCOM equivalent: - 
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 0.1, 0.2, 0.4, 0.3, 0.0, 0.0 /)

    case( 'nearsurface' )
       ! Near surface emissions 
       ! no zonal differences, no dependency on species
       ! --
       ! AR5: solvent production and use, agricultural waste burning, &
       !      ships, grassland fire, mineral dust (AEROCOM or Tegen-Vignati)
       ! EMEP equivalents: SNAP categories [6,8]
       ! AEROCOM equivalents: ships, agricultural waste
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 0.8, 0.2, 0.0, 0.0, 0.0, 0.0 /)

    case( 'surface' )
       ! Surface emissions
       ! no zonal differences, no dependency on species
       ! --
       ! AR5: agriculture, transport, soil, oceanic, biogenic
       ! EMEP equivalents: SNAP categories [7,8,10]
       ! AEROCOM equivalents: surface emissions from road/off-road, ships, agricultural waste
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)

!!$    case ( 'mindust' )
!!$      ! Mineral dust emissions
!!$      ! no zonal differences
!!$      ! -- 
!!$      ! AR5: surface dust emissions from either from AEROCOM or Tegen-Vignati
!!$      !      assumed to be slightly higher than near-surface
!!$      splittype = 'nn' 
!!$      sfglobal  = (/ 0.8, 0.2, 0.0, 0.0, 0.0, 0.0 /)

    case( 'volcanic' )
       ! Volcanic emissions
       ! --
       ! AR5: natural SOx emissions (from MACC repository)
       ! EMEP equivalents: SNAP categories [11]
       ! AEROCOM equivalents: continuous: 2/3 to 1 of height of volcano top
       !                      explosive: 0.5 to 1.5 km above volcano top
       ! Since no data is available about volcano top levels, we assume here, that
       ! a volcano top is usually 200-1000 m higher than the grid's surrounding terrain
       ! height, thus emissions are distributed from 100 - 2000 m almost homogeneously
       ! --
       splittype = 'nn'
       sfglobal  = (/ 0.0, 0.0, 0.1, 0.3, 0.4, 0.2 /)

    case( 'forestfire' )
       ! forest fires: 6 layer model, different splitting for different regions
       !               dont distinguish between boreal Europe and boreal Canada, use Europe globally
       ! Dentener et al., modified in tropics (up to 2 km) with respect to Huijnen et al., GMD 2010, \
       !  who are citing Labonne et al., 2007 ofr new information based on satellite data.
       splittype = 'bb' 
       sfglobal  = (/ 0.1, 0.1, 0.2, 0.2, 0.4, 0.0 /)
       sfintemp  = (/ 0.2, 0.2, 0.2, 0.4, 0.0, 0.0 /)
       sfintrop  = (/ 0.2, 0.2, 0.2, 0.4, 0.0, 0.0 /)

    case default
       write(gol, '("ERROR: wrong specification of emission sector ",a,"in vdist_by_sector called for constituent ",a)') &
            trim(splitname), trim(constituent); call goErr
       status=1; return

    end select

    ! --------------------
    ! Do the splitting -
    !    no need to skip cases of fraction (sfglobal, sfintemp, ..) being zero:
    !    distribute_emis2D is optimized for those and will return right away
    ! --------------------
    select case( splittype )
    case( 'nn' )
       do ilev = 1, vdist_nn_nlev
          call distribute_emis2D( region, emis2D, emis3D, vdist_nn_hlow(ilev), vdist_nn_hhigh(ilev), status, sfglobal(ilev) )
          IF_NOTOK_RETURN(status=1)
       end do
    case( 'bb' )
       do ilev = 1, vdist_bb_nlev
          ! boreal 
          call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, &
               sfglobal(ilev),  -90.,-61.)
          IF_NOTOK_RETURN(status=1)
          call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, &
               sfglobal(ilev),   60., 90.)
          IF_NOTOK_RETURN(status=1)
          ! temperate 
          call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, &
               sfintemp(ilev),  -60.,-31.)
          IF_NOTOK_RETURN(status=1)
          call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, &
               sfintemp(ilev),   30., 59.)
          IF_NOTOK_RETURN(status=1)
          ! tropics
          call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, &
               sfintrop(ilev),  -30., 29.)
          IF_NOTOK_RETURN(status=1)
       end do
    case default
       write(gol, '("ERROR: wrong specification of splitting type in vdist_by_sector.")'); call goErr
       status=1; return
    end select

    ! OK 
    status=0

  END SUBROUTINE EMISSION_VDIST_BY_SECTOR
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    MSG_EMIS
  !
  ! !DESCRIPTION: Provide output to verify the amount of emissions 
  !               entering the calculations
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE MSG_EMIS(year_month, provider, sector, msg_comp, msg_mass, summed_emissions)
    !
    ! !USES:
    !
    use dims, only: idate
    !
    ! !INPUT PARAMETERS:
    !
    integer,          intent(in) :: year_month      ! 1: year, 2: monthly amount
    character(len=*), intent(in) :: provider, sector, msg_comp
    real,             intent(in) :: msg_mass, summed_emissions
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - protex; increase space for output
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=7),dimension(2), parameter:: ym=(/'Yearly ','Monthly'/)

1111 format(' EMISS-INFO -   ',a7,' Emission from ',a8,' -> ',a18,1x,a6,' : mw=',f5.1,', month=',i2.2,', Sum [kg]=',1x,1pe11.4)
1112 format(' EMISS-INFO -   ',a7,' Emission from ',a8,' -> ',a18,1x,a6,' : mw=',f5.1,', --------, Sum [kg]=',1x,1pe11.4)
    
    if (year_month==1) then 
       write (gol,1112) ym(year_month), provider, sector, msg_comp, msg_mass, summed_emissions; call goPr
    else    
       write (gol,1111) ym(year_month), provider, sector, msg_comp, msg_mass, idate(2), summed_emissions; call goPr
    endif
       
  END SUBROUTINE MSG_EMIS
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:     SCALE_CYCLE
  !
  ! !DESCRIPTION:  Evaluates a daily cycle to be used for the biomass 
  !                burning emissions. It is based on the isoprene daily 
  !                cycle, except that the max peak is at 14:00h LT.
  !                Calculates factors to be multiplied with an 
  !                emission field (e.g. co) in order to simulate a 
  !                diurnal cycle in emissions (caused by solar dependency/ 
  !                human activity), e.g.  
  !      rm(i,j,1,isop) = rm(i,j,1,isop) + e_isop(i,j)*scalef(ipos,j)*dt
  !                with ipos depending on time and longitude.
  !                The scalefactor is calculated for -180 longitude and 
  !                the mean value for ntim timesteps during a day is scaled
  !                to 1. 
  !                The routine should be called only once, since 
  !                the position relative to the sun is not relevant here.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SCALE_CYCLE(ntim, scalef)
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)           :: ntim
    !
    ! !OUTPUT PARAMETERS:
    !
    real, dimension(ntim), intent(out) :: scalef      
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - added without_diurnal_emission_cycle
    !                                (to be removed again :-( )
    !   19 Oct 2010 - Narcisa Banda - changed without_diurnal_emission_cycle to with_....
    !   31 Jan 2013 - Ph. Le Sager - get rid of with_diurnal_emission_cycle!
    !
    ! !REMARKS:
    !   - This routine is called only once (see emission.F90), and only if
    !      emis_bb_trop_cycle is .true.. The later is set with the
    !      input.emis.bb.dailycycle in the rc file (chem.input.rc)
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real      :: pi, piby, obliq, deday, delta, dt, lon, hr, lat, houra, smm, w, sig
    integer   :: i

    ! The calling routine (emission_declare) defines NTIM as follows:
    ! dtime=float(ndyn_max)/(2*tref(region)) ! timestep emissions
    ! ntim=86400/nint(dtime)	             ! number of timesteps in 24 hours.

    pi = acos(-1.0)
    piby = pi/180.0
    ! obliq = 23.45 * piby

    w = 0.2
    sig = 2.0

    dt = 24./ntim   !timestep in hours
    lon = -180.0*piby       !calculate for dateline

    hr =  - 0.5*dt   !shift half an interval back
    if (hr .gt. 12) hr = hr - 24
    ! hr =  hr - 2   !shift two hours. This means that Local-day-time maximum will be at 14:00 h LT.
    smm = 0.0
    do i=1,ntim
       hr = hr + dt
       houra = lon - pi + hr * (2.*pi/24.)
       scalef(i) = w + (1-w)*24.0/(sig*sqrt(2*pi))*exp(-0.5*(((hr-1.5)/sig)**2))

       ! scalef(i) =   max(1+(cos(houra)),0.0)
       smm = smm+scalef(i)/ntim
    end do
    if ( smm > 1e-5 ) then 
       scalef(1:ntim) = scalef(1:ntim)/smm 
    else 
       scalef(1:ntim) = 1.0
    end if

  END SUBROUTINE SCALE_CYCLE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DO_ADD_2D
  !
  ! !DESCRIPTION: General routine to add a 2D emission field (given in kg
  !                /month) into tracer mass array. The mass increase is done
  !                at level L_LEVEL, and for tracer I_TRACER.
  !                
  !               Update accordingly the budget.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DO_ADD_2D( region,   i_tracer,   l_level, is, js, emis_field, &
                        mol_mass, mol_mass_e, status,  xfrac               )
    !
    ! !USES:
    !
    use dims,          only : CheckShape, tref, nsrce
    use dims,          only : sec_month, im, jm, okdebug
    use global_data,   only : mass_dat, region_dat
    use chem_param,    only : ntracet, names
    use tm5_distgrid,  only : get_distgrid, dgrid
#ifdef with_budgets
    use budget_global, only : budg_dat, nzon_vg
#endif

    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)          :: region  
    integer, intent(in)          :: i_tracer 
    integer, intent(in)          :: l_level
    integer, intent(in)          :: is, js
    real, intent(in)             :: emis_field(is:,js:)
    real, intent(in)             :: mol_mass    ! mol mass of tracer field
    real, intent(in)             :: mol_mass_e  ! mol mass of emission field (e.g. NOx emission ar in kgN/month)
    real, intent(in), optional   :: xfrac       ! partial addition 
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)         :: status  
    !
    ! !REVISION HISTORY: 
    !   27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer :: i, j, nzone, nzone_v, i1, i2, j1, j2
    real    :: x, efrac, dtime, month2dt
    real    :: mbef, maft, sume  ! mass BEFore, AFTer, SUM Emiss (debug)
    
    real, dimension(:,:,:,:), pointer :: rm, rzm   

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

    call GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
    
    call CheckShape( (/i2-i1+1,j2-j1+1/), shape(emis_field), status )
    IF_NOTOK_RETURN(status=1)
    
    dtime=float(nsrce)/(2*tref(region))  ! timestep emissions
    month2dt=dtime/sec_month             ! conversion to emission per timestep
   
    rm  => mass_dat(region)%rm
    rzm => mass_dat(region)%rzm

    if(present(xfrac)) then
       efrac = xfrac
    else
       efrac = 1.0
    endif

    if (okdebug) then 
       write (gol, '("--------------------------------------------------------------")')            ; call goPr
       write (gol, '(" DO ADD 2D debug in region:",i3) '    ) region                                ; call goPr
       write (gol, '(" itracer                  :",i3,a8)'  ) i_tracer, ' :'//names(i_tracer)       ; call goPr
       write (gol, '(" xfrac                    :",f6.1)'   ) efrac                                 ; call goPr
       write (gol, '(" mol_mass                 :",f6.1 ) ' ) mol_mass                              ; call goPr
       write (gol, '(" mol_mass_e               :",f6.1 )'  ) mol_mass_e                            ; call goPr
       write (gol, '(" emisfield                :",2e12.4 )') minval(emis_field),maxval(emis_field) ; call goPr
       sume =  sum(emis_field)*month2dt*(mol_mass/mol_mass_e)*efrac
       write (gol, '(" sum emmis                :",e12.4 )' ) sume                                  ; call goPr
       mbef = sum(rm(i1:i2,j1:j2,:,i_tracer))
    endif

    do j=j1,j2
       do i=i1,i2
          x = emis_field(i,j)*month2dt*(mol_mass/mol_mass_e)*efrac
          rm(i,j,l_level,i_tracer) = rm(i,j,l_level,i_tracer) + x
#ifdef slopes
          ! injection of emissions at surface, thus vertical slope more negative
          ! (keep concentration at top of layer the same as before emission)
          if ( i_tracer <= ntracet ) then
             rzm(i,j,l_level,i_tracer) = rzm(i,j,l_level,i_tracer) - x
          end if
#endif

#ifdef with_budgets
          nzone   = budg_dat(region)%nzong(i,j)    !global budget
          nzone_v = nzon_vg(l_level) 
          budemi_dat(region)%emi(i,j,nzone_v,i_tracer) = &
               budemi_dat(region)%emi(i,j,nzone_v,i_tracer) + x/mol_mass*1e3
          if(i_tracer == 1) sum_emission(region) = sum_emission(region) + x  !in kg
#endif
       enddo
    enddo

    if (okdebug) then
       maft = sum(rm(i1:i2,j1:j2,:,i_tracer))
       write (gol, '(" Added amount             :",e12.4 )' ) maft-mbef                ; call goPr
       write (gol, '(" Added amount can be different when inner zoom is present!" )' ) ; call goPr
       if (maft-mbef-sume > 1e-8*sume) then
          write (gol, '(" Warning: error in emission ")' )  ; call goErr
       endif
    endif
    nullify(rm)
    nullify(rzm)

    status=0

  END SUBROUTINE DO_ADD_2D
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DO_ADD_3D
  !
  ! !DESCRIPTION: General routine to add a 3D emission field (given in kg
  !                /month) into tracer mass array. The mass increase is done
  !                for tracer I_TRACER.
  !                
  !               Update accordingly the budget.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DO_ADD_3D( region, i_tracer, is, js, emis_field, mol_mass, mol_mass_e, status, xfrac)
    !
    ! !USES:
    !
    use dims,          only : CheckShape, nsrce, tref
    use dims,          only : sec_month, im, jm, lm, okdebug
    use global_data,   only : mass_dat, region_dat
    use chem_param,    only : ntracet, names
    use tm5_distgrid,  only : get_distgrid, dgrid
#ifdef with_budgets
    use budget_global, only : budg_dat, nzon_vg
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)           :: region   ! region #
    integer, intent(in)           :: i_tracer ! id of tracer to increase
    integer, intent(in)           :: is, js     
    real,    intent(in)           :: emis_field(is:,js:,:)
    real,    intent(in)           :: mol_mass, mol_mass_e
    real,    intent(in), optional :: xfrac                 !partial addition 
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          :: status  
    !
    ! !REVISION HISTORY: 
    !   27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer                             :: i, j, l, nzone, nzone_v, i1,j1,i2,j2
    integer, dimension(3)               :: ubound_e
    real                                :: x, efrac, dtime, month2dt
    real, dimension(:,:,:,:), pointer   :: rm
    
    real :: mbef, maft, sume ! debug

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

    call GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    status = 0
    
    ! check arguments
    ubound_e = ubound(emis_field)
    if(ubound_e(1) /= i2) then
       write(gol,'(A)') 'do_add_3d: im emission not consistent' ; call goPr
       status = 1
    endif
    if(ubound_e(2) /= j2) then
       write(gol,'(A)') 'do_add_3d: jm emission not consistent' ; call goPr
       status = 1
    endif
    if(ubound_e(3) > lm(region)) then
       write(gol,'(A)') 'do_add_3d: lm emission not consistent' ; call goPr
       status = 1
    endif
    IF_NOTOK_RETURN(status=1)
    
    rm => mass_dat(region)%rm

    dtime=float(nsrce)/(2*tref(region))  ! timestep emissions
    month2dt=dtime/sec_month             ! conversion to emission per timestep

    if(present(xfrac)) then
       efrac = xfrac
    else
       efrac = 1.0
    endif

!     if (okdebug) then 
!        write (gol, '("--------------------------------------------------------------")')             ; call goPr
!        write (gol, '(" DO ADD 3D debug in region:",i3) '    ) region                                 ; call goPr
!        write (gol, '(" itracer                  :",i3,a8)'  ) i_tracer, ' :'//names(i_tracer)        ; call goPr
!        write (gol, '(" xfrac                    :",f6.1)'   ) efrac                                  ; call goPr
!        write (gol, '(" mol_mass                 :",f6.1 ) ' ) mol_mass                               ; call goPr
!        write (gol, '(" mol_mass_e               :",f6.1 )'  ) mol_mass_e                             ; call goPr
!        write (gol, '(" emisfield                :",2e12.4 )') minval(emis_field), maxval(emis_field) ; call goPr
!        write (gol, '(" ubound(3)                :",i3 )'    ) ubound_e(3)                            ; call goPr
!        sume = sum(emis_field)*month2dt*(mol_mass/mol_mass_e)*efrac
!        write (gol, '(" sume                     :", e12.4 )') sume                                   ; call goPr
!        mbef = sum(rm(i1:i2,j1:j2,:,i_tracer))
!     endif

    do l=1,ubound_e(3)
       do j=j1,j2
          do i=i1,i2
             x = emis_field(i,j,l)*month2dt*(mol_mass/mol_mass_e)*efrac
             rm(i,j,l,i_tracer) = rm(i,j,l,i_tracer) + x

#ifdef with_budgets
             ! budget___
             nzone=budg_dat(region)%nzong(i,j)    !global budget
             nzone_v=nzon_vg(l) 
             budemi_dat(region)%emi(i,j,nzone_v,i_tracer) = &
                  budemi_dat(region)%emi(i,j,nzone_v,i_tracer) + x/mol_mass*1e3
             if(i_tracer ==1) sum_emission(region) = sum_emission(region) + x  !in kg
#endif
          enddo
       enddo
    enddo  !levels

!     if (okdebug) then
!        maft = sum(rm(i1:i2,j1:j2,:,i_tracer))
!        write (gol, '(" after-before             :", e12.4 )' ) maft-mbef ; call goPr
!        write (gol, '(" END debug output          ")' )                   ; call goPr
!     endif

    nullify(rm)

  END SUBROUTINE DO_ADD_3D
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DO_ADD_3D_CYCLE
  !
  ! !DESCRIPTION: Routine to add a 3D emission field (given in kg/month), 
  !               scaled with a daily cycle. To be used for biomass burning 
  !               emissions.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DO_ADD_3D_CYCLE( region, i_tracer, is, js, emis_field, curr_cycle, mol_mass, mol_mass_e, status, xfrac)
    !
    ! !USES:
    !
    use dims,          only : CheckShape, tref, nsrce, itaur
    use dims,          only : sec_month, im,jm,lm, okdebug
    use dims,          only : ndyn_max,dx, xref, xbeg,dy,yref,ybeg
    use global_data,   only : mass_dat, region_dat
    use chem_param,    only : ntracet, names
    use tm5_distgrid,  only : get_distgrid, dgrid
#ifdef with_budgets
    use budget_global, only : budg_dat, nzon_vg
#endif
    use datetime,      only : tau2date
    use toolbox,       only : dumpfield
    !
    ! !INPUT PARAMETERS:
    !
    integer,            intent(in)           :: region  
    integer,            intent(in)           :: i_tracer 
    integer,            intent(in)           :: is, js     
    real,               intent(in)           :: emis_field(is:,js:,:)
    real, dimension(:), intent(in)           :: curr_cycle
    real,               intent(in)           :: mol_mass     ! first is the mass of tracer field
    real,               intent(in)           :: mol_mass_e   ! mass of emission field (e.g. NOx emission ar in kgN/month)
    real,               intent(in), optional :: xfrac        ! partial addition 
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          :: status  
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - v0 based on do_add_3d
    !   20 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer                                    :: i, j, l, ipos, nzone, nzone_v
    integer,dimension(6)                       :: idater
    integer                                    :: sec_in_day, ntim, itim, i1, j1, i2, j2
    integer, dimension(3)                      :: ubound_e
    real                                       :: x, xlon, xlat, efrac, dtime, month2dt, dtime2
    real,dimension(:,:,:,:),pointer            :: rm
    real :: mbef, maft, sume

    ! --- begin --------------------------------------
    
    call GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    status = 0

    ! check arguments
    ubound_e = ubound(emis_field)
    if(ubound_e(1) /= i2) then
       write(gol,'(A)') 'do_add_3d: im emission not consistent' ; call goPr
       status = 1
    endif
    if(ubound_e(2) /= j2) then
       write(gol,'(A)') 'do_add_3d: jm emission not consistent' ; call goPr
       status = 1
    endif
    if(ubound_e(3) > lm(region)) then
       write(gol,'(A)') 'do_add_3d: lm emission not consistent' ; call goPr
       status = 1
    endif
    IF_NOTOK_RETURN(status=1)

    rm => mass_dat(region)%rm

    dtime      = float(nsrce)/(2*tref(region))             !timestep emissions
    month2dt   = dtime/sec_month                           ! conversion to emission per timestep
    dtime2     = float(ndyn_max)/(2*tref(region))          !timestep emissions based on non-CFL interference (CMK 05/2006)
    ntim       = 86400/nint(dtime2)                        !number of timesteps in 24 hours based on non-CFL interference (CMK 05/2006)
    
    call tau2date(itaur(region),idater)
    sec_in_day = idater(4)*3600 + idater(5)*60 + idater(6) !elapsed seconds this day
    itim       = sec_in_day/nint(dtime2)+1                 !time interval

    if(present(xfrac)) then
       efrac = xfrac
    else
       efrac = 1.0
    endif

!     if (okdebug) then 
!        write (gol, '("--------------------------------------------------------------")')                   ; call goPr
!        write (gol, '(" DO ADD 3D CYCLE debug in region:",i3) '     ) region                                ; call goPr
!        write (gol, '(" itracer                        :",i3,a8)'   ) i_tracer, ' :'//names(i_tracer)       ; call goPr
!        write (gol, '(" xfrac                          :",f6.1)'    ) efrac                                 ; call goPr
!        write (gol, '(" mol_mass                       :",f6.1 ) '  ) mol_mass                              ; call goPr
!        write (gol, '(" mol_mass_e                     :", f6.1 )'  ) mol_mass_e                            ; call goPr
!        write (gol, '(" emisfield                      :", 2e12.4 )') minval(emis_field),maxval(emis_field) ; call goPr
!        write (gol, '(" ubound(3)                      :", i3 )'    ) ubound_e(3)                           ; call goPr
!        sume = sum(emis_field)*month2dt*(mol_mass/mol_mass_e)*efrac
!        write (gol, '(" sume                     :", e12.4 )' ) sume                                        ; call goPr
!        mbef = sum(rm(i1:i2,j1:j2,:,i_tracer))
!        write (gol,*) 'lbounds:', is,i1,js,j1 ; call goPr
!     endif

    do l=1,ubound_e(3)
       do j=j1,j2
          do i=i1,i2
             xlat = ybeg(region) + (j-0.5)*dy/yref(region)
             if (xlat .gt. -20 .and. xlat .lt. 20) then
                ! apply emission cycle in tropics only
                ! itim = 1 and lon = -180 --->position 1
                ! itim = ntim ant lon = -180 --->position ntim, etc.
                ! itim = 1 and lon = 0 ---->position ntim/2
                xlon = xbeg(region) + (i-0.5)*dx/xref(region)
                ipos = 1 + mod(int((xlon+180.)*ntim/360.) + (itim-1),ntim)  !position in array depending on time and lon.
                x = emis_field(i,j,l)*month2dt*(mol_mass/mol_mass_e)*efrac*curr_cycle(ipos)
	     else
                x = emis_field(i,j,l)*month2dt*(mol_mass/mol_mass_e)*efrac
             endif

             rm(i,j,l,i_tracer) = rm(i,j,l,i_tracer) + x

#ifdef with_budgets
             ! budget___
             !             if(region==nregions) then   !sub budget finest region
             !                nzone=nzon(i,j)
             !                nzone_v=nzon_v(l)
             !                budemi(nzone,nzone_v,i_tracer)=budemi(nzone,nzone_v,i_tracer)+x/mol_mass*1e3 !mole
             !             endif
             nzone=budg_dat(region)%nzong(i,j)    !global budget
             nzone_v=nzon_vg(l) 
             budemi_dat(region)%emi(i,j,nzone_v,i_tracer) = &
                  budemi_dat(region)%emi(i,j,nzone_v,i_tracer) + x/mol_mass*1e3
             if(i_tracer ==1) sum_emission(region) = sum_emission(region) + x  !in kg
#endif
          enddo
       enddo
    enddo  !levels


!     if (okdebug) then
!        maft = sum(rm(i1:i2,j1:j2,:,i_tracer))
!        write (gol, '(" after-before             :", e12.4 )' ) maft-mbef ; call goPr
!        !if ( sume > tiny(sume) ) then
! ! This is wrong...
! !          if ( (maft-mbef)/sume < 0.5) then
! !             write (gol, '("  spurious do_add_3d_cycle....  ")' )                   ; call goPr
! !             call dumpfield(1, 'spurious.hdf', im(region), jm(region), ubound_e(3), 1, emis_field, names(i_tracer))
! !          endif
!        !endif
!        write (gol, '(" END debug output          ")' )                   ; call goPr
!     endif

    nullify(rm)

  end subroutine do_add_3d_cycle
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DISTRIBUTE_L44
  !
  ! !DESCRIPTION: Scales input field (field2d) such that 
  !               1) NH temperate vegetation fires are distributed over the
  !                  months: july-august-september 
  !               2) SH temperate vegetation fires are distributed over the
  !                  months: January-February-March
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DISTRIBUTE_L44(month, imdim, jmdim, field2d)
    !
    ! !USES:
    !
    use dims,     only : okdebug
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: month         ! month
    integer, intent(in) :: jmdim, imdim  ! dimension of grid
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - protex
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real,dimension(imdim,jmdim) :: field2d
    real               :: bb_nh,bb_sh                 ! multiplication factor for yearly temperate fires
    integer            :: j,i

    bb_nh=0.
    bb_sh=0.
    if (month==7.or.month==8.or.month==9) bb_nh=3./12.
    if (month==1.or.month==2.or.month==3) bb_sh=3./12.
    do i=1, imdim
       do j=1, jmdim/2
          field2d(i,j)=field2d(i,j)*bb_sh
       enddo
       do j=jmdim/2+1,jmdim
          field2d(i,j)=field2d(i,j)*bb_nh
       enddo
    enddo

    if(okdebug) then
       write(gol,*) 'l44 is distributed' ; call goPr
    endif

  END SUBROUTINE DISTRIBUTE_L44
  !EOC



END MODULE EMISSION_DATA