text: ﻿
  # -*- coding: utf-8 -*-
  """
  Created on Tue Feb 18 16:09:30 2014
  
  @author: Greg
  """
  
  
  import numpy as np
  from scipy.interpolate import interp1d
  from scipy.signal import fftconvolve
  from scipy.integrate import cumtrapz
  
  # Radiative efficiencies of each gas, calculated from AR5 & AR5 SM
  co2_re, ch4_re, n2o_re, sf6_re = 1.756E-15, 1.277E-13 * 1.65, 3.845E-13, 2.010E-11
  
  # AR5 2013 IRF values
  a0, a1, a2, a3 = 0.2173, 0.224, 0.2824, 0.2763
  tau1, tau2, tau3 = 394.4, 36.54, 4.304
  
  def f0(t):
      return a0
  def f1(t):
      return a1*np.exp(-t/tau1)
  def f2(t):
      return a2*np.exp(-t/tau2)
  def f3(t):
      return a3*np.exp(-t/tau3)
  def CO2_AR5(t):
      return f0(t) + f1(t) + f2(t) + f3(t)
      
  #Methane response fuction
  CH4tau = 12.4
  def CH4_AR5(t):
      return np.exp(-t/CH4tau)
      
  #N2O response fuction
  N2Otau = 121
  def N2O_AR5(t):
      return np.exp(-t/CH4tau)
  
  #SF6 response fuction   
  SF6tau = 3200
  def SF6_AR5(t):
      return np.exp(-t/CH4tau)
  
  #Temperature response function to radiative forcing
  def AR5_GTP(t):
      c1, c2, d1, d2 = 0.631, 0.429, 8.4, 409.5
      """ The default response function for radiative forcing from AR5. Source is \
      Boucher (2008). ECR is 3.9K, which is on the high side.
      Convolve with radiative forcing to get temperature.
      """
      return c1/d1*np.exp(-t/d1) + c2/d2*np.exp(-t/d2)
  
  def Alt_GTP(t):
      c1, c2, d1, d2 = 0.43, 0.32, 2.57, 82.24
      """ The response function for radiative forcing. Taken from Olivie and Peters (2013),
      Table 4, using the CMIP5 data. This has a slightly lower climate response value than
      Boucher (2008), which is used in AR5.
  
      Convolve with radiative forcing to get temperature.
      """
      return c1/d1*np.exp(-t/d1) + c2/d2*np.exp(-t/d2)
  
  def Alt_low_GTP(t):
      c1, c2, d1, d2 = 0.43 / (1 + 0.29), 0.32 / (1 + 0.59), 2.57 * 1.46, 82.24 * 2.92
      #c1, c2, d1, d2 = 0.48 * (1 - 0.3), 0.20 * (1 - 0.52), 7.15 * 1.35, 105.55 * 1.38
      #c1, c2, d1, d2 = 0.48 * (1 - 0.3), 0.20 * (1 - 0.52), 7.15, 105.55
      #c1, c2, d1, d2 = 0.631 * 0.7, 0.429 * 0.7, 8.4, 409.5
      """ The response function for radiative forcing. Taken from Olivie and Peters (2013),
      Table 4, using the CMIP5 data. This has a lower climate response value than AR5.
      The uncertainty in Table 4 assumes lognormal distributions, which is why values less
      than the median are determined by dividing by (1 + uncertainty).
  
      Convolve with radiative forcing to get temperature.
      """
      return c1/d1*np.exp(-t/d1) + c2/d2*np.exp(-t/d2)
  
  def Alt_high_GTP(t):
      c1, c2, d1, d2 = 0.43 * 1.29, 0.32 * 1.59, 2.57 / (1 + 0.46), 82.24 / (1 + 1.92)
      #c1, c2, d1, d2 = 0.48 * 1.3, 0.20 * 1.52, 7.15 * (1 - 0.35), 105.55 * (1 - 0.38)
      #c1, c2, d1, d2 = 0.48 * 1.2, 0.20 * 1.52, 7.15, 105.55
      #c1, c2, d1, d2 = 0.631, 0.429 * 1.3, 8.4, 409.5    
      """ The response function for radiative forcing. Taken from Olivie and Peters (2013),
      Table 4, using the CMIP5 data. This has a higher climate response value than AR5.
      The uncertainty in Table 4 assumes lognormal distributions, which is why values less
      than the median are determined by dividing by (1 + uncertainty).
  
      Convolve with radiative forcing to get temperature.
      """
      return c1/d1*np.exp(-t/d1) + c2/d2*np.exp(-t/d2)
  
  def CO2_rf(emission, years, tstep=0.01, kind='linear'): 
      """Transforms an array of CO2 emissions into radiative forcing with user-
      defined time-step.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      """
  #emission is a series of emission numbers, years should match up with it
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years) 
      f = interp1d(years, emission, kind=kind)
      time = np.linspace(years[0], end, end/tstep + 1)    
      inter_emissions = f(time)
      atmos = np.resize(fftconvolve(CO2_AR5(time), inter_emissions), time.size) * tstep
      rf = atmos * co2_re
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
          
      return rf[fil]
  
  def CO2_rate(emission, years, tstep=0.01, kind='linear'): 
      """Transforms an array of CO2 emissions into radiative forcing with user-
      defined time-step.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      """
  #emission is a series of emission numbers, years should match up with it
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years) 
      f = interp1d(years, emission, kind=kind)
      time = np.linspace(years[0], end, end/tstep + 1)    
      inter_emissions = f(time)
      atmos = np.resize(fftconvolve(CO2_AR5(time), inter_emissions), time.size) * tstep
      rf = atmos * co2_re
      dx = np.gradient(time)
      rate = np.gradient(rf, dx)
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
          
      return rate[fil]
  
  def CO2_crf(emission, years, tstep=0.01, kind='linear'): 
      """Transforms an array of CO2 emissions into radiative forcing with user-
      defined time-step.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      """
  #emission is a series of emission numbers, years should match up with it
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years) 
      f = interp1d(years, emission, kind=kind)
      time = np.linspace(years[0], end, end/tstep + 1)    
      inter_emissions = f(time)
      atmos = np.resize(fftconvolve(CO2_AR5(time), inter_emissions), time.size) * tstep
      rf = atmos * co2_re
      crf = cumtrapz(rf, dx = tstep, initial = 0)
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
          
      return crf[fil]
  
  
  def CO2_temp(emission, years, tstep=0.01, kind='linear', source='AR5'): 
      """Transforms an array of CO2 emissions into temperature with user-defined
      time-step. Default temperature IRF is from AR5, use 'Alt_low' or 'Alt_high'
      for a sensitivity test.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      source: the source of parameters for the temperature IRF. default is AR5,
      'Alt', 'Alt_low', and 'Alt_high' are also options.
      """
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years)       
      f = interp1d(years, emission, kind=kind, bounds_error=False)
      time = np.linspace(min(years), end, end/tstep + 1) 
      inter_emissions = f(time)
      atmos = np.resize(fftconvolve(CO2_AR5(time), inter_emissions), time.size) * tstep
      rf = atmos * co2_re
      if source == 'AR5':
          temp = np.resize(fftconvolve(AR5_GTP(time), rf), time.size) * tstep
      elif source == 'Alt':
          temp = np.resize(fftconvolve(Alt_GTP(time), rf), time.size) * tstep
      elif source == 'Alt_low':
          temp = np.resize(fftconvolve(Alt_low_GTP(time), rf), time.size) * tstep
      elif source == 'Alt_high':
          temp = np.resize(fftconvolve(Alt_high_GTP(time), rf), time.size) * tstep
  
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
          
      return temp[fil]
  
  
  def ch42co2(t, alpha=0.51):
      """As methane decays some fraction is converted to CO2. This function is 
      from Boucher (2009). By default it converts 51%. The convolution of this 
      function with the methane emission profile gives the CO2 emission profile.
      
      t: time
      alpha: fraction of methane converted to CO2
      """
      ch4tau = 12.4
      return 1/ch4tau * alpha * np.exp(-t/ch4tau)
  
  
  
  def CH4_rf(emission, years, tstep=0.01, kind='linear',
               decay=True):
      """Transforms an array of methane emissions into radiative forcing with user-defined
      time-step.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      """
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years) 
      fch4 = interp1d(years, emission, kind=kind)
      time = np.linspace(years[0], end, end/tstep + 1)    
      ch4_inter_emissions = fch4(time)
      ch4_atmos = np.resize(fftconvolve(CH4_AR5(time), ch4_inter_emissions),
                            time.size) * tstep
      co2 = np.resize(fftconvolve(ch42co2(time), ch4_inter_emissions),
                      time.size) * tstep
      co2_atmos = np.resize(fftconvolve(CO2_AR5(time), co2),
                            time.size) * tstep
      
      if decay == True:
           rf = ch4_atmos * ch4_re + co2_atmos * co2_re
      else:
          rf = ch4_atmos * ch4_re
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
      
      return rf[fil]
      
  def CH4_rf_cc(emission, years, tstep=0.01, kind='linear',
               decay=True):
      """Transforms an array of methane emissions into radiative forcing with user-defined
      time-step, accounting for climate-carbon feedbacks.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      """
      gamma = (44.0/12.0) * 10**12
      
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years) 
      fch4 = interp1d(years, emission, kind=kind)
      time = np.linspace(years[0], end, end/tstep + 1)    
      ch4_inter_emissions = fch4(time)
      ch4_atmos = np.resize(fftconvolve(CH4_AR5(time), ch4_inter_emissions),
                            time.size) * tstep
      co2 = np.resize(fftconvolve(ch42co2(time), ch4_inter_emissions),
                      time.size) * tstep
      co2_atmos = np.resize(fftconvolve(CO2_AR5(time), co2),
                            time.size) * tstep
      cc_co2 = CH4_cc_tempforrf(emission, years) * gamma
      cc_co2_atmos = np.resize(fftconvolve(CO2_AR5(time), cc_co2),
                            time.size) * tstep
      
      if decay == True:
           rf = ch4_atmos * ch4_re + (co2_atmos +cc_co2_atmos) * co2_re
      else:
          rf = ch4_atmos * ch4_re + (cc_co2_atmos) * co2_re
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
      
      return rf[fil]
  
  def CH4_rate(emission, years, tstep=0.01, kind='linear'):
      """Transforms an array of methane emissions into radiative forcing with user-defined
      time-step, accounting for climate-carbon feedbacks.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      """
      gamma = (44.0/12.0) * 10**12
      
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years) 
      fch4 = interp1d(years, emission, kind=kind)
      time = np.linspace(years[0], end, end/tstep + 1)    
      ch4_inter_emissions = fch4(time)
      ch4_atmos = np.resize(fftconvolve(CH4_AR5(time), ch4_inter_emissions),
                            time.size) * tstep
      co2 = np.resize(fftconvolve(ch42co2(time), ch4_inter_emissions),
                      time.size) * tstep
      co2_atmos = np.resize(fftconvolve(CO2_AR5(time), co2),
                            time.size) * tstep
      cc_co2 = CH4_cc_tempforrf(emission, years) * gamma
      cc_co2_atmos = np.resize(fftconvolve(CO2_AR5(time), cc_co2),
                            time.size) * tstep
      
      rf = ch4_atmos * ch4_re + (co2_atmos +cc_co2_atmos) * co2_re
      dx = np.gradient(time)
      rate = np.gradient(rf, dx)
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
      
      return rate[fil]
  
  def CH4_crf(emission, years, tstep=0.01, kind='linear',
               decay=True):
      """Transforms an array of methane emissions into radiative forcing with user-defined
      time-step.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      """
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years) 
      fch4 = interp1d(years, emission, kind=kind)
      time = np.linspace(years[0], end, end/tstep + 1)    
      ch4_inter_emissions = fch4(time)
      ch4_atmos = np.resize(fftconvolve(CH4_AR5(time), ch4_inter_emissions),
                            time.size) * tstep
      co2 = np.resize(fftconvolve(ch42co2(time), ch4_inter_emissions),
                      time.size) * tstep
      co2_atmos = np.resize(fftconvolve(CO2_AR5(time), co2),
                            time.size) * tstep
      
      if decay == True:
           rf = ch4_atmos * ch4_re + co2_atmos * co2_re
      else:
          rf = ch4_atmos * ch4_re
      crf = cumtrapz(rf, dx = 1, initial = 0)
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
      
      return crf[fil]
  
  def CH4_crf_cc(emission, years, tstep=0.01, kind='linear',
               decay=True):
      """Transforms an array of methane emissions into radiative forcing with user-defined
      time-step.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      """
      gamma = (44.0/12.0) * 10**12
  
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years) 
      fch4 = interp1d(years, emission, kind=kind)
      time = np.linspace(years[0], end, end/tstep + 1)    
      ch4_inter_emissions = fch4(time)
      ch4_atmos = np.resize(fftconvolve(CH4_AR5(time), ch4_inter_emissions),
                            time.size) * tstep
      co2 = np.resize(fftconvolve(ch42co2(time), ch4_inter_emissions),
                      time.size) * tstep
      co2_atmos = np.resize(fftconvolve(CO2_AR5(time), co2),
                            time.size) * tstep
      cc_co2 = CH4_cc_tempforrf(emission, years) * gamma
      cc_co2_atmos = np.resize(fftconvolve(CO2_AR5(time), cc_co2),
                            time.size) * tstep
      
      if decay == True:
           rf = ch4_atmos * ch4_re + (co2_atmos +cc_co2_atmos) * co2_re
      else:
          rf = ch4_atmos * ch4_re + (cc_co2_atmos) * co2_re
      crf = cumtrapz(rf, dx = 1, initial = 0)
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
      
      return crf[fil]
  
  def CH4_temp(emission, years, tstep=0.01, kind='linear', source='AR5',
               decay=True): 
      """Transforms an array of methane emissions into temperature with user-defined
      time-step. Default temperature IRF is from AR5, use 'Alt_low' or 'Alt_high'
      for a sensitivity test.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      source: the source of parameters for the temperature IRF. default is AR5,
      'Alt', 'Alt_low', and 'Alt_high' are also options.
      decay: a boolean variable for if methane decay to CO2 should be included
      """
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years) 
      f = interp1d(years, emission, kind=kind)
      time = np.linspace(years[0], end, end/tstep + 1)    
      ch4_inter_emissions = f(time)
      ch4_atmos = np.resize(fftconvolve(CH4_AR5(time), ch4_inter_emissions),
                            time.size) * tstep
      co2 = np.resize(fftconvolve(ch42co2(time), ch4_inter_emissions),
                      time.size) * tstep
      co2_atmos = np.resize(fftconvolve(CO2_AR5(time), co2),
                            time.size) * tstep
      if decay == True:
           rf = ch4_atmos * ch4_re + co2_atmos * co2_re
      else:
          rf = ch4_atmos * ch4_re
      if source == 'AR5':
          temp = np.resize(fftconvolve(AR5_GTP(time), rf), time.size) * tstep
      elif source == 'Alt':
          temp = np.resize(fftconvolve(Alt_GTP(time), rf), time.size) * tstep
      elif source == 'Alt_low':
          temp = np.resize(fftconvolve(Alt_low_GTP(time), rf), time.size) * tstep
      elif source == 'Alt_high':
          temp = np.resize(fftconvolve(Alt_high_GTP(time), rf), time.size) * tstep
  
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
      
      return temp[fil]
  
  def CH4_cc_tempforrf(emission, years, tstep=0.01, kind='linear', source='AR5',
               decay=True): 
      """Transforms an array of methane emissions into temperature with user-defined
      time-step. Default temperature IRF is from AR5, use 'Alt_low' or 'Alt_high'
      for a sensitivity test.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      source: the source of parameters for the temperature IRF. default is AR5,
      'Alt', 'Alt_low', and 'Alt_high' are also options.
      decay: a boolean variable for if methane decay to CO2 should be included
      """
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years) 
      f = interp1d(years, emission, kind=kind)
      time = np.linspace(years[0], end, end/tstep + 1)    
      ch4_inter_emissions = f(time)
      ch4_atmos = np.resize(fftconvolve(CH4_AR5(time), ch4_inter_emissions),
                            time.size) * tstep
      co2 = np.resize(fftconvolve(ch42co2(time), ch4_inter_emissions),
                      time.size) * tstep
      co2_atmos = np.resize(fftconvolve(CO2_AR5(time), co2),
                            time.size) * tstep
      if decay == True:
           rf = ch4_atmos * ch4_re + co2_atmos * co2_re
      else:
          rf = ch4_atmos * ch4_re
      if source == 'AR5':
          temp = np.resize(fftconvolve(AR5_GTP(time), rf), time.size) * tstep
      elif source == 'Alt':
          temp = np.resize(fftconvolve(Alt_GTP(time), rf), time.size) * tstep
      elif source == 'Alt_low':
          temp = np.resize(fftconvolve(Alt_low_GTP(time), rf), time.size) * tstep
      elif source == 'Alt_high':
          temp = np.resize(fftconvolve(Alt_high_GTP(time), rf), time.size) * tstep
  
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
      
      return temp
  
  def CH4_temp_cc(emission, years, tstep=0.01, kind='linear', source='AR5',
               decay=True): 
      """Transforms an array of methane emissions into temperature with user-defined
      time-step. Default temperature IRF is from AR5, use 'Alt_low' or 'Alt_high'
      for a sensitivity test. Accounts for climate-carbon feedbacks.
      
      emission: an array of emissions, should be same size as years
      years: an array of years at which the emissions take place
      tstep: time step to be used in the calculations
      kind: the type of interpolation to use; can be linear or cubic
      source: the source of parameters for the temperature IRF. default is AR5,
      'Alt', 'Alt_low', and 'Alt_high' are also options.
      decay: a boolean variable for if methane decay to CO2 should be included
      """
      gamma = (44.0/12.0) * 10**12
      
      if min(years) > 0:
          years = years - min(years)
      
      end = max(years) 
      f = interp1d(years, emission, kind=kind)
      time = np.linspace(years[0], end, end/tstep + 1)    
      ch4_inter_emissions = f(time)
      ch4_atmos = np.resize(fftconvolve(CH4_AR5(time), ch4_inter_emissions),
                            time.size) * tstep
      co2 = np.resize(fftconvolve(ch42co2(time), ch4_inter_emissions),
                      time.size) * tstep
      co2_atmos = np.resize(fftconvolve(CO2_AR5(time), co2),
                            time.size) * tstep
      cc_co2 = CH4_cc_tempforrf(emission, years) * gamma
      cc_co2_atmos = np.resize(fftconvolve(CO2_AR5(time), cc_co2),
                            time.size) * tstep
  	
      if decay == True:
           rf = ch4_atmos * ch4_re + (co2_atmos + cc_co2_atmos) * co2_re
      else:
          rf = ch4_atmos * ch4_re + cc_co2_atmos * co2_re
      if source == 'AR5':
          temp = np.resize(fftconvolve(AR5_GTP(time), rf), time.size) * tstep
      elif source == 'Alt':
          temp = np.resize(fftconvolve(Alt_GTP(time), rf), time.size) * tstep
      elif source == 'Alt_low':
          temp = np.resize(fftconvolve(Alt_low_GTP(time), rf), time.size) * tstep
      elif source == 'Alt_high':
          temp = np.resize(fftconvolve(Alt_high_GTP(time), rf), time.size) * tstep
  
      fil = np.zeros_like(time, dtype=bool)
      for i in time:
          if i == int(i):
              fil[i/tstep] = True
      
      return temp[fil]