!*****************************************************************
   subroutine obcch3d_vz_exp(uz,vz,boz,i,j,frontiere)
!*****************************************************************
!   objet     : conditions aux limites ouvertes sur vz au point i,j
!               et pour un niveau k
!   appelant  : limauz
!   auteur    : Joel MARIN
!
!   parameters : i,j : position
!                  k : niveau
!*****************************************************************
!
!importation des variables globales
USE parameters
use comvarp2d
use comvars2d
use comvars3d
use comcfo
   USE obccombine, ONLY : obc_vze,obc_vzw,obc_vzn,obc_vzs
implicit none
!
! Declaration des paremetres
integer,intent(in)                                 :: i,j,frontiere
real(kind=rsh),dimension(kmax,imin:imax,jmin:jmax),intent(in) :: uz,boz
real(kind=rsh),dimension(kmax,imin:imax,jmin:jmax),intent(out) :: vz
!
! Declaration locales
integer          :: k,n,nbrModesATraiter
real(kind=rsh)   :: cn,g,lambda1,lambda2,lambda3,h_n,v_xe
real(kind=rlg)   :: v_next,h_next,v_n1,h_n1,w1,w2,w3
real(kind=rsh),dimension(1:kmax) :: vz_xe, p,v_n
!
!*****************************************************************

g=9.80665

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

!! frontiere contient un entier 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 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_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
     v_n(1)=v(i,j)
  else
     v_n(1)=-v(i,j)
  endif


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

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

!Calcul de v_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=(u(i-1,j)+u(i,j))*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 vz_en_xe(vzext,vz_xe,i,j)
      call scalaire(frontiere,vz_xe,v_next,i,j,n) 
      w2=v_next
  else
      call vz_en_xe(vz,vz_xe,i-1,j)
      call scalaire(frontiere,vz_xe,v_n1,i-1,j,n)
      w2=v_n1
  endif
  
!Calcul de v_n et h_n au dernier point en h ...
!
  v_n(n)=w2
  
  enddo
  
!Recalcul de vz sur la frontiere
! Calcul de vz au point de h
  call obcch3d_recompose(frontiere, v_n, vz, i, j, nbrModesATraiter)
! Calcul de vz au point v
! pas fait ... cf rapport


else if (frontiere==4) then

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

  if (barotrope_est_ouest(frontiere,j)) 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_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
     v_n(1)=v(i,j)
  else
     v_n(1)=-v(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=-(u(i+1,j)+u(i,j))*0.5
  
!Calcul des conditions aux limites sur les w
  if (lambda2 < 0) then
!w1 entrant -> on specifie la valeur exterieure
!      write(*,*) 'w2 entrant'
!cval      call vz_en_xe(vzext,vz_xe,i+1,j)
      call scalaire(frontiere,vz_xe,v_next,i+1,j,n)
      w2=v_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

!      write(*,*) 'w2 sortant'
      
      call vz_en_xe(vz,vz_xe,i+2,j)
      call scalaire(frontiere,vz_xe,v_n1,i+2,j,n)
      
      w2=v_n1
  endif
  
!Calcul de v_n et h_n au dernier point en h ...
!
  v_n(n)=w2

  enddo
  
!Recalcul de vz sur la frontiere
  call obcch3d_recompose(frontiere, v_n, vz, i, j, nbrModesATraiter)


   
!Attention ... a changer pour la version finale
!do k=1,kmax
!  vz(k,i,j)=v_n(1)
!enddo
  
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
     v_n(1)=v(i,j)
  else
     v_n(1)=-v(i,j)
  endif


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

    cn=Cn_nord_sud(frontiere,i,n)
    
    v_xe=0.5*(v(i,j)+v(i,j-1))
    
    lambda1=v_xe+cn
    lambda3=v_xe-cn
    
    if (lambda1<0 .or. lambda3<0) then
!cval       call vz_en_xe(vzext,vz_xe,i,j)
       call scalaire(frontiere,vz_xe,v_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 vz_en_xe(vz,vz_xe,i,j-1)
       call scalaire(frontiere,vz_xe,v_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,j-1,h0(i,j)+xext(i,j))
       call scalaire(frontiere,p,h_n1,i,j,n)
       h_n1=h_n1/g

    endif

!Calcul de v_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
  
  v_xe=(v(i,j-1)+v(i,j))*0.5
  
!Calcul des conditions aux limites sur les w
  if (lambda1 < 0) then
!w1 entrant -> on specifie la valeur exterieuret
!      write(*,*) 'w1 entrant'
      w1=v_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=v_n1+(g/cn)*h_n1

  endif
  
  if (lambda3 < 0) then
!w3 entrant -> on specifie la valeur exterieure
!      write(*,*) 'w3 entrant'
      w3=v_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=v_n1-(g/cn)*h_n1

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

  vz(:,i,j)=2*vz(:,i,j)-vz(:,i,j-1)
  
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
     v_n(1)=v(i,j)
  else
     v_n(1)=-v(i,j)
  endif


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

!Calcul de v_n, h_n et cn

  cn=Cn_est_ouest(frontiere,i,n)
  
  v_xe=0.5*(v(i,j)+v(i,j-1))
  
  lambda1=-v_xe-cn
  lambda3=-v_xe+cn
  
  
  if (lambda1<0 .or. lambda3<0) then
  
!cval    call vz_en_xe(vzext,vz_xe,i,j+1)
    call scalaire(frontiere,vz_xe,v_next,i,j+1,n)
  
    call obcch3d_calculep(p,bzext,xext,i,j+1,h0(i,j+1)+xext(i,j+1))
    call scalaire(frontiere,p,h_next,i,j+1,n)
    h_next=h_next/g

  endif

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

    call vz_en_xe(vz,vz_xe,i,j+2)
    call scalaire(frontiere,vz_xe,v_n1,i,j+2,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,j+2,h0(i,j+1)+xext(i+1,j+1))
    call scalaire(frontiere,p,h_n1,i,j+2,n)
    h_n1=h_n1/g
    w1=v_n1+(g/cn)*h_n1

  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
  
  v_xe=(v(i,j+1)+v(i,j))*0.5

!Calcul des conditions aux limites sur les w
  if (lambda1 < 0) then
!w1 entrant -> on specifie la valeur exterieure
!      write(*,*) 'w1 entrant'
      w1=v_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=v_n1+(g/cn)*h_n1
  endif
  
  if (lambda3 < 0) then
!w3 entrant -> on specifie la valeur exterieure
!      write(*,*) 'w3 entrant'
      w3=v_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=v_n1-(g/cn)*h_n1
  endif  
  
!Calcul de u_n et h_n au dernier point en h ...
!
  v_n(n)=(w1+w3)/2
  h_n=cn/g*((w1-w3)/2)
  
  enddo
  
!Recalcul de vz sur la frontiere
! Calcul de vz au point de h
  call obcch3d_recompose(frontiere, v_n, vz, i,j, nbrModesATraiter)
! Calcul de uz au point de v
  vz(:,i,j)=2*vz(:,i,j)-vz(:,i,j+1)

else
   write (*,*) 'routine cfovz ... numero de frontiere invalide'
   stop
endif

end subroutine obcch3d_vz_exp
