  SUBROUTINE dyn3d_stabts(salt,temp,boz)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE dyn3d_stabts  ***
   !&E
   !&E ** Purpose : 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : stepu
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays : 
   !&E
   !&E ** Modified variables : sal, temp, b
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2004-06-13
   !&E       !  2008-12   (V. Garnier) shift of sigma values : -1 < sig < 0
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d
   USE comvars3d
   USE_MPI toolmpi, ONLY: ex_i_rsh,ex_j_rsh

   IMPLICIT NONE

   !! * Arguments
   REAL(kind=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: salt,temp,boz

   !! *  Local declarations
   INTEGER        :: i,j,k,kk
   REAL(kind=rlg) :: pres,sal,tep,c,cc
   REAL(kind=rlg) :: rost,rostp,row,alpha,beta

   !!----------------------------------------------------------------------
   !! * Executable part
   
   alpha=0.16d-03
   beta=0.8d-03

   ! Remark : this subroutine is called from stepu.F90
   ! As stepu.F90 is only called when nbouc is an odd number,
   ! mod(nbouc,24) is never equal to 0.
   ! 24 is changed into 13 so that it is called for each 26th nbouc
   ! ie nbouc= 13, 39, 65, ... (nbouc=26, 52, ...  does not come into stepu)

   IF (MOD(nbouc,13) ==  0) THEN

     DO j=ljminp1,ljmax
       DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
         DO k=kmax,2,-1
           IF( -sig(k)*h0(i,j) <= 500.0_rsh .AND. h0(i,j) > hm &
             .AND. boz(k,i,j) < boz(k-1,i,j)*0.99_rsh) THEN
             temp(k-1:kmax,i,j)=SUM(temp(k-1:kmax,i,j)*dsigu(k-1:kmax))/SUM(dsigu(k-1:kmax))
             salt(k-1:kmax,i,j)=SUM(salt(k-1:kmax,i,j)*dsigu(k-1:kmax))/SUM(dsigu(k-1:kmax))
             !
             IF(.NOT.l_stateeq_lin) THEN
               tep = REAL(temp(k-1,i,j),rlg)
               sal = REAL(salt(k-1,i,j),rlg)
               !    pres=-sig(k-1)*h0(i,j)
               pres=0.0_rlg
               row=999.842594_rlg +6.793952d-02*tep -9.095290d-03*tep**2 &
                                  +1.001685d-04*tep**3 -1.120083d-06*tep**4 &
                                  +6.536332d-09*tep**5 !  densite de l eau pure
               
               !-----------------------------------------------------------
               !  densite de l eau en surface  salinite S (o/oo)
               !  temperature,potentielle T (DEG C)
               !-----------------------------------------------------------
               
               rost=row+sal*(0.824493_rlg -4.0899d-03*tep +7.6438d-05*tep**2           &
                                          +8.2467d-07*tep**3 +5.3875d-09*tep**4)       &
                       +sal**1.5_rsh*(-5.72466d-03 +1.0227d-04*tep -1.6546d-06*tep**2) &
                       +4.8314e-04*sal**2
               c=1449.2_rlg+1.34_rlg*(sal-35.0_rlg) +4.55_rlg*tep-0.045_rlg*tep**2     &
                           +0.00821_rlg*pres +15d-09*pres**2
               cc=c*c
               rostp=rost +10000.0_rlg*(1.0_rlg-0.20_rlg*pres/cc)*pres/cc
               boz(k-1,i,j)=grav*(rhoref-rostp)/rhoref
             ELSE
               boz(k-1,i,j)=grav*(alpha*(temp(k-1,i,j)-tetaref_lin)-beta*(salt(k-1,i,j)-saliref_lin))
             ENDIF
           ENDIF ! fin du premier gros test
         END DO ! fin de boucle sur k
       END DO ! fin de boucle sur i
     END DO ! fin de boucle sur j

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

     
   ENDIF ! fin de test sur nbouc

  END SUBROUTINE dyn3d_stabts
