  SUBROUTINE ts_stateeq_linear(sal,temp,boz,pn2)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE ts_stateeq_linear  ***
   !&E
   !&E ** Purpose : buoyancy estimation from a linear state equation
   !&E
   !&E ** Description : This option is used if l_stateeq_lin=.true.
   !&E                  in the namelist specified for the rank.
   !&E                  
   !&E
   !&E ** Called by : main, stepu, stepv, zoom
   !&E
   !&E ** Used ij-arrays : h0,temp,sal
   !&E 
   !&E ** Modified variables : boz,pn2
   !&E
   !&E ** Reference : 
   !&E
   !&E ** History :
   !&E       !  2004-10-06
   !&E       !  2008-07  (T. Duhaut) different law for baroclinic vortex key_Vortex_Ktest
   !&E       !  2009-04  (A. Mano) MPI2D
   !&E       !  2014-07  (P. Garreau) : Estimate of the brunt-vaisala frequency (pn2)
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d,ONLY : h0,hm
   USE comvars3d,ONLY : saliref_lin,tetaref_lin,rhoref,dsigw
   USE comvarp2d,ONLY : ssh
#ifdef key_profil_back
   USE comvars3d,ONLY : sali_back,temp_back,boz2
#endif
#if defined key_sedim_MUSTANG && defined key_roswat_depend_partcleconc
   USE comsubstance, ONLY : nvpc,cv_wat
   USE comMUSTANG,    ONLY : ros
#endif
#if defined key_siggen || defined key_gencoord
   USE toolgeom, ONLY : f_hzw
#endif
#if defined key_sedim_MUSTANG || defined key_obstructions
   USE comvars3d,  ONLY : roswat_bot
#endif
   USE_MPI toolmpi, ONLY : ex_i_rsh,ex_j_rsh

   IMPLICIT NONE

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in)  :: sal,temp
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(out) :: boz
   REAL(KIND=rsh),DIMENSION(kmax,limin  :limax  ,ljmin  :ljmax  ),INTENT(out) :: pn2

   !! * Local declarations
   INTEGER           :: i,j,k
   REAL(kind=rsh)    :: alpha,beta
   REAL(kind=rsh) :: hz_w
#if defined key_sedim_MUSTANG && defined key_roswat_depend_partcleconc
   INTEGER                         :: iv
   REAL(KIND=rsh),DIMENSION(kmax)  :: sumpart,rostp
#ifdef key_profil_back
   REAL(KIND=rsh),DIMENSION(kmax)  :: sumpart2,rostp2
#endif
#endif


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

!$OMP parallel PRIVATE(alpha,beta,hz_w,k,i,j)
   alpha=0.00016_rsh
   beta=0.0008_rsh

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp1,ljmax
     DO i=liminp1,limax

#if defined key_casbcvortex || defined key_casinnershelf
       IF (h0(i,j)+ssh(i,j) > hm) &
         boz(:,i,j)= -grav*(1030.0_rsh-0.28_rsh*temp(:,i,j)-1024.40_rsh)/1024.40_rsh
#elif key_caslock
       IF (h0(i,j)+ssh(i,j) > hm) &
         boz(:,i,j)= -grav*(1030.0_rsh-0.28_rsh*temp(:,i,j)-rhoref)/rhoref
#elif key_equeta_temp
!FVDM         boz(:,i,j)= temp(:,i,j)
         boz(:,i,j) = -temp(:,i,j)*grav/rhoref
#elif key_tssub_casadvlense
         boz(:,i,j)= boz(:,i,j)
#else
       IF (h0(i,j)+ssh(i,j) > hm) THEN
#if defined key_sedim_MUSTANG && defined key_roswat_depend_partcleconc
         rostp(:)=(1.0_rsh+beta*(sal(:,i,j)-saliref_lin)-alpha*(temp(:,i,j)-tetaref_lin))*rhoref
#ifdef key_profil_back
         rostp2(:)=(1.0_rsh+beta*(sal(:,i,j)-sali_back(:,i,j))-alpha*(temp(:,i,j)-temp_back(:,i,j))*rhoref
#endif
         ! water density function of particles concentration (Eq 7, p235, Le Hir et al, 2001)
         ! "Application of the continuous modeling concept to simulate high-concentration 
         ! suspended sediment in a macrotidal estuary". IN Coastal and Estuarine Fine Sediment Processes
         ! Proceedings in Marine Science, 2001. Ed McAnally & Mehta. Elsevier
         sumpart(:)=0.0_rsh
#ifdef key_profil_back
         sumpart2(:)=0.0_rsh
#endif
         DO iv=1,nvpc
          sumpart(:)=sumpart(:)+cv_wat(iv,:,i,j)*(1.0_rsh-rostp(:)/ros(iv))
#ifdef key_profil_back
          sumpart2(:)=sumpart2(:)+cv_wat(iv,:,i,j)*(1.0_rsh-rostp2(:)/ros(iv))
#endif
         ENDDO
         rostp(:)=rostp(:)+sumpart(:)
#ifdef key_profil_back
         rostp2(:)=rostp2(:)+sumpart2(:)
#endif
         boz(:,i,j)=grav*(rhoref-rostp(:))/rhoref
         roswat_bot(i,j)=REAL(rostp(1),rsh) 
#ifdef key_profil_back
         boz2(:,i,j)=grav*(rhoref-rostp2(:))/rhoref
#endif

#else 
         boz(:,i,j)= grav*(-beta*(sal(:,i,j)-saliref_lin)+alpha*(temp(:,i,j)-tetaref_lin))
#ifdef key_profil_back
         boz2(:,i,j)= grav*(-beta*(sal(:,i,j)-sali_back(:,i,j))+alpha*(temp(:,i,j)-temp_back(:,i,j)))
#endif
#if defined key_sedim_MUSTANG || defined key_obstructions
         roswat_bot(i,j)= (1.0_rsh+beta*(sal(1,i,j)-saliref_lin)-alpha*(temp(1,i,j)-tetaref_lin))*rhoref
#endif         
#endif
       ELSE
#if defined key_sedim_MUSTANG || defined key_obstructions
         roswat_bot(i,j)= (1.0_rsh+beta*(sal(1,i,j)-saliref_lin)-alpha*(temp(1,i,j)-tetaref_lin))*rhoref
#endif         
       ENDIF
#endif

   ! estimation of the Brunt-vaissala frequency (N2)
       IF (h0(i,j)+ssh(i,j) > hm) THEN
       !IF (h0(i,j)+ssh(i,j) > 0.0_rsh) THEN
         DO k=1,kmax-1
#if defined key_siggen || defined key_gencoord
           hz_w=f_hzw(h0(i,j),ssh(i,j),k,i,j)
           pn2(k,i,j)=(boz(k+1,i,j)-boz(k,i,j))/(hz_w*dsigw(k))
#else
           pn2(k,i,j) = (boz(k+1,i,j)-boz(k,i,j))/((h0(i,j)+ssh(i,j))*dsigw(k))
#endif
         END DO      ! fin de la boucle en k
       ELSE
         pn2(:,i,j) = 0.0_rsh
         boz(:,i,j) = 0.0_rsh
#ifdef key_profil_back
         boz2(:,i,j) = 0.0_rsh
#endif
#if defined key_sedim_MUSTANG || defined key_obstructions
         roswat_bot(i,j)= (1.0_rsh+beta*(sal(1,i,j)-saliref_lin)-alpha*(temp(1,i,j)-tetaref_lin))*rhoref
#endif         
       END IF       ! fin du test sur le masque terre 

     END DO
   END DO
!$OMP END DO
!$OMP END PARALLEL

#if ! defined key_profil_back
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,boz(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,boz(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#else
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,boz2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,boz2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif

   PRINT_DBG*, 'exit TS_STATEEQ_LINEAR routine'

  END SUBROUTINE ts_stateeq_linear
