module sulbnd 1,3

   !----------------------------------------------------------------------- 
   ! Purpose: 
   ! Interpolate the GEIA sulfur emissions data (SO2, SO4, DMS) from the
   ! dataset prepared by Mary Barth.  This dataset contains seasonal average
   ! data for 1985.
   !
   ! 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
   !----------------------------------------------------------------------- 
   use shr_kind_mod, only: r8 => shr_kind_r8
   use pmgrid
   use abortutils, only: endrun
!nf90   use netcdf

   implicit none
   save
   private
   public :: &
      sulbndini,      &! initialize sulbnd module
      sulbndint,      &! interpolate sulbnd data to requested date/time
      sulbndget        ! return latitude slice data at current date/time

   ! public module data
   
   public :: sulems_data ! full pathname for time-variant sulfer emissions dataset
   character(len=256) :: sulems_data='sulfur_emissions.nc'

   ! private module data

#include <netcdf.inc>

   integer, parameter ::&
      ndlev=2          ! number of levels in emissions data

   real(r8), allocatable, dimension(:) :: &
      time             ! time coordinate (calander days + frac)

   real(r8), dimension(plon,plat,ndlev,2) :: &
      so2in,          &! SO2 input data.  Upper and lower time bounds.
      so4in,          &! SO4 input data.  Upper and lower time bounds.
      dmsin            ! DMS input data.  Upper and lower time bounds.

   real(r8), dimension(plon,plat,ndlev) :: &
      so2,            &! SO2 interpolated data.
      so4,            &! SO4 interpolated data.
      dms              ! DMS interpolated data.

   integer ::&
      ncid,           &! ID for netCDF file
      nrec,           &! number of records (time samples)
      lotim,          &! time(lotim) .le. current time
      hitim,          &! current time .lt. time(hitim)
      loin,           &! index into input data array containing time(lotim) data
      hiin,           &! index into input data array containing time(hitim) data
      start(4),       &! start vector for netCDF hyperslabs
      count(4)         ! count vector for netCDF hyperslabs

!##############################################################################
contains
!##############################################################################


   subroutine sulbndini( calday ),24

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Open netCDF file containing sulfur emissions data.  Initialize arrays
      ! with the data to be interpolated to the current time.
      !
      ! It is assumed that the time coordinate is increasing and represents
      ! calendar days (range = [1.,366.)).
      ! 
      ! Author: B. Eaton
      !----------------------------------------------------------------------- 

      use error_messages, only: alloc_err, handle_ncerr
!nf90      use netcdf

      implicit none
!-----------------------------------------------------------------------
#include <netcdf.inc>
!-----------------------------------------------------------------------

      real(r8), intent(in) ::&
         calday  ! current time in calendar days + fraction.

      ! Local variables:
      integer ::&
         i, istat, did, nlon, vid, recid
      character(len=*), parameter ::&
         routine_name = 'sulbndini'
      !-----------------------------------------------------------------------

      ! Open file.
      write(6,*)'Sulfer emissions dataset is: ', trim(sulems_data)
      call handle_ncerr( nf_open( trim( sulems_data ), NF_NOWRITE, ncid ), &
         routine_name//': error opening file '//trim( sulems_data ) )

      ! Check that input data is a right resolution.
      call handle_ncerr( nf_inq_dimid( ncid, 'lon', did ), routine_name//': ' )
      call handle_ncerr( nf_inq_dimlen( ncid, did, nlon ), routine_name//': ' )
      if ( nlon .ne. plon ) then
         write(6,*)routine_name//': number of longitudes (', nlon, ')', &
                   ' doesn''t match model resolution.'
         call endrun
      end if

      ! Get size of unlimited dimension.
      call handle_ncerr( nf_inq_unlimdim( ncid, recid ), routine_name//': ' )
      call handle_ncerr( nf_inq_dimlen( ncid, recid, nrec ), routine_name//': ') 

      ! Get time coordinate.
      allocate( time(nrec), stat=istat )
      call alloc_err( istat, routine_name, 'time', nrec )
      call handle_ncerr( nf_inq_varid( ncid, 'time', vid ), &
           routine_name//': cannot find time coordinate variable' )
!nf90      call handle_ncerr( nf90_get_var( ncid, vid, time ), &
!nf90         routine_name//': error getting time coordinate data' )
      call handle_ncerr( nf_get_var_double( ncid, vid, time ), &
         routine_name//': error getting time coordinate data' )

      ! Make sure the time coordinate looks like calander day, and is
      ! increasing.
      call chktime( time, nrec )

      ! Find indices for time samples that bound the current time.
      call findplb( time, nrec, calday, lotim )
      hitim = mod( lotim, nrec ) + 1

      ! Read data.
      loin = 1
      hiin = 2
      start(1) = 1
      start(2) = 1
      start(3) = 1
      count(1) = plon
      count(2) = plat
      count(3) = ndlev
      count(4) = 1

      start(4) = lotim
      call handle_ncerr( nf_inq_varid( ncid, 'SO2', vid ), routine_name//': cannot find variable '//'SO2' )
      call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, so2in(1,1,1,loin) ), &
         routine_name//': cannot read data for '//'SO2' )

      call handle_ncerr( nf_inq_varid( ncid, 'SO4', vid ), routine_name//': cannot find variable '//'SO4' )
      call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, so4in(1,1,1,loin) ), &
         routine_name//': cannot read data for '//'SO4' )
      call handle_ncerr( nf_inq_varid( ncid, 'DMS', vid ), routine_name//': cannot find variable '//'DMS' )
      call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, dmsin(1,1,1,loin) ), &
         routine_name//': cannot read data for '//'DMS' )

      start(4) = hitim
      call handle_ncerr( nf_inq_varid( ncid, 'SO2', vid ), routine_name//': cannot find variable '//'SO2' )
      call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, so2in(1,1,1,hiin) ), &
         routine_name//': cannot read data for '//'SO2' )
      call handle_ncerr( nf_inq_varid( ncid, 'SO4', vid ), routine_name//': cannot find variable '//'SO4' )
      call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, so4in(1,1,1,hiin) ), &
         routine_name//': cannot read data for '//'SO4' )
      call handle_ncerr( nf_inq_varid( ncid, 'DMS', vid ), routine_name//': cannot find variable '//'DMS' )
      call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, dmsin(1,1,1,hiin) ), &
         routine_name//': cannot read data for '//'DMS' )

 

      write(6,*) routine_name//': calendar day = ',calday, ' : read data for days ',time(lotim), &
                 ' and ',time(hitim)

   end subroutine sulbndini

!#######################################################################


   subroutine sulbndint( calday ),17

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Interpolate sulfur data to the current time.  Update the input data
      ! as necessary.
      !
      ! Author: B. Eaton
      !----------------------------------------------------------------------- 

      use error_messages, only: handle_ncerr
!nf90      use netcdf

      implicit none
!-----------------------------------------------------------------------
#include <netcdf.inc>
!-----------------------------------------------------------------------

      real(r8), intent(in) ::&
         calday  ! current calendar day (w/ fraction), range = [1.,366.)

      ! Local variables:
      integer ::&
         oldlotim, oldhitim, &
         vid, ret
      real(r8) ::&
         dt, dt1, tint
      character(len=*), parameter ::&
         routine_name = 'sulbndint'
      !-----------------------------------------------------------------------

      ! Check to see if model time is still bounded by dataset times.
      oldlotim = lotim
      oldhitim = hitim
      call findplb( time, nrec, calday, lotim )
      hitim = mod( lotim, nrec ) + 1

      if ( hitim .ne. oldhitim ) then
         ! Read in new hitim data.  Replace old lotim data.
         loin = hiin
         hiin = mod( loin, 2 ) + 1
         start(4) = hitim
         call handle_ncerr( nf_inq_varid( ncid, 'SO2', vid ), routine_name//': cannot find variable '//'SO2' )
         call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, so2in(1,1,1,hiin) ), &
            routine_name//': cannot read data for '//'SO2' )
         call handle_ncerr( nf_inq_varid( ncid, 'SO4', vid ), routine_name//': cannot find variable '//'SO4' )
         call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, so4in(1,1,1,hiin) ), &
            routine_name//': cannot read data for '//'SO4' )
         call handle_ncerr( nf_inq_varid( ncid, 'DMS', vid ), routine_name//': cannot find variable '//'DMS' )
         call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, dmsin(1,1,1,hiin) ), &
            routine_name//': cannot read data for '//'DMS' )

         write(6,*) routine_name//': read data for day ',time(hitim)

         if ( lotim .ne. oldhitim ) then
            ! Read in new lotim data.  Replace old hitim data.
            start(4) = lotim
            call handle_ncerr( nf_inq_varid( ncid, 'SO2', vid ), routine_name//': cannot find variable '//'SO2' )
            call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, so2in(1,1,1,loin) ), &
               routine_name//': cannot read data for '//'SO2' )
            call handle_ncerr( nf_inq_varid( ncid, 'SO4', vid ), routine_name//': cannot find variable '//'SO4' )
            call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, so4in(1,1,1,loin) ), &
               routine_name//': cannot read data for '//'SO4' )
            call handle_ncerr( nf_inq_varid( ncid, 'DMS', vid ), routine_name//': cannot find variable '//'DMS' )
            call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, dmsin(1,1,1,loin) ), &
               routine_name//': cannot read data for '//'DMS' )

            write(6,*) routine_name//': read data for day ',time(lotim)
         end if

      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( time(hitim) .lt. time(lotim) )then
         dt = 365. - time(lotim) + time(hitim)
         if( calday .le. time(hitim) )then
            dt1 = 365. - time(lotim) + calday
         else
            dt1 = calday - time(lotim)
         end if
      else
         dt = time(hitim) - time(lotim)
         dt1 = calday - time(lotim)
      end if
      tint = dt1/dt
      call linintp( plon*plat*2, 0._r8, 1._r8, tint, so2in(1,1,1,loin), so2in(1,1,1,hiin), so2 )
      call linintp( plon*plat*2, 0._r8, 1._r8, tint, so4in(1,1,1,loin), so4in(1,1,1,hiin), so4 )
      call linintp( plon*plat*2, 0._r8, 1._r8, tint, dmsin(1,1,1,loin), dmsin(1,1,1,hiin), dms )

   end subroutine sulbndint

!#######################################################################


   subroutine sulbndget(pcols, ncol, lat, lon, semis ) 1

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Return sulfur emission data for the requested latitude.
      !
      ! Author: B. Eaton
      !----------------------------------------------------------------------- 

      implicit none

      integer, intent(in) ::     pcols
      integer, intent(in) ::     ncol
      integer, intent(in) ::     lat(pcols)             ! requested latitude indices
      integer, intent(in) ::     lon(pcols)             ! requested longitude

      real(r8), intent(out) ::&
         semis(pcols,6)  ! sulfur emissions

      ! Local variables:
      integer :: i
      !-----------------------------------------------------------------------

      do i = 1, ncol
         semis(i,1) = so2(lon(i),lat(i),1)
         semis(i,2) = so4(lon(i),lat(i),1)
         semis(i,3) = dms(lon(i),lat(i),1)
         semis(i,4) = so2(lon(i),lat(i),2)
         semis(i,5) = so4(lon(i),lat(i),2)
         semis(i,6) = dms(lon(i),lat(i),2)
      end do

   end subroutine sulbndget

!#######################################################################

end module sulbnd