ecearth3/sources/tm5mp/proj/co2/emission_read__co2.F90

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') 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_READ
!
! !DESCRIPTION: This module provides objects and methods related to
!               CMIP6 and IPCC-AR5 emissions.
!               
!    AR5 netCDF files are provided by CMIP5:
!    
!    There are a few keys in the rc-file which control the behaviour of
!    this module and the data used: 
!    # specify the (main) provider of emission sets
!       input.emis.provider        :  AR5
!    #  where to find the emissions (this will be used by install-emis-ar5)
!       input.emis.dir             :  ${TEMP}/EMIS/AR5
!    #  year of emissions (AR5 emissions will be linearly interpolated)
!       input.emis.year            :  2000
!    #  choose RCP out of RCP26, RCP45, RCP60, RCP85
!       input.emis.AR5.RCP : RCP45 
!  
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_READ  
  !
  ! !USES:
  !
  use GO,              only : gol, goErr, goPr, goLabel
  use emission_data,   only : emis_input_dir_cmip6
  use emission_data,   only : emis_input_dir_ar5
  use emission_data,   only : vd_class_name_len
  use dims,            only : nlon360, nlat180, iglbsfc
  use chem_param,      only : xmc, xmco2
  use Dims,            only : okdebug

  USE MDF,             ONLY : MDF_Open,      MDF_NETCDF,  MDF_READ
  USE MDF,             ONLY : MDF_Inq_VarID, MDF_Get_Var, MDF_Close

  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: emission_read_init, emission_read_done
  public :: emission_ar5_regrid_aircraft
  !public :: numb_2dsec, numb_3dsec
  public :: numb_sectors, sectors_def
  public :: numb_providers, providers_def
  public :: sector_name_len
  public :: emission_cmip6_readsector
  public :: emission_ar5_readCO2
  public :: ar5_dim_3ddata

  public :: sector_type, provider_type
  
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter ::  mname = 'emission_read'  

  ! ------------------------------
  ! global characteristics
  ! ------------------------------
  ! AR5: CO2 emissions from land use are provided at 0.5x0.5, from fossil fuel use at 1x1
  integer, parameter           :: nlat360         = 360    ! number of latitudes for CMIP6, AR5 (0.5deg)
  integer, parameter           :: nlon720         = 720    ! number of longitudes for CMIP6, AR5 (0.5deg)
  integer, parameter           :: sector_name_len = 18     ! length of sector descriptor
  integer, parameter           :: categ_name_len  = 14     ! length of category descriptor
  integer, parameter           :: numb_sectors    = 12     ! number of sectors (All providers!)
  !integer, parameter           :: numb_2dsec      = 10      ! number of 2d sectors (all except aircraft)
  !integer, parameter           :: numb_3dsec      = 1      ! number of 3d sectors (aircraft)
  integer, parameter           :: numb_providers  = 2      ! CMIP6, AR5

  ! Since CMIP6 emissions from aviation are provided on the same vertical grid
  ! as for AR5, this variable is also used for CMIP6:
  integer, parameter           :: ar5_dim_3ddata  = 25     ! number of layers for aircraft data

  ! full list of providers
  character(10), dimension(numb_providers), parameter :: all_providers = &
       & (/ 'CMIP6     ', 'AR5       '/)
  ! Once CO2 will be combined with chemical tracers in one model, 
  ! a separate category for AR5 CO2 could be introduced, because the file format is different
  ! The above is an old comment; the code is not working for AR5 CO2 emissions

  ! List of providers effectively used
  character(10), PUBLIC, allocatable :: used_providers(:)       ! CO2

  ! flag for degenerated cases
  logical, PUBLIC :: has_emis      = .true.
  
  ! ------------------------------
  character(len=15), parameter :: filestr_common_pre  = 'IPCC_emissions'
  character(len=25), parameter :: filestr_common_post = '0.5x0.5.nc'

  ! ------------------------------
  ! identifier of RCPs (RCP26, RCP45,...)
  ! ------------------------------
  character(len=5)             :: filestr_rcpiden  

  !------------------------------
  ! SSP scenario name
  !------------------------------
  character(len=14)             :: ssp_name

  ! ------------------------------
  ! available years and related parameters/variables
  ! ------------------------------
  !
  ! CMIP6
  ! availability (min, max years) - Limit MACC and MEGAN to one year for EC-Earth
  integer, dimension(2), parameter :: cmip6_avail = (/1850, 2100/)
  !---------------------------------------------
  ! CMIP5 CO2 emission data
  !---------------------------------------------
  ! historical data
  !---------------------------------------------
  ! emissions from fossil-fuel use are provided as monthly gridded fields for 1751-2007 (1x1 degree)
  ! (http://dods.ipsl.jussieu.fr/cpipsl/ANDRES/) 
  ! emissions from land use are provided as yearly gridded fields for 1850-2005 (0.5x0.5 degree)
  ! (http://www.mpimet.mpg.de/en/wissenschaft/land-im-erdsystem/...
  !         .../wechselwirkung-klima-biogeosphaere/landcover-change-emission-data.html)
  ! we only use the data for the years 1850-2005:
  integer, dimension(2), parameter :: ar5_avail   = (/1850, 2005/)

  ! global totals (Pg C/yr) are provided as well:
  ! the numbers for the reference year 2005 are:
  !real, parameter :: co2_ff_ref = 7.6137 ! Pg C/yr as provided 
  real, parameter :: co2_ff_ref = 7.617692 ! Pg C/yr as calculated from the gridded fields
  !real, parameter :: co2_lu_ref = 1.196  ! Pg C/yr as provided
  real, parameter :: co2_lu_ref = 1.4673 ! Pg C/yr as calculated from the gridded field
  ! for the land-use emissions up to 2001 the totals calculated from the gridded fields agree well 
  ! with the totals given by the data provider
  ! however, for the last four years 2002-2005 the gridded fields give substantially higher totals
  ! this suggests that the emission totals provided for land use have been harmonized with the RCPs 
  ! while the gridded fields have not
  real            :: co2_ref
  
  ! future CO2 emissions for the RCPs (2006-2100) are provided as yearly totals (Pg C/yr)
  ! we currently use the global totals, but regional totals are available as well
  ! values obtained from the IIASA RCP website (http://tntcat.iiasa.ac.at/RcpDb/)
  ! for 2006-2100 we combined these totals with the spatial distribution for 2005

  integer, parameter               :: ar5_nr_avail_yrs = 11
  integer, dimension(ar5_nr_avail_yrs), parameter :: &
       ar5_avail_yrs = (/ 2005, 2010, 2020, 2030, 2040, &
                          2050, 2060, 2070, 2080, 2090, 2100 /)

  real, dimension(ar5_nr_avail_yrs), parameter :: &
       co2ff_rcp26 = (/ 7.971,  8.821,  9.288,  7.157,  4.535,  3.186,  1.419,  0.116, -0.433, -0.870, -0.931/), & 
       co2ff_rcp60 = (/ 7.971,  8.512,  8.950,  9.995, 11.554, 13.044, 14.824, 16.506, 17.281, 14.313, 13.753/), &
       co2ff_rcp45 = (/ 7.971,  8.607,  9.872, 10.953, 11.338, 11.031,  9.401,  7.118,  4.182,  4.193,  4.203/), &
       co2ff_rcp85 = (/ 7.971,  8.926, 11.538, 13.839, 16.787, 20.205, 23.596, 25.962, 27.406, 28.337, 28.740/), &
       ! for 2000 and 2005 the global total fossil-fuel emissions for the RCPs 
       ! are 2.7% resp. 5% higher than the totals given by the provider of the historical dataset
       ! this suggests that this dataset has not been harmonized with the RCPs
       co2lu_rcp26 = (/ 1.196,  1.056,  0.973,  0.789,  0.489,  0.201,  0.615,  0.538,  0.550,  0.602,  0.511/), &
       co2lu_rcp60 = (/ 1.196,  0.877,  0.406, -0.557, -0.714, -0.464, -0.258, -0.029,  0.244,  0.242,  0.181/), &
       co2lu_rcp45 = (/ 1.196,  0.911,  0.341,  0.216,  0.199,  0.249,  0.184,  0.104,  0.008,  0.027,  0.046/), &
       co2lu_rcp85 = (/ 1.196,  1.044,  0.906,  0.715,  0.645,  0.576,  0.501,  0.412,  0.309,  0.194,  0.077/)
  real, dimension(ar5_nr_avail_yrs) :: co2_rcp

  logical, dimension(:), allocatable :: ltimeind

  character(len=7)        :: ar5ff_coverage = 'monthly'
  character(len=7)        :: ar5lu_coverage = 'yearly '

  ! ------------------------------
  ! gridbox area (to be read only once per proc)
  ! ------------------------------
  character(len=25),parameter       :: cmip6_filestr_gridboxarea = 'gridbox_area.nc          '
  character(len=25),parameter       :: ar5_filestr_gridboxarea = 'gridbox_area.nc'
  logical, save                     :: area_found_05   
  real, dimension(:,:), allocatable :: gridbox_area_05  ! gridbox area on 0.5x0.5 deg - used for AR5


  ! -----------------------
  ! data type for sectors
  ! -----------------------
  type sector_type
     sequence
     character(len=sector_name_len)    :: name      ! name of sector
     character(len=categ_name_len)     :: catname   ! name of category to be found in
     logical                           :: f3d       ! 3d-data y/n
     character(len=vd_class_name_len)  :: vdisttype ! vertical distribution type (equal to "classes" still to be defined)
     character(len=8)                  :: prov      ! provider of information (AR5)
     character(len=26), dimension(:), pointer :: species ! list of species available for that sector (use for AR5 only)
  end type sector_type

  type provider_type
     character(len=8)          :: name
     integer                   :: nsect2d, nsect3d
  end type provider_type

  type(sector_type),   dimension(numb_sectors)   :: sectors_def
  type(provider_type), dimension(numb_providers) :: providers_def
  

  !
  ! !REVISION HISTORY: 
  !    1 Oct 2010 - Achim Strunk  - v0 for AR5
  !   19 Jun 2012 - P. Le Sager   - cosmetic for lon-lat MPI domain decomposition
  !                                 (all reading/regridding on root for now)
  !   20 Nov 2012 - Ph. Le Sager  - defined and build lists of used providers
  !                               - deal with inventories years availability
  !                               - switch to MDF interface to read data
  !
  ! !TODO:
  !   - should be renamed something like "emission_inventories" or "emiss_providers"
  !   - and need to get a **SEPARATE** module for each inventories, before it
  !     becomes unmanageable again
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_READ_INIT
  !
  ! !DESCRIPTION: Initialise reading related parameters and 
  !               allocate needed arrays
  !
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_READ_INIT( rcF, status )
    !
    ! !USES:
    !
    use GO,            only : TrcFile, ReadRc
    use partools,      only : isRoot
    use emission_data, only : LCMIP6 
    use emission_data, only : LAR5
    use meteodata,     only : set, gph_dat
    use dims,          only : im, jm, lm, nregions
    !
    ! !INPUT PARAMETERS:
    !
    type(TrcFile)        ::  rcF
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk  - v0 for AR5
    !   20 Nov 2012 - Ph. Le Sager  - build lists of used providers
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter :: rname=mname//'/emission_read_init'
    integer                     :: isect, iprov, nused, region
    logical                     :: mask(numb_providers)

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

    if (LCMIP6) then
      ! The SSP scenario data are provided for the years 2015-2100.
      ! SSP5-8.5 is the official scenario run in C4MIP and is therefore set as the default here
      ! SSP1-1.9 and SSP1-2.6 (with negative emissions) and SSP2-4.5 are now supported
      call ReadRc( rcF, 'input.CMIP6.SSP', ssp_name, status, default = 'SSP5-8.5' )
      IF_ERROR_RETURN(status=1)

      write(gol,'("SSP CMIP6 future scenario for emissions: ",a)') trim(ssp_name); call goPr

    else if (LAR5) then
      call ReadRc( rcF, 'input.emis.AR5.RCP', filestr_rcpiden, status, default='RCP26' )
      IF_ERROR_RETURN(status=1)
    endif

    ! ------------------
    ! build list of used providers
    ! ------------------

    ! CO2
    mask = (/ LCMIP6, LAR5 /)
    nused = count(mask)
    if (nused /= 0) then
       allocate( used_providers(nused) )
       used_providers = pack( all_providers, mask)
    else
       has_emis = .false.
    end if

    ! info
    if ( has_emis ) then 
       write(gol,*) 'EMISS-INFO - Emissions providers used for CO2: ', used_providers      ; call goPr
    else
       write(gol,*) 'EMISS-INFO - Emissions providers used for CO2: NONE'                  ; call goPr
    end if
           
    ! ------------------
    ! initialise sectors
    ! ------------------
    ! Type sequence is (name, category, is_3D_data, vdisttype, providers)
    sectors_def( 1) = sector_type('emiss_ff  ', 'anthropogenic ', .false., 'combenergy  ', 'AR5     ', NULL() ) ! Fossil Fuel
    sectors_def( 2) = sector_type('emiss_lu  ', 'anthropogenic ', .false., 'nearsurface ', 'AR5     ', NULL() ) ! Land Use (assumed near surface for the moment, but that is open for discussion)

    ! CMIP6
    sectors_def( 3) = sector_type('ENE       ', 'anthropogenic ', .false., 'combenergy  ', 'CMIP6   ', NULL() ) ! Energy sector
    sectors_def( 4) = sector_type('RCO       ', 'anthropogenic ', .false., 'combrescom  ', 'CMIP6   ', NULL() ) ! Residential, commercial and other
    sectors_def( 5) = sector_type('IND       ', 'anthropogenic ', .false., 'industry    ', 'CMIP6   ', NULL() ) ! Industrial sector
    sectors_def( 6) = sector_type('WST       ', 'anthropogenic ', .false., 'waste       ', 'CMIP6   ', NULL() ) ! Waste treatment and disposal
    sectors_def( 7) = sector_type('AGR       ', 'anthropogenic ', .false., 'surface     ', 'CMIP6   ', NULL() ) ! Agriculture (excl. agricultural waste burning, which is included in CMIP6 biomass burning emissions)
    ! CO2 emissions from AGR are zero in CMIP6
    sectors_def( 8) = sector_type('SLV       ', 'anthropogenic ', .false., 'nearsurface ', 'CMIP6   ', NULL() ) ! Solvents production and application
    sectors_def( 9) = sector_type('TRA       ', 'anthropogenic ', .false., 'surface     ', 'CMIP6   ', NULL() ) ! Transportation sector (land)
    sectors_def(10) = sector_type('SHP       ', 'ships         ', .false., 'nearsurface ', 'CMIP6   ', NULL() ) ! International shipping
    sectors_def(11) = sector_type('AIR       ', 'aircraft      ', .true. , 'aircraft    ', 'CMIP6   ', NULL() ) ! Aircraft
    sectors_def(12) = sector_type('NEG       ', 'anthropogenic ', .false., 'surface     ', 'CMIP6   ', NULL() ) ! Negative CO2 Emissions

    ! -------------------------
    ! initialise providers info
    ! ------------------------
    do iprov = 1, numb_providers
       providers_def(iprov)%name    = all_providers(iprov)
       providers_def(iprov)%nsect2d = count( (sectors_def%prov == all_providers(iprov)) .and. (sectors_def%f3d .eqv. .false.))
       providers_def(iprov)%nsect3d = count( (sectors_def%prov == all_providers(iprov)) .and. (sectors_def%f3d .eqv. .true.))

       if(okdebug) then
         write(gol,'("EMISS-INFO - Inventory ",a," has ",i3, " 2d-sectors, and ",i3," 3d-sectors")')&
              & all_providers(iprov), providers_def(iprov)%nsect2d, providers_def(iprov)%nsect3d ; call goPr
       endif
    end do

    ! -------------------------
    ! initialise GeopPotential Height on 1x1
    ! ------------------------
    do region=1, nregions
       call Set( gph_dat(region), status, used=.true. )
    end do
    
    ! ----------------------------------------
    ! allocate gridbox_area arrays
    ! ----------------------------------------
    allocate( gridbox_area_05( nlon720, nlat360 ) )
    ! OK 
    status = 0

  END SUBROUTINE EMISSION_READ_INIT
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_READ_DONE
  !
  ! !DESCRIPTION: Free allocated arrays.
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine emission_read_done(  status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - v0
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    deallocate( gridbox_area_05 )
    deallocate( used_providers )

    ! OK
    status = 0

  END SUBROUTINE EMISSION_READ_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:    EMISSION_COARSEN_TO_1X1
  !
  ! !DESCRIPTION: Coarsen the gridded information to 1x1 deg.
  !               (taken from GEMS/MACC repository)
  !\\
  !\\
  ! !INTERFACE:
  !
  function emission_coarsen_to_1x1( emis_in, dim_nlon, dim_nlat, shift_lon, status )
    !
    ! !RETURN VALUE:
    !
    real, dimension(360,180) :: emission_coarsen_to_1x1
    !
    ! !INPUT PARAMETERS:
    !
    integer,  intent(in)    :: dim_nlon
    integer,  intent(in)    :: dim_nlat
    real,     intent(in)    :: emis_in(dim_nlon, dim_nlat)
    logical,  intent(in)    :: shift_lon
    !
    ! OUTPUT PARAMETERS:
    !
    integer          , intent(out)   :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk  - v0 for AR5
    !    1 Dec 2011 - Narcisa Banda - works for any input resolution lower than 1x1 
    !                                 if 1x1 can be divided into exact number of gridcells (no interpolation)
    !    1 Jul 2012 - Narcisa Banda - added the shift_lon logical flag:
    !                                 true if the data is read on longitudes [0,360] (then they need to be shifted on [-180,180])
    !                                 false if the data is read already on [-180,180]
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer                 ::  i, j
    integer                 ::  nri, nrj

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

    ! combine grid cells :
    !   from  [   0,360]x[-90,90]             001:360,361:720    001:360
    !   to    [-180,180]x[-90,90]     001:180,181:360            001:180

    if ((mod(dim_nlon, 360) /= 0 ) .or. (mod(dim_nlat, 180) /= 0)) then
       write(gol,*) 'coarsening of emissions to 1x1 does not work for this resolution'
       status = 1 
       return
    endif

    nri = dim_nlon/360
    nrj = dim_nlat/180

    if (shift_lon) then 
    ! combine grid cells :
    !   from  [   0,360]x[-90,90]             001:360,361:720    001:360
    !   to    [-180,180]x[-90,90]     001:180,181:360            001:180

       do j = 1, 180
          ! west half
          do i = 1, 180
             emission_coarsen_to_1x1(i,j) = sum(emis_in(nri*180+nri*i-nri+1:nri*180+nri*i,nrj*j-nrj+1:nrj*j))
          end do
          ! east half
          do i = 1, 180
             emission_coarsen_to_1x1(180+i,j) = sum(emis_in(nri*i-nri+1:nri*i,nrj*j-nrj+1:nrj*j))
          end do
       end do
    else
       do j=1, 180
          do i=1, 360
             emission_coarsen_to_1x1(i,j) = sum(emis_in(nri*i-nri+1:nri*i,nrj*j-nrj+1:nrj*j))
          end do
       end do
    endif
    ! ok
    status = 0

  end function emission_coarsen_to_1x1
  !EOC

  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:    VALID_YEAR
  !
  ! !DESCRIPTION: return a valid year for an emission inventory, based on
  !               requested year. 
  !\\
  !\\
  ! !INTERFACE:
  !
  FUNCTION VALID_YEAR( iyear, iminmax, provider_name, verbose)
    !
    ! !RETURN VALUE:
    !
    integer :: valid_year
    !
    ! !INPUT PARAMETERS:
    !
    integer,          intent(in) :: iyear, iminmax(2)
    character(len=*), intent(in) :: provider_name
    logical,          intent(in) :: verbose
    !
    ! !REVISION HISTORY: 
    !   26 Nov 2012 - Ph. Le Sager - v0
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    valid_year = MIN(iminmax(2),MAX(iyear,iminmax(1)))

    ! info only once a year, and per inventory
    if (verbose) then
       write(gol,'(a,i4," (avail: ",i4,"-",i4,")")') ' EMISS-INFO - EMISS YEAR for '//trim(provider_name)//' : ', &
            valid_year, iminmax ; call goPr       
    end if

  END FUNCTION VALID_YEAR
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_CMIP6_READSECTOR
  !
  ! !DESCRIPTION: Reading one sector of the files for the requested month and 
  !               returning an interpolated 3d emission field (d3data)
  !               and, for the CMIP6 2-D sectors, an interpolated 2d field
  !               containing emissions from solid biofuel combustion (d2data_bf). 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine emission_cmip6_ReadSector( comp, iyear, imonth, sector, d3data, status, d2data_bf )
    !
    use dims         , only          : icalendo
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*) , intent(in)    :: comp
    integer          , intent(in)    :: iyear 
    integer          , intent(in)    :: imonth 
    integer          , intent(in)    :: sector
    !
    ! !OUTPUT PARAMETERS:
    !
    integer          , intent(out)   :: status
    real, dimension(nlon360,nlat180,ar5_dim_3ddata), intent(out)   :: d3data
    real, dimension(nlon360,nlat180), intent(out), optional :: d2data_bf
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter :: rname = mname//'/emission_cmip6_readsector'
    character(len=256)          :: fname
    character(len=256)          :: fname_gridboxarea
    character(len=256)          :: varfilename, varname, secname
    integer                     :: lt, year, startyear
    character(len=25)           :: sec_str, sec_str2 
    character(len=13)           :: time_str 
    character(len=128)          :: source_str
    character(len=50)           :: version_str
    logical                     :: existfile
    character(len=4)            :: cyear       
    logical                     :: first=.true.

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

    ! only allow CO2 component
    if ( index(comp,'CO2') /= 1 ) then
       write (gol,'(a)') ' CMIP6 emissions of component '//trim(comp)//' not supported in co2 version' ; call goErr
       status=1; TRACEBACK; return 
    endif
    if (present(d2data_bf)) then
       write (gol,'("CMIP6 emissions of solid biofuel combustion (d2data_bf) not supported in co2 version")') ; call goErr
       status=1; TRACEBACK; return 
    endif

    ! initialise target array
    d3data = 0.0
    if (present(d2data_bf)) d2data_bf = 0.0

    ! read in gridbox-area; once per CPU
    if( .not. area_found_05 ) then
       fname_gridboxarea = trim(emis_input_dir_cmip6)//'/'//trim(cmip6_filestr_gridboxarea)
       call emission_ReadGridboxArea(fname_gridboxarea, 'gridbox_area', gridbox_area_05, &
             & nlon720, nlat360, status )
       IF_NOTOK_RETURN(status=1)
       area_found_05=.true.
    endif

    ! deal with out-of-bounds requested years
    year = valid_year( iyear, cmip6_avail, 'CMIP6', first)
    first=.false.

    if ( trim(sectors_def(sector)%catname) == 'aircraft' .and. year < 1920 ) then

      ! CMIP6 aircraft emissions before 1920 are zero and not read anymore
      d3data(:,:,:) = 0. 

    else

      ! cyear will contain strings with the years
      write(cyear,'(I4.4)') year

      ! ------------------------
      ! construct filename
      ! e.g.: <emisdir>/NOx-em-AIR-anthro_input4MIPs_emissions_CMIP_CEDS-v2016-06-18_gr_175001-179912.nc
      ! ------------------------

      if ( index(comp,'CH4') /= 1 ) then
        if (year >= 1750 .and. year < 1800) then
          time_str='175001-179912'
          startyear=1750
        else if (year >= 1800 .and. year < 1850) then
          time_str='180001-184912'
          startyear=1800
        else if (year >= 1850 .and. year < 1851) then
          time_str='185001-185012'
          startyear=1850
        else if (year >= 1851 .and. year < 1900) then
          time_str='185101-189912'
          startyear=1851
        else if (year >= 1900 .and. year < 1950) then
          time_str='190001-194912'
          startyear=1900
        else if (year >= 1950 .and. year < 2000) then
          time_str='195001-199912'
          startyear=1950
        else if (year >= 2000 .and. year < 2015) then
          time_str='200001-201412'
          startyear=2000
        else if (year >= 2015 .and. year <= 2100) then
          time_str='201501-210012'
          startyear=2015
        else
          write (gol,'("CMIP6 emissions beyond 2100 not available")') ; call goErr
          status=1; TRACEBACK; return 
        endif

        if (year >= 1750 .and. year < 2015) then
          if (trim(sectors_def(sector)%catname) == 'aircraft') then
            if ( index(comp,'SO2') /= 1 ) then
              version_str='2017-08-30'
            else
              ! SO2 aicraft emissions have had another update in Oct. 2017
              version_str='2017-10-05'
            endif
          else
            version_str='2017-05-18'
          endif
        else if (year >= 2015 .and. year <=2100) then
          version_str='1-1'
        else
          write (gol,'("CMIP6 emissions beyond 2100 not available")') ; call goErr
          status=1; TRACEBACK; return
        endif
      else
        ! CH4
        if (year >= 1750 .and. year < 1850) then
          write (gol,'("WARNING - Anthropogenic emissions of CH4 not available prior to 1850")') ; call goPr 
          write (gol,'("WARNING - 1850 values are used")') ; call goPr
          year = 1850 
        endif

        if (year >= 1850 .and. year < 1970) then 
          time_str='185001-196012'
          startyear=1850
          version_str='2017-05-18-supplemental-data'
        else if (year >= 1970 .and. year < 2015) then
          time_str='197001-201412'
          startyear=1970
          version_str='2017-05-18'
        else if (year >= 2015 .and. year <= 2100) then
          time_str='201501-210012'
          startyear=2015
          version_str='1-1'
        else
          write (gol,'("CMIP6 emissions beyond 2100 not available")') ; call goErr
          status=1; TRACEBACK; return
        endif

      endif

      if (year <= 2014 ) then
        source_str='input4MIPs_emissions_CMIP_CEDS'
      else
        select case(trim(ssp_name))
        case("SSP1-1.9")
          source_str='input4MIPs_emissions_ScenarioMIP_IAMC-IMAGE-ssp119'
        case("SSP1-2.6")
          source_str='input4MIPs_emissions_ScenarioMIP_IAMC-IMAGE-ssp126'
        case("SSP2-4.5")
          source_str='input4MIPs_emissions_ScenarioMIP_IAMC-MESSAGE-GLOBIOM-ssp245'
        case("SSP5-8.5")
          source_str='input4MIPs_emissions_ScenarioMIP_IAMC-REMIND-MAGPIE-ssp585'
        case default
          write (gol,'("Emissions not implemented for scenario: ", a)') trim(ssp_name); call goErr
          status=1; TRACEBACK; return
        end select
      endif

      ! construct emissions filename
      if ( trim(sectors_def(sector)%catname) == 'anthropogenic' .or. &
           trim(sectors_def(sector)%catname) == 'ships' ) then
        if ( index(comp,'VOC') == 1 ) then
          ! individual VOC species
          sec_str='em-speciated-VOC-anthro'
          sec_str2='em_speciated_VOC_anthro'
          version_str=trim(version_str)//'-supplemental-data'
        else
          sec_str='em-anthro'
          sec_str2='em_anthro'
        endif
        varname=trim(comp)//'_'//trim(sec_str2)
        ! change dash to underscore in few cases
        if ( index(varname, 'VOC') == 1 )        varname(6:6)= '_'
        if ( index(varname, 'hexanes-pl') == 7 ) varname(7:16) = 'hexanes_pl'
        if ( index(varname, 'other-') == 7 )     varname(7:12) = 'other_'
      else if ( trim(sectors_def(sector)%catname) == 'aircraft' ) then
        sec_str='em-AIR-anthro'
        sec_str2='em_AIR_anthro'
        varname=trim(comp)//'_'//trim(sec_str2)
      else
        write (gol,'("Invalid CMIP6 sector")') ; call goErr
        status=1; TRACEBACK; return
      endif
      varfilename=trim(comp)//'-'//trim(sec_str)

      ! For NO, the species name in the file name has to be set to NOx
      fname = trim(emis_input_dir_cmip6)           //'/'//   &
              trim(varfilename)                    //'_'//   &
              trim(source_str)                     //'-'//   &
              trim(version_str)                    //'_'//   &
              'gn'                                 //'_'//   &
              trim(time_str)                            //   &
              '.nc'

      ! test existence of file
      inquire( file=trim(fname), exist=existfile)
      if( .not. existfile ) then 
        write (gol,'(" CMIP6 file `",a,"` not found ")') trim(fname); call goErr
        status = 1; TRACEBACK; return 
      end if

      ! ------------------------------------------------
      ! data record is read by emission_cmip6_Read2/3DRecord
      secname = sectors_def(sector)%name
      ! distinguish 2d/3d sectors
      if( sectors_def(sector)%f3d ) then 
        d3data(:,:,:) = emission_cmip6_Read3DRecord( fname, varname, secname, imonth, year, startyear, status )
      else
        d3data(:,:,1) = emission_cmip6_Read2DRecord( fname, varname, secname, imonth, year, startyear, status )

      endif

    endif

    IF_NOTOK_RETURN(status=1)

  end subroutine emission_cmip6_ReadSector
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:    EMISSION_CMIP6_READ2DRECORD
  !
  ! !DESCRIPTION: Read a single 2d record of a given file and 
  !               return a 2d emission field interpolated on 1x1 grid.
  !\\
  !\\
  ! !INTERFACE:
  !
  function emission_cmip6_Read2DRecord( fname, varname, secname, imonth, year, startyear, status )
    !
    ! !RETURN VALUE:
    !
    real                  :: emission_cmip6_Read2DRecord(nlon360,nlat180)
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*),                intent(in) :: fname, varname
    character(len=sector_name_len),  intent(in) :: secname
    integer,                         intent(in) :: imonth, year, startyear
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,                        intent(out) :: status
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/emission_cmip6_Read2DRecord'
    character(len=256)          ::  fname2

    integer :: fid, varid, isec
    integer :: fid2, varid2, year_first, year_second
    !real    :: emis_in(nlon720, nlat360, 1)
    real    :: emis_in(nlon720, nlat360, 1, 1)
    real, allocatable :: emis_help(:,:,:,:)
    real    :: x

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

    emission_cmip6_Read2DRecord = 0.0

    select case( trim(secname) )
      case ('AGR')
        isec=0
      case ('ENE')
        isec=1
      case ('IND')
        isec=2
      case ('TRA')
        isec=3
      case ('RCO')
        isec=4
      case ('SLV')
        isec=5
      case ('WST')
        isec=6
      case ('SHP')
        isec=7
! Negative emissions are only present in scenario files (year>2014) so make sure we don't read them when reading the historical forcings
      case ('NEG')
        isec=8
        if ( year < 2015 ) then
           if( okdebug ) then
              write (gol,'("EMISS - CMIP6 - no `",a,"` emissions to read in file ",a)') &
                   secname, trim(fname) ; call goErr
           endif
           status = 0
           return
        endif
      case default
        write (gol,'("EMISS - CMIP6 - no `",a,"` emissions in file ",a)') &
            secname, trim(fname) ; call goErr
       status=1; TRACEBACK; return
    end select

    ! initialise

    CALL MDF_Open( TRIM(fname), MDF_NETCDF, MDF_READ, fid, status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Inq_VarID( fid, TRIM(varname), varid, status )
    IF_ERROR_RETURN(status=1)

    if ( varid < 0 ) then

       write (gol,'("EMISS - CMIP6 - no `",a,"` emissions in file ",a)') &
            varname, trim(fname) ; call goErr
       status=1; TRACEBACK; return

    else
       if( okdebug ) then
          write (gol,'("EMISS-INFO - CMIP6 - found `",a,"` emissions in file ",a)') &
               trim(varname), trim(fname) ; call goPr 
       endif

       ! SSP scenario emissions are provided for 2015, 2020, 2030, ... 2100
       if (year >= 2015 .and. year < 2020) then

         ! First year of the period 
         year_first = 2015

         CALL MDF_Get_Var( fid, varid, emis_in, status, start=(/1,1,isec+1,imonth/) )
         IF_NOTOK_RETURN(status=1)

         if (year /= year_first) then

            ! Also read data for 2020
            ! and apply a linear interpolation between the two years
            allocate(emis_help(nlon720, nlat360, 1, 1))
            year_second = 2020

            ! Read data for 2020 by raising the index by 12
            CALL MDF_Get_Var( fid, varid, emis_help, status, start=(/1,1,isec+1,imonth+12/) )

            ! Interpolate montly data between the two provided years
            x = float(year-year_first)/5.
            emis_in(:,:,1,1) = (1.-x) * emis_in(:,:,1,1) + x * emis_help(:,:,1,1)

            deallocate(emis_help)

         endif

       else if (year >= 2020) then

         ! First year of the decade:
         year_first = int(year/10) * 10

         ! Change to reference year to 2020 by raising the index by 12
         CALL MDF_Get_Var( fid, varid, emis_in, status, start=(/1,1,isec+1,imonth+12+12*(year_first-2020)/10/) )
         IF_NOTOK_RETURN(status=1)

         if (year /= year_first) then

            ! Also read data for the first year of the next decade
            ! and apply a linear interpolation between the two years
            allocate(emis_help(nlon720, nlat360, 1, 1))
            year_second = year_first + 10

            CALL MDF_Get_Var( fid, varid, emis_help, status, start=(/1,1,isec+1,imonth+12+12*(year_second-2020)/10/) )

            ! Interpolate monthly data between the two provided years
            x = float(year-year_first)/10.
            emis_in(:,:,1,1) = (1.-x) * emis_in(:,:,1,1) + x * emis_help(:,:,1,1)

            deallocate(emis_help)

         endif

       else 

         ! read yearly emissions directly from file
         CALL MDF_Get_Var( fid, varid, emis_in, status, start=(/1,1,isec+1,imonth+12*(year-startyear)/) )
         IF_NOTOK_RETURN(status=1)

       endif

       ! convert from kg(species)/m^2/s  to  kg(species)/gridbox/s
       emis_in(:,:,1,1) = emis_in(:,:,1,1) * gridbox_area_05
       
       ! now coarsen to nlon360,nlat180
       emission_cmip6_Read2DRecord = emission_coarsen_to_1x1( emis_in(:,:,1,1), nlon720, nlat360,.false., status )
       IF_NOTOK_RETURN(status=1)

    end if

    CALL MDF_Close( fid, status )
    IF_NOTOK_RETURN(status=1)
 
    status = 0
    return

  end function emission_cmip6_Read2DRecord
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:    EMISSION_CMIP6_READ3DRECORD
  !
  ! !DESCRIPTION: read one 3D sector
  ! 
  !\\
  !\\
  ! !INTERFACE:
  !
  function emission_cmip6_Read3DRecord( fname, varname, secname, imonth, year, startyear, status )
    !
    ! !RETURN VALUE:
    !
    real :: emission_cmip6_Read3DRecord(nlon360,nlat180,ar5_dim_3ddata)
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    character(len=*),                   intent(in)    :: fname, varname
    character(32),                      intent(in)    :: secname
    !
    ! !INPUT PARAMETERS:
    !
    integer,                            intent(in)    :: imonth, year, startyear
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,                            intent(out)   :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer                 :: fid, varid, lk
    real, dimension(nlon720,nlat360,ar5_dim_3ddata,1) :: emis_in
    integer :: fid2, varid2, year_first, year_second
    real, allocatable :: emis_help(:,:,:,:)
    real    :: x

    real, parameter      :: layer_depth = 610. ! fixed height (m) of the vertical levels
                                               ! on which the CMIP6 aircraft emissions
                                               ! are provided.
    
    ! --- begin -----------------------------------

    ! initialise
    emission_cmip6_Read3DRecord = 0.0

    CALL MDF_Open( TRIM(fname), MDF_NETCDF, MDF_READ, fid, status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Inq_VarID( fid, TRIM(varname), varid, status )
    IF_ERROR_RETURN(status=1)

    if ( varid < 0 ) then
       write (gol,'("EMISS - CMIP6 - no `",a,"` emissions in file ",a)') &
            secname, trim(fname) ; call goErr
       status=1; TRACEBACK; return
    else
       if( okdebug ) then
          write (gol,'("EMISS-INFO - CMIP6 - found `",a,"` emissions in file ",a)') &
               secname, trim(fname) ; call goPr
       endif

       ! SSP scenario emissions are provided for 2015, 2020, 2030, ... 2100
       if (year >= 2015 .and. year < 2020) then

         ! First year of the period 
         year_first = 2015

         CALL MDF_Get_Var( fid, varid, emis_in, status, start=(/1,1,1,imonth/) )
         IF_NOTOK_RETURN(status=1)

         if (year /= year_first) then

            ! Also read data for 2020
            ! and apply a linear interpolation between the two years
            allocate(emis_help(nlon720, nlat360, ar5_dim_3ddata, 1))
            year_second = 2020

            ! Read data for 2020 by raising the index by 12
            CALL MDF_Get_Var( fid, varid, emis_help, status, start=(/1,1,1,imonth+12/) )

            ! Interpolate montly data between the two provided years
            x = float(year-year_first)/5.
            emis_in(:,:,:,1) = (1.-x) * emis_in(:,:,:,1) + x * emis_help(:,:,:,1)

            deallocate(emis_help)

         endif

       else if (year >= 2020) then

         ! First year of the decade:
         year_first = int(year/10) * 10

         ! Change to reference year to 2020 by raising the index by 12
         CALL MDF_Get_Var( fid, varid, emis_in, status, start=(/1,1,1,imonth+12+12*(year_first-2020)/10/) )
         IF_NOTOK_RETURN(status=1)

         if (year /= year_first) then

            ! Also read data for the first year of the next decade
            ! and apply a linear interpolation between the two years
            allocate(emis_help(nlon720, nlat360, ar5_dim_3ddata, 1))
            year_second = year_first + 10

            CALL MDF_Get_Var( fid, varid, emis_help, status, start=(/1,1,1,imonth+12+12*(year_second-2020)/10/) )

            ! Interpolate monthly data between the two provided years
            x = float(year-year_first)/10.
            emis_in(:,:,:,1) = (1.-x) * emis_in(:,:,:,1) + x * emis_help(:,:,:,1)

            deallocate(emis_help)

         endif

       else

         CALL MDF_Get_Var( fid, varid, emis_in, status, start=(/1,1,1,imonth+12*(year-startyear)/) )
         IF_NOTOK_RETURN(status=1)
        
       endif

       do lk = 1, ar5_dim_3ddata
          ! convert from kg(species)/m^2/s to kg(species)/m/s;
          ! Note that CMIP6 aircraft emissions are given in kg(species)/m^2/s,
          ! while AR5 aircraft emissions are given in kg(species)/m^3/s.
          ! In order to be able to use the same vertical regridding routine lateron, 
          ! we convert to the same unit and include a division by the layer depth.
          emis_in(:,:,lk,1) = emis_in(:,:,lk,1) * gridbox_area_05 / layer_depth
          
          ! now coarsen to nlon360,nlat180
          emission_cmip6_Read3DRecord(:,:,lk) = emission_coarsen_to_1x1( emis_in(:,:,lk,1) ,&
               & nlon720, nlat360, .false., status )
          IF_NOTOK_RETURN(status=1)

       end do

    end if
    
    CALL MDF_Close( fid, status )
    IF_NOTOK_RETURN(status=1)

    status = 0

    return

  end function emission_cmip6_Read3DRecord
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_AR5_READCO2
  !
  ! !DESCRIPTION: Reading one sector of the files to be interpolated and 
  !               returning an interpolated 3d emission field (d3data)
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine emission_ar5_ReadCO2( comp, iyear, imonth, sector, d3data, status )
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*) , intent(in)    :: comp
    integer          , intent(in)    :: iyear 
    integer          , intent(in)    :: imonth 
    integer          , intent(in)    :: sector
    !
    ! !OUTPUT PARAMETERS:
    !
    integer          , intent(out)   :: status
    real, dimension(nlon360,nlat180,ar5_dim_3ddata), intent(out)   :: d3data
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - v0
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter :: rname = mname//'/emission_ar5_readCO2'
    character(len=256)          :: fname
    character(len=256)          :: fname_gridboxarea
    character(32)               :: secname
    integer                     :: lt, year
    logical                     :: existfile
    integer,          dimension(2) :: ltimes         
    character(len=4), dimension(2) :: ar5_cyears       
    real,             dimension(2) :: ar5_ipcoef_years 
    logical                     :: first=.true.
    real                        :: co2_rcp_target, co2_scale 
    
    ! --- begin -----------------------------------

    ! initialise target array
    d3data = 0.0 

    ! read in gridbox-area; once per CPU
    ! For CO2 the area field is read from the CO2 LU file 
    !if( .not. area_found_05 ) then
    !   fname_gridboxarea = trim(emis_input_dir_ar5)//'/'//trim(ar5_filestr_gridboxarea)
    !   call emission_ReadGridboxArea(fname_gridboxarea, 'gridbox_area', gridbox_area_05, &
    !         & nlon720, nlat360, status )
    !   IF_NOTOK_RETURN(status=1)
    !   area_found_05=.true.
    !endif

    ! deal with out-of-bounds requested years
    year = valid_year( iyear, ar5_avail, 'AR5', first)
    first=.false.

    secname = sectors_def(sector)%name

    if ( iyear > year ) then 
      ! ------------------------
      ! data for the year ar5_avail(2)=2005 will be read from file 
      ! and need to be scaled (index=1: earlier year; index=2: later year)
      ! ------------------------

      ! ----------------------------------------
      ! get the right times to interpolate and related coefficients
      ! (ar5_avail_yrs(ltimes))
      !
      ! --> resulting scale factor  will be a linear interpolation between neighbouring values
      !
      ! ----------------------------------------
      allocate( ltimeind( ar5_nr_avail_yrs ) )
      ltimeind = .false.

      where( ar5_avail_yrs < iyear ) ltimeind = .true.

      ! times(1): index representing time instance earlier than current year
      ! times(2):     -"-          -"-               later than current year
      ltimes(2) = count( ltimeind ) + 1
      ltimes(1) = ltimes(2) - 1

      ! check a match with available years 
      ! (in order to use only value instead of two)
      if( ar5_avail_yrs(ltimes(2)) == iyear ) &
           ltimes(1) = ltimes(2)
       deallocate( ltimeind )

      ! ar5_cyears will contain strings with the years
      write(ar5_cyears(1),'(I4.4)') ar5_avail_yrs(ltimes(1))
      write(ar5_cyears(2),'(I4.4)') ar5_avail_yrs(ltimes(2))

      ! ar5_ipcoef_years will contain interpolation coefficients
      ! default: factors 1.0/0.0
      ar5_ipcoef_years(1) = 1.0
      ar5_ipcoef_years(2) = 0.0
      if( ltimes(2) /= ltimes(1) ) then 
         ar5_ipcoef_years(1) = (ar5_avail_yrs(ltimes(2)) - iyear) / &
                          real( ar5_avail_yrs(ltimes(2)) - ar5_avail_yrs(ltimes(1)) )
         ar5_ipcoef_years(2) = 1.0 - ar5_ipcoef_years(1)
      end if

      select case (trim (secname) )
        case ( 'emiss_ff' )
          co2_ref=co2_ff_ref
          select case (trim(filestr_rcpiden) )
            case ('RCP26')
              co2_rcp(:)=co2ff_rcp26(:)
            case ('RCP45')
              co2_rcp(:)=co2ff_rcp45(:)
            case ('RCP60')
              co2_rcp(:)=co2ff_rcp60(:)
            case ('RCP85')
              co2_rcp(:)=co2ff_rcp85(:)
            case default
              write(gol, '("ERROR: no RCP scenario specified for CO2 emissions")') ; call goErr
          end select
        case ( 'emiss_lu')
          co2_ref=co2_lu_ref
          select case (trim(filestr_rcpiden) )
            case ('RCP26')
              co2_rcp(:)=co2lu_rcp26(:)
            case ('RCP45')
              co2_rcp(:)=co2lu_rcp45(:)
            case ('RCP60')
              co2_rcp(:)=co2lu_rcp60(:)
            case ('RCP85')
              co2_rcp(:)=co2lu_rcp85(:)
          end select
      end select
    
      co2_rcp_target=co2_rcp(ltimes(1))*ar5_ipcoef_years(1)+co2_rcp(ltimes(2))*ar5_ipcoef_years(2)
      co2_scale=co2_rcp_target/co2_ref
    else
      ! no scaling for years <= 2005
      co2_scale=1.
    endif

    ! ------------------------
    ! read CO2 emission file
    ! ------------------------ 
    select case ( trim (secname) )
      case ( 'emiss_ff' )
         fname = trim(emis_input_dir_ar5)             //'/'//   &
         'CMIP5_gridcar_CO2_emissions_fossil_fuel_Andres_1751-2007_monthly_SC_mask11.nc'
      case ( 'emiss_lu' )
         fname = trim(emis_input_dir_ar5)             //'/'//   &
         'carbon_emissions_landuse_20person.nc'
      case default
        write(gol, '("ERROR: emission sector ",a,"not available for CO2")') &
            trim(secname); call goErr
       status=1; return
    end select

    ! test existence of file
    inquire( file=trim(fname), exist=existfile)
    if( .not. existfile ) then
       write (gol,'("ERROR: file `",a,"` not found ")') trim(fname); call goErr
       status=1; return
    end if

    select case ( trim (secname) )
      case ( 'emiss_ff' )
        d3data(:,:,1) = d3data(:,:,1) + co2_scale * &
                                        emission_ar5_ReadCO2FF( fname, year, imonth, status )
      case ( 'emiss_lu' )
        d3data(:,:,1) = d3data(:,:,1) + co2_scale * &
                                        emission_ar5_ReadCO2LU( fname, year, status )
      case default
        write(gol, '("ERROR: emission sector ",a,"not available for CO2")') &
            trim(secname); call goErr
       status=1; return
    end select
    IF_NOTOK_RETURN(status=1)


  end subroutine emission_ar5_ReadCO2
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:    EMISSION_AR5_READCO2FF
  !
  ! !DESCRIPTION: Read monthly AR5 fossil-fuel CO2 emissions on a 1x1 grid
  !\\
  !\\
  ! !INTERFACE:
  !
  function emission_ar5_ReadCO2FF( fname, year, imonth, status )
    !
    ! !RETURN VALUE:
    !
    real                  :: emission_ar5_ReadCO2FF(360,180)
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*),                intent(in) :: fname
    integer,                         intent(in) :: year, imonth
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,                        intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !    20 May 2014 - T. van Noije
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer :: fid, varid
    real    :: emis_in(360, 180), area(360,180)

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

    ! initialise
    emission_ar5_ReadCO2FF = 0.0

    CALL MDF_Open( TRIM(fname), MDF_NETCDF, MDF_READ, fid, status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Inq_VarID( fid, 'FF', varid, status )
    IF_ERROR_RETURN(status=1)

    CALL MDF_Get_Var( fid, varid, emis_in, status, start=(/1,1,12*(year-1751)+imonth/) )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Inq_VarID( fid, 'AREA', varid, status )
    IF_ERROR_RETURN(status=1)

    CALL MDF_Get_Var( fid, varid, area, status, start=(/1,1/) )
    IF_NOTOK_RETURN(status=1)
    ! to speed up reading of area could be done only once 

    ! convert from g(C)/m^2/s  to  kg(CO2)/gridbox/s
    emission_ar5_ReadCO2FF(:,:) = emis_in(:,:) * area(:,:) * 1.e-3 * xmco2/xmc

    IF_NOTOK_RETURN(status=1)

    CALL MDF_Close( fid, status )
    IF_NOTOK_RETURN(status=1)
 
    status = 0
    return

  end function emission_ar5_ReadCO2FF
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:    EMISSION_AR5_READCO2LU
  !
  ! !DESCRIPTION: Read yearly AR5 land-use CO2 emissions on a 0.5x0.5 grid 
  !               and convert to a 1x1 grid
  !\\
  !\\
  ! !INTERFACE:
  !
  function emission_ar5_ReadCO2LU( fname, year, status )
    !
    ! !RETURN VALUE:
    !
    real                  :: emission_ar5_ReadCO2LU(nlon360,nlat180)
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*),                intent(in) :: fname
    integer,                         intent(in) :: year
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,                        intent(out) :: status
    !
    ! !REVISION HISTORY:
    !    20 May 2014 - T. van Noije
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer :: fid, varid
    real    :: emis_in(nlon720, nlat360), area(nlon720, nlat360)

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

    ! initialise
    emission_ar5_ReadCO2LU = 0.0

    CALL MDF_Open( TRIM(fname), MDF_NETCDF, MDF_READ, fid, status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Inq_VarID( fid, 'carbon_emission', varid, status )
    IF_ERROR_RETURN(status=1)

    CALL MDF_Get_Var( fid, varid, emis_in, status, start=(/1,1,year-1850+1/) )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Inq_VarID( fid, 'area', varid, status )
    IF_ERROR_RETURN(status=1)

    CALL MDF_Get_Var( fid, varid, area, status, start=(/1,1/) )
    IF_NOTOK_RETURN(status=1)
    ! to speed up reading of area could be done only once

    ! convert from g(C)/m^2/s  to  kg(CO2)/gridbox/s
    !emis_in(:,:) = emis_in(:,:) * gridbox_area_05(:,:) * 1.e-3 * xmco2/xmc
    emis_in(:,:) = emis_in(:,:) * area(:,:) * 1.e-3 * xmco2/xmc

    ! now coarsen to nlon360,nlat180
    emission_ar5_ReadCO2LU = emission_coarsen_to_1x1( emis_in(:,:), nlon720, nlat360,.true., status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Close( fid, status )
    IF_NOTOK_RETURN(status=1)

    status = 0
    return

  end function emission_ar5_ReadCO2LU
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:   EMISSION_READGRIDBOXAREA
  !
  ! !DESCRIPTION: 
  !       reading gridbox surface areas for 0.5 x 0.5 Edgar 4
  !       needed to scale the emissions from mass/m^2 to mass/grid
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine emission_ReadGridboxArea(fname, recname, gridbox_area, dim_nlon, dim_nlat, status )
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*), intent(in)    ::  fname
    character(len=*), intent(in)    ::  recname
    integer,          intent(in)    ::  dim_nlon
    integer,          intent(in)    ::  dim_nlat
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,                             intent(out)  ::  status
    real, dimension(dim_nlon, dim_nlat), intent(out)  ::  gridbox_area
    !
    ! !REVISION HISTORY: 
    !
    !    1 Oct 2010 - Achim Strunk - v0
    !    1 Dec 2011 - Narcisa Banda - generalized it for any gridbox area size
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer :: fid, varid

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

    CALL MDF_Open( TRIM(fname), MDF_NETCDF, MDF_READ, fid, status )
    IF_NOTOK_RETURN(status=1)
    
    CALL MDF_Inq_VarID( fid, TRIM(recname), varid, status )
    IF_ERROR_RETURN(status=1)

    CALL MDF_Get_Var( fid, varid, gridbox_area, status )
    IF_NOTOK_RETURN(status=1)

    CALL MDF_Close( fid, status )
    IF_NOTOK_RETURN(status=1)
    
    status = 0

  end subroutine emission_ReadGridboxArea
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_AR5_REGRID_AIRCRAFT
  !
  ! !DESCRIPTION: Vertical regridding of the AR5 aircraft data. 
  !               The vertical levels of the input data are hard coded. 
  !               (Taken from GFED_daily_AR5 (VH) and left as is)
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine emission_ar5_regrid_aircraft(region, i0, j0, field_in, field_out, status )
    !
    ! !USES:
    !
    use meteodata,       only : gph_dat
    use tm5_distgrid,    only : dgrid, get_distgrid
    use dims,            only : lm
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)                              :: status
    !
    ! !INPUT PARAMETERS:
    !
    integer,                        intent(in) :: region, i0, j0
    real, dimension(i0:, j0:, 1:),  intent(in) :: field_in
    real, dimension(i0:, j0:, 1:), intent(out) :: field_out
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - Taken from GFED_daily_AR5 (VH) and left as is
    !    3 Dec 2012 - Ph. Le Sager - modified for lat-lon mpi decomposition
    !                              - work with zoom regions
    !                              - mass conservation per column     
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer, parameter :: lm_in=25
    
    real, dimension(:,:,:), pointer        :: gph       ! geopotential height (m)
    integer                                :: i,j,l
    real, dimension(lm_in)                 :: height_in_min, height_in_max
    real, allocatable                      :: dz(:), height(:)
    real                                   :: height_min,height_max
    real                                   :: height_out_min,height_out_max
    real, dimension(lm_in), parameter      :: height_in=(/ &  ! Height in meter
         305.,  915., 1525., 2135., 2745., 3355., 3965., 4575., 5185., 5795., &
         6405., 7015., 7625., 8235., 8845., 9455.,10065.,10675.,11285., &
         11895.,12505.,13115.,13725.,14335.,14945./)
    real    :: dz_in(25)
    integer :: l_in, i1, i2, j1, j2, lmr
    real    :: total_in, total_out, total_ratio
    
    ! --- begin --------------------------------------

    call golabel()

    gph => gph_dat(region)%data
    
    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
    lmr = lm(region)

    allocate(dz(lmr), height(lmr+1))
    
    ! sanity check    
    if (okdebug) then
       if (i1/=i0 .or. j1/=j0) then
          status = 1 ; TRACEBACK
          return
       end if
       if (lm_in /= ubound(field_in,3) ) then
          write(gol,*)'wrong vertical input resolution'; call goErr
          status = 1 
          TRACEBACK; return
       endif
       if ((ubound(field_out,3)+1) /= ubound(gph,3)) then
          write(gol,*)'wrong vertical output resolution'; call goErr
          status = 1 
          TRACEBACK; return
       endif
    end if
    
    ! locally flat atmosphere assumed
    ! area linear in i,j 


    ! height above sea level
    height_in_min(1)=0.
    do l_in = 2,lm_in
       height_in_min(l_in)=(height_in(l_in-1)+height_in(l_in))/2.
    enddo

    height_in_max(lm_in)=15555.
    do l_in = 1,lm_in-1
       height_in_max(l_in)=(height_in(l_in)+height_in(l_in+1))/2.
    enddo

    ! init
    field_out = 0.0

    ! regrid
    do i=i1, i2
    do j=j1, j2
          
          total_in = 0.0
       
          ! calculate total input emissions
          do l_in=1, lm_in
             dz_in(l_in) = height_in_max(l_in)-height_in_min(l_in)
             total_in =total_in + field_in(i,j,l_in)*dz_in(l_in)
          enddo

          ! zero based height:
          height(1) = 0.0 
          do l=1, lmr
             dz(l) = gph(i,j,l+1) - gph(i,j,l)
             height(l+1) = height(l) + dz(l)
          enddo
          
          do l=1,lmr-1
             height_out_min=height(l)
             height_out_max=height(l+1)
             !   write(*,*)'DO AR5- regrid - C2',i,j,l,height_out_min,height_out_max
             do l_in=1,lm_in
                if (height_out_max .le. height_in_min(l_in)) exit
                if (height_out_min .lt. height_in_max(l_in)) then
                   height_max=min(height_out_max,height_in_max(l_in))
                   height_min=max(height_out_min,height_in_min(l_in))
                   ! helpfield as field_in is ordered from high to low
                   !                 field_out(i,j,l) = field_out(i,j,l) + helpfield2(i,j,lm_in-l_in+1)* &
                   !                              (height_max-height_min)/(height_in_max(l_in)-height_in_min(l_in)) 
                   ! helpfield as field_in is ordered from low to high
                   !   write(*,*)'DO AR5- regrid - C',i,j,l,l_in,height_max-height_min
                   field_out(i,j,l) = field_out(i,j,l) + field_in(i,j,l_in)* &
                        (height_max-height_min)  ! kg/m -> kg / gridbox
                endif
             enddo
          enddo

          ! conserve total exactly: not possible because units are in kg/m...
          total_out   = sum(field_out(i,j,:))
          if (total_out /= 0) then
             total_ratio = total_in/total_out
             field_out(i,j,:) = field_out(i,j,:)*total_ratio
          end if

    enddo
    enddo

    deallocate(dz, height)

    call golabel()

    ! ok
    status = 0

  end subroutine emission_ar5_regrid_aircraft
  !EOC


END MODULE EMISSION_READ