c $Log$
c =====================================================================
c
c     SUBROUTINE get_atmco2
c
c     PURPOSE: determines atmospheric CO2 concentration at the
c                current time step of the model. Based on code
c                supplied by OCMIP2 (J. Orr)
c
c
c     VARIABLES:
c 
c     IN:  through modules
c           pert_time : perturbation time (real years)
c
c     OUT: atmco2_t
c
c     AUTHOR:  Nicolas Gruber (gruber@splash.princeton.edu)  
c
c     VERSION: MOM3
c
c     REVISIONS:
c
c     date     author   remarks
c     
c     07.12.98  ng      started coding, adapted from update_bc written
c                         by Rick Slater
c     14.12.98  ng      consider using subroutine timeinterp.F
c     17.12.98  ng      corrected error in bc_fract_next_mo calculation
c      
c =====================================================================
c
      subroutine get_atmco2(pert_time,pert_atmco2) 
c
c
c-----------------------------------------------------------------------
c     common blocks
c-----------------------------------------------------------------------
c
c      implicit none
c
c-----------------------------------------------------------------------
c     local variables
c-----------------------------------------------------------------------
c	
      integer maxrec
      parameter (maxrec=1200)
      
      real, intent(in)	:: pert_time
      real, intent(out)	:: pert_atmco2
      
      integer, save		:: nrec_co2
      real, save	:: year_co2(maxrec), atmco2(maxrec)
  
      integer :: n

      real	:: yr
      real	:: x

      integer, save	:: ientry = 0

      real, parameter	:: yrco2min = 1765.0
      real, parameter	:: yrco2max = 2010.5
c
c-----------------------------------------------------------------------
c     Counter "ientry" for number of entries in this SUBROUTINE
c     (Read atmospheric records on 1st entry ONLY)
c     use CIS9 scenario
c-----------------------------------------------------------------------
c	
      ientry = ientry + 1

      if (ientry .eq. 1) then
         call read_co2atm('CIS9',nrec_co2, year_co2, atmco2)
      endif
c
c-----------------------------------------------------------------------
c     assign perturbation time (real years e.g. 1935.67)
c-----------------------------------------------------------------------
c	
      if (pert_time .lt. yrco2min) then 
         yr = yrco2min
      else if((pert_time .ge. yrco2min)  .and.  
     $        (pert_time .le. yrco2max)) then 
         yr = pert_time
      else if (pert_time .gt. yrco2max) then 
         yr = yrco2max
      end if 
c
c-----------------------------------------------------------------------
c     Find relative POSITION n for pert_time in CO2 record
c-----------------------------------------------------------------------
c	
      call locate(year_co2,nrec_co2, yr, n)
c
c-----------------------------------------------------------------------
c     Determine linear interpolation factor "x"
c-----------------------------------------------------------------------
c	
      x = (yr - year_co2(n)) / (year_co2(n+1) - year_co2(n))
c
c-----------------------------------------------------------------------
c     Perform temporal interpolation for atmospheric CO2, subtract
c      pre-industrial concentration to get perturbation CO2
c-----------------------------------------------------------------------
c    
      pert_atmco2 = (atmco2(n) * (1. - x) + atmco2(n+1) * x) -
     $     atmco2(1)
*      print *, '***************************************************'
*      print*, 'Perturbation of ATM CO2', yr, year_co2(n),pert_atmco2
*      print *, '***************************************************'
c     
c-----------------------------------------------------------------------
c     end of subroutine
c-----------------------------------------------------------------------
c
      return
      end subroutine  get_atmco2