!*******************************************************************
!
    subroutine obcch3d_calculep(p,boz,xe,i,j,Heta)
!
!*******************************************************************
!
! objet : calcul de la pression P=p/rhoref sur toute la profondeur
!         en un point i,j du domaine
!
! apppelant : cfouxexp, cfovzexp
!
! auteur : Joel MARIN
!********************************************************************

USE parameters
use comvars2d
use comvars3d
use cometdyn
implicit none

! Declaration des parametres
real(kind=rsh),dimension(1:kmax),intent(out)       :: p
real(kind=rsh),dimension(kmax,imin:imax,jmin:jmax),intent(in) :: boz
real(kind=rsh),dimension(imin:imax,jmin:jmax),intent(in)      :: xe
!jm 12_06_2006
!jm Heta passe en argument pour pouvoir calculer la pression
!jm a la meme profondeur sur la frontiere et sur le point interieur
real(kind=rsh),intent(in)                          :: Heta
real(kind=rsh)                                     :: Heta2 
integer,intent(in) :: i,j
! Declarations locales
integer                                            :: k
real(kind=rsh)                                     :: int

   print*,'STILL WITH 0<sig<1'
   print*,'The simulation is stopped'
   stop

!jm pour la formule de correction du terme grav(1-sigma)
!jm invariant si on a appele la fonction au point (i,j)
Heta2=h0(i,j)+xe(i,j)
!Heta2=Heta



!k=kmax
int=boz(kmax,i,j)*(1.0-sig(kmax))

!jm p(kmax)=Heta2*(grav*(1.0-((Heta2-Heta)/Heta2+Heta/Heta2*sig(kmax)))-int)+pa(i,j)/rhoref
p(kmax)=Heta2*(grav*(1.0-sig(kmax))-int)+pa(i,j)/rhoref


do k=kmax-1,1,-1
! Formule des trapezes
   int=int+0.5*(boz(k,i,j)+boz(k+1,i,j))*dsigw(k)

!jm   p(k)=Heta2*(grav*(1.0-((Heta2-Heta)/Heta2+Heta/Heta2*sig(k)))-int)+pa(i,j)/rhoref
   p(k)=Heta2*(grav*(1.0-sig(k))-int)+pa(i,j)/rhoref

enddo

 end subroutine obcch3d_calculep  
