      subroutine getco2e(co2epak,co2epal,nlon1,nlat,
     &                   irefyr,irefmn,iday,kount,ijpak,nuem,gg)
C

C     * April 28, 2018 - M. Lazare.    Routine built on GETCO2E, but
C                                      for surface CO2 emissions.
C
c***********************************************************************
c Read monthly averaged surface CO2 emissions from a file every mid month day.
c This is a multi year data set (currently for years 1850-2002) so this
c routine also depends on iyear, passed in via common /keeptim/.
c***********************************************************************

      implicit real (a-h,o-z), integer (i-n) 

      !--- Reference year and month.
      integer irefyr, irefmn

      !--- pak is the interpolated value and pal is the target value
      real co2epak(ijpak), co2epal(ijpak)

      real gg(*)

      !--- local
      integer :: year_prev, year_next, mon_prev, mon_next

      !--- keeptime is required for iyear
      COMMON /KEEPTIM/ IYEAR,IYMDH,MYRSSTI,ISAVDTS

      !--- mynode used to control i/o to stdout
      integer*4 mynode
      common /mpinfo/ mynode 

      !--- A list containing mid month days of the year
      integer, parameter, dimension(12) ::
     &  mm_doy = (/ 16,46,75,106,136,167,197,228,259,289,320,350 /)

      !--- A list containing mid month days for each month
      integer, parameter, dimension(12) ::
     &  mm_dom = (/ 16,14,16,15,16,15,16,16,15,16,15,16 /)
C----------------------------------------------------------------------

c.....Determine the year for which data is required

      if (irefyr .gt. 0) then
        !--- When irefyr > 0 read emissions from the data file for
        !--- a repeated annual cycle of year irefyr
        curr_year=irefyr
      else
        !--- Otherwise read time dependent emissions from the data
        !--- file using iyear to determine the current year
        curr_year=iyear
      endif

c.....check on iday
      if (iday.lt.1 .or. iday.gt.365) then
        if (mynode.eq.0) then
          write(6,*)'GETCO2E: IDAY is out of range. IDAY=',iday
        endif
        call xit("GETCO2E",-1)
      endif

c.....Determine previous and next year/month, relative to iday
      mon_prev=0
      if (irefmn .lt. 0) then
         if (iday.ge.1 .and. iday.lt.mm_doy(1)) then
            year_prev=curr_year-1
            year_next=curr_year
            mon_prev=12
            mon_next=1
         else if (iday.ge.mm_doy(12) .and. iday.le.365) then
            year_prev=curr_year
            year_next=curr_year+1
            mon_prev=12
            mon_next=1
         else
            do i=2,12
               if (iday.ge.mm_doy(i-1) .and. iday.lt.mm_doy(i)) then
                  year_prev=curr_year
                  year_next=curr_year
                  mon_prev=i-1
                  mon_next=i
                  exit
               endif
            enddo
         endif
      else
         !--- The month is not changing.
         mon_prev = irefmn
         mon_next = irefmn
      end if

      if (mon_prev.eq.0) then
        if (mynode.eq.0) then
          write(6,*)'GETCO2E: Problem setting year/month on IDAY=',iday
        endif
        call xit("GETCO2E",-2)
      endif

      if (kount.eq.0) then     !--- Starting up

        !--- Get field for previous mid-month day
        if (mynode.eq.0) then
          write(6,*)"GETCO2E: Reading surface CO2 emission values",
     &      " for year,month: ",year_prev,mon_prev," IDAY=",iday,
     &      " KOUNT=",kount
        endif
        ib2p=year_prev*100+mon_prev
        rewind nuem
        call getggbx(co2epak,nc4to8(" CO2"),nuem,
     &               nlon1,nlat,ib2p,1,gg)

        !--- Get field for next mid-month day
        if (mynode.eq.0) then
          write(6,*)"GETCO2E: Reading surface CO2 emission values",
     &      " for year,month: ",year_next,mon_next," IDAY=",iday,
     &      " KOUNT=",kount
        endif
        ib2n=year_next*100+mon_next
        rewind nuem
        call getggbx(co2epal,nc4to8(" CO2"),nuem,
     &               nlon1,nlat,ib2n,1,gg)

        !--- Interpolate to current day
        day1=real(mm_doy(mon_prev)) 
        day2=real(iday)
        day3=real(mm_doy(mon_next))
        if(day2.lt.day1) day2=day2+365. 
        if(day3.lt.day2) day3=day3+365. 
        w1=(day2-day1)/(day3-day1)
        w2=(day3-day2)/(day3-day1)
        if (mynode.eq.0) then
          write(6,*)'GETCO2E: Interpolating between days ',
     &      mm_doy(mon_prev),' and ',mm_doy(mon_next),' to day ',iday,
     &      ' using weights ',w1,w2
        endif

        do i=1,ijpak
          co2epak(i) = w1*co2epal(i) + w2*co2epak(i)
        enddo

      else     !--- Continue a run

        if (mynode.eq.0) then
          write(6,*)"GETCO2E: Reading surface CO2 emission values",
     &      " for year,month: ",year_next,mon_next," IDAY=",iday,
     &      " KOUNT=",kount
        endif
        !--- Get field for next mid-month day
        ib2n=year_next*100+mon_next
        rewind nuem
        call getggbx(co2epal,nc4to8(" CO2"),nuem,
     &               nlon1,nlat,ib2n,1,gg)
      endif

      return
      end
