!*****************************************************************
   subroutine obcch3d_uz_exp(uz,vz,boz,i,j,frontiere)
!*****************************************************************
!   objet     : conditions aux limites ouvertes sur uz au point i,j
!               et pour un niveau k
!   appelant  : limauz
!   auteur    : Joel MARIN
!
!   parametres : i,j : position
!                  k : niveau
!*****************************************************************
!
!importation des variables globales
USE parameters
use comvarp2d
use comvars2d
use comvars3d
use comcfo
implicit none
!
! Declaration des paremetres
integer,intent(in)                                 :: i,j,frontiere
real(kind=rsh),dimension(kmax,imin:imax,jmin:jmax),intent(inout) :: uz,vz
real(kind=rsh),dimension(kmax,imin:imax,jmin:jmax),intent(in) :: boz
!
! Declaration locales
integer          :: k,n,nbrModesATraiter
real(kind=rsh)   :: cn,g,lambda1,lambda2,lambda3,h_n,u_xe
real(kind=rlg)   :: u_next,h_next,u_n1,h_n1,w1,w2,w3
real(kind=rsh),dimension(1:kmax)   :: uz_xe, p, u_n
!
!*****************************************************************

g=9.81 

!On connait les cn   ... et on peut calculer les h_n

!! frontiere contient u_n entieru_n1 representant la frontiere choisie
!! Nord=1,Sud=2,Est=3,Ouest=4;

if (frontiere==3) then 



  if (barotrope_est_ouest(frontiere,j)) then
    nbrModesATraiter=1
  else
    nbrModesATraiter=nModes
!    write(*,*) 'Cas barocline dans cfouzexp non traite ...'
!    stop
  endif

!jm cf annexe : Methode des caracteristiques 3D pour MARS en coordonnees sigma
!jm on n utilise pas le premier mode, mais on le remplace par la valeur barotrope
  if (Mn_est_ouest(frontiere,j,1,1)>=0) then 
!jm attention, il faut tester le signe du premier mode, car le mode barotrope est obtenu
!jm pour un premier mode positif
     u_n(1)=u(i,j)
  else
     u_n(1)=-u(i,j)
  endif


  do n=2,nbrModesATraiter
!  do n=1,nbrModesATraiter

    cn=Cn_est_ouest(frontiere,j,n)

    
    u_xe=0.5*(u(i,j)+u(i-1,j))


    
    lambda1=u_xe+cn
    lambda3=u_xe-cn
    
    if (lambda1<0 .or. lambda3<0) then
!cval       call uz_en_xe(uzext,uz_xe,i,j)
       call scalaire(frontiere,uz_xe,u_next,i,j,n)

    
       call obcch3d_calculep(p,bzext,xext,i,j,h0(i,j)+xext(i,j))

       call scalaire(frontiere,p,h_next,i,j,n)
       h_next=h_next/g

    endif
    
    if (lambda1>=0 .or. lambda3>=0) then
       call uz_en_xe(uz,uz_xe,i-1,j)
       call scalaire(frontiere,uz_xe,u_n1,i,j,n)

!Calcul de h_n1 ...
!jm pour le dernier parametre, il faut que ce soit
!jm le meme appel qu au dessus
       call obcch3d_calculep(p,boz,ssh,i-1,j,h0(i,j)+xext(i,j))
       call scalaire(frontiere,p,h_n1,i,j,n)
       h_n1=h_n1/g

    endif

!Calcul de u_n, h_n et cn
  
!Attention ... la partie precedente est a modifier pour la version finale
  
!Calcul des valeurs propres au dernier point en h
!Pour u, on prend la valeur au point de grille precedent
  
!Calcul des conditions aux limites sur les w
  if (lambda1 < 0) then
!w1 entrant -> on specifie la valeur exterieuret
!      write(*,*) 'w1 entrant'
      w1=u_next+(g/cn)*h_next

  else
!w1 sortant -> on extrapole a partir des valeurs interieures
!extrapolation d ordre 1
! Attention valable seulement dans la cas barotrope
!      write(*,*) 'w1 sortant'
      w1=u_n1+(g/cn)*h_n1

  endif
  
  if (lambda3 < 0) then
!w3 entrant -> on specifie la valeur exterieure
!      write(*,*) 'w3 entrant'
      w3=u_next-(g/cn)*h_next

  else
!w3 sortant -> on extrapole a partir des valeurs interieures
!extrapolation d ordre 1
! Attention valable seulement dans la cas barotrope
!      write(*,*) 'w3 sortant'
      w3=u_n1-(g/cn)*h_n1

  endif  
  
!Calcul de u_n et h_n au dernier point en h ...
!
  u_n(n)=(w1+w3)/2
  h_n=cn/g*((w1-w3)/2) ! inutile ... non utilise pour
  
  enddo
    
!Recalcul de uz sur la frontiere
! Calcul de uz au point de h
  call obcch3d_recompose(frontiere, u_n, uz, i,j, nbrModesATraiter)
! Calcul de uz au point de u

  uz(:,i,j)=2*uz(:,i,j)-uz(:,i-1,j)


else if (frontiere==4) then

!  write(*,*) 'subroutine cfouz ... frontiere ouest'





  if (barotrope_est_ouest(frontiere,j)) then
    nbrModesATraiter=1
  else
    nbrModesATraiter=nModes
!    write(*,*) 'Cas barocline dans cfouzexp non traite ...'
!    stop
  endif

!jm cf annexe : Methode des caracteristiques 3D pour MARS en coordonnees sigma
!jm on n utilise pas le premier mode, mais on le remplace par la valeur barotrope
  if (Mn_est_ouest(frontiere,j,1,1)>=0) then 
!jm attention, il faut tester le signe du premier mode, car le mode barotrope est obtenu
!jm pour un premier mode positif
     u_n(1)=u(i,j)
  else
     u_n(1)=-u(i,j)
  endif


  do n=2,nbrModesATraiter
!  do n=1,nbrModesATraiter

!Calcul de u_n, h_n et cn

  cn=Cn_est_ouest(frontiere,j,n)
  
  u_xe=0.5*(u(i,j)+u(i+1,j))
  
  lambda1=-u_xe-cn
  lambda3=-u_xe+cn


  
  
  if (lambda1<0 .or. lambda3<0) then
  
!cval    call uz_en_xe(uzext,uz_xe,i+1,j)
    call scalaire(frontiere,uz_xe,u_next,i+1,j,n)


  
    call obcch3d_calculep(p,bzext,xext,i+1,j,h0(i+1,j)+xext(i+1,j))

    call scalaire(frontiere,p,h_next,i+1,j,n)
    h_next=h_next/g


  endif

  if (lambda1<0 .or. lambda3<0) then

    call uz_en_xe(uz,uz_xe,i+2,j)
    call scalaire(frontiere,uz_xe,u_n1,i+2,j,n)



!Calcul de h_n1 ...
!jm pour le dernier parametre, il faut que ce soit
!jm le meme appel qu au dessus
    call obcch3d_calculep(p,boz,ssh,i+2,j,h0(i+1,j)+xext(i+1,j))
    call scalaire(frontiere,p,h_n1,i+2,j,n)

    h_n1=h_n1/g

  endif
  
!Attention ... la partie precedente est a modifier pour la version finale
  
!Calcul des valeurs propres au dernier point en h
!Pour u, on prend la valeur au point de grille precedent

!Calcul des conditions aux limites sur les w
  if (lambda1 < 0) then
!w1 entrant -> on specifie la valeur exterieure
!      write(*,*) 'w1 entrant'
      w1=u_next+(g/cn)*h_next
  else
!w1 sortant -> on extrapole a partir des valeurs interieures
!extrapolation d ordre 1
! Attention valable seulement dans la cas barotrope
!      write(*,*) 'w1 sortant'
      w1=u_n1+(g/cn)*h_n1
  endif
  
  if (lambda3 < 0) then
!w3 entrant -> on specifie la valeur exterieure
!      write(*,*) 'w3 entrant'
      w3=u_next-(g/cn)*h_next
  else
!w3 sortant -> on extrapole a partir des valeurs interieures
!extrapolation d ordre 1
! Attention valable seulement dans la cas barotrope
!      write(*,*) 'w3 sortant'
      w3=u_n1-(g/cn)*h_n1
  endif  
  
!Calcul de u_n et h_n au dernier point en h ...
!
  u_n(n)=(w1+w3)/2

  h_n=cn/g*((w1-w3)/2)
  
  enddo
  
!Recalcul de uz sur la frontiere
! Calcul de uz au point de h
  call obcch3d_recompose(frontiere, u_n, uz, i,j, nbrModesATraiter)
! Calcul de uz au point de u
  uz(:,i,j)=2*uz(:,i,j)-uz(:,i+1,j)
 


!   write(*,*) 'uz reconstruit ...'
!   write(*,*) uz(:,i,j)

!jm write(*,*) uz(:,i,j)
  
else if (frontiere==1) then

  if (barotrope_nord_sud(frontiere,i)) then
    nbrModesATraiter=1
  else
    nbrModesATraiter=nModes
!    write(*,*) 'Cas barocline dans cfovzexp non traite ...'
!    stop
  endif

!jm cf annexe : Methode des caracteristiques 3D pour MARS en coordonnees sigma
!jm on n utilise pas le premier mode, mais on le remplace par la valeur barotrope
  if (Mn_nord_sud(frontiere,i,1,1)>=0) then 
!jm attention, il faut tester le signe du premier mode, car le mode barotrope est obtenu
!jm pour un premier mode positif
     u_n(1)=u(i,j)
  else
     u_n(1)=-u(i,j)
  endif


  do n=2,nbrModesATraiter
!  do n=1,nbrModesATraiter

!write(*,*) 'subroutine cfovz ... frontiere est'

!Calcul de u_n, h_n et cn
!Pour la partie barotrope, h_n = xe et cn=sqrt(g*h0), v_n=uz
  
!Valeurs exterieures
!plus tard : produit scalaire ... pour l instant, moyenne

!Attention ... la partie precedente est a modifier pour la version finale
  
!Calcul des valeurs propres au dernier point en h
!Pour u, on prend la valeur au point de grille precedent
  lambda2=(v(i,j)+v(i,j-1))*0.5
  
!Calcul des conditions aux limites sur les w
  if (lambda2 < 0) then
!w2 entrant -> on specifie la valeur exterieure
!      write(*,*) 'w2 entrant'
!cval      call uz_en_xe(uzext,uz_xe,i,j)
      call scalaire(frontiere,uz_xe,u_next,i,j,n) 
      w2=u_next
  else
      call uz_en_xe(uz,uz_xe,i,j-1)
      call scalaire(frontiere,uz_xe,u_n1,i,j-1,n)
      w2=u_n1
  endif
  
!Calcul de u_n et h_n au dernier point en h ...
!
  u_n(n)=w2
  
  enddo
  
!Recalcul de uz sur la frontiere
! Calcul de uz au point de h
  call obcch3d_recompose(frontiere, u_n, uz, i, j, nbrModesATraiter)
! Calcul de uz au point u
! pas fait ... cf rapport
  
else if (frontiere==2) then

  if (barotrope_nord_sud(frontiere,i)) then
    nbrModesATraiter=1
  else
    nbrModesATraiter=nModes
!    write(*,*) 'Cas barocline dans cfovzexp non traite ...'
!    stop
  endif

!jm cf annexe : Methode des caracteristiques 3D pour MARS en coordonnees sigma
!jm on n utilise pas le premier mode, mais on le remplace par la valeur barotrope
  if (Mn_nord_sud(frontiere,i,1,1)>=0) then 
!jm attention, il faut tester le signe du premier mode, car le mode barotrope est obtenu
!jm pour un premier mode positif
     u_n(1)=u(i,j)
  else
     u_n(1)=-u(i,j)
  endif


  do n=2,nbrModesATraiter
!  do n=1,nbrModesATraiter

!Attention ... la partie precedente est a modifier pour la version finale
  
!Calcul des valeurs propres au dernier point en h
!Pour u, on prend la valeur au point de grille precedent
  lambda2=-(v(i,j+1)+v(i,j))*0.5
  
!Calcul des conditions aux limites sur les w
  if (lambda2 < 0) then
!w1 entrant -> on specifie la valeur exterieure
!      write(*,*) 'w1 entrant'
!cval      call uz_en_xe(uzext,uz_xe,i,j+1)
      call scalaire(frontiere,uz_xe,u_next,i,j+1,n)
      w2=u_next
  else
!w1 sortant -> on extrapole a partir des valeurs interieures
!extrapolation d ordre 1
! Attention valable seulement dans la cas barotrope
!      write(*,*) 'w1 sortant'
!plus tard : produit scalaire ... pour l instant moyenne
      
      call uz_en_xe(uz,uz_xe,i,j+2)
      call scalaire(frontiere,uz_xe,u_n1,i,j+2,n)
      
      w2=u_n1
  endif
  
!Calcul de u_n et h_n au dernier point en h ...
!
  u_n(n)=w2

  enddo
  
!Recalcul de uz sur la frontiere
  call obcch3d_recompose(frontiere, u_n, uz, i, j, nbrModesATraiter)
   
else
   write (*,*) 'routine cfouz ... numero de frontiere invalide'
   stop
endif



end subroutine obcch3d_uz_exp
