  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---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,ONLY : h0,hm
   USE comvars3d,ONLY : saliref_lin,tetaref_lin,rhoref,dsigw
   USE comvarp2d,ONLY : ssh




   USE comvars3d,  ONLY : roswat_bot

   USE 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



!!----------------------------------------------------------------------
!! * 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 (h0(i,j)+ssh(i,j) > hm) THEN

         boz(:,i,j)= grav*(-beta*(sal(:,i,j)-saliref_lin)+alpha*(temp(:,i,j)-tetaref_lin))


         roswat_bot(i,j)= (1.0_rsh+beta*(sal(1,i,j)-saliref_lin)-alpha*(temp(1,i,j)-tetaref_lin))*rhoref


       ELSE

         roswat_bot(i,j)= (1.0_rsh+beta*(sal(1,i,j)-saliref_lin)-alpha*(temp(1,i,j)-tetaref_lin))*rhoref

       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

           pn2(k,i,j) = (boz(k+1,i,j)-boz(k,i,j))/((h0(i,j)+ssh(i,j))*dsigw(k))

         END DO      ! fin de la boucle en k
       ELSE
         pn2(:,i,j) = 0.0_rsh
         boz(:,i,j) = 0.0_rsh


         roswat_bot(i,j)= (1.0_rsh+beta*(sal(1,i,j)-saliref_lin)-alpha*(temp(1,i,j)-tetaref_lin))*rhoref

       END IF       ! fin du test sur le masque terre

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


   CALL ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,boz(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,boz(:,liminm1:limaxp2,ljminm1:ljmaxp2))


   !PRINT*, 'exit TS_STATEEQ_LINEAR routine'

  END SUBROUTINE ts_stateeq_linear
