#include "options.inc"


c#define with_topography

c=======================================================================
c     Templates for setup
c
c      spin down of a front
c=======================================================================

c increase this parameter to enhance resolution
#define RESOLVE  1.0


      subroutine set_parameter
      use cpflame_module
      implicit none
      real :: f_0
c ----------------------------------
c       set here main parameter
c ----------------------------------
      imt   = 32*RESOLVE
      jmt   = 3
      km    = 20*RESOLVE

      dx    = 20e3/RESOLVE
      dz    = 50.0 /RESOLVE
      dt_in = 1200.0/RESOLVE

      K_h=000.0 /RESOLVE
      A_h=0000.0/RESOLVE
      Ahbi  = 5e11/RESOLVE**2
      Khbi  = Ahbi

      K_v=0.0e-4/RESOLVE
      A_v=10.e-4/RESOLVE
      cdbot = 0e-5*RESOLVE

      eps2D_sor = 1e-9
      runlen    = 1200./RESOLVE**3
      snap_int  = 5.

      enable_noslip               = .false.
      enable_cyclic_x             = .false.
      enable_cyclic_y             = .true.
      enable_bottom_stress        = .false.
      enable_biharmonic_friction  = .true.
      enable_biharmonic_diffusion = .false.
      enable_hydrostatic          = .true.
      enable_quicker_advection    = .false.  

      gamma=0.01
      lat_ref = 25.
      beta    = 0*2*omega*cos(lat_ref/180.*pi)/radius
      N_0     = 0.004
      M_0     = 0

      if (my_pe==0) then
       f_0 = 2*omega*sin(lat_ref/180.*pi)
       print*,' '
       print*,' Rossby radius '
       print*,'  N h/(f pi)     = ',N_0*(km-2)*dz/(f_0*pi)/1e3,' km'
       print*,' '
      endif

      end subroutine set_parameter



      subroutine initial_conditions
      use cpflame_module
      implicit none
      integer :: i,j,k
      real :: fac,B0,x,y,z
c ----------------------------------
c      add here initial conditions
c ----------------------------------
      fac = 2*0.2/rho_0*g
      do j=1,jmt
       do k=1,km
        do i=1,imt
c         y = (yt(j)-yt(2))/ (yt(jmt-1)-yt(2))
         y = (xt(i)-xt(2))/ (xt(imt-1)-xt(2))
         B0=fac*0.5*(tanh((y-0.5)/0.1)+1.0)
         b(i,j,k,:)  = (B0-N_0**2*zt(k))*maskT(i,j,k)
        enddo
        enddo
      enddo

c      do j=1,jmt
c       do k=km-3,km
c        do i=1,imt
c         y = (yt(j)-yt(2))/ (yt(jmt-1)-yt(2))
c         B0=fac*0.5*(tanh((y-0.5)/0.1)+1.0)
c         b(i,j,k,:)  = (B0-0.0010**2*zt(k))*maskT(i,j,k)
c        enddo
c        enddo
c      enddo

      end subroutine initial_conditions


      subroutine restoring_zones
      use cpflame_module
      implicit none
c ----------------------------------
c      add here restoring zones
c ----------------------------------
      end subroutine restoring_zones


      subroutine momentum_restoring_zones
      use cpflame_module
      implicit none
c ----------------------------------
c      add here restoring zones for momentum
c ----------------------------------
      end subroutine momentum_restoring_zones




      subroutine boundary_conditions
      use cpflame_module
      implicit none
      integer :: i,j
      real :: L_y,y,t_rest,bstar
c ----------------------------------
c    add here boundary conditions
c ----------------------------------
      end subroutine boundary_conditions



      subroutine topography
      use cpflame_module
      implicit none
      integer :: i,j
      real :: L_x,L_y
c ----------------------------------
c      Add here topography  
c ----------------------------------
      end subroutine topography
