#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.: /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