c
c-----------------------------------------------------------------------
c     file contains main module containing all global variables
c
c-----------------------------------------------------------------------
c
      module  barbi_module
c=======================================================================
c     global variables and arrays for BARBI
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use time_manager_module
      implicit none
c
      real (kind=8), parameter :: version = 1.0 ! version number
c
c     physical constants
c
      real(kind=8) :: h_0    ! max. water depth in m
      real(kind=8) :: rho0   ! reference density in kg/m^3
      real(kind=8) :: radius ! mean radius of earth in m
      real(kind=8) :: omega  ! earth rotation frequency in s
      real(kind=8) :: pi     ! 3.14...
c
c     model parameter
c
      integer      :: imt,jmt,n_order !  size of domain, number of moments
      real(kind=8) :: N0,N0sqr        ! stability freq. and square of it
      real(kind=8) :: Am,Ah,Av        ! hor. visc./diff. and vert. visc. in m^2/s
      real(kind=8), allocatable :: gamma(:,:,:) ! truncation parameter
      real(kind=8), allocatable :: gscal(:,:)   ! scaling for barocl. modes
      real(kind=8) :: eps_solver   = 1.0        ! convergence crit. for streamfct.
      real(kind=8) :: eps_solver_p = 1e-6       ! convergence crit for bottom press.
      real(kind=8) :: time_step  ! time step in seconds
      real(kind=8) :: runlen     ! length of experiment in days
      real(kind=8) :: snapint    ! snapshot interval in days
      integer      :: jlat0      ! reference latitude for beta/f plane
c
c     grid and related variables, see MOM 2.1/2.2/3.0 Manual for explanation
c
      real(kind=8),allocatable :: xt(:), yt(:), xu(:), yu(:)
      real(kind=8),allocatable :: dxtdeg(:), dytdeg(:) 
      real(kind=8),allocatable :: dxudeg(:), dyudeg(:)
      real(kind=8),allocatable :: dxt(:), dyt(:), dxu(:), dyu(:)
      real(kind=8),allocatable :: dytr(:),dyt2r(:),dyur(:), dyu2r(:)
      real(kind=8),allocatable :: dxtr(:),dxt2r(:),dxur(:), dxu2r(:)
      real(kind=8),allocatable :: phi(:),phit(:),cst(:),csu(:),sine(:)
      real(kind=8),allocatable :: cstr(:),csur(:),tng(:)
      real(kind=8),allocatable :: cstdytr(:), cstdyt2r(:)
      real(kind=8),allocatable :: csudyur(:),csudyu2r(:)
      real(kind=8),allocatable :: cst_dytr(:),csu_dyur(:)
      real(kind=8),allocatable :: dus(:),dun(:),duw(:),due(:),dxmetr(:)
c
c     land mask
c
      real(kind=8),allocatable :: tmask(:,:),umask(:,:)
c
c     (inverse of) topography
c
      real(kind=8),allocatable :: h(:,:),hu(:,:),hr(:,:),hur(:,:)
c
c     coriolis parameter
c
      real(kind=8),allocatable :: cori(:,:,:)
c
c     short hands for mixing parameter
c
      real(kind=8),allocatable :: am3(:),am4(:,:)
      real(kind=8),allocatable :: ahc_north(:),ahc_south(:)
      real(kind=8),allocatable :: amc_north(:),amc_south(:)
c
c     coefficients to minimize advection and diffusion computation 
c
      real(kind=8),allocatable :: csudxur(:,:),csudxu2r(:,:)
      real(kind=8),allocatable :: cstdxtr(:,:),cstdxt2r(:,:)
      real(kind=8),allocatable :: am_csudxtr(:,:), ah_cstdxur(:,:)
c
c     advection velocity on northern and eastern face of "T" cells 
c
      real(kind=8),allocatable :: adv_vnt(:,:), adv_vet(:,:)
c
c     pressure gradient and so on
c
      real(kind=8),allocatable :: zu(:,:,:), grad_E(:,:,:), press(:,:)
      real(kind=8),allocatable :: psi(:,:,:),dpsi(:,:),guess(:,:)
c
c     forcing of the momentum
c
      real(kind=8),allocatable :: wind(:,:,:)
c
c     the barotropic velocities
c
      real(kind=8),allocatable :: u(:,:,:,:)
c
c     the baroclinic velocities
c
      real(kind=8),allocatable :: ub(:,:,:,:,:)
c
c     the moments of the vertical integrated density
c
      real(kind=8),allocatable :: g(:,:,:,:)
c
c     some switches for configurations
c
      logical :: cyclic           = .false.   ! cyclic boundary conditions if true
      logical :: enable_obc_north = .false.   ! open northern boundary
      logical :: enable_obc_south = .false.   ! south
      logical :: enable_obc_west  = .false.   ! western 
      logical :: enable_obc_east  = .false.   ! eastern
      logical :: enable_beta_plane  = .false. ! use beta plane approx
      logical :: enable_f_plane     = .false. ! f-plane
c
c     diagnostics
c
      logical :: read_restart  = .false. ! read from restart file at begin of run
      logical :: write_restart = .true.  ! write to restart file at end of run
      logical :: write_modes   = .true. ! write baroclinic modes instead of density moments
      character (len=80):: snap_file,restart_file ! filenames for snapshots and restart
      integer :: mscans    ! iteration needed for poisson solver (streamfct)
      integer :: mscans_p  ! same but solving for bottom press
      real(kind=8)    :: ek,ke ! barotrop. kinetic energy, pot energy

      contains

      subroutine  init_barbi_module
c=======================================================================
c     Read in namelist parameter and initialize all arrays in module 
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      implicit none
      namelist /barbi/ 
     &     imt,jmt,n_order,
     &     time_step,
     &     ah,am,Av,N0,runlen,snapint,
     &     snap_file,write_modes,
     &     eps_solver, eps_solver_p,
     &     read_restart,write_restart,restart_file,
     &     cyclic, 
     &     enable_obc_north, enable_obc_south, 
     &     enable_obc_west,  enable_obc_east, 
     &     enable_beta_plane, enable_f_plane,jlat0

      read(5, nml = barbi)
c
c     truncation parameter
c
      allocate( gamma(imt,jmt,n_order) )
      allocate( gscal(n_order/2+1,n_order/2+1) ); gscal = 0.

c     grid
      allocate( xt(imt), yt(jmt), xu(imt), yu(jmt))
      allocate( dxtdeg(imt), dytdeg(jmt), dxudeg(imt), dyudeg(jmt))
      allocate( dxt(imt), dyt(jmt), dxu(imt), dyu(jmt))
      allocate( dytr(jmt),dyt2r(jmt), dyur(jmt),dyu2r(jmt))
      allocate( phi(jmt),phit(jmt), cst(jmt),csu(jmt),sine(jmt))
      allocate( cstr(jmt),csur(jmt),tng(jmt))
      allocate( cstdytr(jmt), cstdyt2r(jmt),csudyur(jmt),csudyu2r(jmt))
      allocate( cst_dytr(jmt),csu_dyur(jmt))
      allocate( dus(jmt),dun(jmt))
      allocate( dxtr(imt),dxt2r(imt),dxur(imt), dxu2r(imt))
      allocate( dxmetr(imt),duw(imt),due(imt))
c     land mask
      allocate( tmask(imt,jmt),umask(imt,jmt) )
      allocate( h(imt,jmt),hu(imt,jmt), hr(imt,jmt),hur(imt,jmt) )
c     coriolis parameter
      allocate( cori(imt,jmt,2) )
c     mixing parameter
      allocate( am3(jmt),am4(jmt,2) )
      allocate( amc_north(jmt),amc_south(jmt) )
      allocate( ahc_north(jmt),ahc_south(jmt) )
c     coefficients to minimize advection and diffusion computation 
      allocate( csudxur(imt,jmt),csudxu2r(imt,jmt),am_csudxtr(imt,jmt) )
      allocate( cstdxtr(imt,jmt),cstdxt2r(imt,jmt),ah_cstdxur(imt,jmt) )
c     advection velocity on northern and eastern face of "T" cells 
      allocate( adv_vet(imt,jmt), adv_vnt(imt,jmt) )
c     pressure gradient and so on
      allocate( zu(imt,jmt,2) )
      allocate( grad_E(imt,jmt,2) )
      allocate( press(imt,jmt))
      allocate( psi(imt,jmt,2),dpsi(imt,jmt),guess(imt,jmt) )
      guess=0.;dpsi=0.
c     forcing of the velocity
      allocate( wind(imt,jmt,2) )
c     the barotropic velocities
      allocate( u(imt,jmt,2,0:2) )
c     the baroclinic velocities
      allocate( ub(imt,jmt,2,0:2,0:n_order) )
c     the moments of the vertical integrated density
      allocate( g(imt,jmt,0:n_order,0:2) )

      end subroutine init_barbi_module

      end module barbi_module
