#include "options.inc"



      module cpflame_module
c=======================================================================
c      Main module for CPFLAME
c=======================================================================
      use time_manager_module
      use timing_module
      use mpp_module
      implicit none
      real, parameter :: version = 0.98
c---------------------------------------------------------------------------------
c      physical parameter
c---------------------------------------------------------------------------------
      real, parameter :: pi     = 3.14159265358979323846264338327950588
      real, parameter :: radius = 6370.0e3        ! Earth radius in m
      real, parameter :: degtom = radius/180.0*pi ! conversion degrees latitude to meters
      real, parameter :: omega  = pi/43082.0      ! earth rotation frequency in 1/s
      real, parameter :: rho_0  = 1024            ! Reference density in kg/m^3
      real, parameter :: g      = 9.806           ! Earth accelaration constant in m/s^2
c---------------------------------------------------------------------------------
c      variable model parameter
c---------------------------------------------------------------------------------
      integer  :: imt=0,km=0,jmt=0      ! Number of grid points in x,y,z direction
      real :: dx=0,dz=0,dt_in=0        ! Delta z, Delta x in m and Delta t in s
      real :: K_h=0, K_v=0             ! horizontal and vertical diffusivitis in m^2/s
      real :: A_h=0, A_v=0             ! horizontal and vertical viscosities in m^2/s
      real :: cdbot=0                  ! Bottom friction in 1/s
      real :: cdint=0                  ! interior friction in 1/s
      real :: Khbi = .0, Kvbi =.0      ! hor/vert biharmonic diffusivities in m^4/s
      real :: Ahbi = .0, Avbi=.0       ! same for biharmonic viscosities
      real :: runlen=0,snap_int=0      ! length of integration and snapshot interval in days 
      real :: eps2D_sor=0,eps3D_sor=0. ! cut off values for Poisson solver
      real :: lat_ref = 45.0           ! Reference values for f/beta plane in 1/s
      real :: beta    = 2e-11          ! beta for beta plane in 1/m/s
      real :: gamma=0.01               ! value for Roberts time filter
      real :: dtex                         ! time step for external mode 
      real :: h_0                          ! maximal possible depth
c
c      parameters for setup and diagnostics (relevant for for hydrostatic version only)
c
      real :: N_0=0.0          ! Square root of constant reference vertical buoyancy gradient
      real :: M_0=0.0          ! Square root of constant meridional reference buoyancy gradient
      real, allocatable :: b_r(:),p_r(:) ! Veetical reference buoyancy and pressure profiles
c---------------------------------------------------------------------------------
c      domain decomposition (see also mpp_module)
c---------------------------------------------------------------------------------
      integer :: n_pes_j = 1 ! number of PEs in meridional (y) direction
      integer :: my_blk_j= 1 ! number of domain in meridional (y) direction
      integer :: js_pe   = 1 ! start j-index of domain of this PE
      integer :: je_pe   = 1 ! end of j_index of domain of this PE
      integer :: j_blk   = 1 ! number of j-rows at this PE
c---------------------------------------------------------------------------------
c      main model variables, units are SI
c---------------------------------------------------------------------------------
      real, allocatable, dimension(:,:,:,:,:) :: u       ! velocicity components for three time levels
      real, allocatable, dimension(:,:,:,:)   :: b       ! (negative) buoyancy for three time levels
      real, allocatable, dimension(:,:,:,:)   :: p_full  ! full pressure for three time levels
      real, allocatable, dimension(:,:,:)     :: p_hydro ! hydrostatic pressure at time level tau
      real, allocatable, dimension(:,:,:)     :: psi     ! non-hdrostatic pressure at time level tau
      real, allocatable, dimension(:,:)       :: p_surf  ! surface pressure
      real, allocatable, dimension(:,:,:)     :: eta     ! free surface elevation
      real, allocatable, dimension(:,:,:)     :: bu     !  barotropic zonal velocity         (only for enable_expl_free_surf)
      real, allocatable, dimension(:,:,:)     :: bv     !  barotriópic meridonal velocity (only for enable_expl_free_surf)
c
c      auxilliary model variables 
c
      real, allocatable, dimension(:,:,:)     :: K_b     ! vertical diffusivity set in case of static instability
      real, allocatable, dimension(:,:,:)     :: maskT   ! mask for T-grid (1 in water, 0 at land)
      real, allocatable, dimension(:,:,:)     :: maskU   ! same for U-grid
      real, allocatable, dimension(:,:,:)     :: maskV   ! same for V-grid
      real, allocatable, dimension(:,:,:)     :: maskW   ! same for W-grid
      real, allocatable, dimension(:,:,:)     :: fu      ! time tendency of u except pressure gradient
      real, allocatable, dimension(:,:,:)     :: fv      ! time tendency of v except pressure gradient
      real, allocatable, dimension(:,:,:)     :: fw      ! time tendency of w except pressure gradient
      real, allocatable, dimension(:)         :: coriolis_t    ! vertical coriolis parameter
      real, allocatable, dimension(:)         :: coriolis_hor  ! horizontal coriolis parameter
      real, allocatable,dimension(:,:)        :: surf_flux     ! surface buoyancy flux
      real, allocatable,dimension(:,:)        :: bottom_flux   ! bottom buoyancy flux
      real, allocatable,dimension(:,:,:)      :: surf_tau      ! surface flux of momentum (u and v component)
      real, allocatable,dimension(:,:,:)      :: bott_tau      ! bottmom flux of momentum (u and v component)
      real, allocatable :: cf2D(:,:,:,:)      ! array containing coefficients to solve 2D Poisson equation
      real, allocatable :: cf3D(:,:,:,:,:,:)  ! array containing coefficients to solve 3D Poisson equation
c
c      model grid variables  
c
      real, allocatable, dimension(:)   :: zt  ! depth of T-grid point  (positive upward) in meters
      real, allocatable, dimension(:)   :: zw  ! depth of W-grid point in meters
      real, allocatable, dimension(:)   :: xt  ! zonal (x) coordinate of T-grid point in meters
      real, allocatable, dimension(:)   :: xu  ! zonal (x) coordinate of U-grid point in meters
      real, allocatable, dimension(:)   :: yt  ! meridional (y) coordinate of T-grid point in meters
      real, allocatable, dimension(:)   :: yu  ! meridional (y) coordinate of V-grid point in meters
      real, allocatable, dimension(:,:) :: hu  ! depth of water column at U-grid points
      real, allocatable, dimension(:,:) :: hv  ! depth of water column at V-grid points
      real, allocatable, dimension(:,:) :: ht  ! depth of water column at T-grid points
      integer,allocatable,dimension(:,:) :: kmt ! number of verticel levels at T grid points
      integer,allocatable,dimension(:,:) :: kmu ! number of verticel levels at U grid points
      integer,allocatable,dimension(:,:) :: kmv ! number of verticel levels at V grid points
c---------------------------------------------------------------------------------
c       switches to configure model setup
c---------------------------------------------------------------------------------
      logical :: enable_noslip              = .false.  ! no lateral slip for momentum
      logical :: enable_bottom_noslip  = .false. ! also for bottom 
      logical :: enable_cyclic_x            = .false.  ! periodic boundary condition in x
      logical :: enable_cyclic_y            = .false.  ! periodic boundary condition in y
      logical :: enable_quicker_advection   = .false.  ! quicker advection scheme for buoyancy
      logical :: enable_4th_advection       = .false.  ! 4. order advection scheme for buoyancy
      logical :: enable_upwind_advection    = .false.  ! upwind advection scheme for buoyancy
      logical :: enable_lax_advection       = .false.  
      logical :: enable_4th_mom_advection = .false. ! momentum advection
      logical :: enable_quicker_mom_advection = .false. ! momentum advection
      logical :: enable_bottom_stress       = .false.  ! Bottom stress
      logical :: enable_interior_stress     = .false.  ! interior stress
      logical :: enable_biharmonic_friction = .false.  ! biharmonic horizontal friction
      logical :: enable_vert_biha_friction  = .false.  ! biharmonic vertical  friction
      logical :: enable_biharmonic_diffusion= .false.  ! same for diffusion
      logical :: enable_vert_biha_diffusion = .false.  ! biharmonic vertical friction
      logical :: enable_free_surface        = .false.  ! implicit free surface
      logical :: enable_hydrostatic         = .true.   ! hydrostatic approximation
      logical :: enable_expl_free_surf     =.false.  ! explicit free surface
c---------------------------------------------------------------------------------
c      other stuff
c---------------------------------------------------------------------------------
      type(time_type) :: restart_time_in   ! Time of last restart
      integer           :: restart_itt=0     ! number of time steps since last restart
      real             :: spval = -9.9e12   ! special value in netcdf files
      integer           :: sor2D_itts        ! number of iterations of 2D Poisson solver
      integer           :: sor3D_itts        ! number of iterations of 3D Poisson solver
      end module cpflame_module



      subroutine allocate_work_space
c=======================================================================
c      Allocate memory for all arrays in main module
c=======================================================================
      use cpflame_module
      implicit none
      allocate(zt(km),zw(km),xt(imt),xu(imt),yt(jmt),yu(jmt) )
      allocate( u(imt,jmt,km,3,0:2) ); u=0.
      allocate( b(imt,jmt,km,0:2)   ); b=0.
      allocate( K_b(imt,jmt,km)  );    K_b=0.
      allocate( p_hydro(imt,jmt,km) ); p_hydro=0.
      allocate( p_full(imt,jmt,km,0:2)  ); p_full=0.
      allocate( p_surf(imt,jmt)     ); p_surf=0.
      allocate( psi(imt,jmt,km)     ); psi=0.
      allocate( eta(imt,jmt,0:2)    ); eta=0.
      allocate( maskT(imt,jmt,km)   ); maskT=0.
      allocate( maskU(imt,jmt,km)   ); maskU=0.
      allocate( maskV(imt,jmt,km)   ); maskV=0.
      allocate( maskW(imt,jmt,km)   ); maskW=0.
      allocate( fu(imt,jmt,km)      ); fu=0.
      allocate( fv(imt,jmt,km)      ); fv=0.
      allocate( fw(imt,jmt,km)      ); fw=0.
      allocate( surf_flux(imt,jmt)  ); surf_flux=0.
      allocate( bottom_flux(imt,jmt)  ); bottom_flux=0.
      allocate( surf_tau(imt,jmt,2) ); surf_tau=0.
      allocate( bott_tau(imt,jmt,2) ); bott_tau=0.
      allocate(coriolis_t(jmt), coriolis_hor(jmt) )
      coriolis_t=0.;coriolis_hor=0.
      allocate( hu(imt,jmt),hv(imt,jmt),ht(imt,jmt))
      hu=0.;hv=0.;ht=0.
      allocate( kmt(imt,jmt),kmu(imt,jmt),kmv(imt,jmt) )
      kmt=0;kmu=0;kmv=0
      allocate( cf2D(imt,jmt,-1:1,-1:1) )
      if (.not.enable_hydrostatic) then
        allocate( cf3D(imt,jmt,km,-1:1,-1:1,-1:1))
      endif
      allocate( b_r(km), p_r(km) ); b_r=0.; p_r=0.
      if (enable_expl_free_surf) then
       allocate( bu(imt,jmt,0:2)    ); bu=0.
       allocate( bv(imt,jmt,0:2)    ); bv=0.
      endif
      end subroutine allocate_work_space
 

