#include "options.inc"

c=======================================================================
c     Templates for setup
c
c      slantwise convection with meridional background flow
c=======================================================================

#define two_dim

      subroutine set_parameter
      use cpflame_module
      implicit none
      real :: fac=4.0
c ----------------------------------
c       set here main parameter
c       T = L^2/A_h 
c ----------------------------------
      imt=32*fac;km=20*fac;jmt=imt
#ifdef two_dim
      jmt=3
#endif
      dx=250/fac; 
      dz=20/fac; 
      dt_in=0.5*30/fac;
      eps2D_sor = 1e-6
      eps3D_sor = 1e-6
      runlen=10 
      snap_int = 0.02
      enable_hydrostatic         =.false.

c      enable_cyclic_x            =.true.
      enable_cyclic_y            =.true.
      enable_noslip  = .true.
c      enable_bottom_stress  = .true.
      enable_4th_advection   =.true.
c      enable_quicker_advection   =.true.
      enable_biharmonic_friction =.true. 
      enable_vert_biha_friction  =.true. 
      enable_biharmonic_diffusion=.true. 
      enable_vert_biha_diffusion =.true. 
      Ahbi = (2*2.0*(2.0/fac)**2)*dx**2
      Khbi = Ahbi
      Avbi = (2*0.005*(2.0/fac)**2)*dz**2
      Kvbi = Avbi

c      A_h=2.0*(2.0/fac)**2
c      K_h=A_h
c      A_v=0.005*(2.0/fac)**2
c      K_v=A_v
      lat_ref = 180/pi*asin(1e-4/2/omega)
      beta   =0.0
      N_0 = 0.007!(8e-4)
      M_0 = 0
      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)
      real :: U_0=-1* 0.8e-4/100*0.2*1000/1e-4 
      real :: hh = 800
      integer :: i,j,k

c      with b=g rho/rho_0!!!
c      f v = p_x,  p_z= -b,  f v_z = - b_x
c      V_z = V_0/hh , b_x = - f V_0/hh 
c
c     thus best to set U=U_0 z/h  and B=f y U_0/h +f(z)

       !M^2=0.2/1e4/100*0.2= B_y = -f U_0/h -> U_0 =0.2e-4 h/f
      do k=1,km ! u0 is meridional velocity
       u0(:,:,k)=(-U_0*zt(1)/hh+U_0*zt(k)/hh)*maskV(:,:,k)
      enddo
      bb=0.
      do k=1,km
       do i=1,imt-1
        bb(i+1,:,k)=bb(i,:,k)-dx*maskT(i,:,k)*U_0/hh*coriolis_t(1)
       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 :: B0,x,fxa,z  ! temp = b/0.2*tho_0/g
c ----------------------------------
c      add here initial conditions
c ----------------------------------
#ifndef enable_back_stratification
      call set_back_stratification(b(:,:,:,0),u(:,:,:,2,0))
      b(:,:,:,1)=b(:,:,:,0)
      b(:,:,:,2)=b(:,:,:,0)
      u(:,:,:,2,1)=u(:,:,:,2,0)
      u(:,:,:,2,2)=u(:,:,:,2,0)
#endif
      do k=1,km
       do j=1,jmt
        do i=1,imt
         b(i,j,k,:) =b(i,j,k,:)+1e-8*sin(xt(i)/xt(imt)*2.5*pi)
     &         *cos(zt(k)/zt(1)*3*pi) *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 ----------------------------------
      integer :: i,j,k,js,je
      real :: T_rest,B0,hh
      js=max(2,js_pe); je = min(je_pe,jmt-1)

      T_rest=1./(5.*dt)
      hh=abs(zt(1))/10.;
      do k=1,km
       B0=t_rest*max(0.,1.-(zt(k)-zt(1))/hh)
       do i=1,imt
        do j=js,je
         b(i,j,k,taup1)=b(i,j,k,taup1)+maskT(i,j,k)*c2dt*
     &                 B0*(0-b(i,j,k,taum1))
        enddo
       enddo
      enddo
      end subroutine restoring_zones


      subroutine momentum_restoring_zones
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je
      real :: T_rest,B0,hh

c ----------------------------------
c      add here restoring zones for momentum
c ----------------------------------
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      T_rest=1./(5.*dt)
      hh=abs(zt(1))/10.;
      do k=1,km
       B0=t_rest*max(0.,1.-(zt(k)-zt(1))/hh)
       do i=1,imt
        do j=js,je
         fu(i,j,k)= fu(i,j,k)+maskU(i,j,k)*
     &            B0*(0.-u(i,j,k,1,taum1))
         fv(i,j,k)= fv(i,j,k)+maskV(i,j,k)*
     &            B0*(0.-u(i,j,k,2,taum1))
         fw(i,j,k)= fw(i,j,k)+maskW(i,j,k)*
     &            B0*(0.-u(i,j,k,3,taum1))
        enddo
       enddo
      enddo
      end subroutine momentum_restoring_zones



      subroutine boundary_conditions
      use cpflame_module
      implicit none
      integer :: i,j
      real :: x,y
      real :: f0  = 1000 *0.2/rho_0*g /4185.5 ,fxa

c     heat = rho cp T, heat flux= rho cp T v = H
c     temp flux T = H /rho_0 /cp [ J/s/m^2 /kg m^3 /J Kg K=m/s K ]
c     density flux = R = T drho/dt [ K m/s kg/m^3 /K ]
c     buoyancy flux = B = -R/\rho_0 g [ kg/m^3 m/s ]
c     B = -H/rho_0 /cp drho/dt /rho_0 g 
c       = 100/1000^2/4185 0.2 *10=5e-8 

c       b=0.2Temp/rho_0*g    ! temp = b/0.2*tho_0/g
c      surf_flux  = 50 *0.2/rho_0*g /4185.5 
      surf_flux=0
      do j=1,jmt
       do i=1,imt
c         surf_flux(i,j)=f0/rho_0
       enddo
      enddo
      end subroutine boundary_conditions



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