

c=======================================================================
c
c     Setup for salinity as tracer
c     linked in the code in ../extra_modules/tracer.F
c=======================================================================


      subroutine tracer_set_number
c ----------------------------------
c     set number of tracers
c ----------------------------------
      use tracer_module
      implicit none
      nt=1
      end subroutine tracer_set_number

      subroutine tracer_initial_conditions
c ----------------------------------
c     initial conditions for tracer
c ----------------------------------
      use cpflame_module
      use tracer_module
      implicit none
      integer :: i,j,k
      real :: fac,B0,x,y,z
      end subroutine tracer_initial_conditions


      subroutine tracer_sources
      use cpflame_module
      use tracer_module
      implicit none
      integer :: i,j,k,js,je
      real :: salt_flux(jmt)
c ----------------------------------
c      add here restoring zones, surface boundary conditions etc.
c ----------------------------------
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      call get_salt_flux(salt_flux)
      k=km-1
      do j=js,je
       tr(:,j,k,taup1,1)=tr(:,j,k,taup1,1)+maskT(:,j,k)*c2dt*
     &      salt_flux(j)/dz
      enddo
      end subroutine tracer_sources




      subroutine salinity_boundary_conditions
c ----------------------------------
c    add here boundary conditions
c ----------------------------------
      use cpflame_module
      use tracer_module
      implicit none
      integer :: i,j,k,js,je
      real :: fxa,fxb,Ly,t_rest,sst(imt,jmt),sss(imt,jmt),sig(imt,jmt)
      real :: salt_flux(jmt),t_star(jmt)
c ----------------------------------
c      add here restoring zones, surface boundary conditions etc.
c ----------------------------------
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      Ly = yt(jmt)-yt(1)
      call get_salt_flux(salt_flux)
      t_rest=30*86400
      do j=js,je
       t_star(j)=-15*exp(-(yt(j)-yt(jmt-1))**2/(Ly/8.)**2 )
      enddo
      k=km-1
      do j=js,je ! sig = 0.8 S - 0.2T ,  T = (0.8S -sig)/0.2
       sig(:,j)=  b(:,j,k,taum1)*1024/g
       sss(:,j)= tr(:,j,k,taum1,1)
       sst(:,j) = (0.8*sss(:,j)-sig(:,j)) /0.2
       surf_flux(:,j)=dz/t_rest*(t_star(j)-sst(:,j)) 
       surf_flux(:,j)=(0.8*salt_flux-0.2*surf_flux(:,j))/1024*g
      enddo

      end subroutine salinity_boundary_conditions



      subroutine get_salt_flux(salt_flux)
      use cpflame_module
      use tracer_module
      implicit none
      real :: salt_flux(jmt)
      integer :: i,j,k,js,je
      real :: fxa,Ly
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      Ly = yt(jmt)-yt(1)
c     surface forcing
      salt_flux(js:je)=0.0
      fxa=0.5/86400.
      do j=js,je
       salt_flux(j)=-fxa*exp(-(yt(j)-yt(jmt-1))**2/(Ly/8.)**2 )
      enddo
      fxa=0.; 
      do j=js,je
       fxa=fxa+salt_flux(j)
      enddo; 
      call global_sum(fxa)
      fxa=fxa/(jmt-2.)
      salt_flux=salt_flux-fxa
      end subroutine get_salt_flux

