From 5c86c952adf9f6effcfb2e3b70d5413b570d1b88 Mon Sep 17 00:00:00 2001 From: Tomas Torsvik Date: Tue, 27 Aug 2024 18:25:53 +0200 Subject: [PATCH 1/4] Add python script for generating a TIPMIP CO2 emission file. --- .../cam/ggas/create_TIPMIP_emission_file.py | 83 +++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py diff --git a/PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py b/PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py new file mode 100644 index 0000000..32d9583 --- /dev/null +++ b/PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py @@ -0,0 +1,83 @@ +import numpy as np +import sys +import xarray as xr +from scipy.io import netcdf +import pandas as pd + + +# load grgid area from fx +dataset = xr.open_dataset('/mnt/bgcdata-ns2980k/ffr043/TipESM/data/NorESM2-LM/1pctCO2/fx/areacella_fx_NorESM2-LM_1pctCO2_r1i1p1f1_gn.nc') +areavar = dataset['areacella'] + +# load example emission file +inpath = '/mnt/bgcdata-ns2980k/ffr043/TipESM/data' +dataset = xr.open_mfdataset(inpath+'/emissions-cmip6_CO2_anthro_surface_175001-201512_fv_1.9x2.5_c20181011.nc', decode_times=False, format="NETCDF3_64BIT") +co2var = dataset['CO2_flux'] + +# calculate emissions based on TCRE (Arora et al., 2020) +TCRE = 1.32 #°C EgC-1 +E = 1000 * 0.02 / TCRE #GtC yr-1 +Edata = np.ones((250*12, 1)) * E + +# distribute emissions over time and space +# translated to python from matlab based on script by J. Schwinger +nyears = 250 +nmonth = nyears * 12 +nlon = 144 +nlat = 96 +start_year=1 +time_bnds = np.zeros((2, nmonth)) +time = np.zeros(nmonth) +date = np.zeros(nmonth, dtype=np.int32) +co2flx = np.zeros((nmonth, nlat, nlon)) +dayim = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] +days = np.cumsum([0] + dayim[:-1]) + 1 +daye = np.cumsum(dayim) + +# (assumed) midpoint of each month in the format MMDD +date_mint = [116, 215, 316, 416, 516, 616, 716, 816, 916, 1016, 1116, 1216] +for iy in range(1, nyears + 1): + for im in range(1, 13): + idx = (iy - 1) * 12 + im - 1 + time_bnds[0, idx] = 365 * (iy - 1) + days[im - 1] - 1 + time_bnds[1, idx] = 365 * (iy - 1) + daye[im - 1] + time[idx] = np.sum(time_bnds[:, idx]) / 2.0 + date[idx] = (1850 + (iy - 1)) * 10000 + date_mint[im - 1] + + if start_year <= iy < start_year + 250: + # unit correction - time + dt = (time_bnds[1, idx] - time_bnds[0, idx]) * 86400.0 + # Spatially, emissions are distributed evenly over the sphere and months + totarea = areavar.sum().values + # unit correction - Gt C to kg CO2 + co2flx[idx, :, :] = Edata[idx - (start_year - 1) *12] * 1e12 * 3.664 / dt / totarea / 12 + else: + co2flx[idx, :, :] = 0.0 + + +# write emissions to example dataset +dataset = dataset.isel(time=slice(0,3000)) +dataset['CO2_flux'].values = co2flx + +# assign attributes +dataset = dataset.assign_attrs({'data_title':'Annual Anthropogenic Emissions of CO2 based on TCRE prepared for TIPMIP'}) +dataset = dataset.assign_attrs({'data_creator':'F. Froeb (friederike.frob@uib.no)'}) +dataset = dataset.assign_attrs({'creation_date':'2024-08-03'}) + +#set encoding for netcdf file +encoding = { +'time':{'_FillValue': None}, +'time_bnds':{'_FillValue': None}, +'lon':{'zlib': True, 'shuffle': False, 'complevel': 1, 'fletcher32': False, 'contiguous': False, + 'dtype': 'float64', '_FillValue':None}, +'lat':{'zlib': True, 'shuffle': False, 'complevel': 1, 'fletcher32': False, 'contiguous': False, + 'dtype': 'float64', '_FillValue':None}, +'CO2_flux':{'zlib': True, 'shuffle': True, 'complevel': 9, 'fletcher32': False, 'contiguous': False, + 'dtype': 'float32', 'missing_value': 1e+20, '_FillValue': 1e+20} +} + +#write netcdf +dataset.to_netcdf('/mnt/bgcdata-ns2980k/ffr043/TipESM/data/emissions-ESM-tipmip_CO2_anthro_surface_185001-209912_fv_1.9x2.5_c20240803.nc', mode="w", format="NETCDF3_64BIT", encoding=encoding, unlimited_dims='time') + + + From aedbe7e81141350d7538d7bddcbcd66bb3f12bd0 Mon Sep 17 00:00:00 2001 From: Tomas Torsvik Date: Wed, 28 Aug 2024 11:02:01 +0200 Subject: [PATCH 2/4] Update python script for generating a TIPMIP CO2 emission file. --- .../cam/ggas/create_TIPMIP_emission_file.py | 27 +++++++++---------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py b/PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py index 32d9583..a0d3a03 100644 --- a/PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py +++ b/PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py @@ -5,27 +5,23 @@ import pandas as pd -# load grgid area from fx +# load grid area from fx dataset = xr.open_dataset('/mnt/bgcdata-ns2980k/ffr043/TipESM/data/NorESM2-LM/1pctCO2/fx/areacella_fx_NorESM2-LM_1pctCO2_r1i1p1f1_gn.nc') areavar = dataset['areacella'] +totarea = areavar.sum().values # load example emission file inpath = '/mnt/bgcdata-ns2980k/ffr043/TipESM/data' dataset = xr.open_mfdataset(inpath+'/emissions-cmip6_CO2_anthro_surface_175001-201512_fv_1.9x2.5_c20181011.nc', decode_times=False, format="NETCDF3_64BIT") co2var = dataset['CO2_flux'] -# calculate emissions based on TCRE (Arora et al., 2020) -TCRE = 1.32 #°C EgC-1 -E = 1000 * 0.02 / TCRE #GtC yr-1 -Edata = np.ones((250*12, 1)) * E - # distribute emissions over time and space # translated to python from matlab based on script by J. Schwinger nyears = 250 nmonth = nyears * 12 nlon = 144 nlat = 96 -start_year=1 +start_year = 1 time_bnds = np.zeros((2, nmonth)) time = np.zeros(nmonth) date = np.zeros(nmonth, dtype=np.int32) @@ -34,6 +30,11 @@ days = np.cumsum([0] + dayim[:-1]) + 1 daye = np.cumsum(dayim) +# calculate emissions based on TCRE (Arora et al., 2020) +TCRE = 1.32 #°C EgC-1 +E = 1000 * 0.02 / TCRE #GtC yr-1 +Edata = np.ones((nmonth, 1)) * E + # (assumed) midpoint of each month in the format MMDD date_mint = [116, 215, 316, 416, 516, 616, 716, 816, 916, 1016, 1116, 1216] for iy in range(1, nyears + 1): @@ -43,20 +44,18 @@ time_bnds[1, idx] = 365 * (iy - 1) + daye[im - 1] time[idx] = np.sum(time_bnds[:, idx]) / 2.0 date[idx] = (1850 + (iy - 1)) * 10000 + date_mint[im - 1] - - if start_year <= iy < start_year + 250: + + # Spatially, emissions are distributed evenly over the sphere and months + if start_year <= iy < start_year + nyears: # unit correction - time dt = (time_bnds[1, idx] - time_bnds[0, idx]) * 86400.0 - # Spatially, emissions are distributed evenly over the sphere and months - totarea = areavar.sum().values # unit correction - Gt C to kg CO2 - co2flx[idx, :, :] = Edata[idx - (start_year - 1) *12] * 1e12 * 3.664 / dt / totarea / 12 + co2flx[idx, :, :] = Edata[idx] * 1e12 * 3.664 / dt / totarea / 12 else: co2flx[idx, :, :] = 0.0 - # write emissions to example dataset -dataset = dataset.isel(time=slice(0,3000)) +dataset = dataset.isel(time=slice(0, nmonth)) dataset['CO2_flux'].values = co2flx # assign attributes From 6dda4f294c7eb5037543843425765fb052d0ac71 Mon Sep 17 00:00:00 2001 From: Tomas Torsvik Date: Fri, 30 Aug 2024 15:23:13 +0200 Subject: [PATCH 3/4] Update TIPMIP esmission script --- .../cam/ggas/create_TIPMIP_emission_file.py | 20 +++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py b/PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py index a0d3a03..d34371d 100644 --- a/PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py +++ b/PreProc/inputdata/atm/cam/ggas/create_TIPMIP_emission_file.py @@ -33,10 +33,14 @@ # calculate emissions based on TCRE (Arora et al., 2020) TCRE = 1.32 #°C EgC-1 E = 1000 * 0.02 / TCRE #GtC yr-1 -Edata = np.ones((nmonth, 1)) * E +Edata = np.ones((nmonth, 1)) * E # (assumed) midpoint of each month in the format MMDD date_mint = [116, 215, 316, 416, 516, 616, 716, 816, 916, 1016, 1116, 1216] + +# Set emissions constant over each month or not +constant_em = True + for iy in range(1, nyears + 1): for im in range(1, 13): idx = (iy - 1) * 12 + im - 1 @@ -47,10 +51,14 @@ # Spatially, emissions are distributed evenly over the sphere and months if start_year <= iy < start_year + nyears: - # unit correction - time - dt = (time_bnds[1, idx] - time_bnds[0, idx]) * 86400.0 - # unit correction - Gt C to kg CO2 - co2flx[idx, :, :] = Edata[idx] * 1e12 * 3.664 / dt / totarea / 12 + if constant_em == True: + # unit correction - Gt C yr-1 to kg CO2 s-1 + co2flx[idx, :, :] = Edata[idx] * 1e12 * 3.664 / totarea / 86400.0 / 365.0 + else: + # unit correction sec month-1 + dt = (time_bnds[1, idx] - time_bnds[0, idx]) * 86400.0 + # unit correction - Gt C month-1 to kg CO2 s-1 + co2flx[idx, :, :] = Edata[idx] * 1e12 * 3.664 / dt / totarea / 12 else: co2flx[idx, :, :] = 0.0 @@ -76,7 +84,7 @@ } #write netcdf -dataset.to_netcdf('/mnt/bgcdata-ns2980k/ffr043/TipESM/data/emissions-ESM-tipmip_CO2_anthro_surface_185001-209912_fv_1.9x2.5_c20240803.nc', mode="w", format="NETCDF3_64BIT", encoding=encoding, unlimited_dims='time') +dataset.to_netcdf('/mnt/bgcdata-ns2980k/ffr043/TipESM/data/emissions-ESM-tipmip_CO2_anthro_surface_185001-209912_fv_1.9x2.5_c20240823_cE.nc', mode="w", format="NETCDF3_64BIT", encoding=encoding, unlimited_dims='time') From eec49f2d31a14cc547f87064f88dfbc8f902a402 Mon Sep 17 00:00:00 2001 From: Tomas Torsvik Date: Fri, 30 Aug 2024 15:39:03 +0200 Subject: [PATCH 4/4] Also update the README.md file --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index 57b5969..70cf157 100644 --- a/README.md +++ b/README.md @@ -5,6 +5,7 @@ Miscellaneous pre- and post-processing scripts for NorESM - mg_budget_plot.py: Script for plotting the the tendency terms that contribute to the total tendencies of the cloud microphysics scheme by Morrison and Gettelman in NorESM2/CAM6 - PostProc: Scripts for post processing model output - rcat_noresm.sh: Concatenate and compress NorESM model output +- PreProc: Scripts for preparing model input datasets - Utilities: Miscellaneous scripts and programs - buildCPRNC_NIRD.sh: File to build the cprnc tool needed by rcat_noresm.sh - test_rcat_noresm.sh: Unit tests for rcat_noresm.sh