#include "options.inc"

c=======================================================================
c     Templates for setup
c
c    internal waves  
c=======================================================================


      subroutine set_parameter
      use cpflame_module
      implicit none
      real :: fac=1.0
c ----------------------------------
c       set here main parameter
c ----------------------------------
      imt=64*fac;km=60*fac;jmt=3
      dx=0.25/fac
      dz=0.25/fac
      dt_in=0.025/fac
      eps2D_sor = 1e-7
      eps3D_sor = 1e-7
      runlen=150./86400.0 !/ fac**2
      snap_int = 0.5/86400.0
      enable_cyclic_x            =.true.
      enable_cyclic_y            =.true.
      enable_quicker_advection   =.true.
      enable_hydrostatic         =.false.
      enable_biharmonic_friction =.true.
      enable_vert_biha_friction  =.true.
      enable_biharmonic_diffusion=.true.
      enable_vert_biha_diffusion =.true.
      Ahbi = (5e-3/fac**2)*dx**2
      Khbi = Ahbi
      Avbi = (5e-3/fac**2)*dx**2
      Kvbi = Avbi
      A_h = 00e-3/fac**2
      A_v = A_h
      K_h = 0
      K_v = 0

      lat_ref=0.
      beta   =0.
      N_0 = 2*pi/10.
      M_0 = 0.0
      end subroutine set_parameter


      subroutine set_back_stratification(bb,u0)
c ----------------------------------
c    overide setting from module back_stratification
c ----------------------------------
      use cpflame_module
      implicit none
      real :: bb(imt,jmt,km),u0(imt,jmt,km),fxa
      integer :: i,j,k
      do k=1,km
       do j=1,jmt
        do i=1,imt
         bb(i,j,k)=-1e-2*tanh( (zt(k)-zt(km/2) ) /zt(1)*100 )
         u0(i,j,k)= tanh( (zt(k)-zt(km/2))/zt(1)*100)
         bb(i,j,k)=bb(i,j,k)*maskT(i,j,k)
         u0(i,j,k)=u0(i,j,k)*maskU(i,j,k)
        enddo
       enddo
      enddo
      end subroutine set_back_stratification


      subroutine initial_conditions
      use cpflame_module
      implicit none
      integer :: i,j,k
      real :: fxa,bb,dh
      do k=1,km
       do j=1,jmt
        do i=1,imt
         bb=-1e-2*tanh( (zt(k)-zt(km/2) ) /zt(1)*100 )
         dh=1e-4*zt(1)*sin(xt(i)/(20*dx)*pi)
         fxa=-1e-2*tanh((dh+zt(k)-zt(km/2) ) /zt(1)*100 )
         fxa=fxa-bb
         b(i,j,k,:)=fxa*maskT(i,j,k)
        enddo
       enddo
      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
c ----------------------------------
c    add here boundary conditions
c ----------------------------------
      end subroutine boundary_conditions


      subroutine topography
      use cpflame_module
      implicit none
c ----------------------------------
c      Add here topography  
c ----------------------------------
      end subroutine topography
