#include "options.inc"


      module diag_float_module

c-----------------------------------------------------------------------
c     Floats are assumed to be of zero mass, neutrally buoyant and move
c     with the local three dimensional flow during each time step.
c     Initially, particles are randomly distributed within
c     the volume given by ptslon,ptelon, ...
c
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c
c     linked in the code in init_diag (initialization)
c     and diag (integration)
c     namelist parameter in spflame_module:
c        enable_diag_float (logical, .false. at default)
c        diag_float_file   (character, 'float.cdf' at default)
c
c     adapted to SPFLAME : mar 2002, C Eden
c-----------------------------------------------------------------------

      implicit none
      integer :: nptraj ! number of float
      real    :: ptslon, ptelon, ptslat, ptelat, ptsdpt, ptedpt
c     ptslon  = starting longitude for initial particle distribution
c     ptelon  = ending longitude for initial particle distribution
c     ptslat  = starting latitude for initial particle distribution
c     ptelat  = ending latitude for initial particle distribution
c     ptsdpt  = starting depth for initial particle distribution
c     ptedpt  = ending depth for initial particle distribution
c
      real, allocatable :: pxyz(:,:) 
c     pxyz    = particle coordinates. index (1,2,3) is for particle 
c               (longitude, latitude, depth).

      integer, allocatable :: pijk(:,:) 
c     pijk    = the particle is bounded by the volume with vertices
c               given by the eight nearest surrounding model grid points
c               on the "xu","yu", and "z" grids. index (1,2,3) locates
c               the (longitude, latitude, depth) index of the deepest 
c               northeast corner of this bounding volume.

      real, allocatable :: pdens(:)  ! potential density of particle
      real, allocatable :: pvolu(:)  ! volume transport of particle
      real, allocatable :: pvolv(:)  ! volume transport of particle

      real :: spval = -1.0e20
      logical, allocatable :: ptlost(:)

      logical :: set_dens_time_step_for_floats = .true.
      integer :: dens_time_step_for_floats = 0
      end module diag_float_module


      subroutine diag_float_init
      use spflame_module
      use diag_float_module
      implicit none
      integer :: io,n,i,j
      character (len=80) :: name,unit
      character (len=32) :: stamp
#ifdef netcdf_diagnostics
#include "netcdf.inc"
      integer :: iret,ncid,timedim,nrdim,xid,yid,zid,timeid,dims(4)
      integer :: rhoid,lat_tdim,lon_tdim,lat_tid,lon_tid,topoid
      integer :: voluid, volvid
#ifdef netcdf_real4
      real (kind=4), allocatable :: v2(:),v4(:,:)
#else
      real, allocatable          :: v2(:),v4(:,:)
#endif
#endif


      write(name,'("floats_",i3,".dta")') sub_domain
      call replace_space_zero(name)
      if (my_pe==0) then
        print*,' --> reading float data from restart file ',
     &        name(1:len_trim(name))
      endif

      call get_free_iounit(io)
      open(io,file=name,form='unformatted',status='old',err=20)
      read (io) stamp
      read (io) nptraj
      if (my_pe==0) print*,' found ',nptraj,' in file (at ',stamp,')'
      allocate( pxyz(3,nptraj) )
      allocate( pijk(3,nptraj) )
      allocate( ptlost(nptraj) ); ptlost = .false.
      allocate( pdens(nptraj) )
      allocate( pvolu(nptraj) )
      allocate( pvolv(nptraj) )

      if (my_pe==0) read (io) pxyz,pijk,pvolu,pvolv
      close (io)

      call bcast_real(pxyz,3*nptraj,0)
      call bcast_integer(pijk,3*nptraj,0)
      call bcast_real(pvolu,nptraj,0)
      call bcast_real(pvolv,nptraj,0)

      goto 30

 20   if (my_pe==0) then
        print*,''
        print*,'---------------------------------------------'
        print*,' WARNING : cannot read file ',name(1:len_trim(name))
        print*,'---------------------------------------------'
        print*,''
        print*,' initialising floats therefore as follows'
      endif
c      
c-----------------------------------------------------------------------
c     distribute particles within volume defined using starting and
c     ending longitudes, latitudes, and depths given by 
c     ptslon, ptelon, ptslat, ptelat. ptsdpt, and ptedpt.
c-----------------------------------------------------------------------
c
c     begin of USER input
c
c      nptraj = 500
c      ptslon = 303.9   ! 
c      ptelon = 304.1   ! 
c      ptslat = 9.0
c      ptelat = 10.0
c      ptsdpt = 1000.e2
c      ptedpt = 3200.e2

      allocate( pxyz(3,nptraj) )
      allocate( pijk(3,nptraj) )
      allocate( pdens(nptraj) )
      allocate( pvolu(nptraj) )
      allocate( pvolv(nptraj) )
      allocate( ptlost(nptraj) ); ptlost = .false.
      call diag_float_dist (1, nptraj)    
      if (ptslat==ptelat) then
        call diag_float_dist_meridional (1,nptraj)
      else
       call diag_float_dist (1, nptraj)    
      endif
      call diag_float_vol     
c
c     end of USER input
c
 30   continue

#ifdef netcdf_diagnostics
c
c     initialize netcdf file
c
      if (my_pe == 0) then
       print*,' initializing NetCDF output file ',
     &          diag_float_file(1:len_trim(diag_float_file))
       ncid = nccre (diag_float_file, NCCLOB, iret)
       iret=nf_set_fill(ncid, NF_NOFILL, iret)
       call store_info_cdf(ncid)
c      dimensions
       Timedim   = ncddef(ncid,'Time', nf_unlimited, iret)
       nrdim     = ncddef(ncid,'Nr  ', nptraj, iret)
       lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
       lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
c      variables
       Lon_tid = ncvdef(ncid,'Longitude_t',NCFLOAT,1,lon_tdim,iret)
       Lat_tid = ncvdef(ncid,'Latitude_t', NCFLOAT,1,lat_tdim,iret)
       Timeid   =  ncvdef(ncid,'Time', NCFLOAT,1,timedim,iret)
       name = 'Time'; unit = 'days'
       call ncaptc(ncid,Timeid, 'long_name', NCCHAR,4, name, iret) 
       call ncaptc(ncid,Timeid, 'units',    NCCHAR,4, unit, iret) 
       call ncaptc(ncid,Timeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
       dims=(/nrdim,timedim,1,1/)
       xid=ncvdef (ncid,'lon', NCFLOAT,2,dims,iret)
       yid=ncvdef (ncid,'lat', NCFLOAT,2,dims,iret)
       zid=ncvdef (ncid,'depth', NCFLOAT,2,dims,iret)
       rhoid=ncvdef (ncid,'rho', NCFLOAT,2,dims,iret)
       voluid=ncvdef (ncid,'volu', NCFLOAT,1,dims,iret)
       volvid=ncvdef (ncid,'volv', NCFLOAT,1,dims,iret)
       name = 'Longitude on T grid'; unit = 'degrees_W'
       call dvcdf(ncid,xid,name,len_trim(name),unit,
     &            len_trim(unit),spval)
       name = 'Latitude on T grid'; unit = 'degrees_N'
       call dvcdf(ncid,yid,name,len_trim(name),unit,
     &            len_trim(unit),spval)
       name = 'Depth on T grid'; unit = 'm'
       call dvcdf(ncid,zid,name,len_trim(name),unit,
     &            len_trim(unit),spval)
       name = 'Potential density'; unit = 'kg/m^3'
       call dvcdf(ncid,rhoid,name,len_trim(name),unit,
     &            len_trim(unit),spval)
       name = 'Zonal volume transport of particle'; unit = 'm^3/s'
       call dvcdf(ncid,voluid,name,len_trim(name),unit,
     &            len_trim(unit),spval)
       name = 'Meridional volume transport of particle'; unit = 'm^3/s'
       call dvcdf(ncid,volvid,name,len_trim(name),unit,
     &            len_trim(unit),spval)
       dims=(/lon_tdim,lat_tdim,1,1/)
       topoid=ncvdef (ncid,'topo', NCFLOAT,2,dims,iret)
       call ncendf(ncid, iret)
       allocate( v4(imt,jmt) ); v4=0.0
       do i=1,imt
        do j=1,jmt
          if (kmt_big(i,j)>0) v4(i,j)=zw(kmt_big(i,j) )/100.0
        enddo
       enddo
       call ncvpt(ncid, topoid, (/1,1/), (/imt,jmt/),v4, iret)
       deallocate(v4)
       allocate( v2(max(imt,jmt,km)) )
       v2(1:imt)=xt
       call ncvpt(ncid, Lon_tid, 1, imt,v2, iret)
       v2(1:jmt)=yt
       call ncvpt(ncid, Lat_tid, 1, jmt,v2, iret)
       deallocate(v2)

       allocate( v2(nptraj) )
       v2=pvolu; call ncvpt(ncid, voluid, 1, nptraj,v2, iret)
       v2=pvolv; call ncvpt(ncid, volvid, 1, nptraj,v2, iret)
       deallocate(v2)

       call ncclos (ncid, iret)
      endif
#endif
      end subroutine diag_float_init




      subroutine diag_float_vol
c
c-----------------------------------------------------------------------
c    relate transport to particles
c-----------------------------------------------------------------------
c
      use spflame_module
      use diag_float_module
      implicit none
      integer :: i,j,k,count,m,n
      logical, dimension(nptraj) :: ptdone
      integer :: is,ie,js,je

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      ptdone(1:nptraj) = .false.
      do n=1,nptraj
       if (pijk(1,n) >= is .and. pijk(1,n) <= ie_pe .and.
     &      pijk(2,n) >= js .and. pijk(2,n) <= je_pe .and.
     &      pijk(3,n) > 1   .and. pijk(3,n) <= km    ) then
          count=0
          i  = pijk(1,n)
          j  = pijk(2,n)
          k  = pijk(3,n)
          do m=1,nptraj
            if ( pijk(1,m) == i .and.
     &            pijk(2,m) == j .and.
     &            pijk(3,m) == k ) count=count+1
          enddo
          pvolv(n)=0.5*( u(i,k,j,2,tau)+u(i,k+1,j,2,tau) )/count
     &         *dxu(i)*csu(j)*dzw(k)/1e6
          pvolu(n)=0.5*( u(i,k,j,1,tau)+u(i,k+1,j,1,tau) )/count
     &         *dyu(j)*dzw(k)/1e6
          ptdone(n)=.true.
       endif
      enddo
c
c-----------------------------------------------------------------------
c     exchange informations around the PEs
c-----------------------------------------------------------------------
c
      do n=1,nptraj
        if (.not. ptdone(n) ) then
          pvolu(n)=0.
          pvolv(n)=0
        endif
        call global_sum_vec(pvolu(n),1)
        call global_sum_vec(pvolv(n),1)
      enddo
      end subroutine diag_float_vol


      subroutine diag_float_dist (ns, ne)
c
c=======================================================================
c     distribute n particles (ne-ns+1) within volume given by
c     ptslon, ptelon, ptslat, ptelat, ptsdpt, ptedpt by uniformly
c     placing approximately n**(1/3) particles along each dimension
c
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c=======================================================================
c
      use spflame_module
      use diag_float_module
      implicit none
      integer :: ns,ne,n,i,k,j,indp
      real :: fxa,cubr,distx,disty,distz

c
      if (ns .gt. nptraj .or. ne .gt. nptraj) then
        write (6,*) ' => Error: ns=',ns,', ne=',ne
	call halt_stop(' in diag_pdist')
      endif
c
c-----------------------------------------------------------------------
c     constrain the volume (containing starting positions of
c     particles) to lie within the model domain.      
c-----------------------------------------------------------------------
c
      if (ptslon .lt. xu(2))     ptslon = xu(2)
      if (ptslon .gt. xu(imt-1)) ptslon = xu(imt-1)
      if (ptslat .lt. yu(1))     ptslat = yu(1)
      if (ptslat .gt. yu(jmt-1)) ptslat = yu(jmt-1)
      if (ptsdpt .lt. zt(1))     ptsdpt = zt(1)
      if (ptsdpt .gt. zt(km))    ptsdpt = zt(km)
c
      if (ptelon .lt. xu(2))     ptelon = xu(2)
      if (ptelon .gt. xu(imt-1)) ptelon = xu(imt-1)
      if (ptelat .lt. yu(1))     ptelat = yu(1)
      if (ptelat .gt. yu(jmt-1)) ptelat = yu(jmt-1)
      if (ptedpt .lt. zt(1))     ptedpt = zt(1)
      if (ptedpt .gt. zt(km))    ptedpt = zt(km)
c
      if (ptslon .gt. ptelon) then
        fxa    = ptslon
        ptslon = ptelon
        ptelon = fxa
      endif
      if (ptslat .gt. ptelat) then
        fxa    = ptslat
        ptslat = ptelat
        ptelat = fxa
      endif
      if (ptsdpt .gt. ptedpt) then
        fxa    = ptsdpt
        ptsdpt = ptedpt
        ptedpt = fxa
      endif
c
c-----------------------------------------------------------------------
c     distribute the particles throughout the volume
c-----------------------------------------------------------------------
c
      cubr       = (float(ne-ns+1))**0.333333
      distx      = (ptelon - ptslon)/cubr
      disty      = (ptelat - ptslat)/cubr
      distz      = (ptedpt - ptsdpt)/cubr
      pxyz(1,ns) = .5*distx + ptslon
      pxyz(2,ns) = .5*disty + ptslat
      pxyz(3,ns) = .5*distz + ptsdpt
      do n=ns,ne
        if (n .gt. ns) then
	  pxyz(1,n) = pxyz(1,n-1)
	  pxyz(2,n) = pxyz(2,n-1)
	  pxyz(3,n) = pxyz(3,n-1)
	endif
	pxyz(1,n)  = pxyz(1,n) + distx
	if (pxyz(1,n) .gt. ptelon) then
	  pxyz(1,n) = ptslon + (pxyz(1,n)-ptelon)
	  pxyz(2,n) = pxyz(2,n) + disty
	  if (pxyz(2,n) .gt. ptelat) then
	    pxyz(2,n) = ptslat + (pxyz(2,n)-ptelat)
	    pxyz(3,n) = pxyz(3,n) + distz
	  endif
	endif
c
        i = indp (pxyz(1,n), xu, imt)
        if (xu(i) .gt. pxyz(1,n)) then
          pijk(1,n) = i
        else
          pijk(1,n) = i+1
        endif
c
        j = indp (pxyz(2,n), yu, jmt)
        if (yu(j) .gt. pxyz(2,n)) then
          pijk(2,n) = j
        else
          pijk(2,n) = j+1
        endif
c
        k = indp (pxyz(3,n), zt, km)
        if (zt(k) .gt. pxyz(3,n)) then
          pijk(3,n) = k
        else
          pijk(3,n) = k+1
        endif
      enddo


      if (my_pe==0) then
       write (6,9000)
       write (6,*) ne-ns+1,' particles were initialized',
     &           ' to lie within the volume described by:'
       write (6,'(1x,f8.2," < lon <",f8.2)') ptslon, ptelon
       write (6,'(1x,f8.2," < lat <",f8.2)') ptslat, ptelat
       write (6,'(1x,e10.3," < dpt < ",e10.3)') ptsdpt, ptedpt
       write (6,*) ' '
      endif

9000  format (/20x,'P A R T I C L E    I N I T I A L I Z A T I O N'/)

      end subroutine diag_float_dist 




      subroutine diag_float_dist_meridional (ns, ne)
c
c=======================================================================
c     distribute n particles (ne-ns+1) within volume given by
c     ptslon, ptelon, ptslat, ptelat, ptsdpt, ptedpt by uniformly
c     placing approximately n**(1/2) particles along each dimension
c
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c=======================================================================
c
      use spflame_module
      use diag_float_module
      implicit none
      integer :: ns,ne,n,i,k,j,indp
      real :: fxa,cubr,distx,disty,distz

c
      if (ns .gt. nptraj .or. ne .gt. nptraj) then
        write (6,*) ' => Error: ns=',ns,', ne=',ne
	call halt_stop(' in diag_pdist')
      endif
c
c-----------------------------------------------------------------------
c     constrain the volume (containing starting positions of
c     particles) to lie within the model domain.      
c-----------------------------------------------------------------------
c
      if (ptslon .lt. xu(2))     ptslon = xu(2)
      if (ptslon .gt. xu(imt-1)) ptslon = xu(imt-1)
      if (ptslat .lt. yu(1))     ptslat = yu(1)
      if (ptslat .gt. yu(jmt-1)) ptslat = yu(jmt-1)
      if (ptsdpt .lt. zt(1))     ptsdpt = zt(1)
      if (ptsdpt .gt. zt(km))    ptsdpt = zt(km)
c
      if (ptelon .lt. xu(2))     ptelon = xu(2)
      if (ptelon .gt. xu(imt-1)) ptelon = xu(imt-1)
      if (ptelat .lt. yu(1))     ptelat = yu(1)
      if (ptelat .gt. yu(jmt-1)) ptelat = yu(jmt-1)
      if (ptedpt .lt. zt(1))     ptedpt = zt(1)
      if (ptedpt .gt. zt(km))    ptedpt = zt(km)
c
      if (ptslon .gt. ptelon) then
        fxa    = ptslon
        ptslon = ptelon
        ptelon = fxa
      endif
      if (ptslat .gt. ptelat) then
        fxa    = ptslat
        ptslat = ptelat
        ptelat = fxa
      endif
      if (ptsdpt .gt. ptedpt) then
        fxa    = ptsdpt
        ptsdpt = ptedpt
        ptedpt = fxa
      endif
c
c-----------------------------------------------------------------------
c     distribute the particles throughout the volume
c-----------------------------------------------------------------------
c
      cubr       = (float(ne-ns+1))**0.5
      distx      = (ptelon - ptslon)/cubr
      disty      = (ptelat - ptslat)/cubr
      distz      = (ptedpt - ptsdpt)/cubr
      pxyz(1,ns) = .5*distx + ptslon
      pxyz(2,ns) = ptslat
      pxyz(3,ns) = .5*distz + ptsdpt
      do n=ns,ne
        if (n .gt. ns) then
	  pxyz(1,n) = pxyz(1,n-1)
	  pxyz(2,n) = pxyz(2,n-1)
	  pxyz(3,n) = pxyz(3,n-1)
	endif
	pxyz(1,n)  = pxyz(1,n) + distx
	if (pxyz(1,n) .gt. ptelon) then
	  pxyz(1,n) = ptslon + (pxyz(1,n)-ptelon)
	  pxyz(3,n) = pxyz(3,n) + distz
	endif
c
        i = indp (pxyz(1,n), xu, imt)
        if (xu(i) .gt. pxyz(1,n)) then
          pijk(1,n) = i
        else
          pijk(1,n) = i+1
        endif
c
        j = indp (pxyz(2,n), yu, jmt)
        if (yu(j) .gt. pxyz(2,n)) then
          pijk(2,n) = j
        else
          pijk(2,n) = j+1
        endif
c
        k = indp (pxyz(3,n), zt, km)
        if (zt(k) .gt. pxyz(3,n)) then
          pijk(3,n) = k
        else
          pijk(3,n) = k+1
        endif
      enddo


      if (my_pe==0) then
       write (6,9000)
       write (6,*) ne-ns+1,' particles were initialized',
     &           ' to lie within the section described by:'
       write (6,'(1x,f8.2," < lon <",f8.2)') ptslon, ptelon
       print*, 'lat=',ptslat
       write (6,'(1x,e10.3," < dpt < ",e10.3)') ptsdpt, ptedpt
       write (6,*) ' '
      endif

9000  format (/20x,'P A R T I C L E    I N I T I A L I Z A T I O N'/)

      end subroutine diag_float_dist_meridional






      subroutine diag_float
c=======================================================================
c     integrate particle trajectories
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c=======================================================================
      use spflame_module
      use diag_float_module
      implicit none
      integer :: n,j,i,k,m
      real :: cmdeg, rcos
      real :: xold,yold,zold
      real :: xe,xw,yn,ys,za,zb,dv
      real :: xeyszb, xwyszb, xeysza, xwysza, xeynzb, xwynzb 
      real :: xeynza, xwynza , uu, vv, ww
      logical, dimension(nptraj) :: ptdone
#ifdef netcdf_diagnostics
#include "netcdf.inc"
      integer :: iret,ncid,timedim,xid,yid,zid,timeid
      integer :: rhoid
      integer :: corner(2),edges(2)
#ifdef netcdf_real4
      real (kind=4) :: v4(nptraj)
#else
      real          :: v4(nptraj)
#endif
#else
      integer :: io
#endif
      integer :: is,ie,js,je

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c     boundary exchange for vertical velocity
c     (this is a bit odd but no way around)
c
      call border_exchg(adv_vbu,km+1,1)
      call set_cyclic(adv_vbu,km+1,1)
c
c
      ptdone(1:nptraj) = .false.
      cmdeg  = 8.982799e-8
c
c-----------------------------------------------------------------------
c    calculate trajectory for all particles 
c-----------------------------------------------------------------------
c
      do n=1,nptraj

        if (.not. ptlost(n) .and.
     &      pijk(1,n) >= is .and. pijk(1,n) <= ie_pe .and.
     &      pijk(2,n) >= js .and. pijk(2,n) <= je_pe .and.
     &      pijk(3,n) > 1   .and. pijk(3,n) <= km    ) then

          ptdone(n) = .true.
c
c-----------------------------------------------------------------------
c         the particle is bounded by the volume with vertices given by
c         the eight nearest surrounding model grid points on the "xu",
c         "yu", and "zt" grids. (i,j,k) is the index of the deepest 
c         northeast corner of this bounding volume.
c-----------------------------------------------------------------------
c         
          i  = pijk(1,n)
          j  = pijk(2,n)
          k  = pijk(3,n)
c
c-----------------------------------------------------------------------
c         compute volume weights for linear interpolation of velocity
c         at vertices of bounding volume to the particle position.
c
c         distances between particle and bounding volume faces
c
c         xe = distance to the east face
c         xw = distance to the west face
c         yn = distance to the north face
c         ys = distance to the south face
c         za = distance above to the top face
c         zb = distance below to the bottom face
c-----------------------------------------------------------------------
c

          xe = (xu(i) - pxyz(1,n))
          xw = (pxyz(1,n) - xu(i-1))
          yn = (yu(j) - pxyz(2,n))
          ys = (pxyz(2,n) - yu(j-1))
          za = (pxyz(3,n) - zt(k-1))
          zb = (zt(k) - pxyz(3,n))
	  dv = 1./((xu(i)-xu(i-1))*(yu(j)-yu(j-1))*(zt(k)-zt(k-1)))
c
c-----------------------------------------------------------------------
c         construct velocity at position of particle by 3-d linear
c         interpolation. 
c-----------------------------------------------------------------------
c
          xeyszb = xe*ys*zb*dv
          xwyszb = xw*ys*zb*dv
          xeysza = xe*ys*za*dv
          xwysza = xw*ys*za*dv
          xeynzb = xe*yn*zb*dv
          xwynzb = xw*yn*zb*dv
          xeynza = xe*yn*za*dv
          xwynza = xw*yn*za*dv
c          
          uu = u(i-1,k-1,j,1,tau)*xeyszb   + u(i,k-1,j,1,tau)*xwyszb
     &        +u(i-1,k  ,j,1,tau)*xeysza   + u(i,k  ,j,1,tau)*xwysza
     &        +u(i-1,k-1,j-1,1,tau)*xeynzb + u(i,k-1,j-1,1,tau)*xwynzb
     &        +u(i-1,k  ,j-1,1,tau)*xeynza + u(i,k  ,j-1,1,tau)*xwynza
c
          vv = u(i-1,k-1,j,2,tau)*xeyszb   + u(i,k-1,j,2,tau)*xwyszb
     &        +u(i-1,k  ,j,2,tau)*xeysza   + u(i,k  ,j,2,tau)*xwysza
     &        +u(i-1,k-1,j-1,2,tau)*xeynzb + u(i,k-1,j-1,2,tau)*xwynzb
     &        +u(i-1,k  ,j-1,2,tau)*xeynza + u(i,k  ,j-1,2,tau)*xwynza
c
c         interpolate vertical velocities at the bases of
c         the "u" cells.
c
          if (pxyz(3,n) .gt. zw(k-1)) then
	    za = pxyz(3,n) - zw(k-1)
	    zb = zw(k) - pxyz(3,n)
	    dv = 1./((xu(i)-xu(i-1))*(yu(j)-yu(j-1))
     &              *(zw(k)-zw(k-1)))
            xeyszb = xe*ys*zb*dv
            xwyszb = xw*ys*zb*dv
            xeysza = xe*ys*za*dv
            xwysza = xw*ys*za*dv
            xeynzb = xe*yn*zb*dv
            xwynzb = xw*yn*zb*dv
            xeynza = xe*yn*za*dv
            xwynza = xw*yn*za*dv
            ww=  adv_vbu(i-1,k-1,j)*xeyszb   + adv_vbu(i,k-1,j)*xwyszb
     &          +adv_vbu(i-1,k  ,j)*xeysza   + adv_vbu(i,k  ,j)*xwysza
     &          +adv_vbu(i-1,k-1,j-1)*xeynzb + adv_vbu(i,k-1,j-1)*xwynzb
     &          +adv_vbu(i-1,k  ,j-1)*xeynza + adv_vbu(i,k  ,j-1)*xwynza
	  else
	    if (k-2 .eq. 0) then
	      za = pxyz(3,n)
	      zb = zw(k-1) - pxyz(3,n)
	      dv = 1./((xu(i)-xu(i-1))*(yu(j)-yu(j-1))
     &              *(zw(k-1)-0.0))
	    else
	      za = pxyz(3,n) - zw(k-2)
	      zb = zw(k-1) - pxyz(3,n)
	      dv = 1./((xu(i)-xu(i-1))*(yu(j)-yu(j-1))
     &              *(zw(k-1)-zw(k-2)))
	    endif
            xeyszb = xe*ys*zb*dv
            xwyszb = xw*ys*zb*dv
            xeysza = xe*ys*za*dv
            xwysza = xw*ys*za*dv
            xeynzb = xe*yn*zb*dv
            xwynzb = xw*yn*zb*dv
            xeynza = xe*yn*za*dv
            xwynza = xw*yn*za*dv
            ww  =adv_vbu(i-1,k-2,j)*xeyszb   + adv_vbu(i,k-2,j)*xwyszb
     &          +adv_vbu(i-1,k-1,j)*xeysza   + adv_vbu(i,k-1,j)*xwysza
     &          +adv_vbu(i-1,k-2,j-1)*xeynzb + adv_vbu(i,k-2,j-1)*xwynzb
     &          +adv_vbu(i-1,k-1,j-1)*xeynza + adv_vbu(i,k-1,j-1)*xwynza
	  endif

c       if (my_pe==0.and.n==1) then
c        print*,' Float #1:',pxyz(1,1)-360.,pxyz(2,1),pxyz(3,1)
c        print*,' Float #1:',uu,vv,ww
c       endif
c
c-----------------------------------------------------------------------
c         remember where the particle was
c-----------------------------------------------------------------------
c
          xold = pxyz(1,n)
          yold = pxyz(2,n)
          zold = pxyz(3,n)
c
c-----------------------------------------------------------------------
c         integrate the particle trajectory forward for one time step
c         taking convergence of meridians into account.
c-----------------------------------------------------------------------
c
          rcos      = cmdeg/cos(pxyz(2,n)/180.*pi)
          pxyz(1,n) = pxyz(1,n) + dt*uu*rcos
          pxyz(2,n) = pxyz(2,n) + dt*vv*cmdeg
          pxyz(3,n) = pxyz(3,n) - dt*ww
c
c-----------------------------------------------------------------------
c         update bottom most northeast index of bounding volume
c-----------------------------------------------------------------------
c
          if (pxyz(1,n) .ge. xu(i)) then
            pijk(1,n) = i + 1
          else if (pxyz(1,n) .lt. xu(i-1)) then
            pijk(1,n) = i - 1
          endif
c
          if (pxyz(2,n) .ge. yu(j)) then
            pijk(2,n) = j + 1
          else if (pxyz(2,n) .lt. yu(j-1)) then
            pijk(2,n) = j - 1
          endif
c
          if (pxyz(3,n) .ge. zt(k)) then
            pijk(3,n) = k + 1
          else if (pxyz(3,n) .lt. zt(k-1)) then
            pijk(3,n) = k - 1
          endif
c
c-----------------------------------------------------------------------
c         do not allow any component of the trajectory to enter
c         land. If it does, reset it to its previous value
c         thereby simulating free slip conditions. hey... not perfect, 
c         but beats loosing particles in land. Also if the grid has
c         isolated "T,S" cells (ones where all eight surrounding
c         velocities are on land), replace references to "tmask(i,k,j)"
c         by umask(i,k,j) + umask(i-1,k,j) + umask(i,k,j-1) + 
c         umask(i-1,k,j-1) to prevent stagnation of particles if this
c         is a problem.
c-----------------------------------------------------------------------
c
          if (cyclic) then

           if (pijk(1,n) .gt. imt) then
c            if (tmask(3,k,j) .ne. 0) then
            if (kmt_big(3,j) >=k ) then
              pijk(1,n) = pijk(1,n) - (imt-2)
              i         = pijk(1,n)
              pxyz(1,n) = xu(i-1) + (pxyz(1,n)-xu(imt))
            else
              pijk(1,n) = i
              pxyz(1,n) = xold
            endif
           endif

           if (pijk(1,n) .lt. 2) then
c            if (tmask(imt-2,k,j) .ne. 0) then
            if (kmt_big(imt-2,j) >= k) then
              pijk(1,n) = pijk(1,n) + (imt-2)
              i         = pijk(1,n)
              pxyz(1,n) = xu(i-1) + (pxyz(1,n)-xu(1))
            else
              pijk(1,n) = i
              pxyz(1,n) = xold
            endif
           endif

          endif
c
c-----------------------------------------------------------------------
c         constrain particles vertically to lie within ocean
c-----------------------------------------------------------------------
c
          if (pijk(3,n) .ne. k) then
            if (pijk(3,n) .eq. 1 .or. pijk(3,n) .gt. km) then
              pxyz(3,n) = zold
              pijk(3,n) = k
            else if (pijk(3,n) .gt. k .and. tmask(i,k+1,j) ==0.) then
              pxyz(3,n) = zold
              pijk(3,n) = k
            endif
          endif
c
c-----------------------------------------------------------------------
c         constrain particles longitudinally to stay within ocean 
c-----------------------------------------------------------------------
c
          if (pijk(1,n) .ne. i) then
            if (pijk(1,n) .gt. i .and. tmask(i+1,k,j) ==0.) then
              pxyz(1,n) = xold
              pijk(1,n) = i
            else if (pijk(1,n) .lt. i .and. tmask(i-1,k,j) ==0.) then
              pxyz(1,n) = xold
              pijk(1,n) = i
            endif 
          endif 
c
c-----------------------------------------------------------------------
c         constrain particles latitudinally to stay within ocean 
c-----------------------------------------------------------------------
c
          if (pijk(2,n) .ne. j) then
            if (pijk(2,n) .gt. j .and. tmask(i,k,j+1) == 0. ) then
              pxyz(2,n) = yold
              pijk(2,n) = j
            else if (pijk(2,n) .lt. j .and. tmask(i,k,j-1) == 0.)
     &        then
              pxyz(2,n) = yold
              pijk(2,n) = j
            endif 
          endif 

c
c    potential density
c
         if (snapshot_time_step) then
          if (set_dens_time_step_for_floats) then
           m=tau
          else
           m=dens_time_step_for_floats
          endif
          i=pijk(1,n)
          j=pijk(2,n)
          k=pijk(3,n)
          pdens(n)=model_dens_scalar(t(i,k,j,1,m),t(i,k,j,2,m),1)
         endif

        endif ! inside model domain

      enddo ! nptraj
c
c-----------------------------------------------------------------------
c     exchange informations around the PEs
c-----------------------------------------------------------------------
c
      do n=1,nptraj

        if (.not. ptdone(n) ) then
          pxyz(:,n)=0.
          pijk(:,n)=0
          pdens(n)=0.
        endif

        call global_lor(ptdone(n))

        if (.not. ptdone(n) ) then

          if (my_pe==0 .and. .not. ptlost(n)  ) 
     &          print*,' Oops, float #',n,' was lost'

          ptlost(n) = .true.
          pxyz(:,n) = spval
          pijk(:,n) = -99

        else

          call global_sum_vec(pxyz(:,n),3)
          call global_sum_int(pijk(1,n))
          call global_sum_int(pijk(2,n))
          call global_sum_int(pijk(3,n))
          call global_sum_vec(pdens(n),1)
c
c         check positions to be safe, however no floats should be lost
c
          if (pxyz(1,n)<xu(1)  -dxudeg(1) .or.
     &        pxyz(1,n)>xu(imt)+dxudeg(imt)) then
           if (my_pe==0) then
            print*,' Oops, float #',n,' dropped outside domain'
            print*,' (zonally, positions : ',pxyz(:,n),' )'
           endif
           ptlost(n)= .true.; pxyz(:,n)=spval; pijk(:,n)=-99
          endif

          if ((pxyz(2,n) < yu(1)  -dyudeg(1) .or.
     &         pxyz(2,n)>yu(jmt)+dyudeg(jmt)).and.
     &        .not. ptlost(n) ) then
           if (my_pe==0) then
            print*,' Oops, float #',n,' dropped outside domain'
            print*,' (meridionally, positions : ',pxyz(:,n),' )'
           endif
           ptlost(n)= .true.; pxyz(:,n)=spval; pijk(:,n)=-99
          endif

          if ((pxyz(3,n) < zt(1) -dzt(1) .or.
     &         pxyz(3,n) > zt(km)+dzt(km))   .and.
     &        .not. ptlost(n) ) then
           if (my_pe==0) then
            print*,' Oops, float #',n,' dropped outside domain'
            print*,' (vertically, positions : ',pxyz(:,n),' )'
           endif
           ptlost(n)= .true.; pxyz(:,n)=spval; pijk(:,n)=-99
          endif

        endif

      enddo
c
c     Output to either NetCDF or binary file
c
      if (snapshot_time_step) then

       za=0.
       do n=1,nptraj
        if (ptlost(n)) za=za+1.
       enddo
       za=za/nptraj*100.

       if (my_pe==0) then
        print*,' --> writing float data to file ',
     &        diag_float_file(1:len_trim(diag_float_file))
        print*,'    (# of floats: ',nptraj,' ',za,'% already lost)'
       endif

#ifdef netcdf_diagnostics
c
c       write trajectories to NetCDF file
c
       if (my_pe==0) then
        iret=nf_open(diag_float_file,NF_WRITE,ncid)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        iret=nf_inq_varid(ncid,'lon',xid)
        iret=nf_inq_varid(ncid,'lat',yid)
        iret=nf_inq_varid(ncid,'depth',zid)
        iret=nf_inq_varid(ncid,'rho',rhoid)
        iret=nf_inq_varid(ncid,'Time',timeid)
        iret=nf_inq_dimid(ncid,'Time',timedim)
        iret=nf_inq_dimlen(ncid, timedim,k)
        k=k+1
        call read_stamp(current_stamp,za)
        print*,' at stamp=',current_stamp,
     &          ' (days since origin : ',za,')',
     &          ' (time steps in file : ',k,')'
        v4(1)=za
        iret= nf_put_vara_real (ncid,timeid, k,1,v4)
        corner=(/1,k/); edges = (/nptraj,1/)
        v4=pxyz(1,:)-360.
        iret= nf_put_vara_real (ncid,xid,corner,edges,v4)
        v4=pxyz(2,:)
        iret= nf_put_vara_real (ncid,yid,corner,edges,v4)
        v4=pxyz(3,:)/(-100.)
        iret= nf_put_vara_real (ncid,zid,corner,edges,v4)
        v4=pdens(:)*1000.0
        iret= nf_put_vara_real (ncid,rhoid,corner,edges,v4)
        call ncclos (ncid, iret)
       endif 

#else
c
c      write trajectories to binary file
c
       if (my_pe==0) then
        call getunit(io,diag_float_file,'usa ieee')
        write (io) current_stamp, itt
        write (io) nptraj
        write (io) pxyz
        close (io)
       endif
#endif
      endif
c
c     write a restart which can be read at the next run
c
      if (last_time_step) call write_float_restart()
      end subroutine diag_float



      subroutine write_float_restart
      use spflame_module
      use diag_float_module
      implicit none
      integer :: io
      character (len=80) :: name
      write(name,'("floats_",i3,".dta")') sub_domain
      call replace_space_zero(name)
      if (my_pe==0) then
        print*,' writing float data to restart file ',
     &        name(1:len_trim(name))
        call get_free_iounit(io)
        open(io,file=name,form='unformatted',status='unknown')
        write (io) current_stamp,itt
        write (io) nptraj
        write (io) pxyz,pijk,pvolu,pvolv
        close (io)
      endif
      end subroutine write_float_restart
