#include "options.inc"

c=======================================================================
c     Templates for setup
c
c      Eady problem
c=======================================================================


      subroutine set_parameter
      use cpflame_module
      implicit none
      real :: fac=1.0
c ----------------------------------
c       set here main parameter
c ----------------------------------
      imt = 32*fac; km  = 20 *fac; jmt = 30*fac

      dx    = 20e3/fac
      dz    = 100.0 /fac
      dt_in = 1200.0 /fac

      Ahbi  = 5e11/fac**2
      Khbi  = Ahbi
      A_v   = 10.e-4 
      K_v   = 1.0e-4 
      cdint = 1e-8 

      eps2D_sor = 1e-9
      runlen    = 50.0 !/ RESOLVE**3
      snap_int  = 0.5

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

!   f_0 = 2 om sin(phi/180*pi),  phi = 180/pi*asin(f_0/2/om)
      lat_ref = 180/pi*asin(1e-4/2/omega)
      beta    = 0*2*omega*cos(lat_ref/180.*pi)/radius
      N_0     = 0.004
      end subroutine set_parameter



      subroutine set_back_stratification(bb,u0)
c ----------------------------------
c    overide setting from module back_stratification
c ----------------------------------
      use cpflame_module
      implicit none
      real :: bb(imt,jmt,km),u0(imt,jmt,km),U_0=0.5
      integer :: i,j,k
c     thus best to set U=U_0 z/h  and B=f y U_0/h +f(z)
c     numerically consistent:  B_{j+1}=B_{j} + dx*U_z (coriolis_t(j  )+coriolis_t(j+1))/2.0

      do k=1,km
       u0(:,:,k)=(U_0/2+U_0*zt(k)/(km*dz))*maskU(:,:,k)
      enddo
      bb=0.
      do k=1,km
       do j=1,jmt-1
        bb(:,j+1,k)=bb(:,j,k)+dx*maskT(:,j,k)*U_0/(km*dz)
     &          *(coriolis_t(j)+coriolis_t(j+1))/2.0
       enddo
      enddo
      do k=1,km
       bb(:,:,k)=(bb(:,:,k)-N_0**2*zt(k))*maskT(:,:,k)
      enddo
      end subroutine set_back_stratification




      subroutine initial_conditions
      use cpflame_module
      implicit none
      integer :: i,j,k
      real :: fxa,kx,ky,L_r

      L_r = N_0*(km-2)*dz/(coriolis_t(jmt/2)*pi)
      kx=1./L_r
      kx=1.6*coriolis_t(jmt/2)/N_0/((km-2)*dz)
      ky=pi/((jmt-2)*dx)
      b=0;u=0
#ifndef enable_back_stratification
      call set_back_stratification(b(:,:,:,0),u(:,:,:,1,0))
      b(:,:,:,1)=b(:,:,:,0)
      b(:,:,:,2)=b(:,:,:,0)
      u(:,:,:,1,1)=u(:,:,:,1,0)
      u(:,:,:,1,2)=u(:,:,:,1,0)
#endif
      do k=1,km
       do j=1,jmt
        do i=1,imt
         fxa=1e-4*sin(kx*xt(i))*sin(ky*yt(j))!*cos(zt(k)/(km*dz)*pi)
         b(i,j,k,:)=b(i,j,k,:)+fxa*maskT(i,j,k)
        enddo
       enddo
      enddo

c      do i=-5,imt+5
c        print*,i,mod(i-1,imt)+1
c      enddo
c      stop
      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 :: T_rest,B0,om0
      type (time_type) :: time
c ----------------------------------
c      add here restoring zones for momentum
c ----------------------------------
#ifdef notdef
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      time = current_time
      om0=2*pi/coriolis_t(jmt/2)
      om0=2*pi/(0.98*om0)

      B0 =  1./(500.*dt)*sin(om0*time%seconds )
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         fv(i,j,k)= fv(i,j,k)-maskv(i,j,k)*B0
     &                *sin(pi*zt(k)/((km-2)*dz))
        enddo
       enddo
      enddo
#endif
      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
