Few lines of code

Python scripts

• Library imports

  import xarray as xr
  import matplotlib.pyplot as plt
  import pandas as pd
  import numpy as np
  

• Opening of datasets

  variables2 = ['vmrn2o', 'Burden', 'lossTot',  'area', 'AIRMASS', 'emin2o']
  
  nb_months = (1859-1850+1)*12
  exp = 'piC.map.ipsl'
  year=1850
  df1 = xr.open_dataset('/ccc/store/cont003/gen2201/laurentk/IGCM_OUT/LMDZORINCA/PROD/GES/'+str(exp)+'/CHM/Output/MO/'+str(exp)+'_'+str(year)+'0101_'+str(year)+'1231_1M_inca_ges.nc', decode_times=False)
  a_df = df1[variables2]
  for year in range(1851,1860):
      df1 = xr.open_dataset('/ccc/store/cont003/gen2201/laurentk/IGCM_OUT/LMDZORINCA/PROD/GES/'+str(exp)+'/CHM/Output/MO/'+str(exp)+'_'+str(year)+'0101_'+str(year)+'1231_1M_inca_ges.nc', decode_times=False)
      df1 = df1[variables2]
      a_df = xr.merge([a_df,df1])
      del df1
  
  alt = -7*np.log(a_df['presnivs']*1e-2/1054)
  

• Diagnostics

  a_df['airvmr'] = a_df['AIRMASS']*a_df['vmrn2o']
  a_df['bur'] = a_df['area']*a_df['Burden']/1e12
  a_df['lot'] = a_df['area']*a_df['lossTot']*3600*24*365/1e12
  a_df['emi'] = a_df['emin2o']*a_df['area']*3600*365*24/1e9*28/44
  a_df['lifetime'] = a_df['bur'].sum(['lat','lon','presnivs']).values/(a_df['lot'].sum(['lat','lon','presnivs']).values)
  a_df = a_df[['vmrn2o', 'lifetime', 'lot', 'bur','airvmr','AIRMASS', 'emi']]
  
  aa_vmrsurf = a_df['airvmr'].isel(presnivs=0).sum(['lat','lon'])/a_df['AIRMASS'].isel(presnivs=0).sum(['lat','lon'])
  aa_vmrcol = a_df['airvmr'].sum(['lat','lon','presnivs']).values/a_df['AIRMASS'].sum(['lat','lon','presnivs']).values
  

• Rolling average

  a_vmrsurf = xr.DataArray(aa_vmrsurf*1e9, coords=[pd.date_range('1850-01-01',freq='M',periods=nb_months)], dims='time')
  a_vmrsurf = a_vmrsurf.rolling(time=12, center=True).mean()
  a_vmrcol = xr.DataArray(aa_vmrcol*1e9, coords=[pd.date_range('1850-01-01',freq='M',periods=nb_months)], dims='time')
  a_vmrcol = a_vmrcol.rolling(time=12, center=True).mean()
  a_lft = xr.DataArray(a_df['lifetime'].values, coords=[pd.date_range('1850-01-01',freq='M',periods=nb_months)], dims='time')
  a_lft = a_lft.rolling(time=16, center=True).mean()
  a_lot = xr.DataArray(a_df['lot'].sum(['lat','lon','presnivs']).values, coords=[pd.date_range('1850-01-01',freq='M',periods=nb_months)], dims='time')
  a_lot = a_lot.rolling(time=16, center=True).mean()
  a_bur = xr.DataArray(a_df['bur'].sum(['lat','lon','presnivs']).values, coords=[pd.date_range('1850-01-01',freq='M',periods=nb_months)], dims='time')
  a_bur = a_bur.rolling(time=12, center=True).mean()
  

• Images

  plt.figure()
  plt.plot(pd.date_range('1850-01-01',freq='M',periods=nb_months), a_lft, color='tab:green', lw=2,label='emissions IPSL')
  plt.title('Annual N2O lifetime (year)')
  plt.xlabel('Year')
  plt.ylabel('N2O lifetime (yr)')
  plt.legend()
  plt.grid()
  plt.tight_layout()
  plt.savefig('name_of_image.png', transparent=False)
  

• Inputs files

  nat = xr.open_dataset('/ccc/store/cont003/gen2201/laurentk/sflx_outputs_N2O/sflx_lmdz_NAT_1850_sep_phy.nc', decode_times=False)
  nat = nat[['flx_n2o_land','flx_n2o_ocean']]
  ant = xr.open_dataset('/ccc/store/cont003/gen2201/laurentk/sflx_outputs_N2O/sflx_lmdz_ANT_1850_dyn.nc', decode_times=False)
  bbg = xr.open_dataset('/ccc/store/cont003/gen2201/laurentk/sflx_outputs_N2O/sflx_lmdz_BBG_1850_dyn.nc', decode_times=False)
  aire = xr.open_dataset('/ccc/store/cont003/gen2201/laurentk/sflx_outputs_N2O/area_phy.nc', decode_times=False)
  airephy = xr.open_dataset('/ccc/store/cont003/gen2201/laurentk/sflx_outputs_N2O/aire_phy_restartlmdz.nc', decode_times=False)
  airephy = airephy.swap_dims({'points_physiques':'vector'})
  
  nat['flx_nat']=(nat['flx_n2o_land']+nat['flx_n2o_ocean'])*airephy['AIRE']*3600*365*24/1e9*28/44
  ant['flx_ant']=ant['flx_n2o']*aire['area']*3600*365*24/1e9*28/44
  bbg['flx_bbg']=bbg['flx_n2o']*aire['area']*3600*365*24/1e9*28/44
  print('nat, ant, bbg =', nat['flx_nat'].sum(['vector']), ant['flx_ant'].sum(['lat','lon']), bbg['flx_bbg'].sum(['lat','lon']), ' TgN/yr')
  
  

• Full script for pisces

  import xarray as xr
  import matplotlib.pyplot as plt
  import pandas as pd
  ipsl = xr.open_dataset('N2Oflx_hist.pisces_18360101_18831231_1M_gastrc.nc')
  cnrm = xr.open_dataset('N2Oflx_HOMOImon_CNRM-ESM2-1_historical_r1i1p1f2_gn_185001-201412.nc')
  areacello = xr.open_dataset('areacello.nc')
  flx_ipsl=ipsl['N2Oflx']*ipsl['cell_area']
  flx_cnrm = cnrm['N2Oflx']*areacello['areacello']
  ipsl = xr.DataArray(flx_ipsl.sum(['y','x'])*3600*365*24*14/1e12*(-1)*44/28, coords=[pd.date_range('1836-01-01',freq='M',periods=576)], dims='time_counter')
  cnrm = xr.DataArray(flx_cnrm.sum(['y','x'])*3600*365*24*14/1e12*(-1)*44/28, coords=[pd.date_range('1850-01-01',freq='M',periods=1980)], dims='time')
  mipsl36 = ipsl.rolling(time_counter=36, center=True).mean()
  mcnrm36 = cnrm.rolling(time=36, center=True).mean()
  plt.figure(figsize=(12,8))
  plt.plot(pd.date_range('1850-01-01',freq='M',periods=1980), mcnrm36,lw=2, color='tab:blue', label='cnrm')
  plt.plot(pd.date_range('1836-01-01',freq='M',periods=576), mipsl36,lw=2, color='orange', label='ipsl')
  plt.legend()
  plt.grid()
  plt.xlabel('Year')
  plt.ylabel('N2Oflx (TgN2O/yr)')
  plt.title('Evolution of N2Oflx from IPSL and CNRM model')
  plt.savefig('ipsl_cnrm.png')
  plt.show()
  

Ferret scripts

• Inputs files

  use sflx_lmdz_NAT_esm_dyn.nc
  let flx_tot=(flx_n2o_ocean[d=1]+flx_n2o_land[d=1])*area
  list flx_tot[i=@sum,j=@sum,l=@ave]*3600*24*365/1e9*28/44
  use sflx_lmdz_bbg 
  use sflx_lmdz_ant
  let ant=area*flx_n2o[d=1,l=@ave]*3600*24*365/1e9*28/44
  list ant[i=@sum,j=@sum]
  let bbg=area*flx_n2o[d=2,l=@ave]*3600*24*365/1e9*28/44
  list bbg[i=@sum,j=@sum] 
  
  use sflx_lmdz_MAP_reduce4.17.nc 
  use aire_phy_restartlmdz.nc 
  let flx=flx_n2o_map[d=1]*aire*3600*365*24/1e9
  list flx[i=@sum,l=@ave]
  list (flx[i=@sum,l=1]*31+ flx[i=@sum,l=2]*28+ flx[i=@sum,l=3]*31+ flx[i=@sum,l=4]*30+ flx[i=@sum,l=5]*31+ flx[i=@sum,l=6]*30+ flx[i=@sum,l=7]*31+ flx[i=@sum,l=8]*31+ flx[i=@sum,l=9]*30+ flx[i=@sum,l=10]*31+ flx[i=@sum,l=11]*30+ flx[i=@sum,l=12]*31)/365
  
  

• Diagnostics

  use inca_ges.nc
  let lot=losstot*area
  let emi=emin2o*area 
  list emi[i=@sum,j=@sum,l=@ave]*3600*24*365/1e9*28/44
  
  let vmr=vmrn2o*airmass
  let surf=vmr[k=79,i=@sum,j=@sum]/airmass[i=@sum,j=@sum,k=79]*1e9
  plot surf[d=1],surf[d=2], surf[d=3],surf[d=4],surf[d=5],surf[d=6]
  let col=vmr[k=@sum,i=@sum,j=@sum]/airmass[i=@sum,j=@sum,k=@sum]*1e9
  plot col[d=1]
  let prof=vmr[l=12,i=@sum,j=@sum]/airmass[i=@sum,j=@sum,l=12]*1e9
  plot/vlog prof[d=1]
  

• Images

  shade/palette=rainbow/levels="(-inf)(0 10 0.5)(inf)" n2oflx[d=1,l=@ave,k=1]*area[d=1], nav_lon, nav_lat