#include "options.inc"


c=======================================================================
c       Integrate particles
c       linked in the code in driver only
c=======================================================================

#ifdef enable_diag_particle

c#define two_dim
c#define  remove_zonal_mean
c
c write out also buoyancy and pressure
#define  particle_buoyancy

c
c
c#define  particle_averages

      module particle_module
       use time_type_module
      implicit none
      integer :: nptraj 
      real, allocatable  :: pxyz(:,:),puvw(:,:),ptr(:,:)
      integer, allocatable :: pijk(:,:)
      real, allocatable  :: mu(:),muu(:)
      real, allocatable  :: mv(:),mvv(:)
      real, allocatable  :: mw(:),mww(:)
      integer :: counter
 
      logical, allocatable :: particle_active(:)
      integer, allocatable :: particle_pe(:)
      integer :: mmax
      type( time_type )          :: float_snapshot_time   
      type( time_type )          :: float_snapshot_inctime
      logical :: float_snapshot_time_step   = .false.
      real :: float_snap_int = 0.05
      real, allocatable  :: ulen(:,:),vlen(:,:),wlen(:,:)
      end module particle_module


      subroutine init_particle
c=======================================================================
c      initialize the particles
c=======================================================================
      use cpflame_module
      use particle_module
      implicit none 
#include "netcdf.inc"
      integer :: n,iret,dims(2),xid,yid,zid,uid,vid,wid,bid,pid
      integer :: ncid,Timedim,Timeid,pdim,uuid,vvid,wwid
      character :: name*24, unit*16
      integer :: i,j,k,indp,io,js,je
      logical :: file_exists = .false.
    
      if (my_pe==0) then
       print*,''
       print*,' Initializing particle module '
       print*,''
       if (n_pes>1) print*,'WARNNING: floats do not work with PE>1'
       print*,''
      endif

      mmax = 20
      if (my_pe==0) 
     &  print*,' Euler forward time step for particles ',dt_in/mmax,' s'

c-----------------------------------------------------------------------
c    read restart file
c-----------------------------------------------------------------------
       if (my_pe==0) then
       inquire( FILE='restart_float.dta', EXIST=file_exists )
       if (file_exists) then 
        print*,' reading data from restart file '
        call get_free_iounit(io)
        open(io,file='restart_float.dta',form='unformatted',
     &         status='old')
        read(io) nptraj
        allocate( pxyz(3,nptraj), pijk(3,nptraj),puvw(3,nptraj) )
        print*,' number of particles = ',nptraj
        do n=1,nptraj
         read(io) pijk(1:3,n),pxyz(1:3,n),puvw(1:3,n)
        enddo
        close(io)
        print*,' done reading restart'
       else
c-----------------------------------------------------------------------
c     initialize new set of particles 
c-----------------------------------------------------------------------
        print*,' cannot read restart file '
        nptraj=500
        allocate( pxyz(3,nptraj), pijk(3,nptraj),puvw(3,nptraj) )
        call seed_particles_randomly
c        call seed_particles_regular
        do n=1,nptraj
          i = indp (pxyz(1,n), xu, imt)
          if (xu(i) > pxyz(1,n)) then
            pijk(1,n) = i
          else
            pijk(1,n) = i+1
          endif
          j = indp (pxyz(2,n), yu, jmt)
          if (yu(j) > pxyz(2,n)) then
            pijk(2,n) = j
          else
           pijk(2,n) = j+1
          endif
          k = indp (pxyz(3,n), zw, km)
          if (zw(k) > pxyz(3,n)) then
           pijk(3,n) = k
          else
           pijk(3,n) = k+1
          endif
        enddo
       endif
      endif
c-----------------------------------------------------------------------
c    communicate particles to other pes
c-----------------------------------------------------------------------
      call barrier
      call bcast_integer(nptraj,1,0)
      if (my_pe/=0) 
     &      allocate( pxyz(3,nptraj), pijk(3,nptraj) ,puvw(3,nptraj))
      allocate(particle_active(nptraj),particle_pe(nptraj) )
      call bcast_real(pxyz,3*nptraj,0)
      call bcast_integer(pijk,3*nptraj,0)
      puvw=0.
#ifdef particle_buoyancy
      allocate( ptr(2,nptraj) ); ptr=0.
#endif

#ifdef  particle_averages
      allocate( mu(nptraj) ); mu=0.
      allocate( muu(nptraj) ); muu=0.
      allocate( mv(nptraj) ); mv=0.
      allocate( mvv(nptraj) ); mvv=0.
      allocate( mw(nptraj) ); mw=0.
      allocate( mww(nptraj) ); mww=0.
      counter=0
#endif

c-----------------------------------------------------------------------
c     define output file
c-----------------------------------------------------------------------
      if (my_pe==0) then
       ncid = nccre ('float.cdf', NCCLOB, iret)
       iret=nf_set_fill(ncid, NF_NOFILL, iret)
       Timedim = ncddef(ncid, 'Time', nf_unlimited, iret)
       pdim    = ncddef(ncid, 'Number', nptraj, iret)
       timeid  = ncvdef (ncid,'Time', NCFLOAT,1,timedim,iret)
       dims = (/pdim,Timedim/)
       xid  = ncvdef (ncid,'x_pos', NCFLOAT,2,dims,iret)
       yid  = ncvdef (ncid,'y_pos', NCFLOAT,2,dims,iret)
       zid  = ncvdef (ncid,'z_pos', NCFLOAT,2,dims,iret)
       uid  = ncvdef (ncid,'u', NCFLOAT,2,dims,iret)
       vid  = ncvdef (ncid,'v', NCFLOAT,2,dims,iret)
       wid  = ncvdef (ncid,'w', NCFLOAT,2,dims,iret)
#ifdef  particle_averages
       uuid  = ncvdef (ncid,'uu', NCFLOAT,2,dims,iret)
       vvid  = ncvdef (ncid,'vv', NCFLOAT,2,dims,iret)
       wwid  = ncvdef (ncid,'ww', NCFLOAT,2,dims,iret)
#endif
#ifdef particle_buoyancy
       pid  = ncvdef (ncid,'p', NCFLOAT,2,dims,iret)
       bid  = ncvdef (ncid,'b', NCFLOAT,2,dims,iret)
#endif
       name = 'Time '; unit = 'days'
       call ncaptc(ncid, timeid, 'long_name', NCCHAR, 24, name, iret) 
       call ncaptc(ncid, timeid, 'units',     NCCHAR, 16, unit, iret) 
       call ncaptc(ncid, Timeid,'time_origin',NCCHAR, 20,
     &  '01-JAN-1900 00:00:00', iret)
       name = 'X position'; unit = 'm'
       call dvcdf(ncid,xid,name,24,unit,16,spval)
       name = 'Y position'; unit = 'm'
       call dvcdf(ncid,yid,name,24,unit,16,spval)
       name = 'Z position'; unit = 'm'
       call dvcdf(ncid,zid,name,24,unit,16,spval)
       name = 'zonal velocity'; unit = 'm/s'
       call dvcdf(ncid,uid,name,24,unit,16,spval)
       name = 'meridional velocity'; unit = 'm/s'
       call dvcdf(ncid,vid,name,24,unit,16,spval)
       name = 'vertical velocity'; unit = 'm/s'
       call dvcdf(ncid,wid,name,24,unit,16,spval)
#ifdef particle_buoyancy
       name = 'buoyancy'; unit = 'm/s^2'
       call dvcdf(ncid,bid,name,24,unit,16,spval)
       name = 'pressure'; unit = 'm^2/s'
       call dvcdf(ncid,pid,name,24,unit,16,spval)
#endif
#ifdef  particle_averages
       name = 'zonal velocity variance'; unit = 'm^2/s^2'
       call dvcdf(ncid,uuid,name,24,unit,16,spval)
       name = 'meridional velocity variance'; unit = 'm^2/s^2'
       call dvcdf(ncid,vvid,name,24,unit,16,spval)
       name = 'vertical velocity variance'; unit = 'm^2/s^2'
       call dvcdf(ncid,wwid,name,24,unit,16,spval)
#endif
       call ncclos (ncid, iret)
      endif

#ifdef remove_zonal_mean
      allocate(ulen(jmt,km),vlen(jmt,km),wlen(jmt,km))
      ulen=0;vlen=0;wlen=0
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do j=js,je
       do i=2,imt-1
        ulen(j,:)=ulen(j,:)+dx*maskU(i,j,:)
        vlen(j,:)=vlen(j,:)+dx*maskV(i,j,:)
        wlen(j,:)=wlen(j,:)+dx*maskW(i,j,:)
       enddo 
      enddo 
#endif

      if (my_pe==0)  then
       print*,' float snapshot interval is ', float_snap_int*86400,'s'
       print*,' Done initializing time averaging module'
       print*,' '
      endif
      end subroutine init_particle


      function ilim(i)
      use cpflame_module
      implicit none 
      integer :: i ,ilim
      if (enable_cyclic_x) then
        ilim=i
        if (i>imt) ilim=i-(imt-2)
        if (i<1) ilim=i+(imt-2)
      else
        ilim=max(1,min(imt,i))
      endif
      end function ilim

      function jlim(j)
      use cpflame_module
      implicit none 
      integer :: j ,jlim
      if (enable_cyclic_y) then
        jlim=j
        if (j>jmt) jlim=j-(jmt-2)
        if (j<1) jlim=j+(jmt-2)
      else
        jlim=max(1,min(jmt,j))
      endif
      end function jlim

      function klim(k)
      use cpflame_module
      implicit none 
      integer :: k ,klim
      klim=max(1,min(km,k))
      end function klim

 
      subroutine integrate_particle
c=======================================================================
c       integrate particles
c=======================================================================
      use cpflame_module
      use particle_module
      implicit none 
      integer :: i,j,k,n,m,js,je
      real :: xe,xw,yn,ys,zu,zl,dv,th,tf
      real :: xeyszu,xeyszl,xeynzu,xeynzl
      real :: xwyszu,xwyszl,xwynzu,xwynzl
      real :: uuh,vvh,wwh,xold,yold,zold
      real :: uuf,vvf,wwf,uu,vv,ww
      real :: bbf,bbh
#include "netcdf.inc"
      integer :: ncid,iret,xid,yid,zid,timedim,ilen,timeid,io
      integer :: uid,vid,wid,bid,pid,uuid,vvid,wwid
      real :: fxa
      type(time_type) :: time
      logical, save :: first = .true.
#ifdef remove_zonal_mean
      real :: um(jmt,km),vm(jmt,km),wm(jmt,km)
      real :: um2(jmt,km),vm2(jmt,km),wm2(jmt,km)
#endif
      integer, external :: ilim,jlim,klim
      if (first) then
c-----------------------------------------------------------------------
c      set next snapshot time
c-----------------------------------------------------------------------
       float_snapshot_time = initial_time
       call set_time(float_snapshot_inctime ,
     &         (float_snap_int-int(float_snap_int))*60*60*24 , 
     &        int(float_snap_int))
       if (float_snap_int>0.) then
 20   float_snapshot_time=float_snapshot_time+float_snapshot_inctime
        if (float_snapshot_time < current_time ) goto 20
       else
        float_snapshot_time=current_time
       endif

      endif
c-----------------------------------------------------------------------
c      is this a float snapshot time step ?
c-----------------------------------------------------------------------
       if ( current_time >= float_snapshot_time .and. 
     &      .not. eulerback1  ) then 
        float_snapshot_time_step = .true.
        float_snapshot_time = float_snapshot_time 
     &        + float_snapshot_inctime
        if (float_snapshot_time <= current_time) then
         float_snapshot_time = current_time
        endif
       else
        float_snapshot_time_step = .false.
       endif

#ifdef remove_zonal_mean
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      um=0.0;vm=0.0;wm=0.0; um2=0.0;vm2=0.0;wm2=0.0
      do j=js,je 
        do i=2,imt-1
         um(j,:)=um(j,:)+u(i,j,:,1,tau)*dx*maskU(i,j,:) 
         vm(j,:)=vm(j,:)+u(i,j,:,2,tau)*dx*maskV(i,j,:) 
         wm(j,:)=wm(j,:)+u(i,j,:,3,tau)*dx*maskW(i,j,:) 
         um2(j,:)=um2(j,:)+u(i,j,:,1,taup1)*dx*maskU(i,j,:) 
         vm2(j,:)=vm2(j,:)+u(i,j,:,2,taup1)*dx*maskV(i,j,:) 
         wm2(j,:)=wm2(j,:)+u(i,j,:,3,taup1)*dx*maskW(i,j,:) 
        enddo 
      enddo 
      where (ulen/=0); um=um/ulen; um2=um2/ulen; else where; 
      um=spval ; um2=spval; end where
      where (vlen/=0); vm=vm/vlen; vm2=vm2/vlen; else where; 
      vm=spval ; vm2=spval; end where
      where (wlen/=0); wm=wm/wlen; wm2=wm2/wlen; else where; 
      wm=spval ; wm2=spval; end where
      call border_exchg_merid(um,1); call border_exchg_merid(vm,1);
      call border_exchg_merid(wm,1); call border_exchg_merid(um2,1);
      call border_exchg_merid(vm2,1); call border_exchg_merid(wm2,1);
#endif


      do m=1,mmax
       call particle_inside_pe_domain()
       do n=1,nptraj
        th= (m-1.)/mmax;  tf = 1.-th;
        if (particle_active(n)) then
#ifdef  particle_buoyancy
c-----------------------------------------------------------------------
c        interpolate buoyancy on particle
c-----------------------------------------------------------------------
          i  = pijk(1,n); j  = pijk(2,n); k  = pijk(3,n)
          if (pxyz(1,n) > xt(i) )  i=i+1
          xe = (xt(i) - pxyz(1,n)); xw = (pxyz(1,n)-xt(i-1))
          if (pxyz(2,n) > yt(j) )  j=j+1
          yn = (yt(j) - pxyz(2,n)); ys = (pxyz(2,n) - yt(j-1))
          if (pxyz(3,n) > zt(k) )  k=k+1
          zu = (zt(k) - pxyz(3,n)); zl = (pxyz(3,n)-zt(k-1))

          if (maskT(ilim(i+1),j,k)==0) then; xe=dx;xw=0; endif
          if (maskT(i-1,j,k)==0) then; xe=0;xw=dx; endif
          if (maskT(i,jlim(j+1),k)==0) then; yn=dx;ys=0; endif
          if (maskT(i,j-1,k)==0) then; yn=0;ys=dx; endif
          if (maskT(i,j,klim(k+1))==0) then; zu=dz;zl=0; endif
          if (maskT(i,j,k-1)==0) then; zu=0;zl=dz; endif

          dv = 1./(dx*dx*dz)
          xeyszu = xe*ys*zu*dv; xwyszu = xw*ys*zu*dv
          xeyszl = xe*ys*zl*dv; xwyszl = xw*ys*zl*dv
          xeynzu = xe*yn*zu*dv; xwynzu = xw*yn*zu*dv
          xeynzl = xe*yn*zl*dv; xwynzl = xw*yn*zl*dv

          bbh=b(i-1,j  ,k-1,tau  )*xeyszu+b(i,j  ,k-1,tau  )*xwyszu
     &       +b(i-1,j  ,k  ,tau  )*xeyszl+b(i,j  ,k  ,tau  )*xwyszl
     &       +b(i-1,j-1,k-1,tau  )*xeynzu+b(i,j-1,k-1,tau  )*xwynzu
     &       +b(i-1,j-1,k  ,tau  )*xeynzl+b(i,j-1,k  ,tau  )*xwynzl
          bbf=b(i-1,j  ,k-1,taup1)*xeyszu+b(i,j  ,k-1,taup1)*xwyszu
     &       +b(i-1,j  ,k  ,taup1)*xeyszl+b(i,j  ,k  ,taup1)*xwyszl
     &       +b(i-1,j-1,k-1,taup1)*xeynzu+b(i,j-1,k-1,taup1)*xwynzu
     &       +b(i-1,j-1,k  ,taup1)*xeynzl+b(i,j-1,k  ,taup1)*xwynzl
          ptr(1,n)=th*bbh+tf*bbf

          ptr(2,n) =p_full(i-1,j  ,k-1,tau  )*xeyszu
     & +p_full(i,j  ,k-1,tau  )*xwyszu
     &       +p_full(i-1,j  ,k  ,tau  )*xeyszl
     & +p_full(i,j  ,k  ,tau  )*xwyszl
     &       +p_full(i-1,j-1,k-1,tau  )*xeynzu
     & +p_full(i,j-1,k-1,tau  )*xwynzu
     &       +p_full(i-1,j-1,k  ,tau  )*xeynzl
     & +p_full(i,j-1,k  ,tau  )*xwynzl

#endif
c-----------------------------------------------------------------------
c            pijk gives tracer box of the particle, 
c            find u-box and distances to borders, account for free slip
c-----------------------------------------------------------------------
          i  = pijk(1,n); j  = pijk(2,n); k  = pijk(3,n)
          xe = (xu(i) - pxyz(1,n)); xw = (pxyz(1,n)-xu(i-1))
          if (pxyz(2,n) > yt(j) )  j=j+1
          yn = (yt(j) - pxyz(2,n)); ys = (pxyz(2,n) - yt(j-1))
          if (pxyz(3,n) > zt(k) )  k=k+1
          zu = (zt(k) - pxyz(3,n)); zl = (pxyz(3,n)-zt(k-1))

          if (maskT(i,jlim(j+1),k)==0) then; yn=dx;ys=0; endif
          if (maskT(i,j-1,k)==0) then; yn=0;ys=dx; endif
          if (maskT(i,j,klim(k+1))==0) then; zu=dz;zl=0; endif
          if (maskT(i,j,k-1)==0) then; zu=0;zl=dz; endif

	  dv = 1./(dx*dx*dz)
          xeyszu = xe*ys*zu*dv; xwyszu = xw*ys*zu*dv
          xeyszl = xe*ys*zl*dv; xwyszl = xw*ys*zl*dv
          xeynzu = xe*yn*zu*dv; xwynzu = xw*yn*zu*dv
          xeynzl = xe*yn*zl*dv; xwynzl = xw*yn*zl*dv
c-----------------------------------------------------------------------
c           interpolate u on particle position
c-----------------------------------------------------------------------
#ifdef remove_zonal_mean
          uuh=
     &  (u(i-1,j  ,k-1,1,tau  )-um(j,k-1))*xeyszu*maskU(i-1,j,k-1)
     & +(u(i,j  ,k-1,1,tau  )-um(j,k-1))*xwyszu*maskU(i,j,k-1)
     & +(u(i-1,j  ,k  ,1,tau  )-um(j,k))*xeyszl*maskU(i-1,j,k)
     & +(u(i,j  ,k  ,1,tau  )-um(j,k))*xwyszl*maskU(i,j,k)
     & +(u(i-1,j-1,k-1,1,tau  )-um(j-1,k-1))*xeynzu*maskU(i-1,j-1,k-1)
     & +(u(i,j-1,k-1,1,tau  )-um(j-1,k-1))*xwynzu*maskU(i,j-1,k-1)
     & +(u(i-1,j-1,k  ,1,tau  )-um(j-1,k))*xeynzl*maskU(i-1,j-1,k)
     & +(u(i,j-1,k  ,1,tau  )-um(j-1,k))*xwynzl*maskU(i,j-1,k)
          uuf=
     &  (u(i-1,j  ,k-1,1,taup1)-um2(j,k-1))*xeyszu*maskU(i-1,j,k-1)
     & +(u(i,j  ,k-1,1,taup1)-um2(j,k-1))*xwyszu*maskU(i,j,k-1)
     & +(u(i-1,j  ,k  ,1,taup1)-um2(j,k))*xeyszl*maskU(i-1,j,k)
     & +(u(i,j  ,k  ,1,taup1)-um2(j,k))*xwyszl*maskU(i,j,k)
     & +(u(i-1,j-1,k-1,1,taup1)-um2(j-1,k-1))*xeynzu*maskU(i-1,j-1,k-1)
     & +(u(i,j-1,k-1,1,taup1)-um2(j-1,k-1))*xwynzu*maskU(i,j-1,k-1)
     & +(u(i-1,j-1,k  ,1,taup1)-um2(j-1,k))*xeynzl*maskU(i-1,j-1,k)
     & +(u(i,j-1,k  ,1,taup1)-um2(j-1,k))*xwynzl*maskU(i,j-1,k)
#else
          uuh=u(i-1,j  ,k-1,1,tau  )*xeyszu+u(i,j  ,k-1,1,tau  )*xwyszu
     &       +u(i-1,j  ,k  ,1,tau  )*xeyszl+u(i,j  ,k  ,1,tau  )*xwyszl
     &       +u(i-1,j-1,k-1,1,tau  )*xeynzu+u(i,j-1,k-1,1,tau  )*xwynzu
     &       +u(i-1,j-1,k  ,1,tau  )*xeynzl+u(i,j-1,k  ,1,tau  )*xwynzl
          uuf=u(i-1,j  ,k-1,1,taup1)*xeyszu+u(i,j  ,k-1,1,taup1)*xwyszu
     &       +u(i-1,j  ,k  ,1,taup1)*xeyszl+u(i,j  ,k  ,1,taup1)*xwyszl
     &       +u(i-1,j-1,k-1,1,taup1)*xeynzu+u(i,j-1,k-1,1,taup1)*xwynzu
     &       +u(i-1,j-1,k  ,1,taup1)*xeynzl+u(i,j-1,k  ,1,taup1)*xwynzl
#endif
          uu=th*uuh+tf*uuf
c-----------------------------------------------------------------------
c           find v-box and distances to borders, account for free slip
c-----------------------------------------------------------------------
          i  = pijk(1,n); j  = pijk(2,n); k  = pijk(3,n)
          if (pxyz(1,n) > xt(i) )  i=i+1
          xe = (xt(i) - pxyz(1,n)); xw = (pxyz(1,n)-xt(i-1))
          yn = (yu(j) - pxyz(2,n)); ys = (pxyz(2,n) - yu(j-1))
          if (pxyz(3,n) > zt(k) )  k=k+1
          zu = (zt(k) - pxyz(3,n)); zl = (pxyz(3,n)-zt(k-1))

          if (maskT(ilim(i+1),j,k)==0) then; xe=dx;xw=0; endif
          if (maskT(i-1,j,k)==0) then; xe=0;xw=dx; endif
          if (maskT(i,j,klim(k+1))==0) then; zu=dz;zl=0; endif
          if (maskT(i,j,k-1)==0) then; zu=0;zl=dz; endif

          dv = 1./(dx*dx*dz)
          xeyszu = xe*ys*zu*dv; xwyszu = xw*ys*zu*dv
          xeyszl = xe*ys*zl*dv; xwyszl = xw*ys*zl*dv
          xeynzu = xe*yn*zu*dv; xwynzu = xw*yn*zu*dv
          xeynzl = xe*yn*zl*dv; xwynzl = xw*yn*zl*dv
c-----------------------------------------------------------------------
c           interpolate v on particle position
c-----------------------------------------------------------------------
#ifdef remove_zonal_mean
          vvh=
     &   (u(i-1,j  ,k-1,2,tau  )-vm(j,k-1))*xeyszu*maskV(i-1,j,k-1)
     & +(u(i,j  ,k-1,2,tau  )-vm(j,k-1))*xwyszu*maskV(i,j,k-1)
     & +(u(i-1,j  ,k  ,2,tau  )-vm(j,k))*xeyszl*maskV(i-1,j,k)
     & +(u(i,j  ,k  ,2,tau  )-vm(j,k))*xwyszl*maskV(i,j,k)
     & +(u(i-1,j-1,k-1,2,tau  )-vm(j-1,k-1))*xeynzu*maskV(i-1,j-1,k-1)
     & +(u(i,j-1,k-1,2,tau  )-vm(j-1,k-1))*xwynzu*maskV(i,j-1,k-1)
     & +(u(i-1,j-1,k  ,2,tau  )-vm(j-1,k))*xeynzl*maskV(i-1,j-1,k)
     & +(u(i,j-1,k  ,2,tau  )-vm(j-1,k))*xwynzl*maskV(i,j-1,k)
          vvf=
     &  (u(i-1,j  ,k-1,2,taup1)-vm2(j,k-1))*xeyszu*maskV(i-1,j,k-1)
     & +(u(i,j  ,k-1,2,taup1)-vm2(j,k-1))*xwyszu*maskV(i,j,k-1)
     & +(u(i-1,j  ,k  ,2,taup1)-vm2(j,k))*xeyszl*maskV(i-1,j,k)
     & +(u(i,j  ,k  ,2,taup1)-vm2(j,k))*xwyszl*maskV(i,j,k)
     & +(u(i-1,j-1,k-1,2,taup1)-vm2(j-1,k-1))*xeynzu*maskV(i-1,j-1,k-1)
     & +(u(i,j-1,k-1,2,taup1)-vm2(j-1,k-1))*xwynzu*maskV(i,j-1,k-1)
     & +(u(i-1,j-1,k  ,2,taup1)-vm2(j-1,k))*xeynzl*maskV(i-1,j-1,k)
     & +(u(i,j-1,k  ,2,taup1)-vm2(j-1,k))*xwynzl*maskV(i,j-1,k)
#else
          vvh=u(i-1,j  ,k-1,2,tau  )*xeyszu+u(i,j  ,k-1,2,tau  )*xwyszu
     &       +u(i-1,j  ,k  ,2,tau  )*xeyszl+u(i,j  ,k  ,2,tau  )*xwyszl
     &       +u(i-1,j-1,k-1,2,tau  )*xeynzu+u(i,j-1,k-1,2,tau  )*xwynzu
     &       +u(i-1,j-1,k  ,2,tau  )*xeynzl+u(i,j-1,k  ,2,tau  )*xwynzl
          vvf=u(i-1,j  ,k-1,2,taup1)*xeyszu+u(i,j  ,k-1,2,taup1)*xwyszu
     &       +u(i-1,j  ,k  ,2,taup1)*xeyszl+u(i,j  ,k  ,2,taup1)*xwyszl
     &       +u(i-1,j-1,k-1,2,taup1)*xeynzu+u(i,j-1,k-1,2,taup1)*xwynzu
     &       +u(i-1,j-1,k  ,2,taup1)*xeynzl+u(i,j-1,k  ,2,taup1)*xwynzl
#endif
          vv=th*vvh+tf*vvf
c-----------------------------------------------------------------------
c           find w-box and distances to borders, account for free slip
c-----------------------------------------------------------------------
          i  = pijk(1,n); j  = pijk(2,n); k  = pijk(3,n)
          if (pxyz(1,n) > xt(i) )  i=i+1
          xe = (xt(i) - pxyz(1,n)); xw = (pxyz(1,n)-xt(i-1))
          if (pxyz(2,n) > yt(j) )  j=j+1
          yn = (yt(j) - pxyz(2,n)); ys = (pxyz(2,n) - yt(j-1))
          zu = (zw(k) - pxyz(3,n)); zl = (pxyz(3,n)-zw(k-1))

          if (maskT(ilim(i+1),j,k)==0) then; xe=dx;xw=0; endif
          if (maskT(i-1,j,k)==0) then; xe=0;xw=dx; endif
          if (maskT(i,jlim(j+1),k)==0) then; yn=dx;ys=0; endif
          if (maskT(i,j-1,k)==0) then; yn=0;ys=dx; endif

	  dv = 1./(dx*dx*dz)
          xeyszu = xe*ys*zu*dv; xwyszu = xw*ys*zu*dv
          xeyszl = xe*ys*zl*dv; xwyszl = xw*ys*zl*dv
          xeynzu = xe*yn*zu*dv; xwynzu = xw*yn*zu*dv
          xeynzl = xe*yn*zl*dv; xwynzl = xw*yn*zl*dv
c-----------------------------------------------------------------------
c           interpolate w on particle position
c-----------------------------------------------------------------------
#ifdef remove_zonal_mean
          wwh=
     &  (u(i-1,j  ,k-1,3,tau  )-wm(j,k-1))*xeyszu*maskW(i-1,j,k-1)
     & +(u(i,j  ,k-1,3,tau  )-wm(j,k-1))*xwyszu*maskW(i,j,k-1)
     & +(u(i-1,j  ,k  ,3,tau  )-wm(j,k))*xeyszl*maskW(i-1,j,k)
     & +(u(i,j  ,k  ,3,tau  )-wm(j,k))*xwyszl*maskW(i,j,k)
     & +(u(i-1,j-1,k-1,3,tau  )-wm(j-1,k-1))*xeynzu*maskW(i-1,j-1,k-1)
     & +(u(i,j-1,k-1,3,tau  )-wm(j-1,k-1))*xwynzu*maskW(i,j-1,k-1)
     & +(u(i-1,j-1,k  ,3,tau  )-wm(j-1,k))*xeynzl*maskW(i-1,j-1,k)
     & +(u(i,j-1,k  ,3,tau  )-wm(j-1,k))*xwynzl*maskW(i,j-1,k)
         if (.not. enable_hydrostatic ) then
          wwf=
     &  (u(i-1,j  ,k-1,3,taup1)-wm2(j,k-1))*xeyszu*maskW(i-1,j,k-1)
     & +(u(i,j  ,k-1,3,taup1)-wm2(j,k-1))*xwyszu*maskW(i,j,k-1)
     & +(u(i-1,j  ,k  ,3,taup1)-wm2(j,k))*xeyszl*maskW(i-1,j,k)
     & +(u(i,j  ,k  ,3,taup1)-wm2(j,k))*xwyszl*maskW(i,j,k)
     & +(u(i-1,j-1,k-1,3,taup1)-wm2(j-1,k-1))*xeynzu*maskW(i-1,j-1,k-1)
     & +(u(i,j-1,k-1,3,taup1)-wm2(j-1,k-1))*xwynzu*maskW(i,j-1,k-1)
     & +(u(i-1,j-1,k  ,3,taup1)-wm2(j-1,k))*xeynzl*maskW(i-1,j-1,k)
     & +(u(i,j-1,k  ,3,taup1)-wm2(j-1,k))*xwynzl*maskW(i,j-1,k)
#else
          wwh=u(i-1,j  ,k-1,3,tau  )*xeyszu+u(i,j  ,k-1,3,tau  )*xwyszu
     &       +u(i-1,j  ,k  ,3,tau  )*xeyszl+u(i,j  ,k  ,3,tau  )*xwyszl
     &       +u(i-1,j-1,k-1,3,tau  )*xeynzu+u(i,j-1,k-1,3,tau  )*xwynzu
     &       +u(i-1,j-1,k  ,3,tau  )*xeynzl+u(i,j-1,k  ,3,tau  )*xwynzl
         if (.not. enable_hydrostatic ) then
          wwf=u(i-1,j  ,k-1,3,taup1)*xeyszu+u(i,j  ,k-1,3,taup1)*xwyszu
     &       +u(i-1,j  ,k  ,3,taup1)*xeyszl+u(i,j  ,k  ,3,taup1)*xwyszl
     &       +u(i-1,j-1,k-1,3,taup1)*xeynzu+u(i,j-1,k-1,3,taup1)*xwynzu
     &       +u(i-1,j-1,k  ,3,taup1)*xeynzl+u(i,j-1,k  ,3,taup1)*xwynzl
#endif          
          ww=th*wwh+tf*wwf
         else
          ww=wwh
         endif
c-----------------------------------------------------------------------
c         integrate the particle trajectory forward for one time step
c-----------------------------------------------------------------------
          xold=pxyz(1,n); yold=pxyz(2,n); zold=pxyz(3,n)
          pxyz(1,n) = pxyz(1,n) + dt_in*uu/mmax
          pxyz(2,n) = pxyz(2,n) + dt_in*vv/mmax
          pxyz(3,n) = pxyz(3,n) + dt_in*ww/mmax
          puvw(1,n)=uu; puvw(2,n)=vv; puvw(3,n)=ww
c-----------------------------------------------------------------------
c           update index of bounding tracer volume
c-----------------------------------------------------------------------
          i  = pijk(1,n); j  = pijk(2,n); k  = pijk(3,n)
          if (pxyz(1,n) >= xu(i)) then
            pijk(1,n) = i + 1
          else if (pxyz(1,n) < xu(i-1)) then
            pijk(1,n) = i - 1
          endif
          if (enable_cyclic_x .and. pijk(1,n)==imt) then
             pijk(1,n)=2; pxyz(1,n)=pxyz(1,n)-dx*(imt-2)
          endif
          if (enable_cyclic_x .and. pijk(1,n)==1) then
            pijk(1,n)=imt-1; pxyz(1,n)=pxyz(1,n)+dx*(imt-2)
          endif

          if (pxyz(2,n) >= yu(j)) then
            pijk(2,n) = j + 1
          else if (pxyz(2,n) < yu(j-1)) then
            pijk(2,n) = j - 1
          endif
          if (enable_cyclic_y .and. pijk(2,n)==jmt) then
             pijk(2,n)=2; pxyz(2,n)=pxyz(2,n)-dx*(jmt-2)
          endif
          if (enable_cyclic_y .and. pijk(2,n)==1) then
            pijk(2,n)=jmt-1; pxyz(2,n)=pxyz(2,n)+dx*(jmt-2)
          endif

          if (pxyz(3,n) >= zw(k)) then
            pijk(3,n) = k + 1
          else if (pxyz(3,n) < zw(k-1)) then
            pijk(3,n) = k - 1
          endif
c-----------------------------------------------------------------------
c         constrain particles vertically to lie within ocean
c-----------------------------------------------------------------------
          if (pijk(3,n) /= k) then
            if (pijk(3,n) > k .and. maskT(i,j,k+1) ==0.) then
              pxyz(3,n) = zold
              pijk(3,n) = k
            endif
            if (pijk(3,n) < k .and. maskT(i,j,k-1) ==0.) then
              pxyz(3,n) = zold
              pijk(3,n) = k
            endif
          endif
c-----------------------------------------------------------------------
c         constrain particles longitudinally to stay within ocean 
c-----------------------------------------------------------------------
          if (pijk(1,n) /= i) then
            if (pijk(1,n) > i .and. maskT(i+1,j,k) ==0.) then
              pxyz(1,n) = xold
              pijk(1,n) = i
            else if (pijk(1,n) < i .and. maskT(i-1,j,k) ==0.) then
              pxyz(1,n) = xold
              pijk(1,n) = i
            endif 
          endif 
c-----------------------------------------------------------------------
c         constrain particles latitudinally to stay within ocean 
c-----------------------------------------------------------------------
          if (pijk(2,n) /= j) then
            if (pijk(2,n) > j .and. maskT(i,j+1,k) == 0. ) then
              pxyz(2,n) = yold
              pijk(2,n) = j
            else if (pijk(2,n) < j .and. maskT(i,j-1,k) == 0.) then
              pxyz(2,n) = yold
              pijk(2,n) = j
            endif 
          endif 
        endif ! particle_active
c-----------------------------------------------------------------------
c      average some quantitites
c-----------------------------------------------------------------------

#ifdef  particle_averages
        mu(n)=mu(n)+uu
        muu(n)=muu(n)+uu**2
        mv(n)=mv(n)+vv
        mvv(n)=mvv(n)+vv**2
        mw(n)=mw(n)+ww
        mww(n)=mww(n)+ww**2
#endif
       enddo ! nptraj
#ifdef  particle_averages
       counter=counter+1
#endif
       call particle_distribute
      enddo ! m
c-----------------------------------------------------------------------
c      write particle positions to netcdf file
c-----------------------------------------------------------------------
      if (float_snapshot_time_step.or.initial_time==current_time) then
#ifdef  particle_averages
         mu=mu/counter; muu=muu/counter-mu**2
         mv=mv/counter; mvv=mvv/counter-mv**2
         mw=mw/counter; mww=mww/counter-mw**2
#endif
        if (my_pe==0) then
         print*,'writing particles to file float.cdf'
         iret=nf_open('float.cdf',NF_WRITE,ncid)
         iret=nf_set_fill(ncid, NF_NOFILL, iret)
         iret=nf_inq_varid(ncid,'x_pos',xid)
         iret=nf_inq_varid(ncid,'y_pos',yid)
         iret=nf_inq_varid(ncid,'z_pos',zid)
         iret=nf_inq_varid(ncid,'u',uid)
         iret=nf_inq_varid(ncid,'v',vid)
         iret=nf_inq_varid(ncid,'w',wid)
         iret=nf_inq_varid(ncid,'uu',uuid)
         iret=nf_inq_varid(ncid,'vv',vvid)
         iret=nf_inq_varid(ncid,'ww',wwid)
#ifdef particle_buoyancy
         iret=nf_inq_varid(ncid,'b',bid)
         iret=nf_inq_varid(ncid,'p',pid)
#endif
         iret=nf_inq_dimid(ncid,'Time',timedim)
         iret=nf_inq_dimlen(ncid, timedim,ilen)
         iret=nf_inq_varid(ncid,'Time',timeid)
         ilen=ilen+1
         time = current_time-initial_time
         fxa = time%days + time%seconds/86400.
         iret= nf_put_vara_double(ncid,timeid,ilen,1,fxa)
         do n=1,nptraj
        iret= nf_put_vara_double(ncid,xid,(/n,ilen/),(/1,1/),pxyz(1,n))
        iret= nf_put_vara_double(ncid,yid,(/n,ilen/),(/1,1/),pxyz(2,n))
        iret= nf_put_vara_double(ncid,zid,(/n,ilen/),(/1,1/),pxyz(3,n))
#ifdef  particle_averages
        iret= nf_put_vara_double(ncid,uid,(/n,ilen/),(/1,1/),mu(n))
        iret= nf_put_vara_double(ncid,vid,(/n,ilen/),(/1,1/),mv(n))
        iret= nf_put_vara_double(ncid,wid,(/n,ilen/),(/1,1/),mw(n))
        iret= nf_put_vara_double(ncid,uuid,(/n,ilen/),(/1,1/),muu(n))
        iret= nf_put_vara_double(ncid,vvid,(/n,ilen/),(/1,1/),mvv(n))
        iret= nf_put_vara_double(ncid,wwid,(/n,ilen/),(/1,1/),mww(n))
#else
        iret= nf_put_vara_double(ncid,uid,(/n,ilen/),(/1,1/),puvw(1,n))
        iret= nf_put_vara_double(ncid,vid,(/n,ilen/),(/1,1/),puvw(2,n))
        iret= nf_put_vara_double(ncid,wid,(/n,ilen/),(/1,1/),puvw(3,n))
#endif
#ifdef particle_buoyancy
        iret= nf_put_vara_double(ncid,bid,(/n,ilen/),(/1,1/),ptr(1,n))
        iret= nf_put_vara_double(ncid,pid,(/n,ilen/),(/1,1/),ptr(2,n))
#endif
         enddo
         call ncclos (ncid, iret)
       endif
       call barrier
#ifdef  particle_averages
       counter=0; mu=0;muu=0;mv=0;mvv=0;mw=0;mww=0
#endif
      endif
c-----------------------------------------------------------------------
c    write restart file
c-----------------------------------------------------------------------
      if (last_time_step.and. my_pe==0) then
         print*,'writing restart file for particles'
         call get_free_iounit(io)
         open(io,file='restart_float.dta',form='unformatted',
     &          status='unknown')
         write(io) nptraj
         do n=1,nptraj
          write(io) pijk(1:3,n),pxyz(1:3,n),puvw(1:3,n)
         enddo
         close(io)
         print*,' done writing restart'
      endif
      first = .false.
      end subroutine integrate_particle


      subroutine particle_inside_pe_domain()
c=======================================================================
c       is particle inside domain of pe?
c=======================================================================
      use cpflame_module
      use particle_module
      implicit none 
      integer :: n

      do n=1,nptraj
       if (pijk(2,n)>=js_pe .and. pijk(2,n)<=je_pe) then
        particle_active(n) = .true. 
        particle_pe(n)=my_pe
       else
        particle_active(n) = .false.
        particle_pe(n)=-1
       endif
      enddo
      call global_max_int(particle_pe,nptraj)
      end subroutine particle_inside_pe_domain


      subroutine particle_distribute
c=======================================================================
c     distribute particles to all pes
c=======================================================================
      use cpflame_module
      use particle_module
      implicit none 
      integer :: n
      call barrier()
      do n=1,nptraj
       if (particle_pe(n)>=0) then
        call bcast_real(    pxyz(3,n),3,particle_pe(n))
        call bcast_integer(pijk(3,n),3,particle_pe(n))
      elseif (my_pe==0) then
        print*,' WARNING : particle #',n,' is out of domain.',
     &          ' pxyz= ',pxyz(:,n),' pijk= ',pijk(:,n),
     &          ' particle PE=',particle_pe(n)
      endif
      enddo
      end subroutine particle_distribute


      subroutine seed_particles_randomly
c=======================================================================
c      seed the particles
c=======================================================================
      use cpflame_module
      use particle_module
      implicit none 
      real :: fxa, xs,xe,ys,ye,zs,ze
      integer :: n
#ifdef two_dim
      xs=xt(2); xe=xt(imt-1)
      ys=yt(2);
      zs=zt(2); ze=zt(km-1)
      do n=1,nptraj
        pxyz(2,n) = ys
        call random_number(fxa)
        pxyz(1,n) = xs+fxa*(xe-xs)
        call random_number(fxa)
        pxyz(3,n) = zs+fxa*(ze-zs)
      enddo
#else
      xs=xt(2); xe=xt(imt-1)
      ys=xt(2); ye=yt(jmt-1)
      zs=zt(2); ze=zt(km-1)
      do n=1,nptraj
        call random_number(fxa)
        pxyz(1,n) = xs+fxa*(xe-xs)
        call random_number(fxa)
        pxyz(2,n) = ys+fxa*(ye-ys)
        call random_number(fxa)
        pxyz(3,n) = zs+fxa*(ze-zs)
      enddo
#endif
      end subroutine seed_particles_randomly


      subroutine seed_particles_regular
c=======================================================================
c      seed the particles
c=======================================================================
      use cpflame_module
      use particle_module
      implicit none 
      real :: cubr,distx,disty,distz
      integer :: i,j,k,n
      real :: xs,xe,ys,ye,zs,ze

#ifdef two_dim
      xs=xt(2); xe=xt(imt-1)
      ys=yt(2);
      zs=zt(2); ze=zt(km-1)
      cubr       = float(nptraj)**0.5
      distx      = (xe - xs)/cubr
      distz      = (ze - zs)/cubr
      pxyz(1,1) = .5*distx + xs
      pxyz(2,1) =  ys
      pxyz(3,1) = .5*distz + zs
      do n=2,nptraj
        pxyz(1,n) = pxyz(1,n-1)
        pxyz(2,n) = pxyz(2,n-1)
        pxyz(3,n) = pxyz(3,n-1)
        pxyz(1,n)  = pxyz(1,n) + distx
        if (pxyz(1,n) > xe) then
         pxyz(1,n) = xs + (pxyz(1,n)-xe)
         pxyz(2,n) = ys
         pxyz(3,n) = pxyz(3,n) + distz
       endif
      enddo
#else

      xs=xt(2); xe=xt(imt-1)
      ys=yt(2); ye=yt(jmt-1)
      zs=zt(2); ze=zt(km-1)
      cubr       = float(nptraj)**0.333333
      distx      = (xe - xs)/cubr
      disty      = (ye - ys)/cubr
      distz      = (ze - zs)/cubr
      pxyz(1,1) = .5*distx + xs
      pxyz(2,1) = .5*disty + ys
      pxyz(3,1) = .5*distz + zs
      do n=2,nptraj
        pxyz(1,n) = pxyz(1,n-1)
        pxyz(2,n) = pxyz(2,n-1)
        pxyz(3,n) = pxyz(3,n-1)
        pxyz(1,n)  = pxyz(1,n) + distx
        if (pxyz(1,n) > xe) then
         pxyz(1,n) = xs + (pxyz(1,n)-xe)
         pxyz(2,n) = pxyz(2,n) + disty
         if (pxyz(2,n) > ye) then
          pxyz(2,n) = ys + (pxyz(2,n)-ye)
          pxyz(3,n) = pxyz(3,n) + distz
         endif
        endif
      enddo
#endif
      end subroutine seed_particles_regular




#else
      subroutine particle_dummy
      end
#endif
