  SUBROUTINE ts_stateeq_blumberg(sal,temp,boz,pn2)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE ts_stateeq_blumberg  ***
   !&E
   !&E ** Purpose : buoyancy estimation from a unlinear state equation
   !&E
   !&E ** Description :
   !&E                  TP potential temperature in DEG C
   !&E                  = in situ temperature for small depths
   !&E                  S salinity in o/oo
   !&E                  p pressure in decibars = depth in meters
   !&E
   !&E ** Called by : main, stepu, stepv
   !&E
   !&E ** Used ij-arrays : h0,temp,sal
   !&E 
   !&E ** Modified variables : boz,pn2
   !&E
   !&E ** Reference : Mellor, 609-611 JAOT 1991
   !&E
   !&E ** History :
   !&E       !  2000-01  (P. Garreau) Original code
   !&E       !  2007-10  (V. Garnier) Cleaning
   !&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
   !&E       !  2009-04  (A. Mano) MPI2D
   !&E       !  2010-12  (P. Garreau) Minimum of S and T to avoid underflow
   !&E                   Later, we should take care of the validity range of this state equation
   !&E       !  2014-07  (P. Garreau) : Estimate of the local brunt-vaisala frequency (pn2) [from NEMO eosnb2.F90 McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987]
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d, ONLY : h0,hm
   USE comvars3d, ONLY : sig,rhoref,dsigw
   USE comvarp2d, ONLY : ssh
   !USE toolgeom, ONLY : f_zu
#ifdef key_profil_back
   USE comvars3d, ONLY : rho2,boz2
#endif
#if defined key_siggen || defined key_gencoord
   USE toolgeom,  ONLY : f_hzw
#endif
#ifdef key_potential_N2
   USE toolgeom,  ONLY : f_dzw,f_zw
#endif
#if defined key_sedim_MUSTANG || defined key_obstructions
   USE comvars3d,  ONLY : roswat_bot,saliref_lin,tetaref_lin
#endif
#if defined key_sedim_MUSTANG && defined key_roswat_depend_partcleconc
   USE comsubstance, ONLY : nvpc,cv_wat
   USE comMUSTANG,    ONLY : ros
#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                  :: k,i,j
   REAL(KIND=rsh)           :: TP,S,p,roref
   REAL(kind=rsh)           :: hz_w
   REAL(KIND=rlg)           :: row,rost,c,pcc,rostp
   REAL(KIND=rlg)           :: zgde3w, zt, zs, zh, zalbet, zbeta
#if defined key_sedim_MUSTANG && defined key_roswat_depend_partcleconc
   INTEGER                  :: iv
   REAL(KIND=rsh)           :: sumpart
#endif

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

   roref=1040.0_rsh
   !roref=rhoref

!$OMP PARALLEL DO PRIVATE(TP,S,P,row,rost,rostp,c,pcc,hz_w,k,i,j, &
!$OMP&                    zgde3w, zt, zs, zh, zalbet, zbeta) &
!$OMP& SCHEDULE(RUNTIME)
   DO j=ljminp1,ljmax
    DO i=liminp1,limax


!     IF (h0(i,j)+ssh(i,j) > 0.0_rsh) THEN
     IF (h0(i,j)+ssh(i,j) > hm) THEN
     !IF (h0(i,j) > -valmanq) THEN

       DO k=1,kmax

         TP = MAX(temp(k,i,j),0.001_rsh)
         S = MAX(sal(k,i,j),0.0001_rsh)
!         P = -f_zu(h0(i,j),0.0_rsh,k,i,j)
         !P = 10.0_rsh   ! = 1000mbar reference atmospherical pressure
         P = 0.0_rsh   ! = reference atmospherical pressure

   ! density of pure water

         row = 999.842594_rlg + (  6.793952d-02 + ( - 9.095290d-03 +                &
               (1.001685d-04 - 1.120083d-06*TP + 6.536332d-09*TP*TP)*TP )*TP  )*TP

   ! density of sea water at the surface

         rost = row + S * (   0.824493_rlg + (  -4.0899d-03 + ( 7.6438d-05 +        &
                          (-8.2467d-07+5.3875d-09*TP)*TP )*TP  )*TP   ) +           &
                S**1.5_rlg * ( -5.72466d-03 + (1.0227d-04-1.6546d-06*TP)*TP ) +     &
                S*S*4.8314d-04

   ! density of sea water at pressure p (decibars)

         c = 1449.2_rlg + 1.34_rlg*(S-35.0_rlg) + 4.55_rlg*TP - 0.045_rlg*TP*TP +   &
             0.00821_rlg*p + 15.0d-09*p*p
         pcc = p/(c*c)
         rostp = rost + 1.0d4*(1.0_rlg-0.2_rlg*pcc)*pcc
#if defined  key_sedim_MUSTANG && defined key_roswat_depend_partcleconc
! 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
         DO iv=1,nvpc
          sumpart=sumpart+cv_wat(iv,k,i,j)*(1.0_rsh-REAL(rostp,rsh)/ros(iv))
         ENDDO
         rostp=rostp+REAL(sumpart,rlg)
#endif
#if defined  key_sedim_MUSTANG || defined key_obstructions
         IF(k==1)roswat_bot(i,j)=REAL(rostp,rsh)
#endif
        
   ! buoyancy

         boz(k,i,j) = grav*( roref - REAL(rostp,rsh) )/roref
#ifdef key_profil_back
         boz2(k,i,j) = grav*( REAL(rho2(k,i,j),riosh) - REAL(rostp,riosh) )/roref
#endif

       END DO      ! fin de la boucle en k

   ! estimation of the Brunt-vaissala frequency (N2)
       DO k=1,kmax-1

#ifdef key_potential_N2
         zgde3w = grav / f_dzw(h0(i,j),ssh(i,j),k,i,j)
         zt = 0.5 * ( temp(k,i,j) + temp(k+1,i,j) )         ! potential temperature at w-pt
         zs = 0.5 * ( sal(k,i,j) + sal(k+1,i,j) ) - 35.0  ! salinity anomaly (s-35) at w-pt
         zh = f_zw(h0(i,j),ssh(i,j),k,i,j)                           !  depth in meters  at w-point
         !
         zalbet = ( ( ( - 0.255019e-07_rlg * zt + 0.298357e-05_rlg ) * zt   &   ! ratio alpha/beta
                     &                                  - 0.203814e-03_rlg ) * zt   &
                     &                                  + 0.170907e-01_rlg ) * zt   &
                     &   +         0.665157e-01_rlg &
                     &   +     ( - 0.678662e-05_rlg * zs &
                     &           - 0.846960e-04_rlg * zt + 0.378110e-02_rlg ) * zs   &
                     &   +   ( ( - 0.302285e-13_rlg * zh &
                     &           - 0.251520e-11_rlg * zs &
                     &           + 0.512857e-12_rlg * zt * zt              ) * zh   &
                     &           - 0.164759e-06_rlg * zs &
                     &        +(   0.791325e-08_rlg * zt - 0.933746e-06_rlg ) * zt   &
                     &                                  + 0.380374e-04_rlg ) * zh
                     !
         zbeta  = ( ( -0.415613e-09_rlg * zt + 0.555579e-07_rlg ) * zt &   ! beta
                     &                               - 0.301985e-05_rlg ) * zt &
                     &   +       0.785567e-03_rlg &
                     &   + (     0.515032e-08_rlg * zs &
                     &         + 0.788212e-08_rlg * zt - 0.356603e-06_rlg ) * zs &
                     &   + ( (   0.121551e-17_rlg * zh &
                     &         - 0.602281e-15_rlg * zs &
                     &         - 0.175379e-14_rlg * zt + 0.176621e-12_rlg ) * zh &
                     &                                + 0.408195e-10_rlg   * zs &
                     &     + ( - 0.213127e-11_rlg * zt + 0.192867e-09_rlg ) * zt &
                     &                                - 0.121555e-07_rlg ) * zh
         !
         pn2(k,i,j) = zgde3w * zbeta            &   ! N^2
                     &          * ( zalbet * ( temp(k+1,i,j) - temp(k,i,j) )   &
                     &                     - ( sal(k+1,i,j) - sal(k,i,j) ) )

#else
#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
#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)=(rhoref*(1.0_rsh+0.0008_rsh*(sal(1,i,j)-saliref_lin)  &
                                     -0.00016_rsh*(temp(1,i,j)-tetaref_lin)))

#endif

     END IF       ! fin du test sur le masque terre 

    END DO        ! fin de la boucle en i
   END DO         ! fin de la boucle en j 

!$OMP END PARALLEL DO
#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_BLUMBERG routine'

  END SUBROUTINE ts_stateeq_blumberg
