#include <misc.h>
module soxbnd 4,5
!-----------------------------------------------------------------------
! Purpose:
! Interpolate the SOX emissions data from the GEIA-SMITH datasets. This
! dataset contains seasonal average data for various years at roughly 15
! year intervals. The interpolation scheme does linear interpolations to
! the requested calendar day in the seasonal cycle for each of
! the two years in the dataset that bound the requested year. Finally a
! linear interpolation is done to the requested year.
!
! It is assumed that the model calling this interface has been
! compiled so that 8 byte real data are being used. On many
! machines this implies compiling with a "-r8" flag.
!
! Author: B. Eaton
!
! Modified 11 Jan 2000 PJR: work with 4 or 8 byte floats
! Modified 15 Jan 2004 D Bundy: can cycle one year of emission data
! using namelist variable rampyear_prognostic_sulfur
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid, only: masterproc, plon, plat
use ppgrid, only: pcols, begchunk, endchunk
use abortutils, only: endrun
use infnan, only: bigint
!nf90 use netcdf
implicit none
save
private
public :: &
soxbndini, &! initialize soxbnd module
soxbndint, &! interpolate soxbnd data to requested date/time
soxbndget ! return latitude slice data at current date/time
! public module data
! Variables set by namelist
character(len=256) :: soxems_data
!++drb
integer, public :: rampyear_prognostic_sulfur = bigint ! if = bigint => unset, else = YYYY to cycle
character(len=16), public :: scenario_prognostic_sulfur = 'RAMPED' ! only option available
logical :: doramp_sox ! at the moment this is the only option
logical :: fixyear_sox ! if true then cycle on year from the emissions dataset
!--drb
#include <netcdf.inc>
integer, parameter ::&
ndlev=2 ! number of levels in emissions data
integer, allocatable, dimension(:) :: &
date, &! date coordinate (yyyymmdd format)
yr ! years of annual cycle data
real(r8), allocatable, dimension(:) :: &
cdays ! mid-season calendar days (fractional)
real(r8), allocatable, dimension(:,:,:,:) :: &
soxyrlo, &! SOX input data for lower bound year (pcols,begchunk:endchunk,ndlev,ntpy)
soxyrhi ! SOX input data for upper bound year
real(r8), allocatable, dimension(:,:,:) :: &
sox ! SOX interpolated data (pcols,begchunk:endchunk,ndlev)
integer ::&
ntpy, &! number of time samples per year in emissions dataset
nyr, &! number of years in emissions dataset
ncid, &! ID for netCDF file
loyri, &! index in yr array for lower bound year
start(4), &! start vector for netCDF hyperslabs
count(4) ! count vector for netCDF hyperslabs
!##############################################################################
contains
!##############################################################################
subroutine soxbndini( inyear ) 1,32
!-----------------------------------------------------------------------
! Purpose:
! Open netCDF file containing SOX emissions data. Initialize arrays
! with the data to be interpolated to the current time.
!
! It is assumed that each year for which data is available contains
! the same number of time samples.
!
! It is assumed that the input data contains a "date" coordinate that is
! increasing and represents dates in the format yyyymmdd.
!
! Author: B. Eaton
!-----------------------------------------------------------------------
use error_messages, only: alloc_err, handle_ncerr, handle_err
use phys_grid, only: scatter_field_to_chunk
use ioFileMod, only: getfil
use filenames, only: bndtvsox
#if ( defined SPMD )
use mpishorthand
#endif
#ifdef MATCH
use calendar
#endif
!-----------------------------------------------------------------------
#include <netcdf.inc>
!-----------------------------------------------------------------------
integer, intent(in) ::&
inyear ! current year
! Local variables:
integer year ! year that is used (may be inyear or rampyear_prognostic_sulfur if fixed)
integer ::&
i, istat, did, nlon, vid, recid, nrec
real(r8), allocatable :: soxlatlonhi(:,:,:,:) ! used for netCDF input
real(r8), allocatable :: soxlatlonlo(:,:,:,:) ! used for netCDF input
character*10 errmsg ! error messages
! Externals
! real(r8) ::&
! caldayr
!-----------------------------------------------------------------------
! Get file name.
call getfil(bndtvsox, soxems_data, 0)
!++drb Set logical flags according to namelist variables
! Currently the only scenario available is RAMPED
if ( scenario_prognostic_sulfur == 'RAMPED' ) then
doramp_sox = .true.
if ( masterproc ) &
write(6,*)'soxbndini: scenario_prognostic_sulfur = RAMPED'
else if ( scenario_prognostic_sulfur == 'FIXED' ) then
call endrun('soxbndini: input namelist SCENARIO_SOX must be set to RAMPED')
else
call endrun('soxbndini: input namelist SCENARIO_SOX must be set to RAMPED')
end if
! Determine if namelist variable to cycle year has been set
if ( rampyear_prognostic_sulfur .ne. bigint ) then
fixyear_sox = .true.
if ( masterproc ) &
write(6,*) 'soxbndini: cycling emissions from year ',rampyear_prognostic_sulfur
else
fixyear_sox = .false.
end if
! Set year to get out of file to either year from model (inyear) or fixed (rampyear_prognostic_sulfur)
if ( .not. fixyear_sox ) then
year = inyear
else
if ( masterproc ) &
write(6,*)'soxbndini: WARNING- overriding model year with rampyear_prognostic_sulfur',rampyear_prognostic_sulfur
year = rampyear_prognostic_sulfur
endif
!--drb fixed year
if (masterproc) then
write(6,*)'soxbndini: called for model year:',year
! Open file.
!nf90 call handle_ncerr( nf90_open( soxems_data, NF_NOWRITE, ncid ), &
!nf90 'soxbndini: error opening file '//trim(soxems_data) )
call handle_ncerr( nf_open( soxems_data, NF_NOWRITE, ncid ), &
'soxbndini: error opening file '//trim(soxems_data) )
! get the record id
!nf90 call handle_ncerr( nf90_inquire( ncid, unlimiteddimid=recid), &
!nf90 'soxbndini: no record variables ' )
call handle_ncerr( nf_inq_unlimdim( ncid, recid), &
'soxbndini: no record variables ' )
! Check that input data is a right resolution.
!nf90 call handle_ncerr( nf90_inq_dimid( ncid, 'lon', did ), 'soxbndini: ' )
!nf90 call handle_ncerr( nf90_inquire_dimension( ncid, did, len=nlon ), 'soxbndini: ' )
call handle_ncerr( nf_inq_dimid( ncid, 'lon', did ), 'soxbndini: ' )
call handle_ncerr( nf_inq_dimlen( ncid, did, nlon ), 'soxbndini: ' )
if ( nlon .ne. plon ) then
write(unit=errmsg,fmt='(i10)') nlon
call endrun('soxbndini: number of longitudes ('//trim(errmsg)//') doesn''t match model resolution.')
end if
! Get size of unlimited dimension.
!nf90 call handle_ncerr( nf90_inquire_dimension( ncid, recid, len=nrec ), 'soxbndini: ' )
call handle_ncerr( nf_inq_dimlen( ncid, recid, nrec ), 'soxbndini: ' )
endif ! end of masterproc
#if ( defined SPMD )
! broadcast nrec to nodes
call mpibcast (nrec, 1, mpiint, 0, mpicom)
#endif
allocate( date(nrec), stat=istat )
call alloc_err( istat, 'soxbndini', 'date', nrec )
! Get date coordinate.
if (masterproc) then
!nf90 call handle_ncerr( nf90_inq_varid( ncid, 'date', vid ), &
!nf90 'soxbndini: cannot find date coordinate variable' )
!nf90 call handle_ncerr( nf90_get_var( ncid, vid, date ), &
!nf90 'soxbndini: error getting date coordinate data' )
call handle_ncerr( nf_inq_varid( ncid, 'date', vid ), &
'soxbndini: cannot find date coordinate variable' )
call handle_ncerr( nf_get_var_int( ncid, vid, date ), &
'soxbndini: error getting date coordinate data' )
endif ! end of masterproc
#if ( defined SPMD )
! broadcast date to nodes
call mpibcast (date, nrec, mpiint, 0, mpicom)
#endif
! Determine number of time samples per year.
ntpy = 1
do i = 2, nrec
if ( date(i)/10000 .eq. date(1)/10000 ) ntpy = ntpy + 1
end do
! Construct the years array.
nyr = nrec/ntpy
allocate( yr(nyr), stat=istat )
call alloc_err( istat, 'soxbndini', 'yr', nyr )
do i = 1, nyr
yr(i) = date((i-1)*ntpy+1)/10000
end do
if ( masterproc ) &
write(6,*)'soxbndini: years in emission dataset:',(yr(i),i=1,nyr)
! Construct array of calendar days for the annual cycle.
allocate( cdays(ntpy), stat=istat )
call alloc_err( istat, 'soxbndini', 'cdays', ntpy )
do i = 1, ntpy
#ifdef MATCH
cdays(i) = caldayr( date(i), 0 ) ! match version
#else
call bnddyi( date(i), 0, cdays(i) ) ! ccm version
#endif
end do
if ( masterproc ) &
write(6,*)'soxbndini: calendar days in annual cycle:', (cdays(i),i=1,ntpy)
if ( masterproc ) &
write (6,*) 'soxbndini: searching for emissions for year ', year
! Check that requested year is contained in the data.
if ( year .lt. yr(1) .or. year .gt. yr(nyr) ) then
call endrun ('SOXBNDINI: requested year outside data limits in '//trim(soxems_data))
end if
! Find index for the data year that is the lower bound of the
! interval containing the input year.
do i = 1, nyr
if ( yr(i) .gt. year ) then
loyri = i - 1
exit
end if
end do
allocate( soxyrlo(pcols,begchunk:endchunk,ndlev,ntpy), stat=istat )
call alloc_err( istat, 'soxbndini', 'soxyrlo', pcols*((endchunk-begchunk)+1)*ndlev*ntpy )
allocate( soxyrhi(pcols,begchunk:endchunk,ndlev,ntpy), stat=istat )
call alloc_err( istat, 'soxbndini', 'soxyrhi', pcols*((endchunk-begchunk)+1)*ndlev*ntpy )
allocate( sox(pcols,begchunk:endchunk,ndlev), stat=istat )
call alloc_err( istat, 'soxbndini', 'sox', pcols*((endchunk-begchunk)+1)*ndlev )
soxyrlo(:,:,:,:) = 0._r8
soxyrhi(:,:,:,:) = 0._r8
sox(:,:,:) = 0._r8
allocate( soxlatlonlo(plon,plat,ndlev,ntpy), &
soxlatlonhi(plon,plat,ndlev,ntpy), stat=istat )
call alloc_err( istat, 'soxbndini', 'soxlatlonlo and soxlatlonhi', &
plon*plat*ndlev*ntpy*2 )
if (masterproc) then
! Read SOx data for years surrounding initial year.
start(1) = 1
start(2) = 1
start(3) = 1
start(4) = 1
count(1) = plon
count(2) = plat
count(3) = ndlev
count(4) = ntpy
!nf90 call handle_ncerr( nf90_inq_varid( ncid, 'SOx', vid ), &
!nf90 'soxbndini: cannot find variable '//'SOx' )
call handle_ncerr( nf_inq_varid( ncid, 'SOx', vid ), &
'soxbndini: cannot find variable '//'SOx' )
start(4) = (loyri-1)*ntpy + 1
!nf90 call handle_ncerr( nf90_get_var( ncid, vid, soxlatlonlo, start, count ), &
!nf90 'soxbndini: cannot read data for '//'SOx' )
call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, soxlatlonlo ), &
'soxbndini: cannot read data for '//'SOx' )
start(4) = start(4) + ntpy
!nf90 call handle_ncerr( nf90_get_var( ncid, vid, soxlatlonhi, start, count ), &
!nf90 'soxbndini: cannot read data for '//'SOx' )
call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, soxlatlonhi ), &
'soxbndini: cannot read data for '//'SOx' )
endif ! end of masterproc
! scatter data to chunked data structures
call scatter_field_to_chunk (1, 1, ndlev*ntpy, plon, soxlatlonlo, &
soxyrlo(1,begchunk,1,1))
call scatter_field_to_chunk (1, 1, ndlev*ntpy, plon, soxlatlonhi, &
soxyrhi(1,begchunk,1,1))
if ( masterproc ) then
write(6,*)'soxbndini: read data for years; ', yr(loyri),' and ',yr(loyri+1)
endif ! end of masterproc
deallocate( soxlatlonlo, soxlatlonhi, stat=istat )
call handle_err( istat, &
'ERROR deallocating memory for soxlatlonlo and soxlatlonhi in routine soxbndini')
end subroutine soxbndini
!#######################################################################
subroutine soxbndint( inyear, calday ) 1,14
!-----------------------------------------------------------------------
! Purpose:
! Interpolate SOX data to the current time. Update the input data
! as necessary.
!
! Author: B. Eaton
!-----------------------------------------------------------------------
use error_messages, only: alloc_err, handle_ncerr, handle_err
use phys_grid, only: scatter_field_to_chunk
#if ( defined SPMD )
use mpishorthand
#endif
!-----------------------------------------------------------------------
#include <netcdf.inc>
!-----------------------------------------------------------------------
integer, intent(in) ::&
inyear ! current year
real(r8), intent(in) ::&
calday ! calendar day (w/ fraction) in current year, range = [1.,366.)
! Local variables:
integer year ! year that is used (may be inyear or rampyear_prognostic_sulfur if fixed)
integer ::&
i, j, k, n, &
vid, &
lotim, hitim
real(r8) :: dt, dt1, tint
real(r8), allocatable :: anncyclo(:,:,:) ! low year data interpolated to calday
real(r8), allocatable :: anncychi(:,:,:) ! high year data interpolated to calday
! TBH: Get rid of all the replicated code for I/O and scattering!
real(r8), allocatable :: soxlatlonhi(:,:,:,:) ! used for netCDF input
integer :: istat
!-----------------------------------------------------------------------
allocate( anncyclo(pcols,begchunk:endchunk,ndlev), &
anncychi(pcols,begchunk:endchunk,ndlev), stat=istat )
call alloc_err( istat, 'soxbndint', 'anncyclo and anncychi', &
pcols*((endchunk-begchunk)+1)*ndlev*2 )
anncyclo(:,:,:) = 0._r8
anncychi(:,:,:) = 0._r8
!++drb fixed year
if ( .not. fixyear_sox ) then
year = inyear
else
year = rampyear_prognostic_sulfur
endif
!--drb fixed year
! Check to see if model year is still bounded by dataset years.
if ( year .gt. yr(nyr) ) then
write(6,*)'soxbndint: requested year = ',year, ' last dataset year = ',yr(nyr)
call endrun
end if
if ( year .gt. yr(loyri+1) ) then
loyri = loyri + 1
soxyrlo = soxyrhi
allocate( soxlatlonhi(plon,plat,ndlev,ntpy), stat=istat )
call alloc_err( istat, 'soxbndini', 'soxlatlonhi', &
plon*plat*ndlev*ntpy )
! Read in new soxyrhi data. Replace old soxyrlo data.
if (masterproc) then
start(4) = start(4) + ntpy
!nf90 call handle_ncerr( nf90_inq_varid( ncid, 'SOx', vid ), &
!nf90 'soxbndint: cannot find variable '//'SOx' )
!nf90 call handle_ncerr( nf90_get_var( ncid, vid, soxlatlonhi, start, count ), &
!nf90 'soxbndint: cannot read data for '//'SOx' )
call handle_ncerr( nf_inq_varid( ncid, 'SOx', vid ), &
'soxbndint: cannot find variable '//'SOx' )
call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, soxlatlonhi ), &
'soxbndint: cannot read data for '//'SOx' )
endif ! end of masterproc
! scatter data to chunked data structures
call scatter_field_to_chunk (1, 1, ndlev*ntpy, plon, soxlatlonhi, &
soxyrhi(1,begchunk,1,1))
if ( masterproc ) then
write(6,*)'soxbndint: read data for year; ',yr(loyri+1)
endif ! end of masterproc
deallocate( soxlatlonhi, stat=istat )
call handle_err( istat, &
'ERROR deallocating memory for soxlatlonhi in routine soxbndint')
end if
! Linear interpolation... Start by computing the number of days between
! the lower and upper bounds, and days between
! the model time and lower bound.
if ( ntpy .gt. 1 ) then
call findplb( cdays, ntpy, calday, lotim )
hitim = mod( lotim, ntpy ) + 1
if( cdays(hitim) .lt. cdays(lotim) )then
dt = 365. - cdays(lotim) + cdays(hitim)
if( calday .le. cdays(hitim) )then
dt1 = 365. - cdays(lotim) + calday
else
dt1 = calday - cdays(lotim)
end if
else
dt = cdays(hitim) - cdays(lotim)
dt1 = calday - cdays(lotim)
end if
tint = dt1/dt
! Annual cycle interpolations.
!TBH call linintp( size(anncyclo), 0._r8, 1._r8, tint, soxyrlo(1,1,1,lotim), &
!TBH soxyrlo(1,1,1,hitim), anncyclo )
!TBH call linintp( size(anncychi), 0._r8, 1._r8, tint, soxyrhi(1,1,1,lotim), &
!TBH soxyrhi(1,1,1,hitim), anncychi )
call linintp( size(anncyclo), 0._r8, 1._r8, tint, &
soxyrlo(1,begchunk,1,lotim), &
soxyrlo(1,begchunk,1,hitim), &
anncyclo )
call linintp( size(anncychi), 0._r8, 1._r8, tint, &
soxyrhi(1,begchunk,1,lotim), &
soxyrhi(1,begchunk,1,hitim), &
anncychi )
else
anncyclo(:,:,:) = soxyrlo(:,:,:,1)
anncychi(:,:,:) = soxyrhi(:,:,:,1)
end if
! Interpolate between years for which annual cycle data is present
dt = yr(loyri+1) - yr(loyri)
dt1 = year - yr(loyri)
tint = dt1/dt
call linintp( size(sox), 0._r8, 1._r8, tint, anncyclo, anncychi, sox )
! deallocate memory
deallocate( anncyclo, anncychi, stat=istat )
call handle_err( istat, &
'ERROR deallocating memory for anncyclo and anncychi in routine soxbndint')
end subroutine soxbndint
!#######################################################################
subroutine soxbndget( ncol, lchnk, x ) 1
!-----------------------------------------------------------------------
! Purpose:
! Return SOX emission data for the requested latitude.
!
! Author: B. Eaton
!-----------------------------------------------------------------------
integer, intent(in) :: ncol ! number of columns used
integer, intent(in) :: lchnk ! chunk index
real(r8), intent(out) :: x(pcols,ndlev) ! SOx emissions in Tg S/m2/s
! Local variables:
integer :: i
!-----------------------------------------------------------------------
do i = 1, ncol
x(i,1) = sox(i,lchnk,1)
x(i,2) = sox(i,lchnk,2)
end do
end subroutine soxbndget
!#######################################################################
end module soxbnd