program kaya

      parameter(imax=19,inow=11)
      parameter(pop_offset=6.,pop_width=35.)

      real year(imax),pop(imax),GDP_cap(imax),W_dol(imax),C_W(imax),
     $     CO2_emit(imax), emit_stab(imax,5),c_free_energy(imax,5),
     $     energy_mod(imax), energy(imax)
c      real pCO2(imax),C350(imax),C450(imax),C550(imax),C650(imax),
c     $     C750(imax),CBAU(imax)
      real pop_mod(imax),GDP_cap_mod(imax),W_dol_mod(imax),
     $     C_W_mod(imax),CO2_emit_mod(imax)
      character*20 inString

      open(7,file="kaya_in.txt")
      read(7,*)
      read(7,*)
      do i=1,inow
         read(7,*) year(i),pop(i),GDP_cap(i),W_dol(i),
     $        C_W(i)
      enddo

      read(7,*)
      read(7,*)
      read(7,*)
      
      do i=inow,imax
         read(7,*) year(i)
c     $        ,C350(i),C450(i),C550(i),
c     $        C650(i),C750(i),CBAU(i)
      enddo

      close(7)

      open(7,file='stab_emit.txt')
      do i=inow,imax
         read(7,*) year(i), (emit_stab(i,j),j=1,5)
      enddo

c      read(5,*) pop_stab, g_GDP_cap, g_W_dol, g_C_W
      call getarg(1, inString)
      read(inString,*) pop_stab

      call getarg(2, inString)
      read(inString,*) g_GDP_cap

      call getarg(3, inString)
      read(inString,*) g_W_dol

      call getarg(4, inString)
      read(inString,*) g_C_W


      pop_preexp = pop_stab - pop_offset

      do i=1,inow
         pop_mod(i) = pop(i)
         energy(i)   = pop(i)*1E9           ! people
     $                   * GDP_cap(i)*1e3   ! $
     $                   * W_dol(i)         ! W
     $                   / 1e12             ! tW
         CO2_emit(i) = pop(i)*1.E9
     $        * GDP_cap(i) * 1.e3
     $        * W_dol(i)
     $        / 1e12
     $        * C_W(i)
      enddo
      do i=inow+1,imax
         year_offset = year(i) - year(inow)
         pop_mod(i) = pop_offset
     $              + pop_preexp
     $                * ( exp( year_offset/pop_width) - 1 )
     $                / ( exp( year_offset/pop_width) + 1 )
      enddo

      do i=1,imax
         year_offset = year(i) - year(inow)
         GDP_cap_mod(i) = GDP_cap(inow) * exp( year_offset
     $                                       * g_GDP_cap
     $                                       / 100.
     $                                       )
         W_dol_mod(i) = W_dol(inow) * exp( year_offset
     $                                       * g_W_dol
     $                                       / 100.
     $                                       )
         C_W_mod(i) = C_W(inow) * exp( year_offset
     $                                       * g_C_W
     $                                       / 100.
     $                                       )
         energy_mod(i)   = pop_mod(i)*1E9   ! people
     $                   * GDP_cap_mod(i)*1e3   ! $
     $                   * W_dol_mod(i)     ! W
     $                   / 1e12             ! tW

         CO2_emit_mod(i) = pop_mod(i)*1E9   ! people
     $                   * GDP_cap_mod(i)*1e3   ! $
     $                   * W_dol_mod(i)     ! W
     $                   / 1e12             ! tW
     $                   * C_W_mod(i)       ! Gt

      enddo

      do i=inow,imax
         do j=1,5
            c_free_energy(i,j) = ( CO2_emit_mod(i)
     $                           - emit_stab(i,j)
     $                           )          ! Gt
     $                           / 0.78     ! tW
         enddo
      enddo

      do i=1,inow
         write(6,10) year(i),pop_mod(i),
     $                      GDP_cap_mod(i),
     $                      W_dol_mod(i),
     $                      C_W_mod(i),
     $                      CO2_emit_mod(i),
     $                      energy_mod(i),
     $                      0., 0., 0., 0., 0.,
     $                      pop(i),
     $                      GDP_cap(i),
     $                      W_dol(i),
     $                      C_W(i),
     $                      CO2_emit(i),
     $                      energy(i)
         
      enddo
      do i=inow+1,imax
         write(6,10) year(i),pop_mod(i),
     $                      GDP_cap_mod(i),
     $                      W_dol_mod(i),
     $                      C_W_mod(i),
     $                      CO2_emit_mod(i),
     $                      energy_mod(i),
     $                      (c_free_energy(i,j),j=1,5)
      enddo
 10   format(F10.0, 18F9.3)
      end