C**** This code calculates changes to the atmospheric methane C**** concentration as a function of anomalous sources of methane into C**** the atmosphere. It uses a look up table with linear interpolation C**** from results using the GISS 2D tropospheric-stratospheric chemistry C**** described in Schmidt and Shindell, Paleoceanography, 2002. C**** C**** Note that these calculations assume modern (or near modern) C**** temperatures and humidity levels. In very different climates this C**** may be important. However, note that this calculation is not C**** directly affected by the base CO2 level. C**** C**** Compilation must be done using a fortran 90 compiler or better, C**** but conversion to F77 or C would not be onerous. C**** C**** The parameters must first be intialised using C**** call init_methane(init_CH4) C**** where init_CH4 is the base (equilibirum) level of methane desired. C**** Then, at every time step the main routine is called: C**** call methane(em_anomaly,CH4,H2O_strat,dCO2_anomaly,dt) C**** This takes the em_anomaly rate (Gt carbon/yr) over the time step C**** dt (years) and uses it to update the CH4 and H2O_strat C**** concentrations, and calculate the dCO2_anomaly rate C**** (in Gt carbon/yr). C**** C**** Please direct any questions, complaints to: C**** Gavin Schmidt: gschmidt@giss.nasa.gov, Tel: 212 678 5627 C**** Drew Shindell: dshindell@giss.nasa.gov, Tel: 212 678 5561 C**** module methane_parameters C**** this module contains the parameters needed to calculate the C**** tropospheric methane sinks, stratospheric water vapour sources C**** and carbon dioxide sources. C**** These are set from routine init_methane. implicit none save C**** Stratosphere assumed to be above 200mb. C**** sfrac: fraction of strat. methane oxidised real, parameter :: sfrac = 0.21 C**** tauS mean residence time for stratospheric air (years) real, parameter :: tauS = 5. C**** tauW mean residence time for stratospheric water (years) C**** nu timescale factor for tauW = nu * tauS real, parameter :: nu = 0.8 real, parameter :: tauW = nu * tauS C**** pre-industrial value of CH4 (DO NOT CHANGE) real, parameter :: pi_CH4 = 0.7 ! ppmv C**** the pre-industiral level of strat. water vapour (in equilibrium C**** with the pi_CH4 concentration), This could be changed C**** (independently of base_CH4) but affects only the strat water C**** calculation. real, parameter :: pi_H2O_strat = 2.6 ! ppmv C**** parameters for main calculation C**** The behaviour of the CH4 depends on the base level (which C**** determines the initial tropospheric sink) real :: base_CH4 ! ppmv C**** tauM: mean methane residence time (or lifetime) for base case C**** This is defined as the total reservoir for atmospheric methane C**** divided by the surface emissions real :: tauM C**** alpha: tropospheric fraction of total sink for base case real :: alpha C**** s1: Scaling needed if base case is not pre-industrial real :: s1 C**** base_H2O_strat: base stratospheric water mixing ratio (in C**** equilibirum with base_CH4), based on the initial pi_H20_strat C**** plus a contribution from base_CH4 change. real :: base_H2O_strat ! ppmv C**** Conversion factors: (do not take the precision here seriously) real, parameter :: matm=5.137 ! mass of atmosphere (10^6 Gt) real, parameter :: mair=28.966 ! mean atomic mass of air (kg/mol) real, parameter :: mC =12.011 ! atomic mass carbon (kg/mol) real, parameter :: Gtcbyppmv = mC*matm/mair end module methane_parameters subroutine init_methane(init_CH4) C**** Sets the base parameters for the main calculation. Needs C**** to be called only once. use methane_parameters implicit none C**** Input: init_CH4 level C**** The init_CH4 value is assumed to be the value at which CH4 is C**** stable, mean sinks and tropospheric oxidation fractions are C**** calculated accordingly. The closure assumption comes from the C**** estimates of pre-industrial sources and sinks. real, intent(in) :: init_CH4 C**** values used for original calculation (with base_CH4=0.7): real :: tauM0 = 8.4 ! years real :: alpha0= 0.88 real :: sink C**** calculate parameters if base_CH4 != pi_CH4 base_CH4=init_CH4 if (base_CH4 .ne. pi_CH4) then s1 = sink(base_CH4) alpha = alpha0*s1/(alpha0*s1+1.-alpha0) tauM = tauM0/(alpha0*s1+1.-alpha0) base_H2O_strat = pi_H2O_strat + 2.*nu*sfrac*(base_CH4 - pi_CH4) else s1 = 1. alpha = alpha0 tauM = tauM0 base_H2O_strat = pi_H2O_strat end if print*,"Methane calculation initialised." print*," Base CH4 level (ppmv): ", base_CH4 print*," Residence time (years): ", tauM print*," Tropospheric fraction of CH4 sink: ", alpha print*," Stratospheric water vapour conc. (ppmv): ", * base_H2O_strat C**** The implied pre-industrial sources are set by the combination C**** of the base lifetime, pre-industrial mixing ratio and the estimate C**** of the number of moles in the atmosphere. C**** pi_source = pi_CH4*Gtcbyppmv/tauM0 = 0.177 Gt carbon/year print*," Multiple of pre-industrial sources (0.177 Gt carbon/yr)" print*," required to produce this base CH4 level: ",base_CH4 * /pi_CH4*(alpha0*s1+1.-alpha0) end subroutine init_methane subroutine methane(em_anomaly,CH4,H2O_strat,dCO2_anomaly,dt) C**** methane: this subroutine calculates the methane, stratospheric C**** water levels and anomalous CO2 production as a function of the C**** emission anomaly, the current CH4 and H2O_strat amounts. C**** The time step for the calculation is set by dt (>= one year). C**** use methane_parameters implicit none C**** Methane mixing ratio (ppmv) real, intent(inout) :: CH4 C**** Stratospheric water mixing ratio (ppmv) real, intent(inout) :: H2O_strat C**** Anomalous amount of carbon dioxide produced (Gt carbon/yr) real, intent(out) :: dCO2_anomaly C**** Anomalous source of methane (Gt carbon/yr) C**** Note Gt of carbon = 3/4 * Gt of methane real, intent(in) :: em_anomaly C**** Time step for methane calculation (years) real, intent(in) :: dt ! should be 1 or greater C**** M: normalised methane conc. C**** lambda: normalised tropospheric sink factor C**** source: normalised anomalous emissions real :: M, lambda, source, sink source = em_anomaly/Gtcbyppmv ! conversion from Gt carbon/yr to ppmv CH4/yr source = source/base_CH4 ! normalise lambda = sink(CH4)/s1 ! relative specific loss term M=CH4/base_CH4 ! non-dimensional methane amount C**** calculate anomalous source of carbon dioxide: C**** All anomalous sinks (soil, tropospheric chemistry, stratospheric C**** oxidation) increase CO2. CO2 is assumed to be stable for the base C**** case dCO2_anomaly=((alpha*lambda*M+(1.-alpha)*M - 1.)/tauM)*base_CH4 dCO2_anomaly=dCO2_anomaly*Gtcbyppmv ! conversion to Gt carbon/yr C**** calculate new strat water amount H2O_strat=H2O_strat + dt*(base_H2O_strat - H2O_strat + * 2.*nu*sfrac*(M -1.)*base_CH4)/tauW C**** Update methane amount M = M + dt*(source + (1.-alpha*lambda*M-(1.-alpha)*M)/tauM) C**** rescale to ppmv CH4=M*base_CH4 return end subroutine methane function sink(CH4) C**** Calculate tropospheric methane sink use methane_parameters, only : pi_CH4 implicit none C**** The look-up data were calculated based on perturbations to the C**** pre-industrial level C**** Look up table data: C**** a = multiple of pre-indust methane conc. (=pi_CH4) C**** b = multiple of pre-indust tropospheric sink for methane C**** = alpha_trop in Schmidt and Shindell, 2001 integer, parameter :: ndat = 8 real,dimension(ndat) :: * a = (/1. ,2.5 ,4. ,7. ,10. ,40. ,100. ,200. /), * b = (/1.0 ,0.808,0.692,0.575,0.435,0.219,0.096,0.088/) real, intent(in) :: CH4 ! current methane mixing ratio (ppmv) real M,sink integer n C**** scale to pre-industrial M = CH4/pi_CH4 do n=1,ndat-1 if (M.ge.a(n).and.M.lt.a(n+1)) then C**** Use linear interpolation btw. data points. C**** Note: linear interpolation does not produce a smooth derivative, C**** therefore second order quantities (such as the equilibrium value C**** as a function of emissions), will not be completely smooth. C**** Feel free to improve the interpolation if this bothers you. sink = b(n) + (M-a(n))*(b(n+1)-b(n))/(a(n+1)-a(n)) return end if end do sink = b(ndat) ! constant above max a(n) return end function sink