#include "options.inc"


c#define with_topography

c=======================================================================
c     Templates for setup
c
c      Wide eddying channel with restoring zones at side walls
#ifdef with_topography
c      reads topography from file
#endif
c=======================================================================

c increase this parameter to enhance resolution
#define RESOLVE  1.0

c 2D version
#define two_dim

#ifdef notdef
c theiss1d_interior:   dt=1200s,  sigma=6.8, L_r=68
c      cdint = 0.5e-7*RESOLVE
      N_0     = 0.004*0.5*2**0.5
      M_0     = sqrt(1e-5*0.1/rho_0*g)*0.5**0.5
c theiss2_interior.cdf dt=2400s, sigma=6.8 L_r=135km     
      N_0     = 0.004*sqrt(2.0)
      M_0     = sqrt(1e-5*0.1/rho_0*g)
c theiss2d_interior         sigma=6.8  L_r= 191
      N_0     = 0.79978099E-02
      M_0     = 0.11667403E-03
c theiss1g_interior         sigma=6.8  L_r=286
      N_0     = 0.004*3.0
      M_0     = sqrt(1e-5*0.1/rho_0*g)*2.15**0.5
c theiss3b_interior         sigma=5.8  L_r= 57
      N_0     = 0.0048*0.5
      M_0     = sqrt(1e-5*0.1/rho_0*g)*0.5**0.5
c theiss3_interior         sigma=5.8  L_r=115
     N_0     = 0.0048
      M_0     = sqrt(1e-5*0.1/rho_0*g)
c theiss3d_interior         sigma=5.8  L_r=191
      N_0     = 0.79978099E-02
      M_0     = 0.12633236E-03
c theiss3c_interior         sigma=5.8  L_r=229
      N_0     = 0.0048*2.0
      M_0     = sqrt(1e-5*0.1/rho_0*g)*2.0**0.5
c theiss1b_interior  dt=1200      sigma=4.8  L_r=48
      N_0     = 0.004*0.5
      M_0     = sqrt(1e-5*0.1/rho_0*g)*0.5**0.5
c theiss1_interior   dt=2400     sigma=4.8  L_r=96
      N_0     = 0.004
      M_0     = sqrt(1e-5*0.1/rho_0*g)
c theiss1c_interior         sigma=4.8  L_r=191
      N_0     = 0.004*2.0
      M_0     = sqrt(1e-5*0.1/rho_0*g)*2**0.5
c theiss1h_interior         sigma=4.8  L_r=239
      N_0     = 0.004*2.5
      M_0     = sqrt(1e-5*0.1/rho_0*g)*2.5**0.5
c theiss2c_interior         sigma=3.2  L_r=72
      N_0     = 0.004*0.75
      M_0     = sqrt(1e-5*0.1/rho_0*g)*sqrt(1.5)*sqrt(0.75)
c theiss2b_interior         sigma=3.2  L_r=96
      N_0     = 0.004
      M_0     = sqrt(1e-5*0.1/rho_0*g)*sqrt(1.5)
c theiss1f_interior    dt=1200     sigma=3.2  L_r=191
      N_0     = 0.004*2.0
      M_0     = sqrt(1e-5*0.1/rho_0*g)*3**0.5
#endif

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

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

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

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


      eps2D_sor = 1e-9
      runlen    = 365*5. !/RESOLVE**3
      snap_int  = 3.0

      enable_noslip               = .true.
      enable_cyclic_x             = .true.
c      enable_bottom_stress        = .true.
      enable_bottom_stress        = .false.
      enable_interior_stress      = .true.
      enable_biharmonic_friction  = .true.
      enable_biharmonic_diffusion = .true.
      enable_hydrostatic          = .true.

      gamma=0.01
      lat_ref = 10.
      beta    = 2*omega*cos(lat_ref/180.*pi)/radius
c      lat_ref = 180/pi*asin(0.7e-4/2/omega)
c      beta    = 0.
      N_0     = 0.004
      M_0     = sqrt(1e-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*,' Theiss setup :'
       print*,' '
       print*,' Prandtl number '
       print*,'   f^2/N^2       = ',f_0**2/N_0**2
       print*,' Isopycnal slope '
       print*,'   s=M^2/N^2     = ',s
       print*,' Aspect ratio '
       print*,'   delta=H/L     = ',delta
       print*,' Rossby number '
       print*,'   H M^2/(L f^2) = ',Ro
       print*,' Burger number '
       print*,'   Ro delta/s    = ',Ro*delta/s
       print*,' '
       print*,' Mean velocity scale'
       print*,'  H M^2/f        = ',M_0**2*(km-2)*dz/f_0,' m/s'
       print*,' Inverse Eady growth rate '
       print*,'   tau=N/M^2     = ',N_0/M_0**2/(60*60*24),' days' 
       print*,' Rossby radius '
       print*,'  N h/(f pi)     = ',N_0*(km-2)*dz/(f_0*pi)/1e3,' km'
       print*,' Rhines scale  '
       print*,'  (H M^2/(beta f))^0.5 = ',
     &               (M_0**2*(km-2)*dz/(f_0*beta))**0.5/1e3,' km'
       print*,' '
      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-3*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
c ----------------------------------
c    add here boundary conditions
c ----------------------------------
      end subroutine boundary_conditions



      subroutine topography
      use cpflame_module
      implicit none

#ifdef with_topography
c ----------------------------------
c      Add here topography  
c ----------------------------------
      logical :: file_exists
      integer :: io=37
      integer :: i,j,k, buf(imt,jmt)

      if (my_pe==0) print*,' reading topography from file'
      inquire( FILE='kmt_theiss2.dta', EXIST=file_exists )
      if (.not.file_exists) then 
        call halt_stop(' file not found in topography ')
      endif
      open(io,file='kmt_theiss2.dta',form='unformatted',status='old')
      read(io) i,j,k
      if (i/=imt.or.j/=jmt.or.k/=km) then
        call halt_stop(' dimensions do not match in topography ')
      endif
      read(io) buf; kmt=buf
      close(io)
      if (my_pe==0) print*,' done reading topography'
      maskT=0.0
      do j=1,jmt
       do i=1,imt
        k=kmt(i,j)
        if (k>0) maskT(i,j,k:km-1)=1.0
       enddo
      enddo
#endif
      end subroutine topography
