module m_ECOSM_initialize_carb
      
      contains

c--------------------------------------------------------------
c---read annual mean atmospheric CO2 from file starting in 1948
c--------------------------------------------------------------      
       subroutine read_pco(atmco2)
       use mod_xc
  
       implicit none

       include 'carbvar.h'
       integer iip, nco
       real, intent(out), dimension(69) :: atmco2
c-------------------------------------------------------------------------------
c read atmospheric co2 19482011 reconstruction; 20802100 A1B future projection
      nco = 2640
c      open(unit=nco,file= pathpc//'atmpco/atmco220482100',
      open(unit=nco,file= 'atmco219482016',
     &form='formatted')
c      do iip=1,53
      do iip=1,69
      read(nco,*) atmco2(iip)
      enddo
      close(nco)
      iyco2=1948
      return
      end
      
c-----------------------------------------------------------
c-- add seasonal cycle and compile monthy mean atmospheric CO2
c-------------------------------------------------------------- 
      subroutine pCO2_month(taup1,pat1,pCO2a)
      use mod_xc
      use mod_year_info, only : rt
      implicit none
      include 'carbvar.h'
      integer yearnr,monthnr, iynum, i,j,l
      real pCO2m,eew, esl
      real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), 
     & intent(out) :: pCO2a
      real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), 
     & intent(in) :: taup1,pat1

c      type(year_info), intent(in) :: rt

      yearnr=rt%iyy
      monthnr=rt%imm
c-----pCO2vec contains annual atmpco from 1948-2016 
      iynum=yearnr-iyco2+1
      pCO2m=pCO2vec(iynum)+seasco(monthnr)
      do j=1-margin,jj+margin
        do l=1,isp(j)
          do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l))
            TT=taup1(i,j)-273.15
            eew=esatw(TT)/9.81*10.**(-4.0)
            esl=pat1(i,j)/9.81*10.**(-2.0)
            pCO2a(i,j)=pCO2m*(esl-eew)*0.997
          enddo 
        enddo
      enddo

 
      return
      end
c-----------------------------------------------------------------------
      function esatw(tt)
c-----------------------------------------------------------------------
c     saturation of vapor pressure at water temperature t in Pa
c     Magnus equation
c-----------------------------------------------------------------------
      ai    = (7.5*tt)/(237.3+tt)
      esatw = 610.7*(10.**ai)
      return
       end
c----------------------------------------------------------------------
        subroutine shumtaup(shum,press,taup)
c----------------------------------------------------------------------
c  specific humidity to dew point temperature ------------------
c-----------------------
      real shum,press,taup
      real ee

        ee=press/100.*shum/(0.622+0.378*shum)
        if (ee.ge.6.176) then
          taup=873.23-27.7124*log(ee)
     &        -sqrt((873.23-27.7124*log(ee))**2 -374871.12)
        else
          taup=6141./(24.3-log(ee))
        endif
        return
        end

      end module