program methanomatic

      implicit none
      
      integer ix
      double precision time,
     $     source_natural, source_chronic, source_spike,
     $     spike_gton, spike_duration,
     $     lifetime_0, lifetime, gton_0, age_exp,
     $     CH4_gton, CH4_src, CH4_sink,
     $     CH4_RF_K, CH4_RF,
     $     pCO2, CO2_tau, CO2_RF, CO2_RF_K
      character*20 inString

      call getarg( 1, inString)
      read( inString, * ) source_chronic
      call getarg( 2, inString)
      read( inString, * ) spike_gton
      call getarg( 3, inString)
      read( inString, * ) spike_duration
c      read(5,*) source_chronic,
c     $     spike_gton, spike_duration

      spike_duration = MAX(spike_duration, 1.1d0)

      time = -101
!      lifetime_0 = 7.5
      gton_0 = 1.6
      source_natural = 0.15  ! gton_0 / lifetime_0
      lifetime_0 = gton_0 / source_natural 
      age_exp = 0.26
      CH4_RF_K = 0.5 * 1.4

      CH4_gton = gton_0
      source_spike = spike_gton / spike_duration
      pCO2 = 280.d0
      CO2_tau = 0.0065
      CO2_RF_K = 3.4

      CH4_gton = 1.e-3

      do ix = 1, 201
         
         time = time + 1.0

         CH4_src = source_natural

         if(time .GT. -50) then
            CH4_src = CH4_src + source_chronic
         endif
         if(time .GT. 0.d0 .AND.
     $      time .LT. spike_duration) then
            CH4_src = CH4_src + source_spike
         endif

         lifetime = lifetime_0
     $        * ( CH4_gton / gton_0 )**age_exp
         CH4_sink = CH4_gton / lifetime

         CH4_gton = CH4_gton + CH4_src - CH4_sink

c         CH4_RF = log(CH4_gton / gton_0) / log(2.d0)
c     $        * CH4_RF_K

c         CH4_RF = 1.93 * ( (CH4_gton / 2)**.338
c     $                   - (gton_0 / 2)**.338
c     $                   )  ! fit to Modtran
c     $          * 1.4  ! efficacy

         CH4_RF = -1.0569 * ( log(0.8)
     $                      - log(CH4_gton/2)
     $                      )
     $            -0.16019* ( log(0.8)**2
     $                      -  log(CH4_gton/2)**2
     $                      )
     $        * 1.4 ! efficacy

         if(time .GT. -50) then
            pCO2 = pCO2 * (1.d0 + CO2_tau)
         endif

c         CO2_RF = log(pco2/280.)/log(2.d0) * CO2_RF_K

         CO2_RF = -2.6222 * ( log(280.d0)
     $                      - log(pCO2)
     $                      )
     $            -0.2960 * ( log(280.d0)**2
     $                      - log(pCO2)**2
     $                      )

         write(6,"(8F12.4)")
     $        time, CH4_gton / 2, CH4_src, CH4_sink, lifetime,
     $        CH4_RF, pCO2, CO2_RF

      enddo

      end