#include "options.inc"


      module passive_tracer_module
c-----------------------------------------------------------------------
c   Module for passive tracer stuff
c-----------------------------------------------------------------------
      implicit none
      logical :: enable_cfc11 = .false.
      integer :: ref_year_cfc11 = 1900   ! model year of CFC11 time series
                                            ! corresponding to 1931 in reality
      integer :: n_cfc11  = 3            ! number of tracer used for CFC11
      logical :: enable_cfc12 = .false.
      integer :: ref_year_cfc12 = 1900   ! reference year
      integer :: n_cfc12  = 4            ! number of tracer used for CFC12

      logical :: enable_npzd = .false. ! activate the ecosystem model (4-7 tracer)
      logical :: enable_doc  = .false. ! activate DOC in the ecosystem model
      logical :: enable_co2  = .false. ! activate CO2 in the ecosystem model
      logical :: enable_o2   = .false. ! activate O2 as well
      integer :: nt_npzd = 4             ! Number of tracers in NPZD model
      integer, allocatable :: n_npzd(:)  ! mapping of tracer number to biological var.
                                            ! (n_npzd lists all tracers in NPZD)
      integer :: n_nutrient,n_phyto      ! again the mapping
      integer :: n_zoo,n_detritus        ! but with explicit naming 
      integer :: n_co2,n_o2,n_doc
      integer :: ref_year_co2 = 1900     ! reference year for atmospheric CO2
      logical :: enable_pco2_series = .false. ! use fixed value for atmos CO2

      logical :: enable_diag_npzd = .false.  ! some diagnostics for npzd model

      end module passive_tracer_module


      subroutine passive_tracer_init
c-----------------------------------------------------------------------
c   setup all passive tracers here
c   called by setup
c-----------------------------------------------------------------------
      use spflame_module
      use passive_tracer_module
      implicit none

      namelist /passive_tracer_nml/ 
     &    enable_cfc11,ref_year_cfc11,
     &    enable_cfc12,ref_year_cfc12,
     &    enable_npzd,enable_co2,enable_o2,enable_doc,
     &    ref_year_co2,enable_pco2_series,enable_diag_npzd
      integer :: k = 3, kk,io

      if (nt <=2) then
       if (my_pe==0) print*,' no passive tracers to be initialized'
       return
      endif

      if (my_pe==0) then
       print*,' '
       print*,' initialization of passive tracers'
       print*,' '
       print*,' reading namelist file namelist.passive_tracer'
      endif
c
c   read in namelist parameter for passive tracers
c
      call getunit(io,'namelist.passive_tracer','fsr')
      read(io, nml = passive_tracer_nml) 
      close(io)

      if (my_pe==0) then
       print*,' read the following namelist parameter: '
       write(6, nml = passive_tracer_nml ) 
      endif

      if (enable_cfc11) then 
        n_cfc11 = k
        call cfc11_init()
        k=k+1
      endif
        
      if (enable_cfc12) then
        n_cfc12 = k
        call cfc12_init()
        k=k+1
      endif

      if (enable_npzd)  then
        nt_npzd=4
        if (enable_co2)  nt_npzd=nt_npzd+1
        if (enable_o2 )  nt_npzd=nt_npzd+1
        if (enable_doc ) nt_npzd=nt_npzd+1
        allocate(n_npzd(nt_npzd))
        do kk=1,nt_npzd
         n_npzd(kk) = kk+k-1    
        enddo
        n_nutrient = n_npzd(1)
        n_phyto    = n_npzd(2)
        n_zoo      = n_npzd(3)
        n_detritus = n_npzd(4)
        kk=n_npzd(4)+1
        if (enable_co2)  then
          n_co2=kk; kk=kk+1
        endif
        if (enable_o2)  then
          n_o2=kk; kk=kk+1
        endif
        if (enable_doc)  then
          n_doc=kk; kk=kk+1
        endif
        call npzd_init()
        k=k+nt_npzd
      endif
        
      if (my_pe==0) then
       k=k-1
       print*,' recognized altogether ',k,' tracers (incl. temp/sal)'
c
c     check number of tracers
c
       if (nt < k) then
        if (my_pe==0) print*,' ERROR: number of tracer is only ',nt
        call halt_stop(' in passive_tracer_init ')
       endif

       if ((nt > k).and.(my_pe==0)) then
         print*,' Note: cannot recognize the rest of the tracers'
       endif

      endif

      if (my_pe==0) print*,' passive tracer initialization done'

      end subroutine passive_tracer_init


      subroutine passive_tracer_add()
c-----------------------------------------------------------------------
c   add to tracer at tau+1 whats necessary for
c   each passive tracer
c-----------------------------------------------------------------------
      use spflame_module
      use passive_tracer_module
      implicit none
      if (nt <=2) return
      if (enable_npzd) call npzd_model()
      end subroutine passive_tracer_add



      subroutine passive_tracer_sflx
c-----------------------------------------------------------------------
c   apply surface fluxes to passive tracers
c   (called by setvbc, icemask is applied afterwards)
c-----------------------------------------------------------------------
      use spflame_module
      use passive_tracer_module
      implicit none
      if (nt <=2) return
      if (enable_npzd)  call npzd_sflx() 
      if (enable_cfc11) call cfc_sflx(n_cfc11 ,'CFC11')
      if (enable_cfc12) call cfc_sflx(n_cfc12 ,'CFC12')
      end subroutine passive_tracer_sflx


      subroutine wanninkhof(tc,u,kw,len)
c-----------------------------------------------------------------------
c   Piston velocity for surface fluxes a la Wanninkhof , 1992
c-----------------------------------------------------------------------
      implicit none
      integer :: len,i
      real :: tc(len),u(len)   ! temperature, wind speed (input)
      real :: kw(len)          ! piston velocity (output)
      do i=1,len
       kw(i) = 0.39*u(i)*u(i)
      enddo
      end subroutine wanninkhof


      subroutine wanninkhof99(tc,u,kw,len)
c-----------------------------------------------------------------------
c   Piston velocity for surface fluxes a la Wanninkhof+McGillis , 1999
c-----------------------------------------------------------------------
      implicit none
      integer :: len,i
      real :: tc(len),u(len)   ! temperature, wind speed (input)
      real :: kw(len)          ! piston velocity (output)
      real :: uu
      do i=1,len
       uu=u(i)*u(i)
       kw(i) = 1.09*u(i)-0.333*uu+0.078*u(i)*uu
      enddo
      end subroutine wanninkhof99



      subroutine flux_delimit_init
c-----------------------------------------------------------------------
c   flux delimiter module
c   linked in the code in tracer and setup
c                               c eden
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      if (my_pe==0) then
       print*,' Initialisation of flux delimiter'
       print*,''
       print*,' limiting advective fluxes for passive tracers'
       print*,' using the multidimensional positive definite'
       print*,' advection scheme (MPDCD) by Lafore, et al 1998 '
       print*,''
       if (nt<=2) then
        print*,' however, there are no passive tracers to limit '
       endif
       if (enable_diffusion_biharmonic) then
        print*,' limiting also biharmonic diffusion scheme '
        print*,' which is also not positiv definite '
       endif
       print*,' done '
      endif
      end subroutine flux_delimit_init


      subroutine delimit_adv_flux(n,adv_fe,adv_fn,adv_fb)
c-----------------------------------------------------------------------
c   delimit advective fluxes for passive tracers
c   using the multidimensional positive definite
c   advection scheme (MPDCD) by Lafore, et al 1998
c   SPFLAME version : c.eden
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      integer, intent(in) :: n
      real adv_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      real adv_fn(is_pe:ie_pe,km,js_pe-1:je_pe)
      real adv_fb(is_pe:ie_pe,0:km,js_pe:je_pe)
c     local array
      real betaf(is_pe-1:ie_pe+1,0:km+1,js_pe-1:je_pe+1)
      integer is,ie,js,je
      integer i,j,k
      real :: fxa,denom
      real, parameter :: small = 1.e-15

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c   compute delimiter for outgoing flux
c
      betaf=2.
      do j=js,je
        do k=1,km
          do i=is,ie
           denom =       max(0.,adv_fe(i  ,k,j)*cstdxt2r(i,j))
           denom = denom-min(0.,adv_fe(i-1,k,j)*cstdxt2r(i,j)) 

           denom = denom+max(0.,adv_fn(i,k,j  )*cstdyt2r(j))
           denom = denom-min(0.,adv_fn(i,k,j-1)*cstdyt2r(j))

           denom = denom-min(0.,adv_fb(i,k,j)*dzt2r(k))
           denom = denom+max(0.,adv_fb(i,k-1,j)*dzt2r(k))

           denom=denom*c2dt
           denom = sign(1.,denom)*max(small,abs(denom)) ! should be positive anyway

           betaf(i,k,j)=t(i,k,j,n,taum1)/denom
          enddo
        enddo
      enddo

      call border_exchg(betaf,km+2,1)
      call set_cyclic(betaf,km+2,1)
c
c   delimit zonal flux
c
      do j=js,je
        do k=1,km
          do i=is-1,ie
           adv_fe(i,k,j)=
     &         min(1.,betaf(i  ,k,j))*max(0.,adv_fe(i,k,j))
     &       + min(1.,betaf(i+1,k,j))*min(0.,adv_fe(i,k,j))
          enddo
        enddo
      enddo
c
c   delimit meridional flux
c
      do j=js-1,je
        do k=1,km
          do i=is,ie
           adv_fn(i,k,j)=
     &         min(1.,betaf(i,k,j  ))*max(0.,adv_fn(i,k,j))
     &       + min(1.,betaf(i,k,j+1))*min(0.,adv_fn(i,k,j))
          enddo
        enddo
      enddo
c
c   delimit vertical flux
c
      do j=js,je
        do k=0,km
          do i=is,ie
           adv_fb(i,k,j)=
     &         min(1.,betaf(i,k  ,j))*min(0.,adv_fb(i,k,j))
     &       + min(1.,betaf(i,k+1,j))*max(0.,adv_fb(i,k,j))
          enddo
        enddo
      enddo
      end subroutine delimit_adv_flux





      subroutine delimit_all_fluxes(n,adv_fe,adv_fn,adv_fb,
     &                                         diff_fe,diff_fn,diff_fb)
c-----------------------------------------------------------------------
c   delimit diffusive fluxes for passive tracers
c   using the multidimensional positive definite
c   advection scheme (MPDCD) by Lafore, et al 1998
c   SPFLAME version : c.eden
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      integer, intent(in) :: n
      real adv_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      real adv_fn(is_pe:ie_pe,km,js_pe-1:je_pe)
      real adv_fb(is_pe:ie_pe,0:km,js_pe:je_pe)
      real diff_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      real diff_fn(is_pe:ie_pe,km,js_pe-1:je_pe)
      real diff_fb(is_pe:ie_pe,0:km,js_pe:je_pe)

      integer is,ie,js,je
      integer i,j,k
      real :: fxa,denom, small = 1.e-15
c     local array
      real betaf(is_pe-1:ie_pe+1,0:km+1,js_pe-1:je_pe+1)

      real adv_tx
      ADV_Tx(i,k,j) = (adv_fe(i,k,j) - adv_fe(i-1,k,j))*cstdxt2r(i,j)
#ifdef partial_cell
     &                             /dht(i,k,j)
#endif
      real adv_ty
      ADV_Ty(i,k,j) = (adv_fn(i,k,j) - adv_fn(i,k,j-1))*cstdyt2r(j)
# ifdef partial_cell
     &                             /dht(i,k,j)
# endif
      real adv_tz
      ADV_Tz(i,k,j) = (adv_fb(i,k-1,j) - adv_fb(i,k,j))
#ifdef partial_cell
     &               *.5/dht(i,k,j)
#else
     &               *dzt2r(k)
#endif


      diff_fn=diff_fn*(-1)
      diff_fe=diff_fe*(-1)
      diff_fb=diff_fb*(-1)

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c   compute delimiter for outgoing flux
c
      betaf=2.
      do j=js,je
        do k=1,km
          do i=is,ie

           denom=      max(0.,
     &                diff_fe(i  ,k,j)*tmask(i+1,k,j)*cstdxtr(i,j))
           denom=denom-min(0.,
     &                diff_fe(i-1,k,j)*tmask(i-1,k,j)*cstdxtr(i,j)) 

           denom=denom+max(0.,
     &                   diff_fn(i,k,j  )*tmask(i,k,j+1)*cstdytr(j))
           denom=denom-min(0.,
     &                   diff_fn(i,k,j-1)*tmask(i,k,j-1)*cstdytr(j))

           denom=denom-min(0.,diff_fb(i,k  ,j)*dztr(k))
           denom=denom+max(0.,diff_fb(i,k-1,j)*dztr(k))

           denom=denom*c2dt
           denom = sign(1.,denom)*max(small,abs(denom)) ! should be positive anyway

           betaf(i,k,j)=(t(i,k,j,n,taum1)
     &  -c2dt*(ADV_Tx(i,k,j)+ADV_Ty(i,k,j)+ADV_Tz(i,k,j))*tmask(i,k,j) )
     &                    /denom
          enddo
        enddo
      enddo

      call border_exchg(betaf,km+2,1)
      call set_cyclic(betaf,km+2,1)
c
c   delimit zonal flux
c
      do j=js,je
        do k=1,km
          do i=is-1,ie
           diff_fe(i,k,j)=
     &         min(1.,betaf(i  ,k,j))*max(0.,diff_fe(i,k,j))
     &       + min(1.,betaf(i+1,k,j))*min(0.,diff_fe(i,k,j))
          enddo
        enddo
      enddo
c
c   delimit meridional flux
c
      do j=js-1,je
        do k=1,km
          do i=is,ie
           diff_fn(i,k,j)=
     &         min(1.,betaf(i,k,j  ))*max(0.,diff_fn(i,k,j))
     &       + min(1.,betaf(i,k,j+1))*min(0.,diff_fn(i,k,j))
          enddo
        enddo
      enddo
c
c   delimit vertical flux
c
      do j=js,je
        do k=0,km
          do i=is,ie
           diff_fb(i,k,j)=
     &         min(1.,betaf(i,k  ,j))*min(0.,diff_fb(i,k,j))
     &       + min(1.,betaf(i,k+1,j))*max(0.,diff_fb(i,k,j))
          enddo
        enddo
      enddo

      diff_fn=diff_fn*(-1)
      diff_fe=diff_fe*(-1)
      diff_fb=diff_fb*(-1)

      end subroutine delimit_all_fluxes



