#include "options.inc"


c#define with_topography

c=======================================================================
c     Templates for setup
c
c      Wide eddying domain with restoring zones at side walls
c=======================================================================

c increase this parameter to enhance resolution
#define RESOLVE  0.25


      subroutine set_parameter
      use cpflame_module
      implicit none
      real :: f_0,Ro,delta,s,fxa
c ----------------------------------
c       set here main parameter
c ----------------------------------
      imt   = 128*RESOLVE 
      km    = 20 !*RESOLVE;
      jmt   = 128*RESOLVE

      dx    = 1./3.*degtom*cos(30./180.*pi)/RESOLVE
      dz    = 100.0!   /RESOLVE
      dt_in = 1800.0/RESOLVE

      K_h=000.0 /RESOLVE
      A_h=0000.0/RESOLVE
      K_v=1.0e-4!  /RESOLVE
      A_v=10.e-4!  /RESOLVE

      Ahbi  = 8e11/RESOLVE**2
      Khbi  = Ahbi
c      cdbot = 1e-6*RESOLVE
      cdbot = 0e-6*RESOLVE
      cdint = 0.1e-7


      eps2D_sor = 1e-9
      runlen    = 365.0/RESOLVE**3
      snap_int  = 3.0

c      enable_noslip               = .true.
      enable_interior_stress      = .true.
      enable_biharmonic_friction  = .true.
c      enable_biharmonic_diffusion = .true.
      enable_quicker_advection    = .true.
      enable_hydrostatic          = .true.

      gamma=0.01
      lat_ref = 20.
      beta    = 2*omega*cos(lat_ref/180.*pi)/radius
      N_0     = 0.004
      M_0     = sqrt(0.0e-5*0.1/rho_0*g)

      if (my_pe==0) then
       f_0 = 2*omega*sin(lat_ref/180.*pi)
       Ro  = (km-2)*dz*M_0**2/((jmt-2)*dx)/f_0**2
       delta = (km-2)*dz/( (jmt-2)*dx )
       s = M_0**2/N_0**2
       print*,' '
       print*,' The setup :'
       print*,' '
       print*,' Prandtl number '
       print*,'   f^2/N^2       = ',f_0**2/N_0**2
       print*,' Aspect ratio '
       print*,'   delta=H/L     = ',delta
       print*,' Rossby radius '
       print*,'  N h/(f pi)     = ',N_0*(km-2)*dz/(f_0*pi)/1e3,' km'
       print*,' Isopycnal slope '
       print*,'   s=M^2/N^2     = ',s
       print*,' '
       print*,' Rossby number '
       print*,'   H M^2/(L f^2) = ',Ro
      endif

      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)+0.5e-4*sin(x*8.5*pi)*exp(-(y-0.5)**2/0.5**2)
         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*RESOLVE
      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
      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) = -2.0e-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 topography
      use cpflame_module
      implicit none

      end subroutine topography
