!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS is distributed in the hope that it will be useful,            !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of    !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     !!
!! GNU General Public License for more details.                      !!
!!                                                                   !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! $Id: drifters_io.F90,v 14.0 2007/03/15 22:38:56 fms Exp $

!!#include <fms_platform.h>

module drifters_io_mod
  implicit none  
  private

  public :: drifters_io_type, drifters_io_new, drifters_io_del, drifters_io_set_time_units
  public :: drifters_io_set_position_names, drifters_io_set_position_units, drifters_io_set_field_names
  public :: drifters_io_set_field_units, drifters_io_write

  ! Globals
  integer, parameter, private   :: MAX_STR_LEN = 128
  character(MAX_STR_LEN), parameter, private :: &
       & version = '$Id: drifters_io.F90,v 14.0 2007/03/15 22:38:56 fms Exp $'

  real :: drfts_eps_t = 10*epsilon(1.)
  

  type drifters_io_type
     real                 :: time
     integer              :: it ! time index
     integer              :: it_id  ! infinite axis index
     integer              :: ncid
     integer              :: nc_positions, nc_fields, nc_ids, nc_time, nc_index_time    
     logical              :: enddef
  end type drifters_io_type

contains

!###############################################################################
  subroutine drifters_io_new(self, filename, nd, nf, ermesg)
    type(drifters_io_type)        :: self
    character(len=*), intent(in)  :: filename
    integer, intent(in)           :: nd  ! number of dims
    integer, intent(in)           :: nf  ! number of fields
    character(len=*), intent(out) :: ermesg

    integer ier, nc_it_id, nc_nd, nc_nf
    integer :: size1(1), size2(2)
    include 'netcdf.inc'

    ermesg=''
    self%enddef = .FALSE.

    ier = nf_create(filename, NF_CLOBBER, self%ncid)
    if(ier/=NF_NOERR) ermesg = 'drifters_io_new::nf_create ('//filename//') '//nf_strerror(ier)

    ! global attributes
    ier = nf_put_att_text(self%ncid, NF_GLOBAL, 'version', len_trim(version), trim(version))
    

    ! dimensions
    ier = nf_def_dim(self%ncid, 'it_id', NF_UNLIMITED, nc_it_id)
    if(ier/=NF_NOERR) ermesg = 'drifters_io_new::nf_def_dim (it_id) '//nf_strerror(ier)

    ier = nf_def_dim(self%ncid, 'nf', nf, nc_nf)
    if(ier/=NF_NOERR) ermesg = 'drifters_io_new::nf_def_dim (nf) '//nf_strerror(ier)

    ier = nf_def_dim(self%ncid, 'nd', nd, nc_nd)
    if(ier/=NF_NOERR) ermesg = 'drifters_io_new::nf_def_dim (nd) '//nf_strerror(ier)

    ! variables
    size1 = (/nc_it_id/)
    ier = nf_def_var(self%ncid, 'index_time', NF_INT, 1, size1, self%nc_index_time)
    if(ier/=NF_NOERR) ermesg = 'drifters_io_new::nf_def_var (index_time)'//nf_strerror(ier)

    ier = nf_def_var(self%ncid, 'time', NF_DOUBLE, 1, size1, self%nc_time)
    if(ier/=NF_NOERR) ermesg = 'drifters_io_new::nf_def_var (time)'//nf_strerror(ier)

    ier = nf_def_var(self%ncid, 'ids', NF_INT, 1, size1, self%nc_ids)
    if(ier/=NF_NOERR) ermesg = 'drifters_io_new::nf_def_var (ids)'//nf_strerror(ier)

    size2 = (/nc_nd, nc_it_id/)
    ier = nf_def_var(self%ncid, 'positions', NF_DOUBLE, 2, size2, self%nc_positions)
    if(ier/=NF_NOERR) ermesg = 'drifters_io_new::nf_def_var (positions)'//nf_strerror(ier)

    size2 = (/nc_nf, nc_it_id/)
    ier = nf_def_var(self%ncid, 'fields', NF_DOUBLE, 2, size2, self%nc_fields)
    if(ier/=NF_NOERR) ermesg = 'drifters_io_new::nf_def_var (fields)'//nf_strerror(ier)
    
    self%time  = -huge(1.)
    self%it    = -1
    self%it_id = 1

  end subroutine drifters_io_new

!###############################################################################
  subroutine drifters_io_del(self, ermesg)
    type(drifters_io_type)        :: self
    character(len=*), intent(out) :: ermesg

    integer ier
    include 'netcdf.inc'

    ermesg = ''
    
    ier = nf_close(self%ncid)
    if(ier/=NF_NOERR) ermesg = 'drifters_io_del::nf_close '//nf_strerror(ier)
    
  end subroutine drifters_io_del

!###############################################################################
  subroutine drifters_io_set_time_units(self, name, ermesg)
    type(drifters_io_type)        :: self
    character(len=*), intent(in)  :: name
    character(len=*), intent(out) :: ermesg

    integer ier
    include 'netcdf.inc'

    ermesg = ''
    ier = nf_put_att_text(self%ncid, NF_GLOBAL, &
         & 'time_units', len_trim(name), trim(name))
    if(ier/=NF_NOERR) &
         & ermesg = 'drifters_io_set_time_units::failed to add time_units attribute ' &
         & //nf_strerror(ier)

  end subroutine drifters_io_set_time_units

!###############################################################################
  subroutine drifters_io_set_position_names(self, names, ermesg)
    type(drifters_io_type)        :: self
    character(len=*), intent(in)  :: names(:)
    character(len=*), intent(out) :: ermesg

    integer n, ier, i
    character(len=128) :: attname
    include 'netcdf.inc'

    n = size(names)
    ermesg = ''

    do i = 1, n
       write(attname, '(i6)' ) i
       attname = 'name_'//adjustl(attname)
       ier = nf_put_att_text(self%ncid, self%nc_positions, &
            & trim(attname), len_trim(names(i)), trim(names(i)))
       if(ier/=NF_NOERR) &
            & ermesg = 'drifters_io_set_position_names::failed to add name attribute to positions '//nf_strerror(ier)
    enddo

  end subroutine drifters_io_set_position_names

!###############################################################################
  subroutine drifters_io_set_position_units(self, names, ermesg)
    type(drifters_io_type)        :: self
    character(len=*), intent(in)  :: names(:)
    character(len=*), intent(out) :: ermesg

    integer n, ier, i
    character(len=128) :: attname
    include 'netcdf.inc'

    n = size(names)
    ermesg = ''

    do i = 1, n
       write(attname, '(i6)' ) i
       attname = 'unit_'//adjustl(attname)
       ier = nf_put_att_text(self%ncid, self%nc_positions, &
            & trim(attname), len_trim(names(i)), trim(names(i)))
       if(ier/=NF_NOERR) &
            & ermesg = 'drifters_io_set_position_names::failed to add unit attribute to positions '//nf_strerror(ier)
    enddo
        
  end subroutine drifters_io_set_position_units

!###############################################################################
  subroutine drifters_io_set_field_names(self, names, ermesg)
    type(drifters_io_type)        :: self
    character(len=*), intent(in)  :: names(:)
    character(len=*), intent(out) :: ermesg

    integer n, ier, i
    character(len=128) :: attname
    include 'netcdf.inc'

    n = size(names)
    ermesg = ''

    do i = 1, n
       write(attname, '(i6)' ) i
       attname = 'name_'//adjustl(attname)
       ier = nf_put_att_text(self%ncid, self%nc_fields, &
            & trim(attname), len_trim(names(i)), trim(names(i)))
       if(ier/=NF_NOERR) &
            & ermesg = 'drifters_io_set_field_names::failed to add name attribute to fields '//nf_strerror(ier)
    enddo

  end subroutine drifters_io_set_field_names

!###############################################################################
  subroutine drifters_io_set_field_units(self, names, ermesg)
    type(drifters_io_type)        :: self
    character(len=*), intent(in)  :: names(:)
    character(len=*), intent(out) :: ermesg

    integer n, ier, i
    character(len=128) :: attname
    include 'netcdf.inc'

    n = size(names)
    ermesg = ''

    do i = 1, n
       write(attname, '(i6)' ) i
       attname = 'unit_'//adjustl(attname)
       ier = nf_put_att_text(self%ncid, self%nc_fields, &
            & trim(attname), len_trim(names(i)), trim(names(i)))
       if(ier/=NF_NOERR) &
            & ermesg = 'drifters_io_set_field_units::failed to add unit attribute to fields '//nf_strerror(ier)
    enddo
    
  end subroutine drifters_io_set_field_units
!###############################################################################

  subroutine drifters_io_write(self, time, np, nd, nf, ids, positions, fields, ermesg)
    type(drifters_io_type)        :: self
    real, intent(in)              :: time
    integer, intent(in)           :: np    ! number of dirfters
    integer, intent(in)           :: nd    ! number of dimensions
    integer, intent(in)           :: nf    ! number of fields
    integer, intent(in)           :: ids(np)          ! of size np
    real, intent(in)              :: positions(nd,np) ! nd times np
    real, intent(in)              :: fields(nf,np)    ! nf times np
    character(len=*), intent(out) :: ermesg

    integer ier, i
    integer :: start1(1), len1(1), start2(2), len2(2)
    integer :: it_indices(np)
    real    :: time_array(np)
    include 'netcdf.inc'

    ermesg = ''
    
    if(.not. self%enddef) then
       ier = nf_enddef(self%ncid)
       if(ier/=NF_NOERR) then 
            ermesg = 'drifters_io_write::nf_enddef failure. No data will be written. '//nf_strerror(ier)
            return
       endif
       self%enddef = .TRUE.
    endif

    if(abs(time - self%time) > drfts_eps_t) then
       self%it = self%it + 1
       self%time = time
    endif

    start1(1) = self%it_id
    len1(1)   = np

    it_indices = (/(self%it,i=1,np)/)
    ier = nf_put_vara_int( self%ncid, self%nc_index_time, start1, len1, it_indices )
    if(ier/=NF_NOERR) &
         & ermesg = 'drifters_io_write::failed to write index_time: ' //nf_strerror(ier)

    time_array = (/(time,i=1,np)/)
    ier = nf_put_vara_double( self%ncid, self%nc_time, start1, len1, time_array )
    if(ier/=NF_NOERR) &
         & ermesg = 'drifters_io_write::failed to write time: ' //nf_strerror(ier)

    ier = nf_put_vara_int(self%ncid, self%nc_ids, start1, len1, ids)
    if(ier/=NF_NOERR) &
         & ermesg = 'drifters_io_write::failed to write ids: '//nf_strerror(ier)

    start2(1) = 1
    start2(2) = self%it_id

    len2(1)   = nd
    len2(2)   = np

    ier = nf_put_vara_double(self%ncid, self%nc_positions, start2, len2, positions)
    if(ier/=NF_NOERR) &
         & ermesg = 'drifters_io_write::failed to write positions: '//nf_strerror(ier)

    len2(1)   = nf
    len2(2)   = np    
    
    ier = nf_put_vara_double(self%ncid, self%nc_fields, start2, len2, fields)
    if(ier/=NF_NOERR) &
         & ermesg = 'drifters_io_write::failed to write fields: '//nf_strerror(ier)

    self%it_id = self%it_id + np
    
  end subroutine drifters_io_write

end module drifters_io_mod
 !###############################################################################
 !###############################################################################
#ifdef _TEST_DRIFTERS_IO
! set FC=pgf95
! set FOPTS='-r8 -g -Mdclchk -Minform=warn'
! set INCS='-I/usr/local/include'
! set LIBS='-L/usr/local/lib -lnetcdf'
! $FC $INCS $FOPTS -D_TEST_DRIFTERS_IO drifters_io.F90 $LIBS
program test
  use drifters_io_mod
  implicit none
  type(drifters_io_type) :: drfts_io
  character(len=128) :: ermesg
  character(len=31) :: filename
  integer :: np, nd, nf, nt, i, j, k, npmax
  real :: dt, time, xmin, xmax, ymin, ymax, u, v, dr, x, y
  integer, allocatable :: ids(:)
  real, allocatable :: positions(:,:), fields(:,:)

  ! number of dimensions
  nd = 3
  ! number of fields 
  nf = 2
  ! max number of dirfters 
  npmax = 20
  ! number of time steps
  nt = 50
  ! starting time
  time = 0.

  ! domain boundary. (drifters outside domain will not be written to file.)
  xmin = 0.
  ymin = 0.
  xmax = 1.
  ymax = 1.

  ! constant velocity
  u = (xmax-xmin)*sqrt(2.)
  v = (ymax-ymin)*sqrt(2.)
  dt = 1/real(nt)

  ! open file
  
  filename = 'test.nc'
  call drifters_io_new(drfts_io, filename, nd, nf, ermesg)
  if(ermesg/='') print *,'ERROR after drifters_io_new: ', ermesg

  ! set attributes

  call drifters_io_set_position_names(drfts_io, (/'x','y','z'/), ermesg)
  if(ermesg/='') print *,'ERROR after drifters_io_position_names: ', ermesg

  ! note the trailing blanks in the first field, which are added here to 
  ! ensure that "salinity" will not be truncated (all names must have the 
  ! same length)
  call drifters_io_set_field_names(drfts_io, (/'temp    ','salinity'/), ermesg)
  if(ermesg/='') print *,'ERROR after drifters_io_field_names: ', ermesg

  call drifters_io_set_position_units(drfts_io, (/'deg east ','deg north','meters'/), ermesg)
  if(ermesg/='') print *,'ERROR after drifters_io_position_units: ', ermesg
  
  call drifters_io_set_field_units(drfts_io, (/'deg K ','ppm'/), ermesg)
  if(ermesg/='') print *,'ERROR after drifters_io_field_units: ', ermesg

  allocate(positions(nd, npmax), ids(npmax), fields(nf, npmax))
  dr = sqrt( (xmax-xmin)**2 + (ymax-ymin)**2 )/real(npmax)


  ! x
  positions(1, :) = +(/ (i*dr,i=0,npmax-1) /)/sqrt(2.)
  ! y
  positions(2, :) = -(/ (i*dr,i=0,npmax-1) /)/sqrt(2.)
  ! z
  positions(3, :) = 0.

  ! drifters' identity array (can be any integer number)
  ids = (/ (i, i=1, npmax) /)
  
  ! set fields as a function of space time
  fields(1, :) = sqrt( (positions(1,:)-xmin)**2 + (positions(2,:)-ymin)**2 )
  fields(2, :) = positions(1,:)-u*time + positions(2,:)-v*time ! invariant

  ! write to disk only drifters inside domain
  do i = 1, npmax
     x = positions(1,i)
     y = positions(2,i)
     if(x>=xmin .and. x<=xmax .and. y>=ymin .and. y<=ymax) then
        call drifters_io_write(drfts_io, time, np=1, nd=nd, nf=nf, &
             & ids=ids(i), positions=positions(:,i), fields=fields(:,i), ermesg=ermesg)
        if(ermesg/='') print *,'ERROR after drifters_io_write: ', ermesg
     endif
  enddo

  ! advect
  
  do j = 1, nt
     time = time + dt
     positions(1, :) = positions(1, :) + u*dt
     positions(2, :) = positions(2, :) + v*dt
     fields(1, :) = sqrt( (positions(1,:)-xmin)**2 + (positions(2,:)-ymin)**2 )
     fields(2, :) = positions(1,:)-u*time + positions(2,:)-v*time ! invariant

     do i = 1, npmax
        x = positions(1,i)
        y = positions(2,i)
        if(x>=xmin .and. x<=xmax .and. y>=ymin .and. y<=ymax) then
           call drifters_io_write(drfts_io, time, np=1, nd=nd, nf=nf, &
                & ids=ids(i), positions=positions(:,i), fields=fields(:,i), ermesg=ermesg)
           if(ermesg/='') print *,'ERROR after drifters_io_write: ', ermesg
        endif
     enddo
     
  enddo

  deallocate(positions, ids, fields)

  call drifters_io_del(drfts_io, ermesg)
  if(ermesg/='') print *,'ERROR after drifters_io_del: ', ermesg

end program test
#endif 
! _TEST_DRIFTERS_IO