#include "options.inc"
    

c=======================================================================
c       integrate passive tracers
c       linked in the code in driver.F and in zonal_averages.F
c       and in bolus_velocity.F and in time_averages.F
c=======================================================================

#ifdef enable_diag_tracer

c=======================================================================
c prevent spurious under shooting by advection of passive tracer
c works only for positive tracer like concentrations
c#define delimit_tracer_fluxes
c=======================================================================

      module tracer_module
      implicit none
      real,allocatable :: tr(:,:,:,:,:)
      integer :: nt=0
      end module tracer_module

c=======================================================================
c      variables setups for tracer in template files
c=======================================================================
c#include "../extra_setup/age_tracer.F"
c#include "../extra_setup/age_theiss1.F"
c#include "../extra_setup/bgc_tracer.F"
c#include "../extra_setup/simple_tracer.F"
c#include "../extra_setup/salinity_THC.F"
#include "../extra_setup/salinity_double_diff1.F"


      subroutine init_tracer
c=======================================================================
c     allocate work space and read restart
c=======================================================================
      use cpflame_module
      use tracer_module
      implicit none
      logical :: file_exists
      if (my_pe==0) then
       print*,' Initializing tracer module '
#ifdef delimit_tracer_fluxes
       print*,' Delimiting passive tracer advective fluxes to prevent'
       print*,' spurios undershooting. Note that this only makes '
       print*,' sense for positive definite tracers like concentrations'
#endif
      endif
      call tracer_set_number
      allocate( tr(imt,jmt,km,0:2,nt) ); tr=0.0
      call tracer_initial_conditions
      if (my_pe==0) call init_tracer_diag
      inquire(FILE='tracer_restart.dta', EXIST=file_exists )
      if (file_exists) call tracer_read_restart
      if (my_pe==0) print*,' done initializing tracer module '
      end subroutine init_tracer


      subroutine tracer_read_restart
c=======================================================================
c     read tracer from restart file
c=======================================================================
      use cpflame_module
      use tracer_module
      implicit none
      integer :: j,io,n
      if (my_pe==0) print*,' reading tracer from restart file '
      call get_free_iounit(io)
      open(io,file='tracer_restart.dta',form='unformatted',status='old')
      do n=1,nt
       do j=1,jmt
        read(io) tr(:,j,:,tau,n),  tr(:,j,:,taum1,n)
       enddo
      enddo
      close(io)
      end subroutine tracer_read_restart


      subroutine tracer_write_restart
c=======================================================================
c     write restart to file
c=======================================================================
      use cpflame_module
      use tracer_module
      implicit none
      integer :: j,io,n,js,je,m
      if (my_pe==0) print*,' writing tracer restart file'
      do m=1,n_pes-1
       if (my_pe==m) then
        call send_integer(js_pe,1,0,0); call send_integer(je_pe,1,0,0)
        do j=js_pe,je_pe
          call send_real(tr(:,j,:,tau  ,:),imt*km*nt,0,0)
          call send_real(tr(:,j,:,taum1,:),imt*km*nt,0,0)
        enddo
       elseif (my_pe==0) then
        call recv_integer(js,1,m,0); call recv_integer(je,1,m,0)
        do j=js,je
          call recv_real(tr(:,j,:,tau  ,:)  ,imt*km*nt,m,0)
          call recv_real(tr(:,j,:,taum1,:)  ,imt*km*nt,m,0)
        enddo
       endif
       call barrier
      enddo
      if (my_pe==0) then
       call get_free_iounit(io)
       open(io,file='tracer_restart.dta',form='unformatted',
     &         status='unknown')
       do n=1,nt
        do j=1,jmt
         write(io) tr(:,j,:,tau,n), tr(:,j,:,taum1,n)
        enddo
       enddo
       close(io)
      endif
      end subroutine tracer_write_restart



      subroutine integrate_tracer
c=======================================================================
c     solve prognostic equation for tracers
c=======================================================================
      use cpflame_module
      use tracer_module
      implicit none
      integer :: i,j,k,js,je,n
      real :: adv_fe(imt,jmt,km), adv_ft(imt,jmt,km)
      real :: adv_fn(imt,jmt,km), diff_fn(imt,jmt,km)
      real :: diff_fe(imt,jmt,km),diff_ft(imt,jmt,km),fxa

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do n=1,nt
c---------------------------------------------------------------------------------
c        advection
c---------------------------------------------------------------------------------
       adv_fe(:,js_pe:je_pe,:)=0
       adv_fn(:,js_pe:je_pe,:)=0
       adv_ft(:,js_pe:je_pe,:)=0
       call adv_flux(adv_fe,adv_fn,adv_ft,tr(:,:,:,:,n))
       call border_exchg3D(adv_fn,1)
       call setcyclic3D(adv_fe)
       call setcyclic3D(adv_fn)
#ifdef delimit_tracer_fluxes
       call delimit_adv_flux(n,adv_fe,adv_fn,adv_ft)
#endif
c---------------------------------------------------------------------------------
c      horizontal diffusion
c---------------------------------------------------------------------------------
       diff_fe(:,js_pe:je_pe,:)=0
       diff_fn(:,js_pe:je_pe,:)=0
       diff_ft(:,js_pe:je_pe,:)=0
       do k=2,km-1
        do j=js,je
         do i=1,imt-1
          fxa=K_h*maskU(i,j,k)
          diff_fe(i,j,k)=fxa*(tr(i+1,j,k,taum1,n)-tr(i,j,k,taum1,n))/dx
         enddo
        enddo
       enddo
       call setcyclic3D(diff_fe)
       do k=2,km-1
        do j=js,je 
         do i=2,imt-1
          fxa=K_h*maskV(i,j,k)
          diff_fn(i,j,k)=fxa*(tr(i,j+1,k,taum1,n)-tr(i,j,k,taum1,n))/dx
         enddo
        enddo
       enddo
       call border_exchg3D(diff_fn,1)
       call setcyclic3D(diff_fn)
c---------------------------------------------------------------------------------
c      vertical diffusion
c---------------------------------------------------------------------------------
       do k=1,km-1
        do j=js,je
         do i=2,imt-1
          fxa=K_v*maskW(i,j,k)
          diff_ft(i,j,k)=fxa*(tr(i,j,k+1,taum1,n)-tr(i,j,k,taum1,n))/dz
         enddo
        enddo
       enddo
c---------------------------------------------------------------------------------
c       time tendency of b
c---------------------------------------------------------------------------------
       do k=2,km-1
        do j=js,je
         do i=2,imt-1
          tr(i,j,k,taup1,n)=tr(i,j,k,taum1,n)+maskT(i,j,k)*c2dt*( 
     &  -(adv_fe(i,j,k)-adv_fe(i-1,j,k))/dx
     &  -(adv_fn(i,j,k)-adv_fn(i,j-1,k))/dx
     &  -(adv_ft(i,j,k)-adv_ft(i,j,k-1))/dz
     &  +(diff_ft(i,j,k)-diff_ft(i,j,k-1))/dz
     &  +(diff_fe(i,j,k)-diff_fe(i-1,j,k))/dx
     &  +(diff_fn(i,j,k)-diff_fn(i,j-1,k))/dx  )
         enddo
        enddo
       enddo
c---------------------------------------------------------------------------------
c       Biharmonic diffusion (horizontal and vertical)
c---------------------------------------------------------------------------------
       if (enable_biharmonic_diffusion) 
     &         call biha_mix(diff_fe,diff_fn,diff_ft,tr(:,:,:,:,n))
c---------------------------------------------------------------------------------
c       remove unstable stratification      
c---------------------------------------------------------------------------------
       if (enable_hydrostatic) then
         call implicit_tracer_mix(tr(:,:,:,:,n),K_b)
       endif
      enddo

c---------------------------------------------------------------------------------
c      Nudging terms, etc
c---------------------------------------------------------------------------------
      call tracer_sources

      do n=1,nt
c---------------------------------------------------------------------------------
c        boundary exchange
c---------------------------------------------------------------------------------
       call border_exchg3D(tr(:,:,:,taup1,n),2)
       call setcyclic3D(tr(:,:,:,taup1,n) )
c---------------------------------------------------------------------------------
c       apply roberts time filter on time levels 
c---------------------------------------------------------------------------------
       do j=js,je
       tr(:,j,:,tau,n) = tr(:,j,:,tau,n) + gamma*
     &     (0.5*(tr(:,j,:,taup1,n)+tr(:,j,:,taum1,n))-tr(:,j,:,tau,n)) 
       enddo
       call border_exchg3D(tr(:,:,:,tau,n),2)
       call setcyclic3D(tr(:,:,:,tau,n) )
      enddo
c---------------------------------------------------------------------------------
c       Diagnostics
c---------------------------------------------------------------------------------
      if (snapshot_time_step.or.initial_time==current_time) then
        call diag_tracer
      endif
      if (last_time_step) call tracer_write_restart
      end subroutine integrate_tracer





      subroutine delimit_adv_flux(n,adv_fe,adv_fn,adv_ft)
c-----------------------------------------------------------------------
c   delimit advective fluxes for passive tracers
c   using the multidimensional positive definite
c   advection scheme (MPDCD) by Lafore, et al 1998
c-----------------------------------------------------------------------
      use cpflame_module
      use tracer_module
      implicit none
      integer, intent(in) :: n
      real :: adv_fe(imt,jmt,km)
      real :: adv_fn(imt,jmt,km)
      real :: adv_ft(imt,jmt,km)
      real :: betaf(imt,jmt,km)
      integer :: js,je,i,j,k
      real :: fxa,denom
      real, parameter :: small = 1.e-15

      js=max(2,js_pe); je=min(je_pe,jmt-1)
c    compute delimiter for outgoing flux
      betaf(:,js_pe:je_pe,:)=2.
      do j=js,je
        do k=2,km-1
          do i=2,imt-1
           denom =       max(0.,adv_fe(i  ,j,k)/dx)
           denom = denom-min(0.,adv_fe(i-1,j,k)/dx) 
           denom = denom+max(0.,adv_fn(i,j,k  )/dx)
           denom = denom-min(0.,adv_fn(i,j-1,k)/dx)
           denom = denom+max(0.,adv_ft(i,j,k)/dz)
           denom = denom-min(0.,adv_ft(i,j,k-1)/dz)
           denom=denom*c2dt
           denom = sign(1.,denom)*max(small,abs(denom)) ! should be positive anyway
           betaf(i,j,k)=tr(i,j,k,taum1,n)/denom
          enddo
        enddo
      enddo
      call border_exchg3D(betaf,1); call setcyclic3D(betaf)
c   delimit zonal flux
      do j=js,je
        do k=2,km-1
          do i=1,imt-1
           adv_fe(i,j,k)=
     &         min(1.,betaf(i ,j ,k))*max(0.,adv_fe(i,j,k))
     &       + min(1.,betaf(i+1,j,k))*min(0.,adv_fe(i,j,k))
          enddo
        enddo
      enddo
c   delimit meridional flux
      do j=js-1,je
        do k=2,km-1
          do i=2,imt-1
           adv_fn(i,j,k)=
     &         min(1.,betaf(i,j,k  ))*max(0.,adv_fn(i,j,k))
     &       + min(1.,betaf(i,j+1,k))*min(0.,adv_fn(i,j,k))
          enddo
        enddo
      enddo
c   delimit vertical flux
      do j=js,je
        do k=1,km-1
          do i=2,imt-1
           adv_ft(i,j,k)=
     &         min(1.,betaf(i,j,k))*max(0.,adv_ft(i,j,k))
     &       + min(1.,betaf(i,j,k+1))*min(0.,adv_ft(i,j,k))
          enddo
        enddo
      enddo
      end subroutine delimit_adv_flux



      subroutine delimit_adv_flux_taup1(n,adv_fe,adv_fn,adv_ft)
c-----------------------------------------------------------------------
c   same as above but cares for tracer at taup1
c-----------------------------------------------------------------------
      use cpflame_module
      use tracer_module
      implicit none
      integer, intent(in) :: n
      real :: adv_fe(imt,jmt,km)
      real :: adv_fn(imt,jmt,km)
      real :: adv_ft(imt,jmt,km)
      real :: betaf(imt,jmt,km)
      integer :: js,je,i,j,k
      real :: fxa,denom
      real, parameter :: small = 1.e-15

      js=max(2,js_pe); je=min(je_pe,jmt-1)
c    compute delimiter for outgoing flux
      betaf(:,js_pe:je_pe,:)=2.
      do j=js,je
        do k=2,km-1
          do i=2,imt-1
           denom =       max(0.,adv_fe(i  ,j,k)/dx)
           denom = denom-min(0.,adv_fe(i-1,j,k)/dx) 
           denom = denom+max(0.,adv_fn(i,j,k  )/dx)
           denom = denom-min(0.,adv_fn(i,j-1,k)/dx)
           denom = denom+max(0.,adv_ft(i,j,k)/dz)
           denom = denom-min(0.,adv_ft(i,j,k-1)/dz)
           denom=denom*c2dt
           denom = sign(1.,denom)*max(small,abs(denom)) ! should be positive anyway
           betaf(i,j,k)=tr(i,j,k,taup1,n)/denom
          enddo
        enddo
      enddo
      call border_exchg3D(betaf,1); call setcyclic3D(betaf)
c   delimit zonal flux
      do j=js,je
        do k=2,km-1
          do i=1,imt-1
           adv_fe(i,j,k)=
     &         min(1.,betaf(i ,j ,k))*max(0.,adv_fe(i,j,k))
     &       + min(1.,betaf(i+1,j,k))*min(0.,adv_fe(i,j,k))
          enddo
        enddo
      enddo
c   delimit meridional flux
      do j=js-1,je
        do k=2,km-1
          do i=2,imt-1
           adv_fn(i,j,k)=
     &         min(1.,betaf(i,j,k  ))*max(0.,adv_fn(i,j,k))
     &       + min(1.,betaf(i,j+1,k))*min(0.,adv_fn(i,j,k))
          enddo
        enddo
      enddo
c   delimit vertical flux
      do j=js,je
        do k=1,km-1
          do i=2,imt-1
           adv_ft(i,j,k)=
     &         min(1.,betaf(i,j,k))*max(0.,adv_ft(i,j,k))
     &       + min(1.,betaf(i,j,k+1))*min(0.,adv_ft(i,j,k))
          enddo
        enddo
      enddo
      end subroutine delimit_adv_flux_taup1




      subroutine init_tracer_diag
c-----------------------------------------------------------------------
c     initialize NetCDF snapshot file
c-----------------------------------------------------------------------
      use cpflame_module
      use tracer_module
      implicit none
#include "netcdf.inc"
      integer :: ncid,iret,i,j,k,n
      integer :: lon_tdim,lon_udim,z_tdim,z_udim,itimedim
      integer :: lat_tdim,lat_udim, tr_id(nt)
      integer :: dims(4), corner(4), edges(4)
      character :: name*24, unit*16

      call def_grid_cdf('tracer.cdf')
      iret=nf_open('tracer.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      call ncredf(ncid, iret)
      iret=nf_inq_dimid(ncid,'xt',lon_tdim)
      iret=nf_inq_dimid(ncid,'yt',lat_tdim)
      iret=nf_inq_dimid(ncid,'zt',z_tdim)
      iret=nf_inq_dimid(ncid,'Time',itimedim)
      dims = (/Lon_tdim,lat_tdim, z_tdim, iTimedim/)
      do n=1,nt
       write(name, '("tr",i2)') n
       call replace_space_zero(name)
       tr_id(n) = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
       name = 'passive tracer'; unit = ' '
       call dvcdf(ncid,tr_id(n),name,24,unit,16,spval)
      enddo
      call ncclos (ncid, iret)
      end subroutine init_tracer_diag


      subroutine diag_tracer
c-----------------------------------------------------------------------
c     write to NetCDF snapshot file
c-----------------------------------------------------------------------
      use cpflame_module
      use tracer_module
      implicit none
#include "netcdf.inc"
      integer :: ncid,iret,n,npe, corner(4), edges(4)
      real :: a(imt,js_pe:je_pe,km),fxa
      integer :: itdimid,ilen,rid,itimeid,fnid
      integer :: tr_id(nt)
      integer :: i,j,is,ie,js,je
      character :: name*24, unit*16
      type(time_type) :: time

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do npe=0,n_pes
       if (my_pe==npe) then
        iret=nf_open('tracer.cdf',NF_WRITE,ncid)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        do n=1,nt
         write(name, '("tr",i2)') n
         call replace_space_zero(name)
         iret=nf_inq_varid(ncid,name,tr_id(n))
        enddo
        iret=nf_inq_dimid(ncid,'Time',itdimid)
        iret=nf_inq_dimlen(ncid, itdimid,ilen)
        iret=nf_inq_varid(ncid,'Time',itimeid)
        if (my_pe==0) then
         ilen=ilen+1
         time = current_time-initial_time
         fxa = time%days + time%seconds/86400.
         iret= nf_put_vara_double(ncid,itimeid,ilen,1,fxa)
        endif
        Corner = (/1,js_pe,1,ilen/); 
        edges  = (/imt,je_pe-js_pe+1,km,1/)
        do n=1,nt
         a=tr(:,js_pe:je_pe,:,tau,n)
         where( maskT(:,js_pe:je_pe,:) == 0.) a = spval
         iret= nf_put_vara_double(ncid,tr_id(n),corner,edges,a)
        enddo
        call ncclos (ncid, iret)
       endif
       call barrier
      enddo
      end subroutine diag_tracer

#else
      subroutine tracer_dummy
      end
#endif
