MODULE metric
  USE genmod
  USE grid_param
  USE domain_param
  
  TYPE t_cell_glo
    INTEGER :: neighbour(0:5)
    INTEGER :: edge(0:5)
    INTEGER :: assign_face
    INTEGER :: assign_i
    INTEGER :: assign_j
    INTEGER :: assign_domain
    REAL(rstd) :: xyz(3)
  END TYPE t_cell_glo
  
  TYPE t_vertex_glo
    REAL(rstd) :: xyz(3)
    INTEGER    :: neighbour(0:5)
    INTEGER    :: ne(0:5)
    INTEGER    :: ind
    INTEGER    :: delta
  END TYPE t_vertex_glo
  
  TYPE t_edge_glo
   INTEGER :: e1
   INTEGER :: e2
   INTEGER :: assign_domain
   INTEGER :: assign_i
   INTEGER :: assign_j
   INTEGER :: assign_pos
   INTEGER :: assign_delta
  END TYPE t_edge_glo
    
 
  TYPE(t_vertex_glo),ALLOCATABLE,SAVE :: vertex_glo(:,:,:)
  TYPE(t_cell_glo),ALLOCATABLE,SAVE :: cell_glo(:)
  TYPE(t_edge_glo),ALLOCATABLE,SAVE :: edge_glo(:)
  INTEGER :: ncell_glo
  
  INTEGER,ALLOCATABLE,SAVE :: Tab_index(:,:,:)  
  INTEGER,PARAMETER :: right=1
  INTEGER,PARAMETER :: rup=2
  INTEGER,PARAMETER :: lup=3
  INTEGER,PARAMETER :: left=4
  INTEGER,PARAMETER :: ldown=5
  INTEGER,PARAMETER :: rdown=6
  
  INTEGER,PARAMETER :: vrup=1
  INTEGER,PARAMETER :: vup=2
  INTEGER,PARAMETER :: vlup=3
  INTEGER,PARAMETER :: vldown=4
  INTEGER,PARAMETER :: vdown=5
  INTEGER,PARAMETER :: vrdown=6
  
CONTAINS

  SUBROUTINE remap_schmidt
    USE spherical_geom_mod
    USE ioipsl
    IMPLICIT NONE
    INTEGER :: nf,i,j
    REAL(rstd) :: schmidt_factor, schmidt_lon, schmidt_lat

    ! Schmidt transform parameters
    schmidt_factor = 1.
    CALL getin('schmidt_factor', schmidt_factor)
    schmidt_factor =  schmidt_factor**2.
    
    schmidt_lon = 0.
    CALL getin('schmidt_lon', schmidt_lon)
    schmidt_lon = schmidt_lon * pi/180.

    schmidt_lat = 45.
    CALL getin('schmidt_lat', schmidt_lat)
    schmidt_lat = schmidt_lat * pi/180.

    DO nf=1,nb_face
       DO j=1,jjm_glo
          DO i=1,iim_glo
             CALL schmidt_transform(vertex_glo(i,j,nf)%xyz, schmidt_factor, schmidt_lon, schmidt_lat)
          END DO
       END DO
    END DO
  END SUBROUTINE remap_schmidt

  SUBROUTINE allocate_metric
  IMPLICIT NONE
    INTEGER :: ind
    
    ALLOCATE(vertex_glo(1-halo:iim_glo+halo,1-halo:jjm_glo+halo,nb_face))
    ncell_glo=10*(iim_glo-2)*(jjm_glo-2)+(iim_glo-2)*20+12
    ALLOCATE(cell_glo(ncell_glo))
    ALLOCATE(tab_index(nb_face,nb_face,0:5))
    ALLOCATE(edge_glo(ncell_glo*3))
    
    DO ind=1,ncell_glo
      cell_glo(ind)%neighbour(:)=0
    ENDDO
    
  END SUBROUTINE allocate_metric
  
  SUBROUTINE Compute_face
  USE spherical_geom_mod
  IMPLICIT NONE
    
    REAL(rstd) :: len_edge
    REAL(rstd) :: rot=0.
    INTEGER :: nf,i,j
    REAL(rstd),DIMENSION(3) :: p1,p2,p3
    REAL(rstd) :: d1,d2,d3
  
  
   len_edge=acos(cos(Pi/5)*cos(2*Pi/5)/(sin(Pi/5)*sin(2*Pi/5)))
   
   DO nf=1,nb_face/2
     CALL lonlat2xyz(0.,Pi/2,vertex_glo(1,1,nf)%xyz)
     CALL lonlat2xyz(rot+(nf-1)*2*Pi/5,Pi/2-len_edge,vertex_glo(iim_glo,1,nf)%xyz)
     CALL lonlat2xyz(rot+2*Pi/5+(nf-1)*2*Pi/5,Pi/2-len_edge,vertex_glo(1,jjm_glo,nf)%xyz)
     CALL lonlat2xyz(rot+Pi/5+(nf-1)*2*Pi/5,-Pi/2+len_edge,vertex_glo(iim_glo,jjm_glo,nf)%xyz)

     CALL lonlat2xyz(0.,-Pi/2,vertex_glo(1,1,nf+nb_face/2)%xyz)
     CALL lonlat2xyz(rot+Pi/5+(nf-1)*2*Pi/5,-(Pi/2-len_edge),vertex_glo(1,jjm_glo,nf+nb_face/2)%xyz)
     CALL lonlat2xyz(rot+Pi/5+2*Pi/5+(nf-1)*2*Pi/5,-(Pi/2-len_edge),vertex_glo(iim_glo,1,nf+nb_face/2)%xyz)
     CALL lonlat2xyz(rot+Pi/5+Pi/5+(nf-1)*2*Pi/5,-(-Pi/2+len_edge),vertex_glo(iim_glo,jjm_glo,nf+nb_face/2)%xyz)
   ENDDO

   DO nf=1,nb_face
     DO i=2,iim_glo-1
       CALL div_arc(vertex_glo(1,1,nf)%xyz,vertex_glo(iim_glo,1,nf)%xyz,(i-1)*1./(iim_glo-1),vertex_glo(i,1,nf)%xyz)
       CALL div_arc(vertex_glo(1,jjm_glo,nf)%xyz,vertex_glo(iim_glo,jjm_glo,nf)%xyz,                                 &
                               (i-1)*1./(iim_glo-1),vertex_glo(i,jjm_glo,nf)%xyz)
     ENDDO
      DO j=2,jjm_glo-1
       CALL div_arc(vertex_glo(1,1,nf)%xyz,vertex_glo(1,jjm_glo,nf)%xyz,(j-1)*1./(jjm_glo-1),vertex_glo(1,j,nf)%xyz)
       CALL div_arc(vertex_glo(iim_glo,1,nf)%xyz,vertex_glo(iim_glo,jjm_glo,nf)%xyz,(j-1)*1./(jjm_glo-1),            &
                    vertex_glo(iim_glo,j,nf)%xyz)
     ENDDO
    
     DO j=2,jjm_glo-1
       DO i=2,iim_glo-1
         IF (i+j-1 > iim_glo) THEN
           CALL div_arc(vertex_glo(iim_glo,i+j-iim_glo,nf)%xyz,vertex_glo(i+j-jjm_glo,jjm_glo,nf)%xyz,    &
                        (iim_glo-i)*1./(2*iim_glo-(i+j)),vertex_glo(i,j,nf)%xyz) 
         ELSE
           CALL div_arc(vertex_glo(i+j-1,1,nf)%xyz,vertex_glo(1,i+j-1,nf)%xyz,(j-1)*1./(i+j-2),vertex_glo(i,j,nf)%xyz)
         ENDIF
       ENDDO
     ENDDO
   ENDDO

  CALL dist_cart(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,d1) 
  CALL dist_cart(vertex_glo(2,1,1)%xyz,vertex_glo(3,1,1)%xyz,d2) 
  CALL dist_cart(vertex_glo(3,1,1)%xyz,vertex_glo(4,1,1)%xyz,d3) 
  CALL div_arc(vertex_glo(1,1,1)%xyz,vertex_glo(3,1,1)%xyz,0.5,p1)

  CALL circumcenter(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,2,1)%xyz,p1)
!  CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,1,1)%xyz,p1)
!  CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,p1)
  CALL dist_cart(vertex_glo(1,1,1)%xyz,p1,d1) 
  CALL dist_cart(vertex_glo(2,1,1)%xyz,p1,d2) 
  CALL dist_cart(vertex_glo(1,2,1)%xyz,p1,d3) 

  END SUBROUTINE compute_face

  SUBROUTINE Compute_face_projection
  USE spherical_geom_mod
  IMPLICIT NONE
    
    REAL(rstd) :: len_edge
    REAL(rstd) :: rot=0.
    INTEGER :: nf,i,j
    REAL(rstd),DIMENSION(3) :: p1,p2,p3
    REAL(rstd) :: d1,d2,d3
  
   len_edge=acos(cos(Pi/5)*cos(2*Pi/5)/(sin(Pi/5)*sin(2*Pi/5)))
   
   DO nf=1,nb_face/2
     CALL lonlat2xyz(0.,Pi/2,vertex_glo(1,1,nf)%xyz)
     CALL lonlat2xyz(rot+(nf-1)*2*Pi/5,Pi/2-len_edge,vertex_glo(iim_glo,1,nf)%xyz)
     CALL lonlat2xyz(rot+2*Pi/5+(nf-1)*2*Pi/5,Pi/2-len_edge,vertex_glo(1,jjm_glo,nf)%xyz)
     CALL lonlat2xyz(rot+Pi/5+(nf-1)*2*Pi/5,-Pi/2+len_edge,vertex_glo(iim_glo,jjm_glo,nf)%xyz)

     CALL lonlat2xyz(0.,-Pi/2,vertex_glo(1,1,nf+nb_face/2)%xyz)
     CALL lonlat2xyz(rot+Pi/5+(nf-1)*2*Pi/5,-(Pi/2-len_edge),vertex_glo(1,jjm_glo,nf+nb_face/2)%xyz)
     CALL lonlat2xyz(rot+Pi/5+2*Pi/5+(nf-1)*2*Pi/5,-(Pi/2-len_edge),vertex_glo(iim_glo,1,nf+nb_face/2)%xyz)
     CALL lonlat2xyz(rot+Pi/5+Pi/5+(nf-1)*2*Pi/5,-(-Pi/2+len_edge),vertex_glo(iim_glo,jjm_glo,nf+nb_face/2)%xyz)
   ENDDO

   DO nf=1,nb_face
     DO i=2,iim_glo-1
       CALL div_arc_bis(vertex_glo(1,1,nf)%xyz,vertex_glo(iim_glo,1,nf)%xyz,(i-1)*1./(iim_glo-1),vertex_glo(i,1,nf)%xyz)
       CALL div_arc_bis(vertex_glo(1,jjm_glo,nf)%xyz,vertex_glo(iim_glo,jjm_glo,nf)%xyz,                                 &
                               (i-1)*1./(iim_glo-1),vertex_glo(i,jjm_glo,nf)%xyz)
     ENDDO
      DO j=2,jjm_glo-1
       CALL div_arc_bis(vertex_glo(1,1,nf)%xyz,vertex_glo(1,jjm_glo,nf)%xyz,(j-1)*1./(jjm_glo-1),vertex_glo(1,j,nf)%xyz)
       CALL div_arc_bis(vertex_glo(iim_glo,1,nf)%xyz,vertex_glo(iim_glo,jjm_glo,nf)%xyz,(j-1)*1./(jjm_glo-1),            &
                    vertex_glo(iim_glo,j,nf)%xyz)
     ENDDO
    
     DO j=2,jjm_glo-1
       DO i=2,iim_glo-1
         IF (i+j-1 > iim_glo) THEN
           CALL div_arc_bis(vertex_glo(iim_glo,i+j-iim_glo,nf)%xyz,vertex_glo(i+j-jjm_glo,jjm_glo,nf)%xyz,    &
                        (iim_glo-i)*1./(2*iim_glo-(i+j)),vertex_glo(i,j,nf)%xyz) 
         ELSE
           CALL div_arc_bis(vertex_glo(i+j-1,1,nf)%xyz,vertex_glo(1,i+j-1,nf)%xyz,(j-1)*1./(i+j-2),vertex_glo(i,j,nf)%xyz)
         ENDIF
       ENDDO
     ENDDO
     
!     DO j=1,jjm_glo
!       DO i=1,iim_glo
!         vertex_glo(i,j,nf)%xyz=vertex_glo(i,j,nf)%xyz/sqrt(sum(vertex_glo(i,j,nf)%xyz**2))
!       ENDDO
!     ENDDO
   ENDDO

!  CALL dist_cart(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,d1) 
!  CALL dist_cart(vertex_glo(2,1,1)%xyz,vertex_glo(3,1,1)%xyz,d2) 
!  CALL dist_cart(vertex_glo(3,1,1)%xyz,vertex_glo(4,1,1)%xyz,d3) 
!  CALL div_arc(vertex_glo(1,1,1)%xyz,vertex_glo(3,1,1)%xyz,0.5,p1)
!  CALL div_arc(vertex_glo(2,1,1)%xyz,vertex_glo(1,3,1)%xyz,1./3,p1)
!  CALL div_arc(vertex_glo(1,2,1)%xyz,vertex_glo(3,1,1)%xyz,1./3,p2)
!  CALL div_arc(vertex_glo(2,2,1)%xyz,vertex_glo(1,1,1)%xyz,1./3,p3)
!  PRINT *, "dist",d1
!  PRINT *, "dist",d2
!  PRINT *, "dist",d3
!  PRINT *,"dist",vertex_glo(2,1,1)%xyz
!  PRINT *,"dist",p1/sqrt(sum(p1**2))
!  CALL Centroide(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,2,1)%xyz,p1)
!  CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,1,1)%xyz,p1)
!  CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,p1)
!  CALL dist_cart(vertex_glo(1,1,1)%xyz,p1,d1) 
!  CALL dist_cart(vertex_glo(2,1,1)%xyz,p1,d2) 
!  CALL dist_cart(vertex_glo(1,2,1)%xyz,p1,d3) 
!  PRINT *, "dist",d1
!  PRINT *, "dist",d2
!  PRINT *, "dist",d3

  END SUBROUTINE compute_face_projection


  SUBROUTINE compute_extended_face_bis
  IMPLICIT NONE
  
    INTEGER :: nf
    INTEGER :: i,j
    INTEGER :: delta
    INTEGER :: ind
    INTEGER :: k,ng

    DO nf=1,nb_face
      DO i=2,iim_glo-1

        CALL set_neighbour(i,1,nf,ldown-1)
        CALL set_neighbour(i,1,nf,rdown-1)

        CALL set_neighbour(i,jjm_glo,nf,lup-1)
        CALL set_neighbour(i,jjm_glo,nf,rup-1)
      
      ENDDO    
     
      DO j=2,jjm_glo-1
        
        CALL set_neighbour(1,j,nf,left-1)
        CALL set_neighbour(1,j,nf,lup-1)

        CALL set_neighbour(iim_glo,j,nf,rdown-1)
        CALL set_neighbour(iim_glo,j,nf,right-1)

      ENDDO
      
      CALL set_neighbour(2,0,nf,left-1)
      CALL set_neighbour(iim_glo,0,nf,right-1)
      CALL set_neighbour(iim_glo+1,jjm_glo-1,nf,rup-1)
      CALL set_neighbour(iim_glo-1,jjm_glo+1,nf,right-1)
      CALL set_neighbour(1,jjm_glo+1,nf,left-1)
      CALL set_neighbour(0,2,nf,ldown-1)
      
      CALL set_neighbour(1,0,nf,left-1)
      CALL set_neighbour(iim_glo,0,nf,right-1)
      CALL set_neighbour(0,jjm_glo,nf,rup-1)
      CALL set_neighbour(iim_glo+1,jjm_glo,nf,rup-1)
    ENDDO

    DO nf=1,nb_face
      DO j=0,jjm_glo+1
        DO i=0,iim_glo+1
          ind=vertex_glo(i,j,nf)%ind
          delta=vertex_glo(i,j,nf)%delta
          DO k=0,5
            ng=MOD(delta+k+6,6)
            vertex_glo(i,j,nf)%neighbour(k)=cell_glo(ind)%neighbour(ng)
          ENDDO
        ENDDO
      ENDDO
 
      i=1 ; j=1
      vertex_glo(i,j,nf)%neighbour(0)=vertex_glo(i+1,j,nf)%ind
      vertex_glo(i,j,nf)%neighbour(1)=vertex_glo(i,j+1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(2)=vertex_glo(i-1,j+1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(3)=vertex_glo(i-1,j,nf)%ind
      vertex_glo(i,j,nf)%neighbour(4)=vertex_glo(i,j-1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(5)=vertex_glo(i+1,j-1,nf)%ind

      i=iim_glo ; j=1
      vertex_glo(i,j,nf)%neighbour(0)=vertex_glo(i+1,j,nf)%ind
      vertex_glo(i,j,nf)%neighbour(1)=vertex_glo(i,j+1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(2)=vertex_glo(i-1,j+1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(3)=vertex_glo(i-1,j,nf)%ind
      vertex_glo(i,j,nf)%neighbour(4)=vertex_glo(i,j-1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(5)=vertex_glo(i+1,j-1,nf)%ind

      i=1 ; j=jjm_glo
      vertex_glo(i,j,nf)%neighbour(0)=vertex_glo(i+1,j,nf)%ind
      vertex_glo(i,j,nf)%neighbour(1)=vertex_glo(i,j+1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(2)=vertex_glo(i-1,j+1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(3)=vertex_glo(i-1,j,nf)%ind
      vertex_glo(i,j,nf)%neighbour(4)=vertex_glo(i,j-1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(5)=vertex_glo(i+1,j-1,nf)%ind

      i=iim_glo ; j=jjm_glo
      vertex_glo(i,j,nf)%neighbour(0)=vertex_glo(i+1,j,nf)%ind
      vertex_glo(i,j,nf)%neighbour(1)=vertex_glo(i,j+1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(2)=vertex_glo(i-1,j+1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(3)=vertex_glo(i-1,j,nf)%ind
      vertex_glo(i,j,nf)%neighbour(4)=vertex_glo(i,j-1,nf)%ind
      vertex_glo(i,j,nf)%neighbour(5)=vertex_glo(i+1,j-1,nf)%ind
   ENDDO      
  
  CONTAINS
  
    SUBROUTINE set_neighbour(i,j,nf,neighbour)
    IMPLICIT NONE
      INTEGER :: i,j,nf,neighbour
      INTEGER :: in,jn
      INTEGER :: k,ng
      INTEGER :: delta
      INTEGER :: ind,ind2
    
       SELECT CASE (neighbour)
         CASE(right-1)
           in=i+1 ; jn=j
         CASE(rup-1)
           in=i ; jn=j+1
         CASE(lup-1)
           in=i-1 ; jn=j+1
         CASE(left-1)
           in=i-1 ; jn=j
         CASE(ldown-1)
           in=i ; jn=j-1
         CASE(rdown-1)
           in=i+1; jn=j-1
        END SELECT
         
        ind=vertex_glo(i,j,nf)%ind
        delta=MOD(vertex_glo(i,j,nf)%delta+neighbour+6,6)
        ind2=cell_glo(ind)%neighbour(delta)
        DO k=0,5
          IF (cell_glo(ind2)%neighbour(k)==ind) ng=k
        ENDDO
        vertex_glo(in,jn,nf)%ind=ind2
        vertex_glo(in,jn,nf)%delta=MOD(ng-neighbour+3+6,6)    
      END SUBROUTINE set_neighbour    
  
  END SUBROUTINE compute_extended_face_bis
    
    
  SUBROUTINE compute_extended_face
  IMPLICIT NONE
  
    INTEGER :: nf,nf2
    INTEGER :: ind,ind2
    INTEGER :: i,j,h
    

    DO h=1,halo
      DO nf=1,nb_face
        DO i=1-h+1,iim_glo+h-1
          j=1-h+1
          ind=vertex_glo(i,j,nf)%ind
          nf2=cell_glo(ind)%assign_face
          ind2=cell_glo(ind)%neighbour(tab_index(nf,nf2,ldown-1))
          vertex_glo(i,j-1,nf)%ind=ind2
          vertex_glo(i,j-1,nf)%xyz=cell_glo(ind2)%xyz
        ENDDO
      
        DO j=1-h,jjm_glo+h-1
          i=iim_glo+h-1
          ind=vertex_glo(i,j,nf)%ind
          nf2=cell_glo(ind)%assign_face
          ind2=cell_glo(ind)%neighbour(tab_index(nf,nf2,right-1))
          vertex_glo(i+1,j,nf)%ind=ind2
          vertex_glo(i+1,j,nf)%xyz=cell_glo(ind2)%xyz
        ENDDO
      
        DO i=iim_glo+h,1-h+1,-1
          j=jjm_glo+h-1
          ind=vertex_glo(i,j,nf)%ind
          nf2=cell_glo(ind)%assign_face
          ind2=cell_glo(ind)%neighbour(tab_index(nf,nf2,rup-1))
          vertex_glo(i,j+1,nf)%ind=ind2
          vertex_glo(i,j+1,nf)%xyz=cell_glo(ind2)%xyz
        ENDDO

         DO j=jjm_glo+h,1-h,-1
          i=1-h+1
          ind=vertex_glo(i,j,nf)%ind
          nf2=cell_glo(ind)%assign_face
          ind2=cell_glo(ind)%neighbour(tab_index(nf,nf2,left-1))
          vertex_glo(i-1,j,nf)%ind=ind2
          vertex_glo(i-1,j,nf)%xyz=cell_glo(ind2)%xyz
        ENDDO
      ENDDO
    ENDDO
    
    
  END SUBROUTINE compute_extended_face
  
  SUBROUTINE Set_index
  IMPLICIT NONE
    INTEGER :: delta
    INTEGER :: n1,n2,i
   
    DO n1=1,5
     DO delta=-2,2
        DO i=0,5
          n2=MOD(n1-1-delta+5,5)+1
          tab_index(n1,n2,i)=MOD(delta+i+6,6)
          n2=MOD(n1-1+delta+5,5)+1
          tab_index(n1+5,n2+5,i)=MOD(delta+i+6,6)
        ENDDO
      ENDDO
    ENDDO
    
    DO n1=1,5
      DO i=0,5
        delta=3

        n2=5+n1
        Tab_index(n1,n2,i)=MOD(delta+i+6,6)
        Tab_index(n2,n1,i)=MOD(delta+i+6,6)
      
        n2=5+n1-1
        IF (n2==5) n2=10
        Tab_index(n1,n2,i)=MOD(delta+i+6,6)
        Tab_index(n2,n1,i)=MOD(delta+i+6,6)
      ENDDO
    ENDDO

  END SUBROUTINE set_index
  
  
  SUBROUTINE set_cell
  IMPLICIT NONE
    INTEGER :: ind,ind2
    INTEGER :: nf,nf2,nfm1,nfp1
    INTEGER :: i,j,k
    INTEGER :: delta
    
    ind=0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Association des cellules !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! interieur des faces
    DO nf=1,nb_face     
      DO i=2,iim_glo-1
        DO j=2,jjm_glo-1
          ind=ind+1
          vertex_glo(i,j,nf)%ind=ind
          cell_glo(ind)%assign_face=nf
          cell_glo(ind)%assign_i=i
          cell_glo(ind)%assign_j=j
        ENDDO
      ENDDO
    ENDDO

! frontieres faces nord
    DO nf=1,nb_face/2

      nfm1=nf-1
      IF (nfm1==0) nfm1=nb_face/2
      
      nf2=nf+nb_face/2

      DO i=2,iim_glo-1
        ind=ind+1
        vertex_glo(i,1,nf)%ind=ind
        vertex_glo(1,i,nfm1)%ind=ind
        cell_glo(ind)%assign_face=nf
        cell_glo(ind)%assign_i=i
        cell_glo(ind)%assign_j=1
      ENDDO

      DO i=2,iim_glo-1
        ind=ind+1
        vertex_glo(i,jjm_glo,nf)%ind=ind
        vertex_glo(iim_glo-i+1,jjm_glo,nf2)%ind=ind
        
        cell_glo(ind)%assign_face=nf
        cell_glo(ind)%assign_i=i
        cell_glo(ind)%assign_j=jjm_glo
      ENDDO
     ENDDO

! frontieres faces sud

    DO nf=nb_face/2+1,nb_face
      
      nfp1=nf+1
      IF (nfp1==nb_face+1) nfp1=nb_face/2+1
      
      nf2=nf-nb_face/2+1
      IF (nf2==nb_face/2+1) nf2=1
      
      DO i=2,iim_glo-1
        ind=ind+1
        vertex_glo(i,1,nf)%ind=ind
        vertex_glo(1,i,nfp1)%ind=ind
        cell_glo(ind)%assign_face=nf
        cell_glo(ind)%assign_i=i
        cell_glo(ind)%assign_j=1
      ENDDO

      DO j=2,jjm_glo-1
        ind=ind+1
        vertex_glo(iim_glo,j,nf)%ind=ind
        vertex_glo(iim_glo,jjm_glo-j+1,nf2)%ind=ind
        
        cell_glo(ind)%assign_face=nf
        cell_glo(ind)%assign_i=iim_glo
        cell_glo(ind)%assign_j=j
      ENDDO
    ENDDO

! vertex nord (sur 3 faces)

    DO nf=1,nb_face/2

      nfp1=nf+1
      IF (nfp1==nb_face/2+1) nfp1=1
      
      nf2=nf+nb_face/2

      ind=ind+1
      vertex_glo(1,jjm_glo,nf)%ind=ind
      vertex_glo(iim_glo,1,nfp1)%ind=ind
      vertex_glo(iim_glo,jjm_glo,nf2)%ind=ind
      cell_glo(ind)%assign_face=nf
      cell_glo(ind)%assign_i=1
      cell_glo(ind)%assign_j=jjm_glo
     ENDDO

! vertex sud (sur 3 faces)

    DO nf=nb_face/2+1,nb_face

      nfp1=nf+1
      IF (nfp1==nb_face+1) nfp1=nb_face/2+1
      
      nf2=nf-nb_face/2+1
      IF (nf2==nb_face/2+1) nf2=1

      ind=ind+1
      vertex_glo(iim_glo,1,nf)%ind=ind
      vertex_glo(1,jjm_glo,nfp1)%ind=ind
      vertex_glo(iim_glo,jjm_glo,nf2)%ind=ind
      cell_glo(ind)%assign_face=nf
      cell_glo(ind)%assign_i=iim_glo
      cell_glo(ind)%assign_j=1
    ENDDO


! vertex nord (sur 5 faces)

    ind=ind+1
    vertex_glo(1,1,1)%ind=ind
    vertex_glo(1,1,2)%ind=ind
    vertex_glo(1,1,3)%ind=ind
    vertex_glo(1,1,4)%ind=ind
    vertex_glo(1,1,5)%ind=ind
    cell_glo(ind)%assign_face=1
    cell_glo(ind)%assign_i=1
    cell_glo(ind)%assign_j=1
      
! vertex sud (sur 5 faces)

    ind=ind+1    
    vertex_glo(1,1,6)%ind=ind
    vertex_glo(1,1,7)%ind=ind
    vertex_glo(1,1,8)%ind=ind
    vertex_glo(1,1,9)%ind=ind
    vertex_glo(1,1,10)%ind=ind
    cell_glo(ind)%assign_face=6
    cell_glo(ind)%assign_i=1
    cell_glo(ind)%assign_j=1

! assignation des coordonnees xyz

  DO ind=1,ncell_glo
    nf=cell_glo(ind)%assign_face
    i=cell_glo(ind)%assign_i
    j=cell_glo(ind)%assign_j
    cell_glo(ind)%xyz=vertex_glo(i,j,nf)%xyz
  ENDDO

!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Association des voisins !
!!!!!!!!!!!!!!!!!!!!!!!!!!!

! interieur des faces
  
    DO nf=1,nb_face
      DO j=2,jjm_glo-1
        DO i=2,iim_glo-1
          ind=vertex_glo(i,j,nf)%ind
! right          
          ind2=vertex_glo(i+1,j,nf)%ind   
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind)%neighbour(right-1)=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,left-1))=ind
          
! rup          
          ind2=vertex_glo(i,j+1,nf)%ind   
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind)%neighbour(rup-1)=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,ldown-1))=ind
              
! lup          
          ind2=vertex_glo(i-1,j+1,nf)%ind  
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind)%neighbour(lup-1)=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,rdown-1))=ind
          
! left          
          ind2=vertex_glo(i-1,j,nf)%ind   
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind)%neighbour(left-1)=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,right-1))=ind
! ldown          
          ind2=vertex_glo(i,j-1,nf)%ind   
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind)%neighbour(ldown-1)=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,rup-1))=ind
          
! rdown          
          ind2=vertex_glo(i+1,j-1,nf)%ind   
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind)%neighbour(rdown-1)=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,lup-1))=ind
        ENDDO
      ENDDO
    ENDDO

! frontiere face nord
    
    DO nf=1,nb_face/2
      DO i=2,iim_glo-1
        j=1
        ind=vertex_glo(i,j,nf)%ind
! right          
        ind2=vertex_glo(i+1,j,nf)%ind   
        nf2=cell_glo(ind2)%assign_face
        cell_glo(ind)%neighbour(right-1)=ind2
        cell_glo(ind2)%neighbour(tab_index(nf,nf2,left-1))=ind
! left          
        ind2=vertex_glo(i-1,j,nf)%ind   
        nf2=cell_glo(ind2)%assign_face
        cell_glo(ind)%neighbour(left-1)=ind2
        cell_glo(ind2)%neighbour(tab_index(nf,nf2,right-1))=ind
! lup          
          ind2=vertex_glo(i-1,j+1,nf)%ind  
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind)%neighbour(lup-1)=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,rdown-1))=ind

      ENDDO

      DO i=2,iim_glo-1
        j=jjm_glo
        ind=vertex_glo(i,j,nf)%ind
! right          
        ind2=vertex_glo(i+1,j,nf)%ind   
        nf2=cell_glo(ind2)%assign_face
        cell_glo(ind)%neighbour(right-1)=ind2
        cell_glo(ind2)%neighbour(tab_index(nf,nf2,left-1))=ind
! left          
        ind2=vertex_glo(i-1,j,nf)%ind   
        nf2=cell_glo(ind2)%assign_face
        cell_glo(ind)%neighbour(left-1)=ind2
        cell_glo(ind2)%neighbour(tab_index(nf,nf2,right-1))=ind
! rdown          
          ind2=vertex_glo(i+1,j-1,nf)%ind   
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind)%neighbour(rdown-1)=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,lup-1))=ind
      ENDDO
    ENDDO
    

! frontiere face sud

    DO nf=nb_face/2+1,nb_face
      DO i=2,iim_glo-1
        j=1
        ind=vertex_glo(i,j,nf)%ind
! right          
        ind2=vertex_glo(i+1,j,nf)%ind   
        nf2=cell_glo(ind2)%assign_face
        cell_glo(ind)%neighbour(right-1)=ind2
        cell_glo(ind2)%neighbour(tab_index(nf,nf2,left-1))=ind
! left          
        ind2=vertex_glo(i-1,j,nf)%ind   
        nf2=cell_glo(ind2)%assign_face
        cell_glo(ind)%neighbour(left-1)=ind2
        cell_glo(ind2)%neighbour(tab_index(nf,nf2,right-1))=ind
! lup          
          ind2=vertex_glo(i-1,j+1,nf)%ind  
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind)%neighbour(lup-1)=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,rdown-1))=ind
      ENDDO

      DO j=2,jjm_glo-1
        i=iim_glo
        ind=vertex_glo(i,j,nf)%ind
! rup          
        ind2=vertex_glo(i,j+1,nf)%ind   
        nf2=cell_glo(ind2)%assign_face
        cell_glo(ind)%neighbour(rup-1)=ind2
        cell_glo(ind2)%neighbour(tab_index(nf,nf2,ldown-1))=ind
! ldown          
        ind2=vertex_glo(i,j-1,nf)%ind   
        nf2=cell_glo(ind2)%assign_face
        cell_glo(ind)%neighbour(ldown-1)=ind2
        cell_glo(ind2)%neighbour(tab_index(nf,nf2,rup-1))=ind
! lup          
          ind2=vertex_glo(i-1,j+1,nf)%ind  
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind)%neighbour(lup-1)=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,rdown-1))=ind

      ENDDO
    ENDDO          

   
! vertex nord (sur 3 faces)

    DO nf=1,nb_face/2
      ind=vertex_glo(1,jjm_glo,nf)%ind
      cell_glo(ind)%neighbour(2)=cell_glo(ind)%neighbour(1)
    ENDDO

! vertex sud (sur 3 faces)

    DO nf=nb_face/2+1,nb_face

      ind=vertex_glo(iim_glo,1,nf)%ind
      cell_glo(ind)%neighbour(5)=cell_glo(ind)%neighbour(0)
    ENDDO


! vertex nord (sur 5 faces)

    ind=vertex_glo(1,1,1)%ind
    cell_glo(ind)%neighbour(3)=cell_glo(ind)%neighbour(4)
      
! vertex sud (sur 5 faces)

    ind=vertex_glo(1,1,6)%ind
    cell_glo(ind)%neighbour(3)=cell_glo(ind)%neighbour(2)
      
    
 !! assignation des delta 
    DO nf=1,nb_face
      DO j=1,jjm_glo
        DO i=1,iim_glo 
          ind=vertex_glo(i,j,nf)%ind
          nf2=cell_glo(ind)%assign_face
          vertex_glo(i,j,nf)%delta=tab_index(nf,nf2,0)
        ENDDO
      ENDDO
    ENDDO    

     
  END SUBROUTINE set_cell
      
          
  SUBROUTINE set_cell_edge
  IMPLICIT NONE
    INTEGER :: ind,k,ind2,ng,ne
    
    DO ind=1,ncell_glo
      cell_glo(ind)%edge(:)=0
    ENDDO
    
    ne=0
    DO ind=1,ncell_glo
      DO k=0,5
        IF (cell_glo(ind)%edge(k)==0) THEN
          ind2=cell_glo(ind)%neighbour(k)
          DO ng=0,5
            IF (cell_glo(ind2)%neighbour(ng)==ind) THEN
              IF (cell_glo(ind2)%edge(ng)==0) THEN
                ne=ne+1
                cell_glo(ind)%edge(k)=ne
                edge_glo(ne)%e1=ind
                edge_glo(ne)%e2=ind2
                cell_glo(ind2)%edge(ng)=ne
              ELSE 
                ne=cell_glo(ind2)%edge(ng)
                cell_glo(ind)%edge(k)=ne
              ENDIF    
              EXIT            
            ENDIF
          ENDDO
        ENDIF
      ENDDO
    ENDDO
    
  END SUBROUTINE  set_cell_edge    
        
  SUBROUTINE set_vertex_edge
  IMPLICIT NONE
    INTEGER :: nf
    INTEGER :: i,j
    INTEGER :: ind,ind2
    INTEGER :: ne
    INTEGER :: k,l     
  
    DO nf=1,nb_face
      DO j=0,jjm_glo+1
        DO i=0,iim_glo+1
          ind=vertex_glo(i,j,nf)%ind
          DO k=0,5
            ind2=vertex_glo(i,j,nf)%neighbour(k)
            DO l=0,5
              ne=cell_glo(ind)%edge(l)
              IF (edge_glo(ne)%e1==ind .AND. edge_glo(ne)%e2==ind2) THEN
                vertex_glo(i,j,nf)%ne(k)=1
              ELSE IF (edge_glo(ne)%e1==ind2 .AND. edge_glo(ne)%e2==ind) THEN
                vertex_glo(i,j,nf)%ne(k)=-1
              ENDIF
            ENDDO
          ENDDO
        ENDDO
      ENDDO
    ENDDO
  END SUBROUTINE set_vertex_edge 
     
  SUBROUTINE compute_metric
  IMPLICIT NONE
  
    CALL allocate_metric
!    CALL compute_face
    CALL compute_face_projection
    CALL remap_schmidt

    CALL set_index
    CALL set_cell
    CALL compute_extended_face_bis
    CALL set_cell_edge
    CALL set_vertex_edge
    
  END SUBROUTINE compute_metric

END MODULE metric