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
!               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_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_HDF4, 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 :: numb_sectors, sectors_def
  public :: numb_providers, providers_def
  public :: sector_name_len
  public :: emission_ar5_readCO2
  public :: ar5_dim_3ddata

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

  ! ------------------------------
  ! global characteristics
  ! ------------------------------
  ! 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 AR5 (0.5deg)
  integer, parameter           :: nlon720         = 720    ! number of longitudes for AR5 data (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    =  2     ! number of sectors (All providers!)
  integer, parameter           :: numb_providers  = 1      ! AR5
  integer, parameter           :: ar5_dim_3ddata  = 25     ! number of layers for aircraft data

  ! full list of providers
  character(8), dimension(numb_providers), parameter :: all_providers = &
       & (/ '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

  ! List of providers effectively used
  character(8), 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  

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

    call ReadRc( rcF, 'input.emis.AR5.RCP', filestr_rcpiden, status, default='RCP26' )
    IF_ERROR_RETURN(status=1)

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

    ! Others gases
    mask = (/ 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     ') ! Fossil Fuel
    sectors_def( 2) = sector_type('emiss_lu  ', 'anthropogenic ', .false., 'nearsurface ', 'AR5     ') ! Land Use (assumed near surface for the moment, but that is open for discussion)

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

       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
    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.
  !\\
  !\\
  ! !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
       print*,'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_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

END MODULE EMISSION_READ