#include "options.inc"

c=======================================================================
c     Templates for setup
c
c      two-dimensional channel with restoring zones at side walls
c=======================================================================

      subroutine set_parameter
      use cpflame_module
      implicit none
c ----------------------------------
c       set here main parameter
c ----------------------------------
      imt=3;km=20;jmt=30
      dx=1./3.*degtom*cos(30./180.*pi)
      dz=100.0
      dt_in=2400.0
      K_h=000.0
      K_v=0.0e-4
      A_h=0000.0
      A_v=2.e-4
      Ahbi = 5e11
      cdbot=5e-5
      cdint=1e-7
      eps2D_sor = 1e-6
      runlen=4*30.0
      snap_int = 3.0

      enable_noslip              = .true.
      enable_cyclic_x            = .true.
      enable_quicker_advection   = .true.
      enable_bottom_stress       = .false.
      enable_interior_stress     = .true.
      enable_biharmonic_friction = .true.
      gamma=0.01
      lat_ref=30.0
      beta    = 2*omega*cos(lat_ref/180.*pi)/radius
      N_0 = 0.004
      M_0 = sqrt(1e-5*0.2/rho_0*g)
      end subroutine set_parameter


      subroutine initial_conditions
      use cpflame_module
      implicit none
      integer :: i,j,k
      real :: B0,x,y
c ----------------------------------
c      add here initial conditions
c ----------------------------------
      do k=1,km
       do j=1,jmt
        do i=1,imt
         x=float(i)/float(imt)
         y = (yt(j)-yt(2))/ (yt(jmt-1)-yt(2))
         B0=M_0**2*yt(j)+1e-6*sin(x*2.5*pi)*cos(y*4*pi)
         b(i,j,k,:)  = (B0-N_0**2*zt(k))*maskT(i,j,k)
        enddo
       enddo
      enddo
      end subroutine initial_conditions



      subroutine restoring_zones
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je, spg_width = 3
      real :: T_rest=1./(5.*86400),B0

c ----------------------------------
c      add here restoring zones
c ----------------------------------
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      if (my_blk_j==1) then
       do k=1,km
        do j=2,spg_width+1
         B0= yt(j)*M_0**2-N_0**2*zt(k)
         b(:,j,k,taup1)=b(:,j,k,taup1)+maskT(:,j,k)*c2dt* 
     &             t_rest/(j-1.)*(B0-b(:,j,k,taum1))
        enddo
       enddo
      endif
      if (my_blk_j==n_pes_j) then
       do k=1,km
        do j=jmt-1,jmt-spg_width,-1
         B0= yt(j)*M_0**2-N_0**2*zt(k)
         b(:,j,k,taup1)=b(:,j,k,taup1)+maskT(:,j,k)*c2dt* 
     &             t_rest/(-1.*(j-jmt))*(B0-b(:,j,k,taum1))
        enddo
       enddo
      endif
      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
