#include <misc.h>


module dmsbnd 3,4

   !----------------------------------------------------------------------- 
   ! Purpose: 
   ! This code does time interpolation for DMS boundary data in a netCDF
   ! file.  Assumptions on the data in the netCDF file are:
   ! 1. Coordinates are ordered (lon,lat,time)
   ! 2. The time coordinate is in days, and the data is assumed to be periodic
   !    annually.
   !
   ! 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 ppgrid, only: pcols, begchunk, endchunk
   use pmgrid, only: masterproc, plon, plat
   use abortutils, only: endrun
!nf90   use netcdf

   implicit none
   save
   private
   public :: &
      dmsbndini,      &! initialize dmsbnd module
      dmsbndint,      &! interpolate dmsbnd data to requested date/time
      dmsbndget        ! return latitude slice data at current date/time

   ! public module data
   
   public :: dmsems_data ! full pathname for time-variant DMS emissions dataset
   character(len=256) :: dmsems_data

   ! private module data

#include <netcdf.inc>

   real(r8), allocatable, dimension(:) :: &
      time            ! time coordinate (calander days + frac)
   real(r8), allocatable :: &
      dmsin(:,:,:)    ! input data (pcols,begchunk:endchunk,2)
   real(r8), allocatable :: &
      dms(:,:)        ! interpolated data (pcols,begchunk:endchunk)

   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(3),      &! start vector for netCDF hyperslabs
      count(3)        ! count vector for netCDF hyperslabs

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


   subroutine dmsbndini( calday ) 1,26

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Open netCDF file containing DMS 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, handle_err
      use phys_grid,      only: scatter_field_to_chunk
      use ioFileMod, only: getfil
      use filenames, only: bndtvdms

#if ( defined SPMD )
      use mpishorthand
#endif

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

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

      ! Local variables:
      integer ::&
         did,   &
         istat, &
         recid, &! record ID
         nlon,  &
         vid
      real(r8), allocatable :: dmslatlon(:,:,:)  ! used for netCDF input
      !-----------------------------------------------------------------------

      start(1) = 1
      start(2) = 1
      start(3) = 1
      count(1) = plon
      count(2) = plat
      count(3) = 1

      ! Get file name.  
      call getfil(bndtvdms, dmsems_data, 0)

      ! Open file.
!nf90      call handle_ncerr( nf90_open( trim(dmsems_data), NF_NOWRITE, ncid ), &
!nf90         'dmsbndini: error opening file '//trim(dmsems_data) )

      if (masterproc) then

        call handle_ncerr( nf_open( trim(dmsems_data), NF_NOWRITE, ncid ), &
           'dmsbndini: error opening file '//trim(dmsems_data) )

      ! get the record id
!nf90      call handle_ncerr( nf90_inquire( ncid, unlimiteddimid=recid), &
!nf90         'dmsbndini: no record variables ' )

        call handle_ncerr( nf_inq_unlimdim( ncid, recid), &
           'dmsbndini: no record variables ' )

      ! Check that input data is a right resolution.
!nf90      call handle_ncerr( nf90_inq_dimid( ncid, 'lon', did ), 'dmsbndini: ' )
!nf90      call handle_ncerr( nf90_inquire_dimension( ncid, did, len=nlon ), 'dmsbndini: ' )

        call handle_ncerr( nf_inq_dimid( ncid, 'lon', did ), 'dmsbndini: ' )
        call handle_ncerr( nf_inq_dimlen( ncid, did, nlon ), 'dmsbndini: ' )
        if ( nlon .ne. plon ) then
           write(6,*)'dmsbndini: model plon = ',plon,', dataset nlon = ',nlon
           call endrun('dmsbndini:  plon != nlon')
        end if

      ! Get size of unlimited dimension.
!nf90      call handle_ncerr( nf90_inquire_dimension( ncid, recid, len=nrec ), 'dmsbndini: ' )
        call handle_ncerr( nf_inq_dimlen( ncid, recid, nrec ), 'dmsbndini: ' )

      endif  ! end of masterproc

#if ( defined SPMD )
      ! broadcast nrec to nodes
      call mpibcast (nrec, 1, mpiint, 0, mpicom)
#endif

      ! Allocate space for time coordinate data.
      allocate( time(nrec), stat=istat )
      call alloc_err( istat, 'dmsbndini', 'time', nrec )

      if (masterproc) then

      ! Get time coordinate.
!nf90      call handle_ncerr( nf90_inq_varid( ncid, 'time', vid ), &
!nf90         'dmsbndini: cannot find time coordinate variable' )
!nf90      call handle_ncerr( nf90_get_var( ncid, vid, time ), &
!nf90         'dmsbndini: error getting time coordinate data' )

        call handle_ncerr( nf_inq_varid( ncid, 'time', vid ), &
           'dmsbndini: cannot find time coordinate variable' )
        call handle_ncerr( nf_get_var_double( ncid, vid, time ), &
           'dmsbndini: error getting time coordinate data' )

      endif  ! end of masterproc

#if ( defined SPMD )
      ! broadcast time to nodes
      call mpibcast (time, nrec, mpir8, 0, mpicom)
#endif

      ! 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

      allocate( dmsin(pcols,begchunk:endchunk,2), stat=istat )
      call alloc_err( istat, 'dmsbndini', 'dmsin', pcols*((endchunk-begchunk)+1)*2 )
      allocate( dms(pcols,begchunk:endchunk), stat=istat )
      call alloc_err( istat, 'dmsbndini', 'dms', pcols*((endchunk-begchunk)+1) )

      allocate( dmslatlon(plon,plat,2), stat=istat )
      call alloc_err( istat, 'dmsbndini', 'dmslatlon', &
                      plon*plat*2 )

      if (masterproc) then

!nf90      call handle_ncerr( nf90_inq_varid( ncid, 'DMS', vid ), &
!nf90         'dmsbndini: cannot find variable '//'DMS' )
        call handle_ncerr( nf_inq_varid( ncid, 'DMS', vid ), &
           'dmsbndini: cannot find variable '//'DMS' )

        start(3) = lotim
!nf90      call handle_ncerr( nf90_get_var( ncid, vid,  dmslatlon(:,:,loin), start, count ), &
!nf90         'dmsbndini: cannot read data for '//'DMS' )
        call handle_ncerr( nf_get_vara_double( ncid, vid, start, count,  dmslatlon(:,:,loin) ), &
           'dmsbndini: cannot read data for '//'DMS' )

        start(3) = hitim
!nf90      call handle_ncerr( nf90_get_var( ncid, vid, dmslatlon(:,:,hiin), start, count ), &
!nf90         'dmsbndini: cannot read data for '//'DMS' )
        call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, dmslatlon(:,:,hiin) ), &
           'dmsbndini: cannot read data for '//'DMS' )

      endif  ! end of masterproc

      ! scatter data to chunked data structures
      call scatter_field_to_chunk (1, 1, 2, plon, dmslatlon, &
                                   dmsin(1,begchunk,1))

      if (masterproc) then
        write(6,*)'dmsbndini: calendar day = ',calday, ' : read data for days ', &
           time(lotim), ' and ',time(hitim)
      endif  ! end of masterproc

      deallocate( dmslatlon, stat=istat )
      call handle_err( istat, &
         'ERROR deallocating memory for dmslatlon in routine dmsbndini')

   end subroutine dmsbndini

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


   subroutine dmsbndint( calday ) 1,12

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Interpolate DMS 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

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

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

      ! Local variables:
      integer ::&
         oldhitim,        &
         vid
      real(r8) ::&
         dt, dt1, tint
      real(r8), allocatable :: dmslatlon(:,:)  ! used for netCDF input
      integer istat
      !-----------------------------------------------------------------------

      ! Check to see if model time is still bounded by dataset times.
      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

         allocate( dmslatlon(plon,plat), stat=istat )
         call alloc_err( istat, 'dmsbndini', 'dmslatlon', &
                         plon*plat )

         if (masterproc) then

           start(3) = hitim
!nf90         call handle_ncerr( nf90_inq_varid( ncid, 'DMS', vid ), &
!nf90            'dmsbndint: cannot find variable '//'DMS' )
!nf90         call handle_ncerr( nf90_get_var( ncid, vid, dmslatlon, start, count ), &
!nf90            'dmsbndint: cannot read data for '//'DMS' )
           call handle_ncerr( nf_inq_varid( ncid, 'DMS', vid ), &
              'dmsbndint: cannot find variable '//'DMS' )
           call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, dmslatlon ), &
              'dmsbndint: cannot read data for '//'DMS' )
           write(6,*)'dmsbndint: read data for day ',time(hitim)

         endif  ! end of masterproc

         ! scatter data to chunked data structures
         call scatter_field_to_chunk (1, 1, 1, plon, dmslatlon, &
                                      dmsin(1,begchunk,hiin))
         if ( lotim .ne. oldhitim ) then
           if (masterproc) then
              ! Read in new lotim data.  Replace old hitim data.
              start(3) = lotim
!nf90            call handle_ncerr( nf90_inq_varid( ncid, 'DMS', vid ), &
!nf90               'dmsbndint: cannot find variable '//'DMS' )
!nf90            call handle_ncerr( nf90_get_var( ncid, vid, dmslatlon, start, count), &
!nf90               'dmsbndint: cannot read data for '//'DMS' )
              call handle_ncerr( nf_inq_varid( ncid, 'DMS', vid ), &
                 'dmsbndint: cannot find variable '//'DMS' )
              call handle_ncerr( nf_get_vara_double( ncid, vid, start, count, dmslatlon), &
                 'dmsbndint: cannot read data for '//'DMS' )
              write(6,*)'dmsbndint: read data for day ',time(lotim)
           endif  ! end of masterproc

           ! scatter data to chunked data structures
           call scatter_field_to_chunk (1, 1, 1, plon, dmslatlon, &
                                        dmsin(1,begchunk,loin))

         end if

         deallocate( dmslatlon, stat=istat )
         call handle_err( istat, &
                 'ERROR deallocating memory for dmslatlon in routine dmsbndint')

      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( size(dms), 0._r8, 1._r8, tint, &
                    dmsin(1,begchunk,loin),        &
                    dmsin(1,begchunk,hiin),        &
                    dms )

   end subroutine dmsbndint

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


   subroutine dmsbndget( ncol, lchnk, x ) 1

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Return DMS 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)     ! DMS emissions (kg DMS/m2/s)

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

      do i = 1, ncol
         x(i) = dms(i,lchnk)
      end do

   end subroutine dmsbndget

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

end module dmsbnd