#include <misc.h>


module ghg_surfvals 11,3

!-----------------------------------------------------------------------------------
! Purpose: Provides greenhouse gas (ghg) values at the Earth's surface.
!          These values may be time dependent.
!
! Author: Brian Eaton (assembled module from existing scattered code pieces)
!-----------------------------------------------------------------------------------

   use shr_kind_mod, only: r8=>shr_kind_r8
   use pmgrid,       only: masterproc
   use abortutils,   only: endrun

!-----------------------------------------------------------------------
!- module boilerplate --------------------------------------------------
!-----------------------------------------------------------------------
   implicit none
   private                   ! Make default access private
   save

! Public methods
   public ::&
      ghg_surfvals_init,        &! initialize options that depend on namelist input
      ghg_surfvals_set,         &! set ghg surface values when ramp option is on
      ghg_surfvals_ramp,        &! returns true when ramp option is on
      ghg_surfvals_get_co2mmr    ! computes and returns the co2 mass mixing ratio

! Public data
   real(r8), public :: co2vmr = 3.550e-4               ! co2   volume mixing ratio 
   real(r8), public :: n2ovmr = 0.311e-6               ! n2o   volume mixing ratio 
   real(r8), public :: ch4vmr = 1.714e-6               ! ch4   volume mixing ratio 
   real(r8), public :: f11vmr = 0.280e-9               ! cfc11 volume mixing ratio 
   real(r8), public :: f12vmr = 0.503e-9               ! cfc12 volume mixing ratio 
   character(len=16), public :: scenario_ghg = 'FIXED' ! 'FIXED','RAMPED' or 'RAMP_CO2_ONLY'
   integer, public  :: rampYear_ghg = 0                ! ramped gases fixed at this year (if > 0)
   character(len=256), public  :: bndtvghg             ! filename for ramped data
   integer, public  :: ramp_co2_start_ymd = 0          ! start date for co2 ramping (yyyymmdd)
   real(r8), public :: ramp_co2_annual_rate = 1.0      ! % amount of co2 ramping per yr; default is 1% 
   real(r8), public :: ramp_co2_cap = -9999.0          ! co2 ramp cap if rate>0, floor otherwise 
                                                       ! as multiple or fraction of inital value
                                                       ! ex. 4.0 => cap at 4x initial co2 setting 

! Private methods
   private ::&
      ghg_surfvals_set_all,     &! set all ghg surface values when ramp option is on
      ghg_surfvals_set_co2       ! set just co2 surface values when ramp option is on

! Private module data
   logical :: doRamp_ghg    ! true => turn on ramping for ghg
   logical :: ramp_just_co2 ! true => ramping to be done just for co2 and not other ghg's
   integer :: fixYear_ghg   ! year at which Ramped gases are fixed
   integer :: co2_start     ! date at which co2 begins ramping
   real(r8) :: co2_daily_factor    ! daily multiplier to achieve annual rate of co2 ramp
   real(r8) :: co2_limit    ! value of co2vmr where ramping ends
   real(r8) :: co2_base     ! initial co2 volume mixing ratio, before any ramping
   integer :: ntim = -1               ! number of yearly data values
   integer,  allocatable :: yrdata(:) ! yearly data values
   real(r8), allocatable :: co2(:)    ! co2 mixing ratios in ppmv 
   real(r8), allocatable :: ch4(:)    ! ppbv
   real(r8), allocatable :: n2o(:)    ! ppbv
   real(r8), allocatable :: f11(:)    ! pptv
   real(r8), allocatable :: f12(:)    ! pptv
   real(r8), allocatable :: adj(:)    ! unitless adjustment factor for f11 & f12

!=========================================================================================
contains
!=========================================================================================


subroutine ghg_surfvals_init() 2,8

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Initialize the ramp options that are controlled by namelist input.
! Set surface values at initial time.
! N.B. This routine must be called after the time manager has been initialized
!      since ghg_surfvals_set calls time manager methods.
! 
! Author: B. Eaton - merged code from parse_namelist and rampnl_ghg.
! 
!-----------------------------------------------------------------------
   use infnan,       only: inf

!---------------------------Local variables-----------------------------


   if (scenario_ghg == 'FIXED') then
      doRamp_ghg = .false.
      ramp_just_co2 = .false.
      if (masterproc) &
         write(6,*)'ghg_surfvals_init: ghg surface values are fixed as follows'

   else if (scenario_ghg == 'RAMPED') then
      doRamp_ghg = .true.
      ramp_just_co2 = .false.
      call ghg_ramp_read

      fixYear_ghg = rampYear_ghg     ! set private member to namelist var
      if (masterproc) then
         if ( fixYear_ghg > 0 ) then
            write(6,*) '  FIXED values from year ',fixYear_ghg
         else
            write(6,*) '  RAMPED values initialized to'
         end if
      end if
      call ghg_surfvals_set()

   else if (scenario_ghg == 'RAMP_CO2_ONLY') then
      if(ramp_co2_start_ymd == 0) then
         call endrun ('ghg_surfvals_init: RAMP_CO2_START_YMD must be set for SCENARIO_GHG=RAMP_CO2_ONLY')
      else
         co2_start = ramp_co2_start_ymd
      end if

      if(ramp_co2_annual_rate <= -100.0) then
         write(6,*) 'RAMP_CO2:  invalid ramp_co2_annual_rate= ',ramp_co2_annual_rate
         call endrun ('ghg_surfvals_init: RAMP_CO2_ANNUAL_RATE must be greater than -100.0')
      end if

      doRamp_ghg = .true.
      ramp_just_co2 = .true.
      co2_base = co2vmr        ! save initial setting 
      if (masterproc) &
           write(6,*) '  RAMPED values initialized to'

      co2_daily_factor = (ramp_co2_annual_rate*0.01_r8+1.0)**(1.0_r8/365.0_r8)

      if(ramp_co2_cap > 0.0) then  
         co2_limit = ramp_co2_cap * co2_base
      else                                  ! if no cap/floor specified, provide default
         if(ramp_co2_annual_rate < 0.0) then
            co2_limit = 0.0
         else
            co2_limit = inf
         end if
      end if
      if((ramp_co2_annual_rate<0.0 .and. co2_limit>co2_base) .or. &
         (ramp_co2_annual_rate>0.0 .and. co2_limit<co2_base)) then
         write(6,*) 'RAMP_CO2: ramp_co2_cap is unreachable'
         write(6,*) 'RAMP_CO2: ramp_co2_annual_rate= ',ramp_co2_annual_rate,' ramp_co2_cap= ',ramp_co2_cap
         call endrun('GHG_SURFVALS_INIT:  ramp_co2_annual_rate and ramp_co2_cap incompatible')
      end if

      call ghg_surfvals_set()
   else
      call endrun ('ghg_surfvals_init: input namelist SCENARIO_GHG must be set to either FIXED, RAMPED or RAMP_CO2_ONLY')
   endif

   if (masterproc) then
      write(6,*) '  co2 volume mixing ratio = ',co2vmr
      write(6,*) '  ch4 volume mixing ratio = ',ch4vmr
      write(6,*) '  n2o volume mixing ratio = ',n2ovmr
      write(6,*) '  f11 volume mixing ratio = ',f11vmr
      write(6,*) '  f12 volume mixing ratio = ',f12vmr
   end if

end subroutine ghg_surfvals_init

!=========================================================================================


subroutine ghg_ramp_read() 1,29

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Read ramped greenhouse gas surface data.  
! 
! Author: T. Henderson
! 
!-----------------------------------------------------------------------

   use ioFileMod, only: getfil
#if ( defined SPMD )
   use mpishorthand, only: mpicom, mpiint, mpir8
#endif

   include 'netcdf.inc'

!---------------------------Local variables-----------------------------
   integer :: ncid
   integer :: co2_id
   integer :: ch4_id
   integer :: n2o_id
   integer :: f11_id
   integer :: f12_id
   integer :: adj_id
   integer :: date_id
   integer :: time_id
   integer :: ierror
   character(len=256) :: locfn          ! netcdf local filename to open

   if (masterproc) then
     call getfil (bndtvghg, locfn, 0)
     call wrap_open (trim(locfn), NF_NOWRITE, ncid)
     write(6,*)'GHG_RAMP_READ:  reading ramped greenhouse gas surface data from file ',trim(locfn)
     call wrap_inq_varid( ncid, 'date', date_id )
     call wrap_inq_varid( ncid, 'CO2', co2_id )
     call wrap_inq_varid( ncid, 'CH4', ch4_id )
     call wrap_inq_varid( ncid, 'N2O', n2o_id )
     call wrap_inq_varid( ncid, 'f11', f11_id )
     call wrap_inq_varid( ncid, 'f12', f12_id )
     call wrap_inq_varid( ncid, 'adj', adj_id )
     call wrap_inq_dimid( ncid, 'time', time_id )
     call wrap_inq_dimlen( ncid, time_id, ntim )
   endif
#if (defined SPMD )
   call mpibcast (ntim, 1, mpiint, 0, mpicom)
#endif
   ! these arrays are never deallocated
   allocate ( yrdata(ntim), co2(ntim), ch4(ntim), n2o(ntim),    &
                 f11(ntim), f12(ntim), adj(ntim), stat=ierror )
   if (ierror /= 0) then
     write(6,*)'GHG_RAMP_READ:  ERROR, allocate() failed!'
     call endrun
   endif
   if (masterproc) then
     call wrap_get_var_int   (ncid, date_id, yrdata )
     yrdata = yrdata / 10000
     call wrap_get_var_realx (ncid, co2_id, co2 )
     call wrap_get_var_realx (ncid, ch4_id, ch4 )
     call wrap_get_var_realx (ncid, n2o_id, n2o )
     call wrap_get_var_realx (ncid, f11_id, f11 )
     call wrap_get_var_realx (ncid, f12_id, f12 )
     call wrap_get_var_realx (ncid, adj_id, adj )
     call wrap_close (ncid)
     write(6,*)'GHG_RAMP_READ:  successfully read ramped greenhouse gas surface data from years ',yrdata(1),' through ',yrdata(ntim)
   endif
#if (defined SPMD )
   call mpibcast (co2, ntim, mpir8, 0, mpicom)
   call mpibcast (ch4, ntim, mpir8, 0, mpicom)
   call mpibcast (n2o, ntim, mpir8, 0, mpicom)
   call mpibcast (f11, ntim, mpir8, 0, mpicom)
   call mpibcast (f12, ntim, mpir8, 0, mpicom)
   call mpibcast (adj, ntim, mpir8, 0, mpicom)
   call mpibcast (yrdata, ntim, mpiint, 0, mpicom)
#endif

   return

end subroutine ghg_ramp_read

!=========================================================================================


function ghg_surfvals_ramp()
   logical :: ghg_surfvals_ramp
   ghg_surfvals_ramp = doRamp_ghg
end function ghg_surfvals_ramp

!=========================================================================================


function ghg_surfvals_get_co2mmr() 1,1
  use physconst,    only: mwdry, mwco2

  real(r8), parameter :: rmwco2 = mwco2/mwdry    ! ratio of molecular weights of co2 to dry air
  real(r8) :: ghg_surfvals_get_co2mmr            ! co2 mass mixing ratio

  ghg_surfvals_get_co2mmr = rmwco2 * co2vmr  

end function ghg_surfvals_get_co2mmr

!=========================================================================================



subroutine ghg_surfvals_set() 3,4

   use time_manager, only: get_curr_date, is_end_curr_day

!---------------------------Local variables-----------------------------

   integer  :: yr, mon, day, ncsec ! components of a date
   integer  :: ncdate              ! current date in integer format [yyyymmdd]

   if(ramp_just_co2) then
      call ghg_surfvals_set_co2()
   else
      call ghg_surfvals_set_all()
   end if

   if (masterproc .and. is_end_curr_day()) then
      call get_curr_date(yr, mon, day, ncsec)
      ncdate = yr*10000 + mon*100 + day
      write(6,*) 'ghg_surfvals_set: ncdate= ',ncdate,' co2vmr=',co2vmr
   end if

   return
end subroutine ghg_surfvals_set

!=========================================================================================


subroutine ghg_surfvals_set_all() 1,6
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Computes greenhouse gas volume mixing ratios via interpolation of
! yearly input data.
! 
! Author: B. Eaton - updated ramp_ghg for use in ghg_surfvals module
! 
!-----------------------------------------------------------------------
   use time_manager, only: get_curr_date, get_curr_calday
   use timeinterp,   only: validfactors

!---------------------------Local variables-----------------------------

   integer yrmodel           ! model year
   integer nyrm              ! year index
   integer nyrp              ! year index
   integer :: yr, mon, day   ! components of a date
   integer :: ncdate         ! current date in integer format [yyyymmdd]
   integer :: ncsec          ! current time of day [seconds]

   real(r8) :: calday            ! current calendar day
   real(r8) doymodel             ! model day of year
   real(r8) doydatam             ! day of year for input data yrdata(nyrm)
   real(r8) doydatap             ! day or year for input data yrdata(nyrp)
   real(r8) deltat               ! delta time
   real(r8) fact1, fact2         ! time interpolation factors
   real(r8) cfcscl               ! cfc scale factor for f11

!
! ---------------------------------------------------------------------
!
   calday = get_curr_calday()
   call get_curr_date(yr, mon, day, ncsec)
   ncdate = yr*10000 + mon*100 + day
!
! determine index into input data
!
   if ( fixYear_ghg > 0) then
      yrmodel  = fixYear_ghg
   else
      yrmodel  = yr   
   end if

   nyrm       = yrmodel - yrdata(1) + 1
   nyrp       = nyrm + 1
!
! if current date is before yrdata(1), quit
!
   if (nyrm < 1) then
      write(6,*)'ghg_surfvals_set_all: data time index is out of bounds'
      write(6,*)'nyrm = ',nyrm,' nyrp= ',nyrp, ' ncdate= ', ncdate
      call endrun
   endif
!
! if current date later than yrdata(ntim), call endrun.
! if want to use ntim values - uncomment the following lines
! below and comment the call to endrun and previous write
!
   if (nyrp > ntim) then
      call endrun ('ghg_surfvals_set_all: error - current date is past the end of valid data')
!         write(6,*)'ghg_surfvals_set_all: using ghg data for ',yrdata(ntim)
!         co2vmr = co2(ntim)*1.e-06
!         ch4vmr = ch4(ntim)*1.e-09
!         n2ovmr = n2o(ntim)*1.e-09
!         f11vmr = f11(ntim)*1.e-12*(1.+cfcscl)
!         f12vmr = f12(ntim)*1.e-12
!         co2mmr = rmwco2 * co2vmr
!         return
   endif
!
! determine time interpolation factors, check sanity
! of interpolation factors to within 32-bit roundoff
! assume that day of year is 1 for all input data
!
   doymodel = yrmodel*365.    + calday
   doydatam = yrdata(nyrm)*365. + 1.
   doydatap = yrdata(nyrp)*365. + 1.
   deltat   = doydatap - doydatam
   fact1    = (doydatap - doymodel)/deltat
   fact2    = (doymodel - doydatam)/deltat

   if (.not. validfactors (fact1, fact2)) then
      write(6,*)'ghg_surfvals_set_all: Bad fact1 and/or fact2=',fact1,fact2
   end if
!
! do time interpolation:
!   co2     in ppmv
!   n2o,ch4 in ppbv
!   f11,f12 in pptv
!
   co2vmr = (co2(nyrm)*fact1 + co2(nyrp)*fact2)*1.e-06
   ch4vmr = (ch4(nyrm)*fact1 + ch4(nyrp)*fact2)*1.e-09
   n2ovmr = (n2o(nyrm)*fact1 + n2o(nyrp)*fact2)*1.e-09

   cfcscl = (adj(nyrm)*fact1 + adj(nyrp)*fact2)
   f11vmr = (f11(nyrm)*fact1 + f11(nyrp)*fact2)*1.e-12*(1.+cfcscl)
   f12vmr = (f12(nyrm)*fact1 + f12(nyrp)*fact2)*1.e-12

   return
end subroutine ghg_surfvals_set_all

!=========================================================================================


subroutine ghg_surfvals_set_co2() 1,4
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Computes co2 greenhouse gas volume mixing ratio via ramping info 
! provided in namelist var's
! 
! Author: B. Eaton - updated ramp_ghg for use in ghg_surfvals module
! 
!-----------------------------------------------------------------------
   use shr_kind_mod, only: r8 => shr_kind_r8
   use time_manager, only: get_curr_date, timemgr_datediff

!---------------------------Local variables-----------------------------

   real(r8) :: daydiff             ! number of days of co2 ramping
   integer  :: yr, mon, day, ncsec ! components of a date
   integer  :: ncdate              ! current date in integer format [yyyymmdd]
!-----------------------------------------------------------------------

   call get_curr_date(yr, mon, day, ncsec)
   ncdate = yr*10000 + mon*100 + day

   call timemgr_datediff(co2_start, 0, ncdate, ncsec, daydiff)

   if (daydiff > 0.0) then

      co2vmr = co2_base*(co2_daily_factor)**daydiff

      if(co2_daily_factor < 1.0) then
         co2vmr = max(co2vmr,co2_limit)
      else
         co2vmr = min(co2vmr,co2_limit)
      end if
   end if

   return
end subroutine ghg_surfvals_set_co2


!=========================================================================================

end module ghg_surfvals