#include "options.inc"

c=======================================================================
c     Templates for setup
c
c        double periodic domain 
c=======================================================================


      subroutine set_parameter
      use cpflame_module
      implicit none
c ----------------------------------
c       set here main parameter
c ----------------------------------
      imt=32;km=10;jmt=imt
      dx=5e3
      dz=100.0
      dt_in=1200.0
      eps2D_sor = 1e-6
      eps3D_sor = 1e-6
      runlen=10*365.0
      snap_int = 0.05
      lat_ref=30
      beta   =0.

c      enable_noslip=.true.
      enable_hydrostatic         = .true.
      enable_cyclic_x            = .true.
      enable_cyclic_y            = .true.
      enable_biharmonic_friction = .true.
      enable_biharmonic_diffusion= .true.
      Ahbi = 2e8
      Khbi = Ahbi
      A_v=10e-4
      K_v=10e-4
      N_0 = 0.005
!      M_0 = sqrt(1e-5*0.2/rho_0*g)
      M_0 = 0.0
      end subroutine set_parameter



      subroutine initial_conditions
      use cpflame_module
      implicit none
      integer :: i,j,k
      real :: z,r,x,y,B0,fxa
c ----------------------------------
c      add here initial conditions
c ----------------------------------
      do k=1,km
       do j=1,jmt
        do i=1,imt

c         x = 8.*pi*xt(i)/(imt-2.)/dx
c         y = (yt(j)-yt(2))/ (yt(jmt-1)-yt(2))
c         B0=M_0**2*yt(j)+0.2*0.1/100*sin(x)

         x = (xt(i)-xt(imt/2))/20e3
         y = (yt(j)-yt(jmt/2))/20e3
         B0 = 2.0*exp(-x**2-y**2)*0.2/100
     &           *exp(-(zt(km)-zt(k))/200. )

c         y = 8.*pi*yt(j)/(jmt-2.)/dx
c         z = pi*zt(k)/(km-2.)/dz
c         B0 = sin(x)*sin(y)*cos(z)*2*0.2/100

         b(i,j,k,tau) = (B0-N_0**2*zt(k))*maskT(i,j,k)
        enddo
       enddo
      enddo
      call setcyclic3D(b(:,:,:,tau))
      b(:,:,:,taum1)=b(:,:,:,tau)

#ifdef notdef
c---------------------------------------------------------------------------------
c      hydrostatic pressure
c---------------------------------------------------------------------------------
      p_hydro = 0
      do j=2,jmt-1
       do k=km-1,1,-1
        p_hydro(:,j,k)= maskT(:,j,k)*
     &    (p_hydro(:,j,k+1)+(b(:,j,k+1,tau)+b(:,j,k,tau))/2.0*dz)
       enddo
      enddo
      call setcyclic3D(p_hydro)

      do k=2,km-1
       do j=2,jmt-1
        do i=2,imt-1
         fxa = 0.5*(p_hydro(i  ,j+1,k)-p_hydro(i  ,j,k))/dx
     &        +0.5*(p_hydro(i+1,j+1,k)-p_hydro(i+1,j,k))/dx
         u(i,j,k,1,tau) = fxa*maskU(i,j,k)/coriolis_t(j)
         fxa = 0.5*(p_hydro(i+1,j  ,k)-p_hydro(i,j  ,k))/dx
     &        +0.5*(p_hydro(i+1,j+1,k)-p_hydro(i,j+1,k))/dx
         u(i,j,k,2,tau) =-fxa*maskV(i,j,k)/coriolis_t(j)
        enddo
       enddo
      enddo
c
c      remove vertical mean from u and v
c
      p_hydro=0;
      do k=2,km-1
       do j=2,jmt-1
        do i=2,imt-1
          p_hydro(i,j,1)=p_hydro(i,j,1)+u(i,j,k,1,tau)*dz*maskU(i,j,k)
          p_hydro(i,j,2)=p_hydro(i,j,2)+u(i,j,k,2,tau)*dz*maskV(i,j,k)
          p_hydro(i,j,3)=p_hydro(i,j,3)+dz*maskU(i,j,k)
          p_hydro(i,j,4)=p_hydro(i,j,4)+dz*maskV(i,j,k)
        enddo
       enddo
      enddo

      p_hydro(:,:,3)=p_hydro(:,:,3)+1e-16
      p_hydro(:,:,4)=p_hydro(:,:,4)+1e-16

      do k=2,km-1
       do j=2,jmt-1
        do i=2,imt-1
         u(i,j,k,1,tau)=
     &       (u(i,j,k,1,tau)-p_hydro(i,j,1)/p_hydro(i,j,3))*maskU(i,j,k)
         u(i,j,k,2,tau)=
     &       (u(i,j,k,2,tau)-p_hydro(i,j,2)/p_hydro(i,j,4))*maskV(i,j,k)
        enddo
       enddo
      enddo
      p_hydro=0;
      call setcyclic3D(u(:,:,:,1,tau))
      call setcyclic3D(u(:,:,:,2,tau))
c
c      w from continuity
c
      u(:,:,1,3,:) = 0.0
      do k=2,km
       do j=2,jmt-1
        do i=2,imt-1
          u(i,j,k,3,tau) = u(i,j,k-1,3,tau)-maskW(i,j,k)*dz*( 
     &            (u(i,j,k,1,tau)-u(i-1,j,k,1,tau))/dx
     &           +(u(i,j,k,2,tau)-u(i,j-1,k,2,tau))/dx)
        enddo
       enddo
      enddo
      do j=2,jmt-1
       u(:,j,:,3,tau) = u(:,j,:,3,tau)*maskW(:,j,:)
      enddo
      call setcyclic3D(u(:,:,:,3,tau))

      u(:,:,:,:,taum1)=u(:,:,:,:,tau)
#endif
      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
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 ----------------------------------
      end subroutine boundary_conditions



      subroutine topography
      use cpflame_module
      implicit none
      integer :: i,k
      real :: z
c ----------------------------------
c      Add here topography  
c ----------------------------------
      end subroutine topography
