#include "options.inc"
c
c-----------------------------------------------------------------------
c    inertial instability
c-----------------------------------------------------------------------
c
      module jets_module
c
c    u(t=0) = u_0 + u_y y + delta y**2 + amp noise
c    T(t=0) = 20.0+rho0/(g alpha)*N**2
c    T*     = 20.0+rho0/(g alpha)*(N**2- u_z f y)
c    N**2   = Bv0**2 int_(-h)^0  exp(-z/h_W)  dz
c
      real,parameter :: g   = 9.81e2 ! cm/s^2
      real,parameter :: Bv0 = 0.0005 ! 1/s
      real,parameter :: h_W = 300.0e12
      real,parameter :: delta = 0.0e-13 ! 1/s /cm
      real,parameter :: u_y   = 10.0e-6  ! 1/s
      real,parameter :: u_0   = 0.0     ! cm/s
      real,parameter :: u_z   = 0.0e-3  ! 1/s
      real,parameter :: amp   = 0.1     ! cm/s
      real,parameter :: tscl   = 1.0/(3.0*86400.0)   ! 1/s
      end module jets_module



      subroutine grid_template()
c-----------------------------------------------------------------------
c     set up a grid here
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      real :: dx=0.,dy=0.,dz=0.,x0=0.,y0=0.,sumt,sumu
      integer i,k
c     equidistant vertical grid
      logical :: flame_grid = .false.
      integer,parameter :: km_in = 200
      real ddz(km_in)
      data ddz/ 200*5.e2/
c
      x0=280.; y0=-3.0; dx = 1./6.; dy = dx
      if (my_pe==0) then
        print*,''
        print*,'setting up a simple grid '
        print*,' dx= ',dx,' dy= ',dy
        print*,''
      endif
      dxtdeg =  dx  ; dytdeg =  dy  ; dxudeg =  dx  ; dyudeg =  dy  
      xt(1)=x0; xu(1)=x0+dxtdeg(1)/2.
      do i=2,imt
        xt(i)=xt(i-1)+dxtdeg(i); xu(i)=xu(i-1)+dxudeg(i)
      enddo
      yt(1)=y0; yu(1)=y0+dytdeg(1)/2.
      do i=2,jmt
        yt(i)=yt(i-1)+dytdeg(i); yu(i)=yu(i-1)+dyudeg(i)
      enddo
      if (km /= km_in ) then
        print*,' km must be ',km_in
        stop
      endif
      dzt=ddz; dzw(1) = dzt(1); sumt = dzt(1)
      sumu = 0.5*dzt(1) + dzw(1)
      do i=2,km
       sumt = sumt + dzt(i)
       dzw(i) = 2*(sumt-sumu)
       sumu = sumu + dzw(i)
      enddo
      zt(1) = 0.5*dzt(1); zw(1) = dzt(1)
      do i=2,km
        zw(i) = zw(i-1) + dzt(i); zt(i) = zt(i-1) + dzw(i-1)
      enddo
      dzw(0)  = zt(1); dzw(km) = zw(km) - zt(km)
      end subroutine grid_template



      subroutine kmt_template
c-----------------------------------------------------------------------
c     set up a simple topography
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      kmt_big=km
      end subroutine kmt_template




      subroutine htp_template
c-----------------------------------------------------------------------
c     set up a simple topography for partial cells
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      integer :: i,j
      htp=0.
      do i=1,imt
        do j=1,jmt
           if (kmt_big(i,j)>0)  htp(i,j)=zw(kmt_big(i,j))
        enddo
      enddo
      end subroutine htp_template




      subroutine init_cond_template
c-----------------------------------------------------------------------
c     setup simple initial conditions
c     for temperature, salinity and velocity
c     Note that streamfunction is diagnosed from full velocities.
c-----------------------------------------------------------------------
      use spflame_module
      use freesurf_module
      use jets_module
      implicit none
      integer :: i,j,k,io,n
      real    :: baru(is_pe:ie_pe)
      real    :: x,y,z,Bvsqr(km)

      Bvsqr(1)= Bv0**2*dzw(0)
      do k=2,km
        Bvsqr(k)= Bvsqr(k-1)+Bv0**2*exp(-zw(k-1)/h_W) *dzw(k-1)
      enddo

      do j=js_pe,je_pe
        do k=1,km
         do i=is_pe,ie_pe
           y=yu(j)* radius/180.0*pi  ! in centimeters
           u(i,k,j,1,:) = u_0 + u_y*y+ delta*y*y/2.0 
     &       + amp*(sin(  zt(k)/120.0e2  *2*pi)
     &       +      sin(  zt(k)/60.0e2  *2*pi)
     &       +      sin(  zt(k)/30.0e2  *2*pi)
     &       +      sin(  zt(k)/15.0e2  *2*pi)
     &       +      sin(  zt(k)/7.5e2  *2*pi)
     &       +      sin(  zt(k)/220.0e2  *2*pi)
     &       +      sin(  zt(k)/320.0e2  *2*pi)
     &       +      sin(  zt(k)/420.0e2  *2*pi)
     &       +      sin(  zt(k)/520.0e2  *2*pi)
     &       +      sin(  zt(k)/620.0e2  *2*pi)  )
c           u(i,k,j,1,:) = min(u(i,k,j,1,:),20.0)
c           u(i,k,j,1,:) = max(u(i,k,j,1,:),-20.0)
           y=y/100./1000. ! in km
           u(i,k,j,1,:) = u(i,k,j,1,:)*exp(-(y/200.)**2)

           u(i,k,j,2,:) = 0.   

           y     =  yt(j)* radius/180.0*pi  ! in centimeters
           t(i,k,j,1,:) = 20.0+rho0/(g*alpha_lin)*Bvsqr(k)
           t(i,k,j,2,:) = 0.   ! salinity equals 35 PSU

           t(i,k,j,1,:) = t(i,k,j,1,:)*tmask(i,k,j)
           t(i,k,j,2,:) = t(i,k,j,2,:)*tmask(i,k,j)
           u(i,k,j,1,:) = u(i,k,j,1,:)*umask(i,k,j)
           u(i,k,j,2,:) = u(i,k,j,2,:)*umask(i,k,j)

         enddo 
        enddo
      enddo

      do n=1,nt
       call set_cyclic_4th_order(t(:,:,:,n,0),km,2)
       call set_cyclic_4th_order(t(:,:,:,n,1),km,2)
      enddo
      do n=1,2
       call set_cyclic_4th_order(u(:,:,:,n,0),km,2)
       call set_cyclic_4th_order(u(:,:,:,n,1),km,2)
      enddo

c-----------------------------------------------------------------------
c        barotropic mode
c-----------------------------------------------------------------------

         if (enable_freesurf) then
c
c         set the free surface
c

         else ! enable_freesurf
c
c        diagnose streamfunction from full velocities
c        (works only for du/dx=0)
c
         psi = 0.0
         do j=js_pe,je_pe
c         depth integrated velocity
          baru=0.
          do k=1,km
           do i=is_pe,ie_pe
            baru(i) = baru(i) + u(i,k,j,1,0)*umask(i,k,j)
#ifdef partial_cell
     &                                    *dhu(i,k,j)
#else
     &                                    *dzt(k)
#endif
           enddo
          enddo
c         local integrant
          do i=is_pe,ie_pe
            psi(i,j+1,1) = -dyu(j)*baru(i)*umask(i,1,j)
          enddo
         enddo
         call border_exchg(psi(:,:,1),1,1)
         call set_cyclic(  psi(:,:,1),1,1)
c        Meridionally integrate psiu.
         do j=js_pe,je_pe
          do i=is_pe,ie_pe
            psi(i,j,1) = psi(i,j,1) + psi(i,j-1,1)
            psi(i,j,2) = psi(i,j,1) 
          enddo
         enddo
        endif ! enable_freesurf
      end subroutine init_cond_template



      subroutine sbc_template
c-----------------------------------------------------------------------
c     setup simple surface boundary conditions
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      smf=0.
      stf_clim=0.
      stf_rest=0.
      if (enable_ktmix.or.enable_tkemix) ustar=0.     ! no wind stirring
      end subroutine sbc_template



      subroutine blue_template
c-----------------------------------------------------------------------
c     setup simple temp / salinity arrays for BLUE
c-----------------------------------------------------------------------
      use spflame_module
      use blue_module
      implicit none
      end subroutine blue_template


      subroutine sponge_template(n)
c-----------------------------------------------------------------------
c     setup simple sponge layers. 
c-----------------------------------------------------------------------
      use spflame_module
      use jets_module
      implicit none
      integer, intent(in) :: n ! number of the tracer
      integer :: i,j,k,is,ie,js,je
      real :: y,fx,tstar,Bvsqr(km)

c-----------------------------------------------------------------------
c     inertial instability
c-----------------------------------------------------------------------

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

      if (n==1) then

       Bvsqr(1)= Bv0**2*dzw(0)
       do k=2,km
        Bvsqr(k)= Bvsqr(k-1)+Bv0**2*exp(-zw(k-1)/h_W) *dzw(k-1)
       enddo

       do j=js,je
        do k=1,km
         do i=is,ie
           y     =  yt(j)* radius/180.0*pi  ! in centimeters
           fx    =  (cori(i,j-1,1)+cori(i,j,1))/2.0
           tstar = 20.0+rho0/(g*alpha_lin)*(Bvsqr(k)-fx*u_z*y )
           source(i,k,j)=source(i,k,j)+
     &        tscl*(tstar-t(i,k,j,n,taum1))*tmask(i,k,j)
         enddo 
        enddo
       enddo
      endif

      end subroutine sponge_template




