!********************************************************************
!     
#ifdef key_obc_char
      subroutine cal_modes(frontiere,nbr_modes,tailleFrontiere)
!     
!********************************************************************
!     
!     objet    : calcul des differents modes pour la methode des caracteristiques
!                Appele en debut de calcul
!                
!
!     appelant         : main 
!                      :
!     appeles          : 
!     canal d ecriture :
!     canal de lecture :
!     cle cpp          : 
!                      :
!     variables in     : frontiere,nbr_modes,tailleFrontiere
!                      :
!     variables out    : Mn,cn
!                      :
!     variables inout  :  

!fortcom: Error: smallf/cal_modes.f90, line 335: Subscript #2 of the array N_2 has value 2 which is greater than the upper bound of 1
!    N_2(i,1)=N_2(i,2)
!-------------^
!fortcom: Error: smallf/cal_modes.f90, line 336: Subscript #2 of the array N_2 has value 0 which is less than the lower bound of 1
!    N_2(i,kmax)=N_2(i,kmax-1)
!----------------^
!fortcom: Error: smallf/cal_modes.f90, line 359: Subscript #2 of the array N_2 has value 0 which is less than the lower bound of 1
!     a= 1.0_rlg/((N_2(i,kmax-1)*(sig(kmax)-sig(kmax-1)))&
!------------------^
!compilation aborted for smallf/cal_modes.f90 (code 1) 
!
!***********************************************************

USE parameters
use comvars2d
use comvars3d
use comvarp2d
use comvarp3d
!use comdiag
!use comdem
!use comfluxdo
!use comtraceur
!use comcol2d
use comcfo
implicit none

!****************************************************************************!
! Variables importees ...
! -----------------------
! sig,sigw : subdivision ..
! kmax : taille de la subdivision
! temp : temperature de l eau en degres !!!
!****************************************************************************!


!nbr_modes designe le nombre de modes, si nbr_modes < 2 seul le mode barotrope
!est calcule dans coef_ligne et coef_colonne

!**********************************************************************
!! * Local declarations
!**********************************************************************
	
! Parametre generaux
integer, intent(in) :: frontiere,nbr_modes,tailleFrontiere
	
!i parcourt les points de la frontiere
integer    :: jj,jk,i,ik,indmin,indmax 
real(kind=rlg)                  ::temp_cn
real(kind=rlg),dimension(1:kmax)::temp_mn
real(kind=rlg), dimension(1,kmax) :: N_2 

real(kind=rlg),dimension(kmax,imin:imax,jmin:jmax)  :: mrho
real(kind=rlg)                    :: rho_temp
logical                           :: rho_cst
real(kind=rlg), dimension(0:tailleFrontiere,1:kmax)       :: rho_front
real(kind=rlg)   :: rho0

! Variables pour la fonction DGEEV de LAPACK
character :: JOBVL, JOBVR

integer :: info, LDVL, LDVR
integer :: LWORK, LDA

real(kind=rlg), allocatable, dimension(:,:) :: MAT, VL, VR
real(kind=rlg), dimension(10*kmax) :: WORK
real(kind=rlg), dimension(kmax) :: 	WR, WI

!jm des compteurs pour les stqtistiques
integer                         :: nbr_rho_cst,nbr_rho_mauvais

! Declarations rajoutees par J.MARIN (provisoires)
        	
real(kind=rlg) :: a,b,H,alpha,norm

! Declarations provisoires pour l impact du nombre de modes
! sur la qualite de la recomposition ...

real(kind=rsh),dimension(kmax,imin:imax,jmin:jmax) :: prov
real(kind=rsh)                        :: ZZu0
real(kind=rsh),dimension(1:kmax)      :: ZZun
integer                               :: k
real(kind=rsh)                        :: res

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

rho0=rhoref
nbr_rho_cst=0
nbr_rho_mauvais=0

!! Lecture du profil moyen de densite choisi pour la periode
!

!On n a besoin de rho qu a proximite des frontieres

!Frontieres Nord/Sud

do i=imin,imax
   do ik=1,kmax
      mrho(ik,i,jb(i))=-rho0/grav*boz(ik,i,jb(i))+rho0
      mrho(ik,i,jb(i)+1)=-rho0/grav*boz(ik,i,jb(i)+1)+rho0
      mrho(ik,i,jh(i))=-rho0/grav*boz(ik,i,jh(i))+rho0
      mrho(ik,i,jh(i)-1)=-rho0/grav*boz(ik,i,jh(i)-1)+rho0
   enddo
enddo

!Frontieres Est/Ouest

do jj=jmin,jmax
   do ik=1,kmax
      mrho(ik,ig(jj),jj)=-rho0/grav*boz(ik,ig(jj),jj)+rho0
      mrho(ik,ig(jj)+1,jj)=-rho0/grav*boz(ik,ig(jj)+1,jj)+rho0
      mrho(ik,id(jj),jj)=-rho0/grav*boz(ik,id(jj),jj)+rho0
      mrho(ik,id(jj)-1,jj)=-rho0/grav*boz(ik,id(jj)-1,jj)+rho0
   enddo
enddo

!
!! frontiere contient un entier representant la frontiere choisie
!! Nord=1,Sud=2,Est=3,Ouest=4; on regle les indices en consequence

 indmin=-1
 indmax=0

 if (frontiere==1) then ! Nord
    
    do i=imin,imax
!jm       if (id(jh(i))>=i .and. i>=ig(jh(i))) then
       if (h0(i,jmax)/=-valmanq) then  !jm test pour savoir si on est sur la terre
          if (indmin==-1) then
	     indmin=i
	  endif
	  if (i>=indmax) then
	     indmax=i
	  endif
          do jk=1,kmax
	     rho_front(i,jk)=mrho(jk,i,jh(i))
	  enddo
       endif
    enddo
    
 else if (frontiere==2) then ! Sud
    
    do i=imin,imax
!jm       if (id(jb(i))>=i .and. i>=ig(jb(i))) then
       if (h0(i,jmin+1)/=-valmanq) then  !jm test pour savoir si on est sur la terre
          if (indmin==-1) then
	     indmin=i
	  endif
	  if (i>=indmax) then
	     indmax=i
	  endif
          do jk=1,kmax
	     rho_front(i,jk)=mrho(jk,i,jb(i))
	  enddo
       endif
    enddo
     
 else if (frontiere==3) then ! Est
    do i=jmin,jmax
!jm       if (jh(id(i))>=i .and. i>=jb(id(i))) then
       if (h0(imax,i)/=-valmanq) then  !jm test pour savoir si on est sur la terre
          if (indmin==-1) then
	     indmin=i
	  endif
	  if (i>=indmax) then
	     indmax=i
	  endif
          do jk=1,kmax
	     rho_front(i,jk)=mrho(jk,id(i),i)
	  enddo
       endif
    enddo
    
 else if (frontiere==4) then ! Ouest
    do i=jmin,jmax
!jm       if (jh(ig(i))>=i .and. i>=jb(ig(i))) then !jm test pour savoir si on est sur la terre
       if (h0(imin+1,i)/=-valmanq) then  !jm test pour savoir si on est sur la terre
          if (indmin==-1) then
	     indmin=i
	  endif
	  if (i>=indmax) then
	     indmax=i
	  endif
          do jk=1,kmax
#ifdef debug_cal_modes
!jm             write(*,*) 'rho_front(',i,',',jk,')=',mrho(jk,ig(i),i)
#endif
	     rho_front(i,jk)=mrho(jk,ig(i),i)
	  enddo
       endif
    enddo
    
 else
    write(*,*) 'cal_modes ... ERREUR ... numero de frontiere invalide !'
    stop
 endif

#ifdef debug_cal_modes
    write(*,*) 'indmin=',indmin,'indmax=',indmax
#endif
!write(*,*) 'rho_front = '
!do k=1,kmax
!   write(*,*) 'NIVEAU ',k,' :'
!   write(*,*) rho_front(:,k)
!enddo
!stop
	
!
!! Boucle sur tous les points de la frontiere ouverte
!

!


DO i=indmin,indmax

#ifdef debug_cal_modes      
  write(*,*)'cal_modes ... tour numero ',i,'indmin',indmin,'indmax',indmax
#endif

!! Pour chaque point de la frontiere ouverte on determine 
!! le nombre de niveaux non masques pour avoir le nombre de
! modes

    ! allocation dynamique des matrices

    allocate(MAT(kmax,kmax))
    allocate(VL(kmax,kmax)) 
    allocate(VR(kmax,kmax)) 

    VL(:,:)=0.0_rlg
    VR(:,:)=0.0_rlg

    WORK(:)=0.0_rlg
    WR(:)=0.0_rlg
    WI(:)=0.0_rlg

    MAT(:,:)=0.0_rlg




!
!! On doit traiter le cas rho_front constant (pas de mode barocline)
!
      
!*****************************************************************
!! test pour savoir si rho est constant ... 
!! ATTENTION : a remplacer plus tard par un test statistique ...
!*****************************************************************
  rho_cst=.false.
  rho_temp=rho_front(i,1)
  
  if (i==5) then
     write(*,*) rho_front(i,:)
  endif
	
  do jk=2,kmax
    !jm on teste si l ecart est suffisamment important
    !jm et si rho est bien une fonction decroissante de sigma
    if (abs((rho_front(i,jk)-rho_temp)/rho_temp) < 0.000001) then
       rho_cst=.true.
       nbr_rho_cst=nbr_rho_cst+1
    endif
    if ((rho_front(i,jk)-rho_temp)>0) then
       write(*,*) '*****************************************************************'
       write(*,*) 'Attention, dans cal_modes ... la densite n est pas une fonction'
       write(*,*) 'decroissante de sigma -> desactivation du calcul des modes'
       write(*,*) '*****************************************************************'
       rho_cst=.true.    
       nbr_rho_cst=nbr_rho_cst+1
       nbr_rho_mauvais=nbr_rho_mauvais+1
    endif
    rho_temp=rho_front(i,jk)
  enddo

!jm on enregistre si rho est constant ou non pour cfouzexp et cfovzexp ...	
  if (frontiere <3) then
    barotrope_nord_sud(frontiere,i)=rho_cst;
  else   
    barotrope_est_ouest(frontiere,i)=rho_cst; 	
  endif

#ifdef debug_cal_modes
if (.not.rho_cst) then
   write(*,*) rho_front(i,:)
   write(*,*) 'On a un point ou on calcule les modes'
endif
#endif

#ifdef debug_cal_modes  
  write(*,*) 'rho_constant ?',rho_cst
  write(*,*) rho_front(i,:)
#endif
	
  if (.not.rho_cst) then ! rho_front pas constant 

!*********************************************************************
!        provisoire ... il faut calculer N_2 a partir de rho ..      !
!*********************************************************************

    if (frontiere==1) then
       H=h0(i,jh(i))    
    else if (frontiere==2) then
       H=h0(i,jh(i))
    else if (frontiere==3) then
       H=h0(id(i),i)    
    else if (frontiere==4) then
       H=h0(ig(i),i)    
    else    
       write(*,*) 'cal_modes ... numero de frontiere non valide : ',frontiere
       stop    
    endif


    alpha=1026.0_rlg
    a=-4.2_rlg
    b=1000.0_rlg

!    do jk=1,kmax
!! Attention ... N doit
!      N_2(i,jk)=-grav/(rho0)*(alpha*a/(-H*sigw(jk)+H+b)**2*exp(a/(-H*sigw(jk)+H+b)))
!!      N_2(i,jk)=0.000001_rlg
!    enddo

!*********************************************************************
!                  CE QU IL FAUDRA FAIRE PLUS TARD ...               !
!*********************************************************************
!
    do jk=2,kmax-1
       ! schema centre (Ordre 1 car les sigma ne sont pas equidistants)
       N_2(i,jk)=-grav/(rho0*H)*(rho_front(i,jk+1)-rho_front(i,jk))/(sigw(jk+1)-sigw(jk))
       if (N_2(i,jk)==0.0_rlg) then
          write(*,*) 'cal_modes : ATTENTION N_2=0 ......... ARRET !'
	  stop
       endif
    enddo

! Pour les valeurs aux bords on extrapole ...
! Pour l instant, extrapolation d ordre 1 -> on recopie les valeurs interieures

    N_2(i,1)=N_2(i,2)
    
    N_2(i,kmax)=N_2(i,kmax-1)
#ifdef debug_cal_modes
   write(*,*) 'Rho sur la frontiere consideree ...'
   write(*,*) rho_front(i,:)
   write(*,*) 'N_2 calcule ... '
   write(*,*) N_2(i,:)
#endif

#ifdef debug_cal_modes            
    open(201,FILE='defN.m',FORM='formatted',access='sequential')


!    do jk=1,kmax
!      write(201,*) 'N(',jk,')=',sqrt(N_2(i,jk)),';'
!   enddo

    close(201)
!jm debug_cal_modes
#endif

!! Initialisation : en un point de la frontiere : 
!! boucle sur tous les niveaux pour calculer (u,v,t,s) 
!

!*************************************************
!! Calcul des coefficients de la matrice a inverser 
!! pour obtenir Mn et Cn 
!*************************************************

#ifdef debug_cal_modes
   write(*,*) 'Sigma et sigmaw ...'
   write(*,*) sig
   write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
   write(*,*) sigw
#endif

!     write(*,*) 'b',N_2(i,1),sig(1+1)
	
     b=1.0_rlg/((N_2(i,1)*(sig(1+1)-sig(1)))&
      		* (sigw(1)-sigw(1-1)))

		
     MAT(1,1)=b
     MAT(1,2)=-b
	
     do jk=2,kmax-1
       a= 1.0_rlg/((N_2(i,jk-1)*(sig(jk)-sig(jk-1)))&
       	        * (sigw(jk)-sigw(jk-1)))
       b=1.0_rlg/((N_2(i,jk)*(sig(jk+1)-sig(jk)))&
      		* (sigw(jk)-sigw(jk-1)))
		
       MAT(jk,jk)=a+b
       MAT(jk,jk-1)=-a
       MAT(jk,jk+1)=-b
       
#ifdef debug_cal_modes
   write(*,*) 'M(',jk,',',jk,')=',MAT(jk,jk)
   write(*,*) 'M(',jk,',',jk-1,')=',MAT(jk,jk-1)
   write(*,*) 'M(',jk,',',jk+1,')=',MAT(jk,jk+1)
#endif
!       write(*,*)'Diag de MAT',jk,a+b
     enddo

     a= 1.0_rlg/((N_2(i,kmax-1)*(sig(kmax)-sig(kmax-1)))&
       	* (sigw(kmax)-sigw(kmax-1)))

     b= 1.0_rlg/(grav/H+N_2(i,kmax)*(sigw(kmax)-sig(kmax)))&
       	* 1.0_rlg/(sigw(kmax)-sigw(kmax-1))

     MAT(kmax,kmax-1)=-a
     MAT(kmax,kmax)=a+b

#ifdef debug_cal_modes
   write(*,*) 'M(',kmax,',',kmax,')=',MAT(kmax,kmax)
   write(*,*) 'M(',kmax,',',kmax-1,')=',MAT(kmax,kmax-1)
#endif
	
!     write(*,*) 'MATRICE !!!!'
!     write(*,*) MAT(:,:)

!! Diagonalisation LAPACK, les vecteurs propres sont normes 
     LDA=kmax
! Modification J.MARIN : avant, c etait 10*nbr_modes
! Ca plantait si nbr_modes<<kmax
     LWORK=10*kmax   
     JOBVL='N'
     JOBVR='V'
     LDVR=kmax
     LDVL=kmax  

#ifdef key_obc_char_barocline
     call DGEEV( JOBVL, JOBVR, kmax, MAT, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )  
#endif

!     write(*,*) 'Renormalisation des vecteurs propres ...'

#ifdef debug_cal_modes
   write(*,*) 'COUCOU ... la matrice est inversee ...'
#endif


! Calcul des cn :
    do jj=1,kmax
      WR(jj)=H*sqrt(1.0_rlg/WR(jj))
    enddo	

#ifdef debug_cal_modes
   write(*,*) 'COUCOU ... les valeurs propres sont calculees ...'
#endif

!!On trie les valeurs propres par ordre decroissant (par bulle)
    ik=kmax
    do while (ik>1)
      jk=2
      do while (jk<=ik) 
        if (WR(jk-1) < WR(jk)) then 
	  temp_cn=WR(jk-1) 
	  temp_mn(1:kmax)=VR(:,jk-1)
!
	  WR(jk-1)=WR(jk)
	  VR(1:kmax,jk-1)=VR(1:kmax,jk)
!
	  WR(jk)=temp_cn
	  VR(1:kmax,jk)=temp_mn(1:kmax)
	endif
	jk=jk+1
      enddo
      ik=ik-1
    enddo

!write(*,*) 'Renormalisation des vecteurs propres ...'

!write(*,*) VR(:,1)
!write(*,*) VR(:,2)

! Renormalisation des vecteurs propres ...
#ifdef debug_cal_modes  
   write(*,*) 'Renormalisation des vecteurs propres ...'
#endif

     do jk=1,kmax
       call obcch3d_scalaire_modes(norm,VR,jk,jk)
       do jj=1,kmax
          VR(jj,jk)=VR(jj,jk)/sqrt(norm)
       enddo
     enddo	

!write(*,*) VR(:,1)
!write(*,*) VR(:,2)


! call obcch3d_scalaire_modes(norm,VR,1,2)
! write(*,*) norm
 
#ifdef debug_cal_modes  
    write(*,*) 'Calcul des cn ...'

    write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
    write(*,*) 'Vecteurs propres'
    write(*,*) 'Dans le fichier defVR.m'

    open(201,FILE='defVR.m',FORM='formatted',access='sequential')

    do jj=1,nbr_modes
      do jk=1,kmax
	write(201,*) 'VR(',jk,',',jj,')=',VR(jk,jj),';'
      enddo
    enddo
 
    close(201)

    write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
    write(*,*) 'Valeurs propres'

    write(*,*) WR

    write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' 
!jm debug_cal_modes    
#endif
  ! partie a ameliorer en utilisant du FORTRAN 90 :-)
  
  
  if (frontiere<3) then
    do jk=1,kmax
       Cn_nord_sud(frontiere,i,jk)=WR(jk)
       do jj=1,kmax
          Mn_nord_sud(frontiere,i,jk,jj)=VR(jk,jj)
       enddo
    enddo
  else
    do jk=1,kmax
       Cn_est_ouest(frontiere,i,jk)=WR(jk)
       do jj=1,kmax
          Mn_est_ouest(frontiere,i,jk,jj)=VR(jk,jj)
!	  write(*,*) 'Mn_est_ouest(',frontiere,',',i,',',jk,',',jj,')=',Mn_est_ouest(frontiere,i,jk,jj)
       enddo
    enddo
  endif
   

  else ! rho_cst == .true.
    !Dans le cas barotrope, on a Mn(k)=0 pour k<>0
    !Mn(0)=cst=1/(eta+H)
    !Cn(0)=sqrt(grav*h0)
    if (frontiere==1) then !frontiere nord
       if (hy(i,jh(i)-1) /= -valmanq-h0fond) then
          Mn_nord_sud(frontiere,i,:,1)=1.0_rsh
          Cn_nord_sud(frontiere,i,1)=sqrt(grav*hx(i,jh(i)-1))
       endif
!       Mn_nord_sud(frontiere,i,jk,1)=1.0_rsh
!       Cn_nord_sud(frontiere,i,jk)=sqrt(grav*hx(i-1,j))
    else if (frontiere==2) then !frontiere sud
       if (hy(i,jb(i)+1) /= -valmanq-h0fond) then
          Mn_nord_sud(frontiere,i,:,1)=1.0_rsh
          Cn_nord_sud(frontiere,i,1)=sqrt(grav*hx(i,jb(i)+1))
       endif
!       Mn_nord_sud(frontiere,i,jk,1)=1.0_rsh
!       Cn_nord_sud(frontiere,i,jk)=sqrt(grav*hx(i-1,j))
    else if (frontiere==3) then !frontiere est
       if (hx(id(i)-1,i) /= -valmanq-h0fond) then
          Mn_est_ouest(frontiere,i,:,1)=1.0_rsh
          Cn_est_ouest(frontiere,i,1)=sqrt(grav*hx(id(i)-1,i))
       endif
    else !frontiere ouest
       if (hx(ig(i)+1,i) /= -valmanq-h0fond) then
          Mn_est_ouest(frontiere,i,:,1)=1.0_rsh
          Cn_est_ouest(frontiere,i,1)=sqrt(grav*hx(ig(i)+1,i))
       endif
    endif

!    write(*,*) 'Cas rho constant non traite ...'
!    stop
  endif ! rho_cst == .false.
  !stop



  !write(*,*) 'Calcul des modes au point ',i,' termine !'

    deallocate(MAT)
    deallocate(VL) 
    deallocate(VR)
    
 !---------------------------------------------------
 ! ATTENTION PROVISOIRE : POUR VOIR LA DEGRADATION 
 ! EN PROJETANT SUR UN NOMBRE LIMITE DE MODES
 !---------------------------------------------------
! 
! prov(:,:,:)=0.0_rsh
! 
! do k=1,kmax
!    prov(k,1,1)=sqrt(N_2(1,k))
! enddo
! ZZun(:)=0.0_rsh
! 
!!test de la recomposition dans le cas de la frontiere est ...
! 
!write(*,*) 'Nombre de modes pour la recomposition ',nModes 
! do k=1,nModes
!   call scalaire(frontiere,prov(:,1,1),ZZu0,1,1,k)
!   ZZun(k)=ZZu0
! enddo
!
!write(*,*) ZZun
!
! prov(:,:,:)=0.0_rsh
!
! call obcch3d_recompose(frontiere, ZZun, prov, 1,1,nModes) 
! 
! write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
! write(*,*) N_2(1,:)
! write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' 
! write(*,*) prov(:,1,1)
! write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
! 
! !---------------------------------------------------
! ! CALCUL DU DEGRE DE RECOMPOSITION (RMS)
! !---------------------------------------------------
! res=0
! do k=1,kmax
!    res=res+dsigw(k)*((prov(k,1,1)-sqrt(N_2(1,k)))/sqrt(N_2(1,k)))**2
! enddo
! res=100*sqrt(res)
! 
! write(*,*) 'pourcentage de recomposition de N^2 : ',100-res
! 
! 
! stop
    

enddo


!jm affichage des statistiques
write(*,*) 'Nombre de profils avec variation de rho insuffisante :',nbr_rho_cst
write(*,*) 'Nombre de profils avec rho non croissant avec la profondeur :',nbr_rho_mauvais

write(*,*)'Fin routine cal_modes'

#else
subroutine cal_modes
print*,'should not see this routine cal_modes (not ready yet)'
#endif
end subroutine cal_modes

