module chem_surfvals 8,6
!-----------------------------------------------------------------------------------
! 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 spmd_utils, only: masterproc
use time_manager, only: get_curr_date, get_start_date, is_end_curr_day, &
timemgr_datediff, get_curr_calday
use abortutils, only: endrun
use netcdf
use error_messages, only: handle_ncerr
use cam_logfile, only: iulog
!-----------------------------------------------------------------------
!- module boilerplate --------------------------------------------------
!-----------------------------------------------------------------------
implicit none
private ! Make default access private
save
! Public methods
public ::&
chem_surfvals_readnl, &! read namelist input
chem_surfvals_init, &! initialize options that depend on namelist input
chem_surfvals_set, &! set ghg surface values when ramp option is on
chem_surfvals_get, &! return surface values for: CO2VMR, CO2MMR, CH4VMR
! N2OVMR, F11VMR, and F12VMR
chem_surfvals_co2_rad ! return co2 for radiation
! Private module data
! Default values for namelist variables -- now set by build-namelist
real(r8) :: o2mmr = .23143_r8 ! o2 mass mixing ratio
real(r8) :: co2vmr_rad = -1.0_r8 ! co2 vmr override for radiation
real(r8) :: co2vmr = -1.0_r8 ! co2 volume mixing ratio
real(r8) :: n2ovmr = -1.0_r8 ! n2o volume mixing ratio
real(r8) :: ch4vmr = -1.0_r8 ! ch4 volume mixing ratio
real(r8) :: f11vmr = -1.0_r8 ! cfc11 volume mixing ratio
real(r8) :: f12vmr = -1.0_r8 ! cfc12 volume mixing ratio
character(len=16) :: scenario_ghg = 'FIXED' ! 'FIXED','RAMPED' or 'RAMP_CO2_ONLY'
integer :: rampYear_ghg = 0 ! ramped gases fixed at this year (if > 0)
character(len=256) :: bndtvghg ! filename for ramped data
integer :: ramp_co2_start_ymd = 0 ! start date for co2 ramping (yyyymmdd)
real(r8) :: ramp_co2_annual_rate = 1.0_r8 ! % amount of co2 ramping per yr; default is 1%
real(r8) :: ramp_co2_cap = -9999.0_r8 ! co2 ramp cap if rate>0, floor otherwise
! as multiple or fraction of inital value
! ex. 4.0 => cap at 4x initial co2 setting
integer :: ghg_yearStart_model = 0 ! model start year
integer :: ghg_yearStart_data = 0 ! data start year
logical :: ghg_use_calendar ! true => data year = model year
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 chem_surfvals_readnl(nlfile) 1,21
! Read chem_surfvals_nl namelist group.
use namelist_utils, only: find_group_name
use units, only: getunit, freeunit
use mpishorthand
character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
! Local variables
integer :: unitn, ierr, i
character(len=*), parameter :: subname = 'chem_surfvals_readnl'
namelist /chem_surfvals_nl/ co2vmr, n2ovmr, ch4vmr, f11vmr, f12vmr, &
co2vmr_rad, scenario_ghg, rampyear_ghg, bndtvghg, &
ramp_co2_start_ymd, ramp_co2_annual_rate, ramp_co2_cap, &
ghg_yearStart_model, ghg_yearStart_data
!-----------------------------------------------------------------------------
if (masterproc) then
unitn = getunit()
open( unitn, file=trim(nlfile), status='old' )
call find_group_name(unitn, 'chem_surfvals_nl', status=ierr)
if (ierr == 0) then
read(unitn, chem_surfvals_nl, iostat=ierr)
if (ierr /= 0) then
call endrun(subname // ':: ERROR reading namelist')
end if
end if
close(unitn)
call freeunit(unitn)
end if
#ifdef SPMD
! Broadcast namelist variables
call mpibcast (co2vmr, 1, mpir8, 0, mpicom)
call mpibcast (n2ovmr, 1, mpir8, 0, mpicom)
call mpibcast (ch4vmr, 1, mpir8, 0, mpicom)
call mpibcast (f11vmr, 1, mpir8, 0, mpicom)
call mpibcast (f12vmr, 1, mpir8, 0, mpicom)
call mpibcast (co2vmr_rad, 1, mpir8, 0, mpicom)
call mpibcast (scenario_ghg, len(scenario_ghg), mpichar, 0, mpicom)
call mpibcast (rampyear_ghg, 1, mpiint, 0, mpicom)
call mpibcast (bndtvghg, len(bndtvghg), mpichar, 0, mpicom)
call mpibcast (ramp_co2_start_ymd, 1, mpiint, 0, mpicom)
call mpibcast (ramp_co2_annual_rate, 1, mpir8, 0, mpicom)
call mpibcast (ramp_co2_cap, 1, mpir8, 0, mpicom)
call mpibcast (ghg_yearstart_model, 1, mpiint, 0, mpicom)
call mpibcast (ghg_yearstart_data, 1, mpiint, 0, mpicom)
#endif
if (co2vmr_rad > 0._r8) then
if (masterproc) &
write(iulog,*) trim(subname)//': co2vmr_rad override is set to ', co2vmr_rad
end if
end subroutine chem_surfvals_readnl
!================================================================================================
subroutine chem_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 chem_surfvals_set calls time manager methods.
!
! Author: B. Eaton - merged code from parse_namelist and rampnl_ghg.
!
!-----------------------------------------------------------------------
use infnan, only: inf
!---------------------------Local variables-----------------------------
integer :: yr, mon, day, ncsec
!-----------------------------------------------------------------------
if (scenario_ghg == 'FIXED') then
doRamp_ghg = .false.
ramp_just_co2 = .false.
if (masterproc) &
write(iulog,*)'chem_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(iulog,*) ' FIXED values from year ',fixYear_ghg
else
write(iulog,*) ' RAMPED values initialized to'
end if
end if
call chem_surfvals_set()
else if (scenario_ghg == 'RAMP_CO2_ONLY') then
if(ramp_co2_start_ymd == 0) then
! by default start the ramp at the initial run time
call get_start_date(yr, mon, day, ncsec)
ramp_co2_start_ymd = yr*10000 + mon*100 + day
end if
co2_start = ramp_co2_start_ymd
if(ramp_co2_annual_rate <= -100.0_r8) then
write(iulog,*) 'RAMP_CO2: invalid ramp_co2_annual_rate= ',ramp_co2_annual_rate
call endrun ('chem_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(iulog,*) ' RAMPED values initialized to'
co2_daily_factor = (ramp_co2_annual_rate*0.01_r8+1.0_r8)**(1.0_r8/365.0_r8)
if(ramp_co2_cap > 0.0_r8) then
co2_limit = ramp_co2_cap * co2_base
else ! if no cap/floor specified, provide default
if(ramp_co2_annual_rate < 0.0_r8) then
co2_limit = 0.0_r8
else
co2_limit = inf
end if
end if
if((ramp_co2_annual_rate<0.0_r8 .and. co2_limit>co2_base) .or. &
(ramp_co2_annual_rate>0.0_r8 .and. co2_limit<co2_base)) then
write(iulog,*) 'RAMP_CO2: ramp_co2_cap is unreachable'
write(iulog,*) 'RAMP_CO2: ramp_co2_annual_rate= ',ramp_co2_annual_rate,' ramp_co2_cap= ',ramp_co2_cap
call endrun('chem_surfvals_init: ramp_co2_annual_rate and ramp_co2_cap incompatible')
end if
call chem_surfvals_set()
else
call endrun ('chem_surfvals_init: input namelist SCENARIO_GHG must be set to either FIXED, RAMPED or RAMP_CO2_ONLY')
endif
if (masterproc) then
write(iulog,*) ' co2 volume mixing ratio = ',co2vmr
write(iulog,*) ' ch4 volume mixing ratio = ',ch4vmr
write(iulog,*) ' n2o volume mixing ratio = ',n2ovmr
write(iulog,*) ' f11 volume mixing ratio = ',f11vmr
write(iulog,*) ' f12 volume mixing ratio = ',f12vmr
end if
end subroutine chem_surfvals_init
!=========================================================================================
subroutine ghg_ramp_read() 1,30
!-----------------------------------------------------------------------
!
! Purpose:
! Read ramped greenhouse gas surface data.
!
! Author: T. Henderson
!
!-----------------------------------------------------------------------
use ioFileMod, only: getfil
#if ( defined SPMD )
use mpishorthand, only: mpicom, mpiint, mpir8
#endif
!---------------------------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 handle_ncerr( nf90_open (trim(locfn), NF90_NOWRITE, ncid),&
'chem_surfvals.F90:351')
write(iulog,*)'GHG_RAMP_READ: reading ramped greenhouse gas surface data from file ',trim(locfn)
call handle_ncerr( nf90_inq_varid( ncid, 'date', date_id ),&
'chem_surfvals.F90:354')
call handle_ncerr( nf90_inq_varid( ncid, 'CO2', co2_id ),&
'chem_surfvals.F90:356')
call handle_ncerr( nf90_inq_varid( ncid, 'CH4', ch4_id ),&
'chem_surfvals.F90:358')
call handle_ncerr( nf90_inq_varid( ncid, 'N2O', n2o_id ),&
'chem_surfvals.F90:360')
call handle_ncerr( nf90_inq_varid( ncid, 'f11', f11_id ),&
'chem_surfvals.F90:362')
call handle_ncerr( nf90_inq_varid( ncid, 'f12', f12_id ),&
'chem_surfvals.F90:364')
call handle_ncerr( nf90_inq_varid( ncid, 'adj', adj_id ),&
'chem_surfvals.F90:366')
call handle_ncerr( nf90_inq_dimid( ncid, 'time', time_id ),&
'chem_surfvals.F90:368')
call handle_ncerr( nf90_inquire_dimension( ncid, time_id, len=ntim ),&
'chem_surfvals.F90:370')
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(iulog,*)'GHG_RAMP_READ: ERROR, allocate() failed!'
call endrun
endif
if (masterproc) then
call handle_ncerr( nf90_get_var (ncid, date_id, yrdata ),&
'chem_surfvals.F90:384')
yrdata = yrdata / 10000
call handle_ncerr( nf90_get_var (ncid, co2_id, co2 ),&
'chem_surfvals.F90:387')
call handle_ncerr( nf90_get_var (ncid, ch4_id, ch4 ),&
'chem_surfvals.F90:389')
call handle_ncerr( nf90_get_var (ncid, n2o_id, n2o ),&
'chem_surfvals.F90:391')
call handle_ncerr( nf90_get_var (ncid, f11_id, f11 ),&
'chem_surfvals.F90:393')
call handle_ncerr( nf90_get_var (ncid, f12_id, f12 ),&
'chem_surfvals.F90:395')
call handle_ncerr( nf90_get_var (ncid, adj_id, adj ),&
'chem_surfvals.F90:397')
call handle_ncerr( nf90_close (ncid),&
'chem_surfvals.F90:399')
write(iulog,*)'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 chem_surfvals_get(name) 14,2
use physconst, only: mwdry, mwco2
character(len=*), intent(in) :: name
real(r8) :: rmwco2
real(r8) :: chem_surfvals_get
rmwco2 = mwco2/mwdry ! ratio of molecular weights of co2 to dry air
select case (name)
case ('CO2VMR')
chem_surfvals_get = co2vmr
case ('CO2MMR')
chem_surfvals_get = rmwco2 * co2vmr
case ('N2OVMR')
chem_surfvals_get = n2ovmr
case ('CH4VMR')
chem_surfvals_get = ch4vmr
case ('F11VMR')
chem_surfvals_get = f11vmr
case ('F12VMR')
chem_surfvals_get = f12vmr
case ('O2MMR')
chem_surfvals_get = o2mmr
case default
call endrun('chem_surfvals_get does not know name')
end select
end function chem_surfvals_get
!=========================================================================================
function chem_surfvals_co2_rad(vmr_in) 2,1
! Return the value of CO2 (as mmr) that is radiatively active.
! This method is used by ghg_data to set the prescribed value of CO2 in
! the physics buffer. If the user has set the co2vmr_rad namelist
! variable then that value will override either the value set by the
! co2vmr namelist variable, or the values time interpolated from a
! dataset.
! This method is also used by cam_history to write the radiatively active
! CO2 to the history file. The optional argument allows returning the
! value as vmr.
use physconst, only: mwdry, mwco2
! Arguments
logical, intent(in), optional :: vmr_in ! return CO2 as vmr
! Return value
real(r8) :: chem_surfvals_co2_rad
! Local variables
real(r8) :: convert_vmr ! convert vmr to desired output
!-----------------------------------------------------------------------
! by default convert vmr to mmr
convert_vmr = mwco2/mwdry ! ratio of molecular weights of co2 to dry air
if (present(vmr_in)) then
! if request return vmr
if (vmr_in) convert_vmr = 1.0_r8
end if
if (co2vmr_rad > 0._r8) then
chem_surfvals_co2_rad = convert_vmr * co2vmr_rad
else
chem_surfvals_co2_rad = convert_vmr * co2vmr
end if
end function chem_surfvals_co2_rad
!=========================================================================================
subroutine chem_surfvals_set( phys_state ) 3,5
! phys_state argument is unused in this version
use ppgrid, only: begchunk, endchunk
use physics_types, only: physics_state
type(physics_state), intent(inout), dimension(begchunk:endchunk), optional :: phys_state
!---------------------------Local variables-----------------------------
integer :: yr, mon, day, ncsec ! components of a date
integer :: ncdate ! current date in integer format [yyyymmdd]
if (.not. doRamp_ghg) return
if(ramp_just_co2) then
call chem_surfvals_set_co2()
else
call chem_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(iulog,*) 'chem_surfvals_set: ncdate= ',ncdate,' co2vmr=',co2vmr
if (.not. ramp_just_co2 .and. mon==1 .and. day==1) then
write(iulog,*) 'chem_surfvals_set: ch4vmr=', ch4vmr, ' n2ovmr=', n2ovmr, &
' f11vmr=', f11vmr, ' f12vmr=', f12vmr
end if
end if
return
end subroutine chem_surfvals_set
!=========================================================================================
subroutine chem_surfvals_set_all() 1,7
!-----------------------------------------------------------------------
!
! Purpose:
! Computes greenhouse gas volume mixing ratios via interpolation of
! yearly input data.
!
! Author: B. Eaton - updated ramp_ghg for use in chem_surfvals module
!
!-----------------------------------------------------------------------
use interpolate_data, only: get_timeinterp_factors
!---------------------------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
integer yearRan_model ! model ran year
!
! ---------------------------------------------------------------------
!
calday = get_curr_calday()
call get_curr_date(yr, mon, day, ncsec)
ncdate = yr*10000 + mon*100 + day
!
! determine ghg_use_calendar
!
if ( ghg_yearStart_model > 0 .and. ghg_yearStart_data > 0 ) then
ghg_use_calendar = .false.
else
ghg_use_calendar = .true.
end if
!
! determine index into input data
!
if ( fixYear_ghg > 0) then
yrmodel = fixYear_ghg
nyrm = fixYear_ghg - yrdata(1) + 1
else
if ( ghg_use_calendar) then
yrmodel = yr
nyrm = yr - yrdata(1) + 1
else
yearRan_model = yr - ghg_yearStart_model
if ( yearRan_model < 0 ) then
call endrun('chem_surfvals_set_all: incorrect ghg_yearStart_model')
endif
yrmodel = yearRan_model + ghg_yearStart_data
nyrm = ghg_yearStart_data + yearRan_model - yrdata(1) + 1
end if
end if
nyrp = nyrm + 1
!
! if current date is before yrdata(1), quit
!
if (nyrm < 1) then
write(iulog,*)'chem_surfvals_set_all: data time index is out of bounds'
write(iulog,*)'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 ('chem_surfvals_set_all: error - current date is past the end of valid data')
! write(iulog,*)'chem_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._r8 + calday
doydatam = yrdata(nyrm)*365._r8 + 1._r8
doydatap = yrdata(nyrp)*365._r8 + 1._r8
call get_timeinterp_factors(.false.,2,doydatam,doydatap, doymodel, &
fact1, fact2,'chem_surfvals')
!
! do time interpolation:
! co2 in ppmv
! n2o,ch4 in ppbv
! f11,f12 in pptv
!
co2vmr = (co2(nyrm)*fact1 + co2(nyrp)*fact2)*1.e-06_r8
ch4vmr = (ch4(nyrm)*fact1 + ch4(nyrp)*fact2)*1.e-09_r8
n2ovmr = (n2o(nyrm)*fact1 + n2o(nyrp)*fact2)*1.e-09_r8
cfcscl = (adj(nyrm)*fact1 + adj(nyrp)*fact2)
f11vmr = (f11(nyrm)*fact1 + f11(nyrp)*fact2)*1.e-12_r8*(1._r8+cfcscl)
f12vmr = (f12(nyrm)*fact1 + f12(nyrp)*fact2)*1.e-12_r8
return
end subroutine chem_surfvals_set_all
!=========================================================================================
subroutine chem_surfvals_set_co2() 1,3
!-----------------------------------------------------------------------
!
! 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 chem_surfvals module
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
!---------------------------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_r8) then
co2vmr = co2_base*(co2_daily_factor)**daydiff
if(co2_daily_factor < 1.0_r8) then
co2vmr = max(co2vmr,co2_limit)
else
co2vmr = min(co2vmr,co2_limit)
end if
end if
return
end subroutine chem_surfvals_set_co2
!=========================================================================================
end module chem_surfvals