 MODULE comvars3d
   !!======================================================================
   !!                   ***  MODULE comvars3d  ***
   !! 
   !! ** Description :
   !! 
   !! ** Reference :
   !! 
   !! ** History :
   !!
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters

   !! * Shared or public module variables
   
   !-------------------------------------------------------------------
   ! l_equisig : logique positionne a .true. pour repartir uniformement les
   !             niveaux sigma sur la verticale.
   ! sig       : tableau contenant la valeur de sigma sur les differents
   !             niveau de calcul des points (uz, vz, sal, temp...). Cette
   !             valeur est interpretee comme un pourcentage de la colonne
   !             d eau sous le niveau considere (valeur comprise entre 0 et 1).
   ! dsigu     : epaisseur sigma de la couche encadrant le niveau de calcul uz.
   ! dsigw     : epaisseur sigma de la couche encadrant le niveau de calcul wz
   !             (ou nz, kz, ect).
   LOGICAL                                 :: l_equisig
   REAL(kind=rsh),DIMENSION(0:kmax+1)      :: sig
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:) :: dsigu,dsigw

   !-------------------------------------------------------------------
   ! l_stability : logique positionne a .true. pour melanger les colonnes d eau
   !             trop instables.
   ! turb_nbeq : nombre d equations dans le modele de fermeture turbulente (0,1 ou 2)
   ! turb_0eq_option : choix de l option pour le modele a zero equation
   !              1- Kz et Nz constants
   !              2- Csanady
   !              3- Quetin
   !              4- Pacanowski et Philander
   ! kzinit    : valeur initiale du coefficient vertical de melange en m2/s.
   ! nzinit    : valeur initiale du coefficient vertical de viscosite en m2/s.
   ! kz        : tableau du coefficient turbulent vertical de melange en m2/s
   !             estime au cours du calcul en 3d.
   ! kz        : tableau du coefficient turbulent vertical de viscosite en m2/s
   !             estime au cours du calcul en 3d.
   ! bgdiff    : background diffusivity : valeur en deca de laquelle ne descendent
   !             ni les kz ni les nz. en m2/s
   LOGICAL        :: l_stability
   INTEGER        :: turb_nbeq,turb_0eq_option,turb_2eq_option,OptionSF
   REAL(kind=rsh) :: kzinit,nzinit,bgdiff
   REAL(kind=rsh),DIMENSION(:,:,:),ALLOCATABLE        :: kz,nz

   !-------------------------------------------------------------------
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: wz
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:)     :: sigw
#ifdef  key_diag_physwz
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: wz_phys
#endif
   !-------------------------------------------------------------------
   LOGICAL        :: l_stateeq_lin
   REAL(kind=rsh) :: saliref,tetaref,rhoref,saliref_lin,tetaref_lin
#ifdef key_profil_back
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: rho2,boz2,sali_back,temp_back
#endif
#if defined key_sedim_MUSTANG || defined key_obstructions
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:) :: roswat_bot
#endif

   !-------------------------------------------------------------------
   ! seuil pour choix schema centre ou upstream sur la verticale (key_tssub_wquickupwind uniquement)
   REAL(kind=rsh) :: qmaxz

   !-------------------------------------------------------------------
   ! Variables used for turbulence advection
   ! uzsw       : zonal component of total (barotropic+baroclinic) velocity (in m/s)
   !              at vertical half-grid cells (sigw)
   ! vzsw       : meridional component of total (barotropic+baroclinic) velocity (in m/s)
   !              at vertical half-grid cells (sigw)
   ! wzsw       : vertical component of velocity at vertical half-grid cells
   ! fxsw       : fxeau re-sampled at vertical half-grid cells (sigw)
   ! fysw       : fyeau re-sampled at vertical half-grid cells (sigw)
   ! phisw      : phieau re-sampled at vertical half-grid cells (sigw)
   ! fxect, fyect
   ! fxpsi, fypsi
#ifdef key_turbadv
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: uzsw,vzsw,wzsw,phisw
   REAL(kind=rlg),ALLOCATABLE,DIMENSION(:,:,:) :: fxsw,fysw,fxect,fyect,fzect,fxpsi,fypsi,fzpsi
#endif

#ifdef key_adapt_imp_vert_adv
   ! Adaptive implicit vertical advection
   REAL(kind=rsh),PARAMETER :: theta_wz=0.55_rsh
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: f_alpha_T  ! w dT/dz * 1 / f(alpha)
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: wmo_exp_x  ! w du/dz * 1 / f(alpha)
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: wmo_exp_y  ! w dv/dz * 1 / f(alpha)
#ifdef key_adapt_imp_vert_adv_suivi
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: adapt_imp_uz , adapt_imp_vz , adapt_imp_ts 
#endif
#endif

 CONTAINS

   !!==============================================================================

  SUBROUTINE alloc_comvars3d

   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE alloc_comvars3d  ***
   !&E  
   !&E--------------------------------------------------------------------------
   !! * Modules used
   IMPLICIT NONE

   !!----------------------------------------------------------------------
   !! * Executable part

   ALLOCATE(nz(0:kmax,limin  :limaxp1,ljmin  :ljmaxp1) )
   !ALLOCATE(kz(kmax,limin:limax  ,ljmin:ljmax  ) )
   ! larger extension because of traj_3d routine
   ! might be not sufficient because of the medium location
   ALLOCATE(kz(0:kmax,liminm1:limaxp1,ljminm1:ljmaxp1) )
   ALLOCATE(sigw(0:kmax+1))
   ALLOCATE(dsigu(kmax),dsigw(kmax))

#ifdef key_turbadv
   ALLOCATE(uzsw(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2))
   ALLOCATE(vzsw(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2))
   ALLOCATE(wzsw(1:kmax,liminm1:limaxp1,ljminm1:ljmaxp1))

   ALLOCATE(fxsw(1:kmax,liminm1:limax,ljmin:ljmax))
   ALLOCATE(fysw(1:kmax,limin:limax,ljminm1:ljmax))
   ALLOCATE(phisw(1:kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(fxect(1:kmax,liminm1:limax,ljmin:ljmax))
   ALLOCATE(fyect(1:kmax,limin:limax,ljminm1:ljmax))
   ALLOCATE(fzect(1:kmax,limin:limax,ljminm1:ljmax))
   ALLOCATE(fxpsi(1:kmax,liminm1:limax,ljmin:ljmax))
   ALLOCATE(fypsi(1:kmax,limin:limax,ljminm1:ljmax))
   ALLOCATE(fzpsi(1:kmax,limin:limax,ljminm1:ljmax))
   uzsw(:,:,:) = 0.0_rsh
   vzsw(:,:,:) = 0.0_rsh
   wzsw(:,:,:) = 0.0_rsh
   fxsw(:,:,:) = 0.0_rlg
   fysw(:,:,:) = 0.0_rlg
   phisw(:,:,:) = 0.0_rsh
   fxect(:,:,:) = 0.0_rlg
   fyect(:,:,:) = 0.0_rlg
   fzect(:,:,:) = 0.0_rlg
   fxpsi(:,:,:) = 0.0_rlg
   fypsi(:,:,:) = 0.0_rlg
   fzpsi(:,:,:) = 0.0_rlg
#endif

#ifdef key_adapt_imp_vert_adv
   ALLOCATE( wmo_exp_x(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2) )
   ALLOCATE( wmo_exp_y(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2) )
   ALLOCATE( f_alpha_T(1:kmax,liminm1:limaxp1,ljminm1:ljmaxp1) )
   wmo_exp_x = 1.0_rsh
   wmo_exp_y = 1.0_rsh
   f_alpha_T = 1.0_rsh
#ifdef key_adapt_imp_vert_adv_suivi
   ALLOCATE( adapt_imp_uz(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2) )
   ALLOCATE( adapt_imp_vz(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2) )
   ALLOCATE( adapt_imp_ts(1:kmax,liminm1:limaxp1,ljminm1:ljmaxp1) )
   adapt_imp_uz = 0.
   adapt_imp_vz = 0.
   adapt_imp_ts = 0.
#endif
#endif
#if defined key_sedim_MUSTANG || defined key_obstructions
   ALLOCATE( roswat_bot(limin:limax,ljmin:ljmax) )
   roswat_bot(:,:)=0.0_rsh
#endif

  END SUBROUTINE alloc_comvars3d

   !!==============================================================================

  SUBROUTINE dealloc_comvars3d

   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE dealloc_comvars3d  ***
   !&E  
   !&E--------------------------------------------------------------------------
   !! * Modules used
   IMPLICIT NONE

   !!----------------------------------------------------------------------
   !! * Executable part

   DEALLOCATE(kz,nz)
   DEALLOCATE(sigw,dsigu,dsigw)

#ifdef key_turbadv
   DEALLOCATE(uzsw,vzsw,wzsw,fxsw,fysw,phisw)
   DEALLOCATE(fxect,fyect,fzect,fxpsi,fypsi,fzpsi)
#endif

#ifdef key_adapt_imp_vert_adv
   DEALLOCATE(wmo_exp_x,wmo_exp_y,f_alpha_T)
#ifdef key_adapt_imp_vert_adv_suivi
   DEALLOCATE(adapt_imp_uz,adapt_imp_vz,adapt_imp_ts)
#endif
#endif

  END SUBROUTINE dealloc_comvars3d

   !!==============================================================================


 END MODULE comvars3d
