#include "options.inc"



      module spflame_module

c==============================================================
c     Main module for SPFLAME
c                               Jul 2001   c.eden
c==============================================================

      use time_manager_module
      use dens_module
      use timing_module
      implicit none

      real, parameter :: version = 1.0

!     set some physical and mathematical constants       
      real, parameter :: rho0     = 1.035
      real, parameter :: rho0r    = 1./rho0
      real, parameter :: grav     = 980.6
      real, parameter :: radius   = 6370.0e5
      real, parameter :: pi       = 
     &        3.14159265358979323846264338327950588
      real, parameter :: omega    = pi/43082.0

!     extent of model domain, number of tracers
      integer            :: imt=0,jmt=0,km=0,nt=0 
!     time step in sec., length of run in days
      real               :: time_step=0.,runlen=0.
!     restart stuff, time stamp of last restart
      character (len=32) :: restart_stamp    
      integer            :: itt_restart  ! time step of last restart
      logical :: read_restart = .false.  ! read from restart file

!     PE number and number of PEs in sub domain
      integer :: my_pe,n_pes      
!     number of all PEs in i/j-direction
      integer :: n_pes_i=0,n_pes_j=0
!     hor. extent of PE domains in i/j-direction
      integer :: i_blk,j_blk      
!     hor. start/end index of this PE domain
      integer :: is_pe,ie_pe,js_pe,je_pe  
!     number of this PE in i/j-dir.
      integer :: my_blk_j,my_blk_i        
!     max number of sub domain and actual number
      integer :: domains   
      integer, parameter :: max_domains=100
      integer :: sub_domain_pe0(max_domains) ! first PE of each sub domain
!     mpi communicator for sub domain and number of this sub domain
!                                         (only used for IO file name)
      integer :: my_comm,sub_domain       
      integer :: mother = -1, childs(max_domains) ! number of the mother domain 
      integer :: nr_childs = 0                    ! and numbers of the child domains
!     factor for zooming of sub domain and indices of sub domain in mother domain
      integer :: zoom_fac=0,is_zoom=0,ie_zoom=0,js_zoom=0,je_zoom=0
!     zoom factor in k-dir. and number of levels in mother domain
      integer :: zoom_fac_k = 1, km_mother 
!     smooth mother-topography for sub domain nn times
      integer :: topo_smooth_iterations = 1 


!     general setup
      logical :: cyclic              = .false. ! cyclic boundary conditions
      logical :: enable_freesurf     = .false. ! free surface or rigid lid
      logical :: enable_beta_plane   = .false. ! beta plane
      integer :: beta_plane_j0       = 1       ! reference grid point 
      logical :: enable_f_plane      = .false. ! f plane
      logical :: enable_rotated_grid = .false. ! use a rotated grid a binary file of
                                               ! unrotated latitudes have to be provided
      logical :: enable_simple_grid = .false. ! generate a simple grid, otherwise 
                                              ! binary file will be read, 
                                              ! see also file setup_template.F
      logical :: enable_simple_topo = .false. ! same for topography
      logical :: enable_simple_initial_cond = .false. 
      logical :: enable_simple_sbc = .false. 
      logical :: enable_simple_obc = .false. 
      logical :: enable_simple_spg = .false. 
      logical :: enable_sponge     = .false.   ! try to read sponge data from file

!   SBC setup
      logical :: enable_salt_flux_sbc            = .false.  ! salt flux instead of restoring
      logical :: enable_shortwave_sbc            = .false.  ! penetrative solar radiation

!     open boundary setup
      logical :: enable_obc_north = .false.  ! open northern boundary 
      logical :: enable_obc_south = .false.
      logical :: enable_obc_west  = .false.
      logical :: enable_obc_east  = .false.

      logical :: restore_TS_obc_north = .false. ! restore to TS field in addition
      logical :: restore_TS_obc_south = .false. ! to the radiation condition
      logical :: restore_TS_obc_west  = .false.
      logical :: restore_TS_obc_east  = .false.

      logical :: prescribe_psi_obc_north = .false.! prescribe streamfct. do not use
      logical :: prescribe_psi_obc_south = .false. ! an orlanski radiation cond.
      logical :: prescribe_psi_obc_west  = .false.
      logical :: prescribe_psi_obc_east  = .false.

      logical :: enable_obc_adjust_islands = .false.

      ! the following parameters define a frictional sponge for obc south
      ! for mother domain only
      logical :: enable_obc_south_sponge       = .false.
      logical :: enable_obc_south_sponge_harm  = .false.
      logical :: enable_obc_south_sponge_diff  = .false.
      integer :: obc_south_sponge_width   = 45 
      real    :: obc_south_sponge_scale   = 6.
      real    :: obc_south_sponge_fac     = 10.
      ! same for obc_north
      logical :: enable_obc_north_sponge       = .false.
      logical :: enable_obc_north_sponge_harm  = .false.
      logical :: enable_obc_north_sponge_diff  = .false.
      integer :: obc_north_sponge_width   = 45 
      real    :: obc_north_sponge_scale   = 6.
      real    :: obc_north_sponge_fac     = 10.
      ! same for obc_east
      logical :: enable_obc_east_sponge       = .false.
      logical :: enable_obc_east_sponge_harm  = .false.
      logical :: enable_obc_east_sponge_diff  = .false.
      integer :: obc_east_sponge_width   = 45 
      real    :: obc_east_sponge_scale   = 6.
      real    :: obc_east_sponge_fac     = 10.
      ! same for obc_west
      logical :: enable_obc_west_sponge       = .false.
      logical :: enable_obc_west_sponge_harm  = .false.
      logical :: enable_obc_west_sponge_diff  = .false.
      integer :: obc_west_sponge_width   = 45 
      real    :: obc_west_sponge_scale   = 6.
      real    :: obc_west_sponge_fac     = 10.

!     basic horizontal mixing setup
      logical :: enable_diffusion_harmonic   = .false.! harmonic (z-level) diffusion
      logical :: enable_diffusion_biharmonic = .false.
      logical :: enable_diffusion_isoneutral = .false. ! new scheme by Griffith
      logical :: enable_diffusion_isopycnic  = .false. ! old Redi/Cox scheme
      logical :: enable_friction_harmonic    = .false.! harmonic friction
      logical :: enable_friction_biharmonic  = .false.! biharmonic ( del^4 u )

      logical :: enable_isopycnic_horvar     = .false. ! horizontal varying thickness diff 
                                                        ! in Redi/Cox scheme (for eddy-perm. appl.)
      logical :: enable_no_isopyc_medwater   = .false. ! no isopycnic diffusion in Gibraltar area (only HS_4)


      logical :: enable_friction_cosine_scal  = .false. ! scale all diff/fric.parameters
      logical :: enable_diffusion_cosine_scal = .false. ! with cosine of latitude

!     basic vertical mixing setup
      logical :: enable_implicit_vert_diff   = .false. ! semi-implicit vert diffusion
      logical :: enable_implicit_vert_fric   = .false. ! semi-implicit vert friction
      logical :: enable_const_vert_mixing        = .false.! constant vertical mix.
      logical :: enable_cgh_vert_mixing          = .false.! Cummins et al mix.
      logical :: enable_cgh_vert_momentum_mixing = .false.! also for momentum
      logical :: enable_cgh_impl_convection      = .false.! implicit convection scheme
      logical :: enable_kpp_vert_mixing          = .false.! not yet implemented
      logical :: enable_expl_convection          = .false.! explicit convection scheme
      logical :: enable_impl_convection          = .false.! implicit convection scheme
      logical :: enable_ktmix                    = .false.! Kraus Turner wind stirring
      logical :: enable_tkemix                   = .false.! Gaspar et al 90 TKE mixing 


!     BLUE setup
      logical :: enable_blue         = .false. ! enables the BLUE approach
      logical :: enable_blue_smooth  = .false. ! smoothed version of BLUE
      logical :: enable_blue_const   = .false. ! constant forcing version of BLUE
      logical :: enable_blue_tapered = .false. ! boundary tapering
      integer :: ismooth_blue  = 2             ! index for smoothing running mean 
      integer :: ilook_blue = 1                ! index for boundary tapering
      logical :: enable_blue_mean    = .false. ! enables the mean BLUE approach
      logical :: enable_blue_simple  = .false. ! setup BLUE for simple configuration
                                               ! see also file setup_template.F

!     basic horizontal mixing parameter
      real :: ah=2e7              ! constant harmonic diffusivity (cm^2/s)
      real :: am=1e8              ! constant harmonic viscosity (cm^2/s)
      real :: ahbi=2e19           ! constant biharmonic diffusivity (cm^4/s)
      real :: ambi=2e19           ! constant biharmonic viscosity (cm^4/s)
      real :: slmx    = 1.0/100.0 ! Maximum isoneutral slope for -Dgkw_taper
      real :: del_dm = 4.0/1000.0 ! transition for scaling isoneutral 
                                  ! diffusivities with -Ddm_taper
      real :: s_dm = 1.0/1000.0   ! half width scaling for diffusivity with -Ddm_taper
      real :: ahisop  = 1.e7      ! Redi isoneutral diffusivity (cm2/sec)
      real :: athkdf  = 1.e7      ! Gent and McWilliams diffusivity (cm2/sec)
      real :: athkdf_ep = 1.e7    ! same for module epflux
      real :: athkdf_ep2= 1.e7    ! same for module epflux
      real :: ahsteep = 1.e7      ! Horizontal diffusivity used in regions 
                                    ! of steep slopes (cm2/sec)
      real :: isohv_a1  = 0.5     ! for horizontal varying ahisop/ahthkdf coeff.
      real :: isohv_a2  = 2.0     ! for horizontal varying ahisop/ahthkdf coeff.
     
      real :: aidif_tracer   = 0.0 ! parameter for implicit vertical mixing
      real :: aidif_momentum = 0.0 ! 0.5 is Crank-Nicholsen, they are set automatically

!     constant vertical mixing parameters
      real :: kappa_h=0.3         ! constant vertical diffusivity (cm^2/s)
                                  ! NOTE: always used at open boundaries
      real :: kappa_m=10          ! constant vertical viscosity (cm^2/s)

!     parameter for Cummins, et al (JPO, 20, 1990, 817-830) vertical mixing
      real :: cgh_vdcfac(2) = (/1.e-3,1.e-2/)
      real :: diff_cbt_cut  = 4.0
      real :: diff_cbt_back = 0.1
      real :: wndmix        = 50.0
      real :: visc_cbu_cut  = 10.0
      real :: visc_cbu_back = 1.0

!     Kraus Turner mixed layer model, see vmixc.F
      real    :: effwind=0.8, z_scale=5000.0

!     penetration scale for shortwave radiation
      real ::  rpart_shortwave  = 0.58
      real ::  efold1_shortwave = 35.0e0
      real ::  efold2_shortwave = 23.0e2

!     TKE closure scheme a la Gaspar et al
!     standard parameters were given by A.Oschlies
      real :: C_eps        =0.7   ! coef. of the Kolmogoroff dissipation
      real :: C_kappa      =0.1   ! coeff. for vertical eddy coef.
      real :: C_tke        =30.   ! coef. for the tke vert. diff. coeff.
      real :: tke_min      =4.e-2 ! minimum value of tke (cm2/s2)
      real :: tke_min_surf =1.    ! surface minimum value of tke (cm2/s2)
      real :: tke_surf_fac =3.0   ! coef. of the surface input of tke
      integer :: nmxl      =2     ! = 0/1/2/3 flag for mixing length used
!     NOTE:  also used for tkemix are parameters diff_cbt_cut, diff_cbt_back,
!            visc_cbu_cut, visc_cbu_back from CGH mixing

      real :: visc_cnu,visc_ceu    ! do not touch these
      real :: diff_cet,diff_cnt    ! do not know why

!     parameter for bottom drag, see setvbc.F
      real :: cdbot=0.,tidaloff=0.

!     bottom boundary layer parametrisation
      logical :: enable_bbl            = .false. ! enable BBL
      logical :: enable_bbl_advection  = .false. ! enable adv. transport in BBL
      logical :: enable_bbl_diffusion  = .false. ! enable diffusive transport
      real :: ah_sigma=1.0e8           ! diffusivity in sigma layer (cm^2/s)
      real :: ah_sigma_min=1e1         ! minimal diffusivity (cm^2/s)
      logical :: enable_bbl_only_north = .false. ! enable BBl only north of 45N and south of 45S
                                                  ! i.e. not for the Med Sea overflow

!     advection schemes
      logical :: enable_quicker_advection = .false.  ! quicker advection for tracers
      logical :: enable_upstream_advection= .false.  ! upstream advection for tracers
      logical :: enable_4th_advection     = .false.  ! 4.th ord. advection for tracers
      logical :: enable_fct_advection     = .false.  ! FCT advection for tracers
      logical :: enable_flux_delimiter    = .false.  ! limit flux for passive tracers
!     using the multidim. pos. definite adv. scheme (MPDCD) by Lafore, et al 1998

!     icemask shots off the surface heat and salt fluxes
!     if SST falls below the freezing point
      logical :: enable_icemask =  .false.
!     account for effect of ocean currents on wind stress
!     (should be small)
      logical :: enable_windstress_feedback =  .false.

!     hor. diagnostic fields for convective adjustment
c     lctot = total of number of levels involved in convection
c     lcven = number of levels ventilated (convection to surface)
      integer, allocatable :: lctot(:,:), lcven(:,:)

!     parameter for the poisson solver
      real    :: eps_solver=1.e8 
      integer :: max_itts_solver=500

!     Diagnostics setup
      real    :: snap_int=0. ! snaphost interval, also used for all other diagnostics
      logical :: enable_snapshots   = .false.    ! NetCDF snapshot
      logical :: enable_timeseries  = .false.    ! Timeseries output
      logical :: enable_ts_monitor  = .true.     ! Time step monitor
      logical :: enable_overturning = .false.    ! meridional overturning
      logical :: enable_heat_tr     = .false.    ! meridional tracer transports
      logical :: enable_diag_press     = .false. ! diagnostic surface/bottom pressure
      logical :: enable_diag_float     = .false. ! float diagnostics
      logical :: enable_diag_blue      = .false. ! BLUE diagnostics
      real    :: eps_surf_press = 1.e-4          ! epsilon for surf. press. in solver 

      character (len=80) :: snap_file='snap.cdf'  ! snapshot netcdf filename
      character (len=80) :: over_file      ='vsf.cdf' 
      character (len=80) :: heat_tr_file   ='heat_tr.cdf'  
      character (len=80) :: diag_press_file='diag_press.cdf'  
      character (len=80) :: diag_float_file='float.cdf'  
 
      logical :: enable_annual_averages   = .false.  ! annual averaged quantities
      logical :: enable_seasonal_averages = .false.
      logical :: enable_monthly_averages  = .false.
      logical :: enable_daily_averages    = .false.
      logical :: enable_seasonal_variances= .false.  ! same for second order moments
      logical :: enable_monthly_variances = .false.
      logical :: enable_daily_variances   = .false.

      logical :: enable_show_island_map   = .false.  ! show a map of islands in stdout
      logical :: enable_stability_tests   = .false.  ! CFL, etc. However very slow


!     basinwide integrated area and volume of T-cell and U-cell grid boxes
      real, allocatable :: tcella(:),ucella(:)
      real :: tcellv,ucellv
      integer :: itts_solver ! last number of iteration needed in the poisson solver
        
!     main working arrays
      real, allocatable :: rho(:,:,:)
      real, allocatable :: u(:,:,:,:,:),t(:,:,:,:,:)

!     advection and diffusion
      real, allocatable :: adv_vnt(:,:,:)
      real, allocatable :: adv_vet(:,:,:)
      real, allocatable :: adv_vbt(:,:,:)

      real, allocatable :: adv_vnu(:,:,:)
      real, allocatable :: adv_veu(:,:,:)
      real, allocatable :: adv_vbu(:,:,:)

      real, allocatable :: visc_cbu(:,:,:)
      real, allocatable :: diff_cbt(:,:,:)

!     barotropic mode and forcing for it
      real, allocatable :: zu(:,:,:),psi(:,:,:),ptd(:,:),guess(:,:)

!     topography related arrays
      integer, allocatable :: kmt(:,:),kmu(:,:),map(:,:)
      integer, allocatable :: kmt_big(:,:) ! topography for entire model
      real,    allocatable :: hr(:,:)
      real,    allocatable :: tmask(:,:,:)
      real,    allocatable :: umask(:,:,:)

!     other arrays
      real, allocatable :: source(:,:,:)
      real, allocatable :: am_scale(:,:)
      real, allocatable :: ambi_scale(:,:)
      real, allocatable :: ah_scale(:,:)
      real, allocatable :: ahbi_scale(:,:)
      real, allocatable :: adv_scale(:,:)

!     grid variables
      real,allocatable :: xt(:), yt(:), xu(:), yu(:),
     & dxtdeg(:), dytdeg(:), dxudeg(:), dyudeg(:),
     & dxt(:), dyt(:), dxu(:), dyu(:),
     & dytr(:),dyt2r(:), dyt4r(:), dyur(:),dyu2r(:), dyu4r(:),
     & phi(:),phit(:), cst(:),csu(:),sine(:),
     & cstr(:),csur(:),tng(:),
     & cstdytr(:), cstdyt2r(:),csudyur(:),csudyu2r(:),
     & cst_dytr(:),csu_dyur(:),
     & dus(:),dun(:),
     & dxtr(:),dxt2r(:),dxur(:), dxu2r(:),
     & dxt4r(:),dxu4r(:),
     & dxmetr(:),duw(:),due(:),
     & zt(:),zw(:),dzt(:),dzw(:),
     & c2dzt(:), dztr(:), dzt2r(:),
     & rho0dztr(:), dzwr(:), dzw2r(:),
     & rho0csudyur(:), dztur(:), dztlr(:)

!     grid related variables
      real, allocatable :: rho0csudxur(:,:), csudxu2r(:,:),
     &                     am_csudxtr(:,:)  
      real, allocatable :: cstdxtr(:,:), cstdxt2r(:,:),
     &                     ah_cstdxur(:,:) 

      real, allocatable :: cori(:,:,:)  ! +/-coriolis factor, always 2 dimensional
      real, allocatable :: am3_biha(:),am4_biha(:,:),advmet(:,:)
      real, allocatable :: am3_ha(:),am4_ha(:,:)

      real, allocatable :: amc_north_ha(:),amc_south_ha(:)
      real, allocatable :: ahc_north_ha(:),ahc_south_ha(:)

!     surface  boundary conditions
      real, allocatable :: smf(:,:,:),bmf(:,:,:)
      real, allocatable :: stf(:,:,:),btf(:,:,:)
      real, allocatable :: stf_rest(:,:,:),stf_clim(:,:,:)
      real, allocatable :: icemask(:,:),ustar(:,:),dml(:,:)
      real, allocatable :: qsol(:,:),divpen_shortwave(:)

!     open boundary conditions
      real, allocatable :: c1ps(:),c1pn(:),c1pe(:),c1pw(:)
      real, allocatable :: psi_wall_north(:)
      real, allocatable :: psi_wall_east(:)
      real, allocatable :: psi_wall_west(:)
      real, allocatable :: psi_wall_south(:)
      real, allocatable :: ts_obc_south(:,:,:)
      real, allocatable :: ts_obc_north(:,:,:)
      real, allocatable :: ts_obc_east(:,:,:)
      real, allocatable :: ts_obc_west(:,:,:)

!     TKE stuff
      real, allocatable :: eke(:,:,:,:)

!     partial bottom cells
!
!     htp   = depth (cm) from surface to the bottom of partial T-cells
!     hup   = depth (cm) from surface to the bottom of partial U-cells
!     dh = delta thickness, dhwt = vertical distance between T points 
!     fracdz(k,0) = fractional distance between grid point and cell top
!     fracdz(k,1) = fractional distance between grid point and cell bot
!     ztp : corresponding to zt
!
      real, dimension(:,:)  , allocatable :: htp, hup
      real, dimension(:,:,:), allocatable :: dht, dhu, dhwt
      real, dimension(:,:,:), allocatable :: pc_sink_ha
      real, dimension(:,:,:), allocatable :: pc_sink_biha
      real, dimension(:,:,:), allocatable :: ztp
      real, dimension(:,:)  , allocatable :: fracdz

      end module spflame_module

      subroutine init_spflame_module
c
c-----------------------------------------------------------------------
c     Initialization of the main module
c-----------------------------------------------------------------------
c
      use spflame_module
      implicit none
c
c     the main namelist 
c

      namelist /spflame_basic/ 
c
c                 domain, number of tracers and PE decomposition
c
     &            imt,jmt,km,nt,n_pes_i,n_pes_j,
c
c                 time stepping
c
     &            time_step,runlen

      namelist /spflame_nesting_setup/ 
c
c                 nesting of sub domains
c
     &            mother, childs, nr_childs,
     &            zoom_fac,is_zoom,ie_zoom,js_zoom,je_zoom, zoom_fac_k,
     &            topo_smooth_iterations

      namelist /spflame_general_setup/ 
c
c                 general setup
c
     &            cyclic, read_restart,  enable_freesurf,
     &            enable_beta_plane, beta_plane_j0, enable_f_plane,
     &            enable_rotated_grid,
     &            enable_simple_grid, enable_simple_topo ,
     &            enable_simple_initial_cond  ,
     &            enable_simple_sbc, enable_simple_obc ,
     &            enable_simple_spg , enable_sponge

      namelist /spflame_hor_mixing_setup/ 
c
c                 basic horizontal mixing setup
c
     &            enable_diffusion_harmonic   ,
     &            enable_diffusion_biharmonic ,
     &            enable_diffusion_isoneutral ,
     &            enable_diffusion_isopycnic ,
     &            enable_isopycnic_horvar ,
     &            enable_no_isopyc_medwater ,
     &            enable_friction_harmonic    ,
     &            enable_friction_biharmonic  ,
     &            enable_friction_cosine_scal  ,
     &            enable_diffusion_cosine_scal  

      namelist /spflame_vert_mixing_setup/ 
c
c                 basic vertical mixing setup
c
     &            enable_implicit_vert_diff,
     &            enable_implicit_vert_fric,
     &            enable_const_vert_mixing        ,        
     &            enable_cgh_vert_mixing          ,        
     &            enable_cgh_vert_momentum_mixing ,        
     &            enable_cgh_impl_convection      ,  
     &            enable_kpp_vert_mixing          ,  
     &            enable_expl_convection          ,  
     &            enable_impl_convection          ,  
     &            enable_ktmix,               
     &            enable_tkemix                     

      namelist /spflame_advection_scheme/ 
c
c                 advection schemes
c
     &            enable_quicker_advection,
     &            enable_upstream_advection,
     &            enable_4th_advection,
     &            enable_fct_advection,
     &            enable_flux_delimiter


      namelist /spflame_obc_setup/ 
c
c                 Open boundaries set up
c
     &            enable_obc_north, enable_obc_south,
     &            enable_obc_west,  enable_obc_east,
     &            restore_TS_obc_north, restore_TS_obc_south,
     &            restore_TS_obc_west,  restore_TS_obc_east,
     &            prescribe_psi_obc_north, prescribe_psi_obc_south,
     &            prescribe_psi_obc_west,  prescribe_psi_obc_east,
     &            enable_obc_south_sponge ,  obc_south_sponge_fac   ,
     &            enable_obc_south_sponge_harm ,
     &            enable_obc_south_sponge_diff,
     &            obc_south_sponge_width,obc_south_sponge_scale ,
     &            enable_obc_north_sponge ,  obc_north_sponge_fac  ,
     &            enable_obc_north_sponge_harm ,
     &            enable_obc_north_sponge_diff ,
     &            obc_north_sponge_width , obc_north_sponge_scale ,
     &            enable_obc_east_sponge ,  obc_east_sponge_fac  ,
     &            enable_obc_east_sponge_harm ,
     &            enable_obc_east_sponge_diff ,
     &            obc_east_sponge_width , obc_east_sponge_scale ,
     &            enable_obc_west_sponge ,  obc_west_sponge_fac  ,
     &            enable_obc_west_sponge_harm ,
     &            enable_obc_west_sponge_diff ,
     &            obc_west_sponge_width , obc_west_sponge_scale ,
     &            enable_obc_adjust_islands

      namelist /spflame_sbc_setup/ 
c
c                   surface boundary setup
c
     &            enable_salt_flux_sbc,     ! fixed salt flux
     &            enable_shortwave_sbc,     ! shortwave penetration
     &            enable_icemask,           !  icemask, zero order ice model
     &            enable_windstress_feedback,! effect of ocean currents on wind stress
     &            efold1_shortwave,efold2_shortwave,rpart_shortwave
c                                   Shortwave radiation penetration scale

      namelist /spflame_blue_setup/ 
c
c                 Blue stuff
c
     &            enable_blue,enable_blue_simple,
     &            enable_blue_smooth,enable_blue_const,ismooth_blue,
     &            enable_blue_tapered, ilook_blue, enable_blue_mean

      namelist /spflame_bbl_setup/
c
c                 bottom boundary layer parametrisation
c
     &            enable_bbl,  
     &            enable_bbl_advection, enable_bbl_diffusion,
     &            ah_sigma, ah_sigma_min,enable_bbl_only_north


      namelist /spflame_mixing_parameters/ 
c
c                 horizontal mixing parameters
c
     &            ambi,ahbi,ah,am, isohv_a1,isohv_a2,
     &            ahisop,athkdf,slmx,ahsteep,del_dm,s_dm,
c
c                 vertical mixing parameters
c
     &            kappa_m,kappa_h,
     &            cgh_vdcfac, diff_cbt_cut, diff_cbt_back, 
     &            wndmix,     visc_cbu_cut, visc_cbu_back,
c
c                 bottom friction
c
     &            tidaloff,cdbot,
c
c                 Kraus Turner model
c
     &            effwind, z_scale,
c
c                 TKE model
c
     &            C_eps,C_kappa,C_tke,
     &            tke_min,tke_min_surf,tke_surf_fac,nmxl

      namelist /spflame_solver/ 
c
c                 poisson solver
c
     &            eps_solver,
     &            max_itts_solver


      namelist /spflame_diagnostics/ 
c
c                 Diagnostics
c
     &            snap_int,
     &            enable_snapshots,
     &            enable_timeseries,
     &            enable_ts_monitor,
     &            enable_overturning,over_file,
     &            enable_heat_tr,heat_tr_file,
     &            enable_diag_press,diag_press_file,eps_surf_press,
     &            enable_annual_averages,
     &            enable_seasonal_averages,
     &            enable_monthly_averages,
     &            enable_daily_averages,
     &            enable_seasonal_variances,
     &            enable_monthly_variances,
     &            enable_daily_variances,
     &            enable_show_island_map,
     &            enable_stability_tests,
     &            enable_diag_float,diag_float_file,
     &            enable_diag_blue


      integer k,ierr,n
      character (len=80) :: name

      call mpi_comm_rank(my_comm,my_pe,ierr)
      call mpi_comm_size(my_comm,n_pes,ierr)

      if (my_pe == 0) then
       print*,''
       print*,' ------------------------------------------'
       print*,''
       print'("      SPFLAME version ",f5.3)',version
       print*,''
       print*,' ------------------------------------------'
       print*,''
       print*,'       reading namelist input'
      endif

c     read in namelist for this specific sub_domain 

      write(name, '("namelist.spflame_",i3)') sub_domain
      do k=1,len_trim(name); if (name(k:k)==' ')name(k:k)='0'
      enddo
      call get_free_iounit(n)
      open(n,file=name,form='formatted',status='old')
      read(n,nml = spflame_basic) 
      read(n,nml = spflame_nesting_setup) 
      read(n,nml = spflame_general_setup) 
      read(n,nml = spflame_hor_mixing_setup) 
      read(n,nml = spflame_vert_mixing_setup) 
      read(n,nml = spflame_advection_scheme) 
      read(n,nml = spflame_obc_setup) 
      read(n,nml = spflame_sbc_setup) 
      read(n,nml = spflame_blue_setup) 
      read(n,nml = spflame_bbl_setup) 
      read(n,nml = spflame_mixing_parameters) 
      read(n,nml = spflame_solver) 
      read(n,nml = spflame_diagnostics) 
      close(n)
c
c     do the MPP setup
c
      if (mother<0) then
       if (my_pe==0) then
        print*,''
        print*,' sub_domain #',sub_domain
        call getunit(n,'document_000.dta','fsr')
        write(n,'("SPFLAME Version ",f5.3)') version
        write(n,*) 'spflame namelist input:'
        write(n,nml = spflame_basic) 
        write(n,nml = spflame_hor_mixing_setup) 
        write(n,nml = spflame_vert_mixing_setup) 
        write(n,nml = spflame_blue_setup) 
        write(n,nml = spflame_bbl_setup) 
        write(n,nml = spflame_mixing_parameters) 
        write(n,nml = spflame_advection_scheme) 
        write(n,nml = spflame_solver) 
        write(n,nml = spflame_nesting_setup) 
        write(n,nml = spflame_general_setup) 
        write(n,nml = spflame_obc_setup) 
        write(n,nml = spflame_sbc_setup) 
        write(n,nml = spflame_diagnostics) 
        close(n)
        print*,''
        print*,' namelist input written to file document_000.dta '
        print*,''
       endif
       km_mother=km
      else
       if (my_pe==0) then
        print*,''
        print*,' sub_domain #',sub_domain
        write(name, '("document_",i3,".dta")') sub_domain
        do k=1,len_trim(name); if (name(k:k)==' ')name(k:k)='0'
        enddo
        call getunit(n,name,'fsr')
        write(n,'("SPFLAME Version ",f5.3)') version
        write(n,*) 'spflame namelist input:'
        write(n,nml = spflame_basic) 
        write(n,nml = spflame_hor_mixing_setup) 
        write(n,nml = spflame_vert_mixing_setup) 
        write(n,nml = spflame_blue_setup) 
        write(n,nml = spflame_bbl_setup) 
        write(n,nml = spflame_mixing_parameters) 
        write(n,nml = spflame_advection_scheme) 
        write(n,nml = spflame_solver) 
        write(n,nml = spflame_nesting_setup) 
        write(n,nml = spflame_general_setup) 
        write(n,nml = spflame_obc_setup) 
        write(n,nml = spflame_sbc_setup) 
        write(n,nml = spflame_diagnostics) 
        close(n)
        print*,''
        print*,' namelist input written to file ',name(1:len_trim(name))
        print*,''
       endif
c
c      checks for zoom_fac
c      zoom_fac must be odd
c
       if ( mod(zoom_fac,2) /= 1 .or. zoom_fac < 3) 
     &     call halt_stop(' zoom_fac must be odd and greater than 3')
c       if ( mod(zoom_fac,2) /= 1 .or. zoom_fac < 1) 
c     &     call halt_stop(' zoom_fac must be odd and greater than 0')
       if ( mod(zoom_fac_k,2) /= 1 .or. zoom_fac_k < 1) 
     &     call halt_stop(' zoom_fac_k must be odd and greater than 0')
c
c      checks for is_zoom and ie_zoom
       if ( is_zoom < 1 .or. js_zoom < 1 ) 
     &   call halt_stop(' sub domain chosen by is_zoom,'//
     &      ' etc must be inside mother domain ')
c
       imt=(ie_zoom-is_zoom)*zoom_fac+1
       jmt=(je_zoom-js_zoom)*zoom_fac+1
       km_mother=km
       km = km*zoom_fac_k
      endif ! mother <0

      if (n_pes /= n_pes_i*n_pes_j) then
       if (my_pe == 0) then
        print*,' Number of PEs (=',n_pes,') do not '
        print*,' match the chosen domain decomposition '
        print*,' n_pes_i=',n_pes_i,' n_pes_j=',n_pes_j
       endif
       call halt_stop('error ')
      endif
c
c     domain decomposition for each PE
c
      i_blk = (imt-1)/n_pes_i + 1    ! i-extent of each block
      j_blk = (jmt-1)/n_pes_j + 1    ! j-extent of each block
      my_blk_i = mod(my_pe,n_pes_i)+1! number of PE in i-dir.
      my_blk_j = (my_pe)/n_pes_i + 1 ! number of PE in j-dir.
      is_pe = (my_blk_i-1)*i_blk + 1 ! start index in i-dir of this PE
      ie_pe = min(my_blk_i*i_blk,imt)
      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     it happens that weird values are computed be the above algorythm
c
      if (my_blk_j==n_pes_j .and. js_pe>=je_pe-2) then
       print*,' ERROR:'
       print*,' domain decompositon impossible in j-direction'
       print*,' choose other number of PEs in j-direction'
       call halt_stop(' in init_spflame_module')
      endif

      if (my_blk_i==n_pes_i .and. is_pe>=ie_pe-2) then
       print*,' ERROR:'
       print*,' domain decompositon impossible in i-direction'
       print*,' choose other number of PEs in i-direction'
       call halt_stop(' in init_spflame_module')
      endif
c
c     print out the PE decomposition, let all PEs talk
c
      if (my_pe==0) then
       print*,''
       print*,' SPFLAME sub domain setup: '
       print*,''
       print*,' imt=',imt,' jmt=',jmt,' km=',km
       print*,' n_pes_i=',n_pes_i,' n_pes_j=',n_pes_j
       print*,' i_blk=',i_blk,' j_blk=',j_blk
      endif
      call barrier
      do n=0,n_pes-1
       if (my_pe==n) then
        print*,''
        print*,' sub domain PE #',n
        print*,' my_blk_i=',my_blk_i,' my_blk_j=',my_blk_j
        print*,' is_pe=',is_pe,' ie_pe=',ie_pe
        print*,' js_pe=',js_pe,' je_pe=',je_pe
        print*,''
       endif
       call barrier
      enddo

      call allocate_spflame_module

      end subroutine init_spflame_module


      subroutine  allocate_spflame_module
c
c-----------------------------------------------------------------------
c     Allocation of all arrays in the main module
c     dimensions are restricted to the PE domain
c-----------------------------------------------------------------------
c
      use spflame_module
      implicit none
      integer  :: is,ie,js,je,k
      real, parameter    :: rini = 0.0
      integer, parameter :: iini = 0

      is=is_pe; ie=ie_pe; js=js_pe; je=je_pe

      allocate( u(is-2:ie+2,km,js-2:je+2,2,0:2) ); u=rini
      allocate( t(is-2:ie+2,km,js-2:je+2,nt,0:2) ); t=rini
      allocate( rho(is-1:ie+1,km,js-1:je+1)); rho=rini
      allocate( adv_vnt(is-1:ie+1,  km,js-1:je+1) );adv_vnt=rini
      allocate( adv_vet(is-1:ie+1,  km,js-1:je+1) );adv_vet=rini
      allocate( adv_vbt(is-1:ie+1,0:km,js-1:je+1) );adv_vbt=rini
      allocate( adv_vnu(is-1:ie+1,  km,js-1:je+1) );adv_vnu=rini
      allocate( adv_veu(is-1:ie+1,  km,js-1:je+1) );adv_veu=rini
      allocate( adv_vbu(is-1:ie+1,0:km,js-1:je+1) );adv_vbu=rini
      allocate( visc_cbu(is:ie,km,js:je) ); visc_cbu=rini
      allocate( diff_cbt(is:ie,km,js:je) ); diff_cbt=rini
      allocate( smf(is:ie,js:je,2) ); smf=rini
      allocate( bmf(is:ie,js:je,2) ); bmf=rini
      allocate( stf(is:ie,js:je,nt), btf(is:ie,js:je,nt) )
      stf=rini; btf=rini
      allocate( stf_rest(is:ie,js:je,nt), stf_clim(is:ie,js:je,nt) )
      stf_rest=rini; stf_clim=rini
      if (enable_icemask) then
       allocate( icemask(is:ie,js:je) ); icemask=rini
      endif
      if (enable_ktmix) then
       allocate( dml(is:ie,js:je) );   dml=rini
      endif
      if (enable_tkemix .or. enable_ktmix) then
       allocate( ustar(is:ie,js:je) ); ustar=rini
      endif
      if (enable_shortwave_sbc) then
       allocate( qsol(is:ie,js:je) ); qsol=rini
       allocate( divpen_shortwave(km) ); divpen_shortwave=rini
      endif
      if (enable_tkemix) then
       allocate( eke(is:ie,km,js:je,0:2) ); eke= rini
      endif
      allocate( zu(is-1:ie+1,js-1:je+1,2) ); zu=rini
      allocate( psi(is-1:ie+1,js-1:je+1,2) ); psi=rini
      allocate( ptd(is-1:ie+1,js-1:je+1) ); ptd=0.
      allocate( guess(is-1:ie+1,js-1:je+1) ); guess=0.
      allocate( kmt(is:ie,js:je) ); kmt=iini
      allocate( kmt_big(imt,jmt) ); kmt_big=iini
      allocate( kmu(is:ie,js:je) ); kmu=iini
      allocate( map(is:ie,js:je) ); map=iini
      allocate( hr(is-1:ie+1,js-1:je+1) ); hr=rini
      if (mother<0) then
       allocate( umask(is-2:ie+2,km,js-2:je+2) )  
       allocate( tmask(is-2:ie+2,km,js-2:je+2) )
      else
c      these boundaries are needed for smoothing the forcing from mother domain
       k=max(2,zoom_fac/2)
       allocate( umask(is-k:ie+k,km,js-k:je+k) )  
       allocate( tmask(is-k:ie+k,km,js-k:je+k) )
      endif
      umask=rini; tmask=rini
      allocate( source(is:ie,km,js:je) )
!     scaling of ah,ahbi am and ambi, eg. near open boundaries
      allocate( am_scale(is-2:ie+2,js-2:je+2) )
      am_scale=1.
      allocate( ambi_scale(is-2:ie+2,js-2:je+2) )
      ambi_scale=1.
      allocate( ah_scale(is-2:ie+2,js-2:je+2) )
      ah_scale=1.
      allocate( ahbi_scale(is-2:ie+2,js-2:je+2) )
      ahbi_scale=1.
!     grid variables are dimensioned over the whole model domain
      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),dyt4r(jmt) )
      allocate( dyur(jmt),dyu2r(jmt), dyu4r(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( dxt4r(imt),dxu4r(imt))
      allocate( dxmetr(imt),duw(imt),due(imt))
      allocate( zt(km),zw(km),dzt(km),dzw(0:km) )
      allocate( c2dzt(km), dztr(km), dzt2r(km) )
      allocate( rho0dztr(km), dzwr(0:km), dzw2r(0:km) )
      allocate( rho0csudyur(jmt) )
      allocate( dztur(km), dztlr(km) )

      allocate(rho0csudxur(is-1:ie+1,js-1:je+1) )
      allocate(csudxu2r   (is-1:ie+1,js-1:je+1) )
      allocate(am_csudxtr (is-1:ie+1,js-1:je+1)   )
      allocate(cstdxtr    (is-1:ie+1,js-1:je+1) )
      allocate(cstdxt2r   (is-1:ie+1,js-1:je+1) )
      allocate(ah_cstdxur (is-1:ie+1,js-1:je+1)  )
      allocate( cori(is:ie,js:je,2) )
      allocate( am3_biha(jmt),am4_biha(jmt,2),advmet(jmt,2) )
      allocate( am3_ha(jmt)  ,am4_ha(jmt,2) )
      allocate( amc_north_ha(jmt), amc_south_ha(jmt) )
      allocate( ahc_north_ha(jmt), ahc_south_ha(jmt) )
!     needed for diagnostics
      allocate( tcella(km), ucella(km) )
      if (enable_expl_convection) then
       allocate( lcven(is:ie,js:je), lctot(is:ie,js:je) )
      endif
!     open boundary conditions
      allocate( c1ps(is:ie),c1pn(is:ie),c1pe(js:je),c1pw(js:je) )
      allocate( psi_wall_north(imt) ); psi_wall_north=0.
      allocate( psi_wall_east(jmt) );  psi_wall_east=0.
      allocate( psi_wall_west(jmt) );  psi_wall_west=0.
      allocate( psi_wall_south(imt) ); psi_wall_south=0.
      if (restore_TS_obc_south) then
       allocate( ts_obc_south(is:ie,km,nt) ); ts_obc_south =0.
      endif
      if (restore_TS_obc_north) then
       allocate( ts_obc_north(is:ie,km,nt) ); ts_obc_north =0.
      endif
      if (restore_TS_obc_east) then
       allocate( ts_obc_east(js:je,km,nt) ); ts_obc_east =0.
      endif
      if (restore_TS_obc_west) then
       allocate( ts_obc_west(js:je,km,nt) ); ts_obc_west =0.
      endif
#ifdef partial_cell
      allocate( htp(imt,jmt), hup(imt,jmt) )
      htp=rini;hup=rini
      allocate( dht(is_pe-2:ie_pe+2,km,js_pe-2:je_pe+2) )
      allocate( dhu(is_pe-2:ie_pe+2,km,js_pe-2:je_pe+2) )
      allocate( dhwt(is_pe-2:ie_pe+2,0:km,js_pe-2:je_pe+2) )
      allocate( ztp(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) )
      allocate( pc_sink_biha(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) )
      allocate( pc_sink_ha(is_pe:ie_pe,km,js_pe:je_pe) )
      allocate( fracdz(km,0:1) )
      dht=rini;dhu=rini;dhwt=rini;ztp=rini;pc_sink_biha=rini
      pc_sink_ha=rini;fracdz=rini
#endif
      end subroutine allocate_spflame_module

