#include "options.inc"

c=======================================================================
c choose one of the following template setups
c 
c approved non-hydrostatic setups
c#include "../setup/kelv_helm1.F"
c#include "../setup/convection1.F"
c#include "../setup/rayleigh1.F"
c#include "../setup/mixedlayer1.F"
c#include "../setup/wave1.F"
c#include "../setup/shear1.F"
c
c approved hydrostatic setups
c
#include "../setup/barbi_exp1.F"
c#include "../setup/barbi_exp2.F"
c#include "../setup/wardle+marshall1.F"
c#include "../setup/wardle+marshall2.F"
c#include "../setup/edenetal1.F"
c#include "../setup/eady1.F"

c
c test setups use at your own risk
c#include "../setup/advection_test2.F"
c#include "../setup/small_channel1.F"
c#include "../setup/eddy1.F"
c#include "../setup/channel1.F"
c#include "../setup/channel2.F"
c#include "../setup/schwap1.F"
c#include "../setup/convection2.F"
c#include "../setup/afterconvection1.F"
c#include "../setup/theiss1.F"
c#include "../setup/balu1.F"
c#include "../setup/ACC1.F"
c#include "../setup/ACC_coarse1.F"
c#include "../setup/bolus1.F"

c#include "../setup/slantwise1.F"
c#include "../setup/slantwise2.F"
c#include "../setup/THC_hysterese1.F"
c#include "../setup/double_diff1.F"
c
c=======================================================================



      subroutine setup
c=======================================================================
c      setup cpflame
c=======================================================================
      use cpflame_module
      implicit none
      integer :: i,j,k,n,ierr,indp
      character(len=32) :: stamp
      logical :: file_exists = .false.

      if (my_pe==0) then
       print*,''
       print*,'--------------------------------------------------'
       print'(" setting up CPFLAME version ",f5.3)',version
       print*,'--------------------------------------------------'
       print*,''
      endif
c ----------------------------------
c       Set main parameter
c ----------------------------------
      call set_parameter
c ----------------------------------
c       Some more or less useful output
c ----------------------------------
      if (my_pe==0) then
       print*,''

       print*,''
       if (enable_hydrostatic) then
        print*,' hydrostatic model version'
       else
        print*,' nonhydrostatic model version'
       endif
       if (enable_free_surface) then
        print*,' using implicit linear free surface formulation'
       elseif (enable_expl_free_surf) then
        print*,' using explicit linear free surface formulation'
       else
        print*,' using rigid lid formulation'
       endif
       if (enable_quicker_advection) then
        print*,' using quicker advection scheme'
       else if (enable_4th_advection) then
        print*,' using 4th order advection scheme'
       else
        print*,' using second order advection scheme'
       endif
       if (enable_4th_mom_advection) then
        print*,' using 4th order advection scheme for momentum'
       elseif (enable_quicker_mom_advection) then
        print*,' using quicker advection scheme for momentum'
       else
        print*,' using second order advection scheme for momentum'
       endif
       if (enable_cyclic_x) then
        print*,' using cyclic zonal boundary conditions'
       else
        print*,' using solid zonal boundary conditions'
       endif
       if (enable_cyclic_y) then
        print*,' using cyclic meridional boundary conditions'
       else
        print*,' using solid meridional boundary conditions'
       endif
       if (enable_noslip) then
        print*,' using no slip lateral boundary conditions'
       else
        print*,' using free slip lateral boundary conditions'
       endif
       if (enable_bottom_noslip) then
        print*,' using no slip bottom boundary conditions'
       else
        print*,' using free slip bottom boundary conditions'
       endif

       print*,' grid size    : imt=',imt,' jmt=',jmt,' km=',km
       print*,' grid spacing : Delta x =',dx,'m  Delta z =',dz,'m'
       print*,' domain size  : ',imt*dx,'m X ',jmt*dx,'m X ',km*dz,'m'
       print*,' time step    : ',dt_in,'s'
       print*,' run length   : ',runlen*86400,'s, i.e. ',
     &        runlen*86400/dt_in,' time steps'
       print*,' snap interval: ',snap_int*86400,'s, i.e. any ',
     &        snap_int*86400/dt_in,' time steps'
       print*,' lateral  diffusivity : K_h=',K_h,' m^2/s'
       print*,' vertical diffusivity : K_v=',K_v,' m^2/s'
       print*,' lateral  viscosity   : A_h=',A_h,' m^2/s'
       print*,' vertical viscosity   : A_v=',A_v,' m^2/s'
       print*,' epsilon for 2D solver : ',eps2D_sor 
       print*,' f_0  = ',2*omega*sin(lat_ref/180.*pi),' 1/s'
       print*,' beta = ',beta,' 1/(ms)'
      endif

c ----------------------------------
c      time step splitting for explicit free surface
c ----------------------------------
       h_0 = (km-2)*dz
       dtex=dx/sqrt(g*h_0) /4.
       n=nint(dt_in/dtex)
       dtex=dt_in/n
       if (my_pe==0.and.enable_expl_free_surf) then
        print*,' free surface time step : ',dtex,' s'
       endif
c ----------------------------------
c      domain decomposition for each PE
c ----------------------------------
      if (n_pes>1) then
       n_pes_j = n_pes
       j_blk = (jmt-1)/n_pes_j + 1    ! j-extent of each block
       my_blk_j = my_pe + 1 ! number of PE in j-dir.
       js_pe = (my_blk_j-1)*j_blk + 1
       je_pe = min(my_blk_j*j_blk,jmt)
       call barrier
c ----------------------------------
c      check for incorrect domain decomposition
c ----------------------------------
       if (my_blk_j==n_pes_j .and. js_pe>=je_pe-2) then
        print*,' ERROR: on PE: ', my_pe
        print*,' domain decompositon impossible in j-direction'
        print*,' choose other number of PEs in j-direction'
        call halt_stop(' ... ')
       endif
      else
       n_pes_j = n_pes
       j_blk = jmt
       my_blk_j = 1 
       js_pe = 1
       je_pe = jmt
      endif
c ----------------------------------
c      print out the PE decomposition, let all PEs talk
c ----------------------------------
      if (my_pe==0) print*,' Domain decomposition:'
      do n=0,n_pes-1
       if (my_pe==n) then
         print*,''
         print*,' sub domain for PE #',n
         print*,' my_blk_j=',my_blk_j
         print*,' js_pe=',js_pe,' je_pe=',je_pe
         print*,''
       endif
       call barrier
      enddo
c ----------------------------------
c      allocate work space
c ----------------------------------
      call allocate_work_space
c ----------------------------------
c      setup grid
c ----------------------------------
      xt(1)=0.0
      xu(1)=xt(1)+dx/2.0
      do i=2,imt
        xt(i)=xt(i-1)+dx
        xu(i)=xu(i-1)+dx
      enddo
      yt(1)=0.0
      yu(1)=yt(1)+dx/2.0
      do i=2,jmt
        yt(i)=yt(i-1)+dx
        yu(i)=yu(i-1)+dx
      enddo
      zw(1)=0
      zt(1)=zw(1)-dz/2.0
      do k=2,km
        zt(k)=zt(k-1)+dz
        zw(k)=zw(k-1)+dz
      enddo
      zt=zt-dz*(km-2)
      zw=zw-dz*(km-2)
c ----------------------------------
c      set coriolis parameter
c ----------------------------------
      do j=1,jmt
       coriolis_t(j)  = 2*omega*sin(lat_ref/180.*pi)+beta*yt(j)
       coriolis_hor(j)= 2*omega*cos(lat_ref/180.*pi)
c     &             -2*omega*sin(lat_ref/180.*pi)/radius*yt(j)
      enddo
c
c       Reference buoyancy and pressure profiles for diagnostics
c
      b_r = yt(jmt/2)*M_0**2-zt*N_0**2
      p_r(km) = 0
      do k=km-1,1,-1
       p_r(k)=p_r(k+1)+0.5*(b_r(k+1)+b_r(k))*dz
      enddo
c ----------------------------------
c      Land mask
c ----------------------------------
      maskT=0.0
      maskT(2:imt-1,2:jmt-1,2:km-1)=1.0
      call topography
      if (enable_cyclic_x) maskT(1,:,:)=0.0
      if (enable_cyclic_x) maskT(imt,:,:)=0.0
      if (enable_cyclic_y) maskT(:,1,:)=0.0
      if (enable_cyclic_y) maskT(:,jmt,:)=0.0
      call setcyclic3D(maskT)
      maskU=maskT
      do i=1,imt-1
       maskU(i,:,:)=min(maskT(i,:,:),maskT(i+1,:,:))
      enddo
      call setcyclic3D(maskU)
      maskV=maskT
      do j=1,jmt-1
       maskV(:,j,:)=min(maskT(:,j,:),maskT(:,j+1,:))
      enddo
      call setcyclic3D(maskV)
      maskW=maskT
      do k=1,km-1
       maskW(:,:,k)=min(maskT(:,:,k),maskT(:,:,k+1))
      enddo
      call setcyclic3D(maskW)
c ----------------------------------
c     depth 
c ----------------------------------
      ht=0.0;hu=0.0;hv=0.0 
      do k=1,km
       ht=ht+maskT(:,:,k)*dz
       hu=hu+masku(:,:,k)*dz
       hv=hv+maskv(:,:,k)*dz
      enddo
      kmt=0;kmu=0;kmv=0
      do k=km,1,-1
       where (maskT(:,:,k)>0) kmt=k
       where (maskU(:,:,k)>0) kmu=k
       where (maskV(:,:,k)>0) kmv=k
      enddo
c ----------------------------------
c   initial conditions for rho
c ----------------------------------
      inquire( FILE='restart.dta', EXIST=file_exists )
      if (file_exists) then 
        call read_restart
      else
       if (my_pe==0) print*,' using initial conditions '
       call set_calendar_type(no_leap)
       restart_time_in = set_date(1900,1,1,0,0,0.0)
       u=0.0; b=0.0; p_surf=0.0;p_hydro=0.0;p_full=0.0;eta=0.0;psi=0.0
       if (enable_expl_free_surf) then
           bu=0; bv=0.
       endif
       call initial_conditions
      endif
      call setcyclic3D(b(:,:,:,0) )
      call setcyclic3D(b(:,:,:,1) )
      call setcyclic3D(b(:,:,:,2) )
c ----------------------------------
c   initialize diagnostic       
c ----------------------------------
      if (my_pe==0) call init_snap_cdf()

      end subroutine setup



