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





























































!! * Modules used
   USE parameters
   USE comvars2d, ONLY : h0,hm
   USE comvars3d, ONLY : sig,rhoref,dsigw
   USE comvarp2d, ONLY : ssh
!USE toolgeom, ONLY : f_zu




   USE comvars3d,  ONLY : roswat_bot,saliref_lin,tetaref_lin


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


!!----------------------------------------------------------------------
!! * 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(k==1)roswat_bot(i,j)=REAL(rostp,rsh)

        
! buoyancy

         boz(k,i,j) = grav*( roref - REAL(rostp,rsh) )/roref


       END DO      ! fin de la boucle en k

! estimation of the Brunt-vaissala frequency (N2)
       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)=(rhoref*(1.0_rsh+0.0008_rsh*(sal(1,i,j)-saliref_lin)  &
                                     -0.00016_rsh*(temp(1,i,j)-tetaref_lin)))



     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

   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_BLUMBERG routine'

  END SUBROUTINE ts_stateeq_blumberg
