#include "options.inc"

c=======================================================================
c     Templates for setup
c=======================================================================


      subroutine set_parameter
      use cpflame_module
      implicit none
c ----------------------------------
c       set here main parameter
c ----------------------------------
c      integer, parameter :: imt=32,km=10,jmt=32

      imt=32;km=10;jmt=32
      dx=1.0*degtom
      dz=100.0
      dt_in=3600.0
      K_h=1000.0
      K_v=1.e-4
      A_h=40000.0
      A_v=10.e-4
      eps2D_sor = 1e-6
      runlen=10*365.0
      snap_int = 5.0
      enable_quicker_advection=.false.
      lat_ref = 45.0      
      beta    = 2*omega*cos(lat_ref/180.*pi)/radius
      N_0     = 0.005
      M_0     = 0.0
      end subroutine set_parameter


      subroutine initial_conditions
      use cpflame_module
      implicit none
      integer :: i,j,k
c ----------------------------------
c      add here initial conditions
c ----------------------------------
      do k=1,km
       do j=1,jmt
        do i=1,imt
         b(i,j,k,:) = (0.2-N_0**2*zt(k))*maskT(i,j,k)
        enddo
       enddo
      enddo
      end subroutine initial_conditions



      subroutine boundary_conditions
      use cpflame_module
      implicit none
      integer :: k,i,j
      real :: x,y
c ----------------------------------
c    add here boundary conditions
c ----------------------------------
! units are m^2/s^3  
! since  b_z = N^2   b[m/s^2]   b_t = F_z   [m/s^3]=[ m^2/s^3]_z

      do j=1,jmt
       do i=1,imt
        surf_flux(i,j)  = dz/dt_in*0.0*maskT(i,j,km-1)
        y = (yu(j)-yu(1))/ (yu(jmt)-yu(1))
        surf_tau(i,j,1) = -0.5e-4*cos(2*pi*y)*maskU(i,j,km-1)
        surf_tau(i,j,2) = 0.0*maskV(i,j,km-1)
       enddo
      enddo

      end subroutine boundary_conditions



      subroutine restoring_zones
      use cpflame_module
      implicit none
      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 topography
      use cpflame_module
      implicit none
      integer :: i,k
      real :: z,y,hcm,h0
c ----------------------------------
c      Add here topography  
c ---------------------------------- 
      h0 = -zt(1)
      do k=2,km-1
       do i=1,imt
        z=H0*(1-0.3*exp(-(xt(i)-xt(imt/2))**2/800e3**2 ) )
        if (zt(k).lt.-z) maskT(i,:,k)=0.0
       enddo
      enddo
      end subroutine topography
