! Calculate the total emission by collecting output from CABLE, ocean model, 
! published Fossil Fuel emission and published Landuse emission.  
! Then calculate the new CO2 concentration for the simulation of next year.  
! Written by Bernard.Pak@csiro.au 
! 3 May 2011 
!  
! compile on vayu using
! ifort -o atmCO2conc  atmCO2conc.f90  -L /apps/netcdf/3.6.3/lib/Intel/ -lnetcdf -I /apps/netcdf/3.6.3/include/Intel/
!ifort -o atmCO2conc atmCO2conc.f90 -L/opt/netcdf/lib -lnetcdf -I/opt/netcdf/include
!
PROGRAM atmCO2conc
   use getopt_m
  IMPLICIT NONE
  INTEGER :: iyear, fyear, fmon
  INTEGER :: ii, jj, yy    ! counter
  INTEGER :: ioerror
  REAL    :: LC, LU, FF, OC, monFF
  REAL    :: oldconc, newconc
  INTEGER, DIMENSION(11) :: futureYr
  REAL,    DIMENSION(11) :: FFfuture, LUfuture
  CHARACTER(LEN=150) :: dummy
  DATA futureYr /  2005,  2010,  2020,  2030,  2040,  2050, &
                   2060,  2070,  2080,  2090,  2100 /
  DATA FFfuture / 7.971, 8.926,11.538,13.839,16.787,20.205, &
                 23.596,25.962,27.406,28.337,28.740 /
  DATA LUfuture / 1.196, 1.044, 0.906, 0.715, 0.645, 0.576, &
                  0.501, 0.412, 0.309, 0.194, 0.077 /
   integer :: nopt, opt, ios
   character(len=80) :: optarg
  LOGICAL :: PISPIN  ! remember to switch on/off (Q.Zhang 13/07/2011)
  PISPIN = .True.

! get actual year to use for emissions files
   do 
      call getopt("c:",nopt,opt,optarg)
      if ( opt == -1 ) then
         exit ! End of options
      end if
      select case ( char(opt) )
      case ( "c" )
         read(optarg,*) iyear
	  print *, iyear
      case default
         print*, " Error: Unknown option "
         print*, "Usage: atmCO2conc -c year"
         stop
      end select
   end do
        
! Obtain the former CO2 concentration from the radiation file
  OPEN(UNIT=9,FILE='co2_data.18l')
  READ(9,*)
  DO ii = 1, 3
    READ(9,'(a150)') dummy
  enddo
   READ(9,*) oldconc
  oldconc = oldconc *1e6  ! convert to ppm
  CLOSE(9)
        
!  OPEN(UNIT=9,FILE='co2_data')
!  READ(9,*)
!  DO ii = 1, 50
!    READ(9,'(a150)') dummy
!    IF (dummy(3:7) == 'fixed') THEN
!      READ(dummy(14:19),'(f6.2)') oldconc
!      EXIT
!    ENDIF
!  ENDDO
!  CLOSE(9)


! Read NEE from CABLE output file
! and sum up to global total
  LC = 0.0
  CALL getNEE(LC)

! read ocean co2 fluxes from output file
! and sum up to give global annual total
! by default the ocean uptake is positive
  OC=0
  call getOCN(OC)
  OC = OC*(-1.)


 print *,iyear
! Read landuse emissions
  IF (iyear <= 2005) THEN
    OPEN(UNIT=11,FILE='carbon_emissions_landuse_20person_global.txt')
    READ(11,*)
    DO yy = 1850, 2005
      READ(11,*) fyear, LU
      IF (yy /= fyear) THEN
        PRINT *, 'yy, fyear = ', yy, fyear, ' not matching.'
        STOP
      ENDIF
      IF (iyear == fyear) EXIT
    ENDDO
    CLOSE(11)
  ELSE     ! interpolate linearly
    CALL interp(iyear,futureYr,LUfuture,LU)
  ENDIF

  ! Read fossil fuel emission
  FF = 0.0
  IF (iyear <= 2007) THEN
    OPEN(UNIT=12,FILE='CMIP5_gridcar_CO2_emissions_fossil_fuel_Andres_1751-2007_monthly_SC.txt')
    DO yy = 1, 1189
      READ(12,*)
    ENDDO
    DO yy = 1850, 2007
      DO jj = 1, 12
        READ(12,*) fyear, fmon, monFF
        IF (yy /= fyear .OR. jj /= fmon) THEN
          PRINT *, 'yy, fyear, jj, fmon not matching: '
          PRINT *,  yy, fyear, jj, fmon
          STOP
        ENDIF
        IF (iyear == fyear) THEN
          FF = FF + monFF
        ENDIF
      ENDDO
      IF (iyear == fyear) EXIT
    ENDDO
    CLOSE(12)
  ELSE     ! interpolate linearly
    CALL interp(iyear,futureYr,FFfuture,FF)
  ENDIF

  ! Q.Zhang (13/07/2011) for PI spinup, assums no LU and FF.
  IF(PISPIN) THEN
    LU = 0.
    FF = 0.
  ENDIF

  ! calculation
  newconc = oldconc + (LU + FF + LC + OC) * 0.47079

  ! output to text file
  OPEN(UNIT=13,FILE='co2.txt')
  WRITE(13,'(2f9.2,i5,4f10.5)') newconc,oldconc,iyear,FF,LU,LC,OC

  ! overwrite namelist file
  ! comment out by Q.Zhang (12/7/2011), in the coupled model, co2 read by cable
  ! is transfered by atmosphere model by default. "fixCO2" from the namelist 
  ! is not used anymore 
!  OPEN(UNIT=9,FILE='cable_mk3l.nml')
!  OPEN(UNIT=14,FILE='temp.nml')
!  READ(9,*) dummy
!  WRITE(14,'(a150)') dummy
!  DO ii = 1, 50
!    READ(9,'(a150)',IOSTAT=ioerror) dummy
!    IF (ioerror < 0) EXIT     ! reached end of file
!    IF (dummy(3:7) == 'fixed') THEN
!      WRITE(dummy(14:19),'(f6.2)') newconc
!    ENDIF
!    WRITE(14,'(a150)') dummy
!  ENDDO
!  CLOSE(9)
!  CLOSE(14)
  
END PROGRAM atmCO2conc


SUBROUTINE interp(iyr,futureYr,predict,outVar)
  IMPLICIT NONE
  INTEGER,                INTENT(IN)    :: iyr
  INTEGER, DIMENSION(11), INTENT(IN)    :: futureYr
  REAL,    DIMENSION(11), INTENT(IN)    :: predict
  REAL,                   INTENT(INOUT) :: outVar
  INTEGER :: jj, yy    ! counter

  DO yy = 2, 11
    IF (iyr <= futureYr(yy)) EXIT
  ENDDO
  outVar = (predict(yy) - predict(yy-1)) / (futureYr(yy) - futureYr(yy-1))
  outVar = outVar * (iyr - futureYr(yy-1)) + predict(yy-1)

END SUBROUTINE interp


SUBROUTINE getNEE(annNEE)
  USE netcdf
  IMPLICIT NONE
  REAL, INTENT(INOUT) :: annNEE
  REAL*8              :: totNEE
  REAL,   DIMENSION(:,:,:,:), ALLOCATABLE :: NEE
  REAL*8, DIMENSION(:,:,:),   ALLOCATABLE :: NEE3d
  REAL*8, DIMENSION(:,:),     ALLOCATABLE :: NEE2d
  REAL,   DIMENSION(:,:,:),   ALLOCATABLE :: patchfrac
  REAL,   DIMENSION(:,:),     ALLOCATABLE :: area
  INTEGER  :: i, j, k, ok, ncid
  INTEGER  :: xdimID, ydimID, pdimID, tdimID, neeID, pfracID, areaID
  INTEGER  :: xdimsize, ydimsize, pdimsize, tdimsize
  INTEGER, DIMENSION(12) :: monLength
  DATA monLength / 31,28,31,30,31,30,31,31,30,31,30,31 /

  ! Read CABLE output file
  ok = NF90_OPEN('out_cable.nc',0,ncid)
! no land output file, land carbon uptake=0
  IF(ok/=NF90_NOERR) return 

  ok = NF90_INQ_DIMID(ncid,'x', xdimID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring xdimID'
  ok = NF90_INQUIRE_DIMENSION(ncid,xdimID,LEN=xdimsize)
  IF(ok/=NF90_NOERR) STOP 'error inquiring xdimsize'
  IF(xdimsize /= 64) STOP 'error: xdimsize not 64'

  ok = NF90_INQ_DIMID(ncid,'y', ydimID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring ydimID'
  ok = NF90_INQUIRE_DIMENSION(ncid,ydimID,LEN=ydimsize)
  IF(ok/=NF90_NOERR) STOP 'error inquiring ydimsize'
  IF(ydimsize /= 56) STOP 'error: ydimsize not 56'

  ok = NF90_INQ_DIMID(ncid,'patch', pdimID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring pdimID'
  ok = NF90_INQUIRE_DIMENSION(ncid,pdimID,LEN=pdimsize)
  IF(ok/=NF90_NOERR) STOP 'error inquiring pdimsize'

  ok = NF90_INQ_DIMID(ncid,'time', tdimID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring tdimID'
  ok = NF90_INQUIRE_DIMENSION(ncid,tdimID,LEN=tdimsize)
  IF(ok/=NF90_NOERR) STOP 'error inquiring tdimsize'

  ALLOCATE(NEE(xdimsize,ydimsize,pdimsize,tdimsize))
  ok = NF90_INQ_VARID(ncid,'NEE',neeID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring neeID'
  ok= NF90_GET_VAR(ncid,neeID,NEE)
  IF(ok/=NF90_NOERR) STOP 'error getting NEE'

  ALLOCATE(patchfrac(xdimsize,ydimsize,pdimsize))
  ok = NF90_INQ_VARID(ncid,'patchfrac',pfracID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring pfracID'
  ok= NF90_GET_VAR(ncid,pfracID,patchfrac)
  IF(ok/=NF90_NOERR) STOP 'error getting patchfrac'

  ok = NF90_CLOSE(ncid)

  ! Obtain grid area
  ALLOCATE(area(xdimsize,ydimsize))
  ok = NF90_OPEN('csiro_agcm_area.nc',0,ncid)
  IF(ok/=NF90_NOERR) STOP 'error opening nc file for area'
  ok = NF90_INQ_VARID(ncid,'area',areaID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring areaID'
  ok= NF90_GET_VAR(ncid,areaID,area)
  IF(ok/=NF90_NOERR) STOP 'error getting area'
  ok = NF90_CLOSE(ncid)

  ! calculate total
  ALLOCATE(NEE3d(xdimsize,ydimsize,pdimsize),NEE2d(xdimsize,ydimsize))
  NEE3d = 0.0
  IF (tdimsize == 12) THEN ! monthly output
    DO i = 1, 12
      NEE3d(:,:,:) = NEE3d(:,:,:) &
          + NEE(:,:,:,i) * monLength(i) * 24.0 * 3600.0 * 1.201E-5  ! gC/m2/yr
    ENDDO
  ELSE IF (tdimsize == 365 .OR. tdimsize == 366) THEN ! daily output
    NEE3d = SUM(NEE,4) * 24.0 * 3600.0 * 1.201E-5  ! gC/m2/yr
  ELSE 
    PRINT *, 'tdimsize not expected: ', tdimsize
    PRINT *, 'need more codes in subroutine getNEE'
    STOP
  ENDIF

  NEE2d = 0.0
  IF (pdimsize == 1) THEN
    WHERE(patchfrac(:,:,1) > 0.005) NEE2d(:,:) = NEE3d(:,:,1)
  ELSE IF (pdimsize > 1) THEN
    DO k = 1, pdimsize
    DO j = 1, ydimsize
    DO i = 1, xdimsize
      IF (patchfrac(i,j,k) > 0.005) &
        NEE2d(i,j) = NEE2d(i,j) + NEE3d(i,j,k) * patchfrac(i,j,k)
    ENDDO
    ENDDO
    ENDDO
  ELSE
    PRINT *, 'pdimsize is negative: ', pdimsize
    STOP
  ENDIF

  totNEE = 0.0
  DO j = 1, ydimsize
  DO i = 1, xdimsize
    totNEE = totNEE + NEE2d(i,j) * area(i,j) / 1.0E15   ! PgC/yr
  ENDDO
  ENDDO
  annNEE = REAL(totNEE)

END SUBROUTINE getNEE

SUBROUTINE getOCN(annOCN)
  USE netcdf
  IMPLICIT NONE
  REAL, INTENT(INOUT) :: annOCN
  REAL,   DIMENSION(:,:,:), ALLOCATABLE :: NEE
  REAL,   DIMENSION(:,:),     ALLOCATABLE :: area
  INTEGER  :: i, j, k, ok, ncid
  INTEGER  :: xdimID, ydimID, pdimID, tdimID, neeID, pfracID, areaID
  INTEGER  :: xdimsize, ydimsize, pdimsize, tdimsize
  INTEGER, DIMENSION(12) :: monLength
  REAL*8              :: totNEE
  DATA monLength / 31,28,31,30,31,30,31,31,30,31,30,31 /

  ! Read output file
  ok = NF90_OPEN('annual.nc',0,ncid)
! no ocean output file, ocean carbon uptake =0
  IF(ok/=NF90_NOERR) return 

  ok = NF90_INQ_DIMID(ncid,'lonts', xdimID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring xdimID'
  ok = NF90_INQUIRE_DIMENSION(ncid,xdimID,LEN=xdimsize)
  IF(ok/=NF90_NOERR) STOP 'error inquiring xdimsize'

  ok = NF90_INQ_DIMID(ncid,'latts', ydimID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring ydimID'
  ok = NF90_INQUIRE_DIMENSION(ncid,ydimID,LEN=ydimsize)
  IF(ok/=NF90_NOERR) STOP 'error inquiring ydimsize'

  ok = NF90_INQ_DIMID(ncid,'time', tdimID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring tdimID'
  ok = NF90_INQUIRE_DIMENSION(ncid,tdimID,LEN=tdimsize)
  IF(ok/=NF90_NOERR) STOP 'error inquiring tdimsize'

  ALLOCATE(NEE(xdimsize,ydimsize,tdimsize))
  ok = NF90_INQ_VARID(ncid,'bflx04',neeID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring neeID'
  ok= NF90_GET_VAR(ncid,neeID,NEE)
  IF(ok/=NF90_NOERR) STOP 'error getting NEE'

  ok = NF90_CLOSE(ncid)

  ! Obtain grid area
  ALLOCATE(area(xdimsize,ydimsize))
  ok = NF90_OPEN('grid_spec_mk3l_128_112_21_v2.nc',0,ncid)
  IF(ok/=NF90_NOERR) STOP 'error opening nc file for area'
  ok = NF90_INQ_VARID(ncid,'dats',areaID)
  IF(ok/=NF90_NOERR) STOP 'error inquiring areaID'
  ok= NF90_GET_VAR(ncid,areaID,area)
  IF(ok/=NF90_NOERR) STOP 'error getting area'
  ok = NF90_CLOSE(ncid)


  ! calculate total

  totNEE = 0.0
  DO j = 1, ydimsize
  DO i = 1, xdimsize
   if (NEE(i,j,1) .le. -9999) NEE(i,j,1)=0
   totNEE = totNEE + NEE(i,j,1) *area(i,j) *86400.*365.*1e-5*12.0E-15 ! PgC/yr
  ENDDO
  ENDDO
  annOCN = REAL(totNEE)

END SUBROUTINE getOCN