#include "options.inc"


c#define with_topography

c=======================================================================
c     Templates for setup
c
c      zonally integrated idealised ACC
c=======================================================================

c increase this parameter to enhance resolution
#define RESOLVE  0.125


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

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

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

      K_v=1.0e-4!/RESOLVE
      A_v=10.e-4!/RESOLVE
      cdbot = 1e-5!*RESOLVE

      eps2D_sor = 1e-9
      runlen    = 365*100
      snap_int  = 50.0

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

      gamma=0.01
      lat_ref = -25.
      beta    = 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 :: bstar,B0,x,y
c ----------------------------------
c      add here initial conditions
c ----------------------------------
      do j=1,jmt
       call b_surface(bstar,j)
       do k=1,km
        do i=1,imt
         x=float(i)/float(imt)
         y = (yt(j)-yt(2))/ (yt(jmt-1)-yt(2))
         B0=0.0e-3*sin(x*8.5*pi)*exp(-(y-0.5)**2/0.5**2)
         b(i,j,k,:)  = 0*(bstar+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
c ----------------------------------
c      add here restoring zones
c ----------------------------------
      end subroutine restoring_zones


      subroutine momentum_restoring_zones
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je
      real :: fxa,fric(jmt)
      real, parameter :: r=1./(1*86400.)
      real, parameter :: L_x = 1280e3
      real, parameter :: epsilon = 0.5
c ----------------------------------
c      add zonal pressure gradient to 
c      int_{x1}^{x2}  u_t = - p(x2)+p(x1) 
c
c    u_t = f v - dp/dl - r u
c    v_t  = -fu - p_y -r v + ...
c ----------------------------------

#ifdef notdef
      js=max(2,js_pe); je = min(je_pe,jmt-1); i=2
      do k=2,km-1
        fric=0
        do j=je,max(jmt/2+1,js),-1
         fric(j)=fric(j+1)+dx*u(i,j,k,2,taum1)*maskV(i,j,k)
        enddo
        fric=fric*r/L_x
        do j=max(jmt/2,js+1),je
          fu(i,j,k) = fu(i,j,k)+maskU(i,j,k)*(fric(j)+fric(j+1))/2.0
        enddo
      enddo
#endif
#ifndef notdef
      js=max(2,js_pe); je = min(je_pe,jmt-2); i=2
      do k=2,km-1
        do j=max(jmt/2+1,js),je
          fxa=epsilon*(p_full(i,j+1,k,tau)-p_full(i,j,k,tau))/dx
          fu(i,j,k) = fu(i,j,k)-maskU(i,j,k)*fxa
        enddo
      enddo
#endif
      end subroutine momentum_restoring_zones


      subroutine b_surface(bstar,j)
      use cpflame_module
      implicit none
      integer :: j
      real :: db,bstar,y2,L_y
      L_y = yt(jmt)-yt(1)
      db = -30e-3
      bstar=db
      if (yt(j)<L_y/2.0) then
       bstar=db*yt(j)/(L_y/2.0)
      endif
      y2=L_y*0.75
      if (yt(j)>y2) then
       bstar=db*(1-(yt(j)-y2)/(yt(jmt)-y2) )
      endif
      end subroutine b_surface


      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 ----------------------------------
      L_y = yt(jmt)-yt(1)
      t_rest=30*86400
      do j=1,jmt
       call b_surface(bstar,j)
       do i=1,imt
        if (yt(j)<L_y/2.0) then
         y = yu(j)/L_y
         surf_tau(i,j,1) = .1e-3*sin(2*pi*y)*maskU(i,j,km-1)
        endif
        surf_flux(i,j)=dz/t_rest*(bstar-b(i,j,km-1,tau)) 
       enddo
      enddo
      end subroutine boundary_conditions


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