MODULE metric_icosa_area_old_mod
  USE genmod
  USE grid_param, ONLY : iim_glo, jjm_glo
  USE domain_param
  USE metric_var_mod  
  

  
CONTAINS

  SUBROUTINE remap_schmidt_glo ! FIXME : not called any more
    USE spherical_geom_mod
    USE getin_mod
    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_glo

  SUBROUTINE allocate_metric
  IMPLICIT NONE
    INTEGER :: ind
    
    nb_face=3
    ALLOCATE(vertex_glo(0:iim_glo+3,-2:jjm_glo+1,nb_face))
    ncell_glo=3*(iim_glo+3)*(iim_glo+3)-3*(iim_glo+2)-2
    ALLOCATE(cell_glo(ncell_glo))
    ALLOCATE(tab_index(nb_face,nb_face,0:5))
    ALLOCATE(edge_glo(ncell_glo*3))
    ALLOCATE(vertices_glo(ncell_glo*2))
    
    DO ind=1,ncell_glo
      cell_glo(ind)%neighbour(:)=0
    ENDDO
    
  END SUBROUTINE allocate_metric
  

  SUBROUTINE compute_face_area
  USE spherical_geom_mod
  USE vector
  USE getin_mod
  use earth_const
  IMPLICIT NONE
    REAL(rstd) :: lonA,lonB,lonC,lonD,lonE,lonF,lonO, testLonO, testLatO
    REAL(rstd) :: latA,latB,latC,latD,latE,latF,latO
    REAL(rstd) :: width,length,d, theta
    REAL(rstd) :: A_xyz(3), B_xyz(3), C_xyz(3), D_xyz(3), E_xyz(3), F_xyz(3), O_xyz(3)
    REAL(rstd) :: A_xyz1(3), B_xyz1(3), C_xyz1(3), D_xyz1(3), E_xyz1(3), F_xyz1(3), O_xyz1(3)
    REAL(rstd) :: A_xyz2(3), B_xyz2(3), C_xyz2(3), D_xyz2(3), E_xyz2(3), F_xyz2(3), O_xyz2(3)
    REAL(rstd) :: A_xyz3(3), B_xyz3(3), C_xyz3(3), D_xyz3(3), E_xyz3(3), F_xyz3(3), O_xyz3(3)
    INTEGER :: i,j, nf
    
    ! A point of the rectangle area
    CALL getin("area_center_lon", lonO)
    CALL getin("area_center_lat", latO)
    CALL getin("area_radius", length)
    lonO=lonO*Pi/180. ; latO=latO*Pi/180 
    
    TestLonO =0
    TestLatO = Pi/2
    lonA = 0 ;
    latA = Pi/2-length*1000/radius 
    lonB = lonA + Pi/3
    latB = latA   
    lonC = lonA + 2*Pi/3 
    latC = latA 
    lonD = lonA + Pi
    latD = latA
    lonE = lonA - 2*Pi/3
    latE = latA 
    lonF = lonA - Pi/3
    latF = latA   
                                                                            
    
    CALL lonlat2xyz(lonA, latA, A_xyz2) 
    CALL lonlat2xyz(lonB, latB, B_xyz2) 
    CALL lonlat2xyz(lonC, latC, C_xyz2) 
    CALL lonlat2xyz(lonD, latD, D_xyz2) 
    CALL lonlat2xyz(lonE, latE, E_xyz2) 
    CALL lonlat2xyz(lonF, latF, F_xyz2) 
    CALL lonlat2xyz(testLonO, testLatO, O_xyz2) 

    CALL rotate_Oy(A_xyz2,Pi/2-latO,A_xyz1)
    CALL rotate_Oy(B_xyz2,Pi/2-latO,B_xyz1)
    CALL rotate_Oy(C_xyz2,Pi/2-latO,C_xyz1)
    CALL rotate_Oy(D_xyz2,Pi/2-latO,D_xyz1)
    CALL rotate_Oy(E_xyz2,Pi/2-latO,E_xyz1)
    CALL rotate_Oy(F_xyz2,Pi/2-latO,F_xyz1)
    CALL rotate_Oy(O_xyz2,Pi/2-latO,O_xyz1)

    CALL rotate_Oz(A_xyz1,lonO,A_xyz)
    CALL rotate_Oz(B_xyz1,lonO,B_xyz)
    CALL rotate_Oz(C_xyz1,lonO,C_xyz)
    CALL rotate_Oz(D_xyz1,lonO,D_xyz)
    CALL rotate_Oz(E_xyz1,lonO,E_xyz)
    CALL rotate_Oz(F_xyz1,lonO,F_xyz)
    CALL rotate_Oz(O_xyz1,lonO,O_xyz)
  
    !face 1
    vertex_glo(1,-2,1)%xyz                   =  A_xyz
    vertex_glo(iim_glo+3,-2, 1)%xyz          =  B_xyz
    vertex_glo(iim_glo+3,jjm_glo,1)%xyz      =  C_xyz
    vertex_glo(1,jjm_glo,1)%xyz              =  O_xyz
    !face 2
    vertex_glo(1,-2,2)%xyz                   =  C_xyz
    vertex_glo(iim_glo+3,-2, 2)%xyz          =  D_xyz
    vertex_glo(iim_glo+3,jjm_glo,2)%xyz      =  E_xyz
    vertex_glo(1,jjm_glo,2)%xyz              =  O_xyz
    !face 3
    vertex_glo(1,-2,3)%xyz                   =  E_xyz
    vertex_glo(iim_glo+3,-2, 3)%xyz          =  F_xyz
    vertex_glo(iim_glo+3,jjm_glo,3)%xyz      =  A_xyz
    vertex_glo(1,jjm_glo,3)%xyz              =  O_xyz
    
     
    DO nf=1,nb_face
      DO i=2,iim_glo+2
       !AB
        CALL div_arc_bis(vertex_glo(1,-2,nf)%xyz, vertex_glo(iim_glo+3,-2,nf)%xyz,(i-1)*1./(iim_glo+2),vertex_glo(i,-2,nf)%xyz)
       !OC
        CALL div_arc_bis(vertex_glo(1,jjm_glo,nf)%xyz, vertex_glo(iim_glo+3,jjm_glo,nf)%xyz,(i-1)*1./(iim_glo+2),vertex_glo(i,jjm_glo,nf)%xyz)

      ENDDO
    
      DO j=-1,jjm_glo-1
       !BC
       CALL div_arc_bis(vertex_glo(iim_glo+3,-2,nf)%xyz,vertex_glo(iim_glo+3,jjm_glo,nf)%xyz,(j+2)*1./(jjm_glo+2),vertex_glo(iim_glo+3,j,nf)%xyz)
       !AO      
       CALL div_arc_bis(vertex_glo(1,-2,nf)%xyz,vertex_glo(1,jjm_glo,nf)%xyz,(j+2)*1./(jjm_glo+2),vertex_glo(1,j,nf)%xyz)
      ENDDO

      DO j=-1,jjm_glo-1
        DO i=2,iim_glo+2
          IF ( i+j > iim_glo+1) THEN
            CALL div_arc_bis(vertex_glo(iim_glo+3,i+j-(iim_glo+3),nf)%xyz,vertex_glo(i+j-jjm_glo,jjm_glo,nf)%xyz,    &
                             (iim_glo+3-i)*1./(2*iim_glo+3-(i+j)),vertex_glo(i,j,nf)%xyz) 
          ELSE
            CALL div_arc_bis(vertex_glo(i+j+2,-2,nf)%xyz,vertex_glo(1,i+j-1,nf)%xyz,(j+2)*1./(i+j+1),vertex_glo(i,j,nf)%xyz)
          ENDIF
        ENDDO
      ENDDO
   
    ENDDO    
   
  END SUBROUTINE compute_face_area
 
  SUBROUTINE Set_index_area
  IMPLICIT NONE
    INTEGER :: delta
    INTEGER :: nf1,nf2,i
   
    DO nf1=1,nb_face
      DO nf2=1,nb_face
        delta=-2*(nf2-nf1)
        DO i=0,5
         Tab_index(nf1,nf2,i) = MOD(delta+i+6,6)
        ENDDO
      ENDDO
    ENDDO

  END SUBROUTINE set_index_area
  
  SUBROUTINE set_cell_area
  IMPLICIT NONE
    INTEGER :: ind,ind2
    INTEGER :: nf,nf2,nfm1,nfp1
    INTEGER :: i,j
    
    ind=0

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

! interieur des faces
    DO nf=1,nb_face     
      DO j=-2,jjm_glo-1
        DO i=2,iim_glo+3
          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 intrieurs (sauf vertex central)

!   face 1->2
   DO i=2,iim_glo+3   
     ind=ind+1
     vertex_glo(i,jjm_glo,1)%ind=ind
     vertex_glo(1,jjm_glo+1-i,2)%ind=ind
     cell_glo(ind)%assign_face=1
     cell_glo(ind)%assign_i=i
     cell_glo(ind)%assign_j=jjm_glo
   ENDDO
   
!   face 2->3
   DO i=2,iim_glo+3   
     ind=ind+1
     vertex_glo(i,jjm_glo,2)%ind=ind
     vertex_glo(1,jjm_glo+1-i,3)%ind=ind
     cell_glo(ind)%assign_face=2
     cell_glo(ind)%assign_i=i
     cell_glo(ind)%assign_j=jjm_glo
   ENDDO
     
!   face 3->1
   DO i=2,iim_glo+3   
     ind=ind+1
     vertex_glo(i,jjm_glo,3)%ind=ind
     vertex_glo(1,jjm_glo+1-i,1)%ind=ind
     cell_glo(ind)%assign_face=3
     cell_glo(ind)%assign_i=i
     cell_glo(ind)%assign_j=jjm_glo
   ENDDO

! vertex central
    ind=ind+1
    vertex_glo(1,jjm_glo,1)%ind = ind  
    vertex_glo(1,jjm_glo,2)%ind = ind  
    vertex_glo(1,jjm_glo,3)%ind = ind  
    cell_glo(ind)%assign_face=1
    cell_glo(ind)%assign_i=1
    cell_glo(ind)%assign_j=jjm_glo

! cellule externe
   vertex_glo(0,-2,1)%ind =-1
   vertex_glo(0,-2,2)%ind =-1
   vertex_glo(0,-2,3)%ind =-1
   vertex_glo(iim_glo+3,jjm_glo+1,1)%ind =-1
   vertex_glo(iim_glo+3,jjm_glo+1,2)%ind =-1
   vertex_glo(iim_glo+3,jjm_glo+1,3)%ind =-1

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

!interieur des faces

    DO nf=1,nb_face     
      DO j=-2,jjm_glo-1
        DO i=2,iim_glo+3
          ind=vertex_glo(i,j,nf)%ind
! right   
          IF (i==iim_glo+3) THEN
            cell_glo(ind)%neighbour(right-1)=-1
          ELSE    
            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
          ENDIF
! 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          
          IF (j==-2) THEN
           cell_glo(ind)%neighbour(ldown-1)=-1
          ELSE
            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
          ENDIF 
! rdown          
          IF (i==iim_glo+3 .OR. j==-2) THEN
            cell_glo(ind)%neighbour(rdown-1)=-1
          ELSE
            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
          ENDIF
        ENDDO
      ENDDO
    ENDDO

! frontiere des faces

    DO nf=1,nb_face     
      DO i=2,iim_glo+2
        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
! 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
      ind=vertex_glo(iim_glo+3,jjm_glo,nf)%ind
      cell_glo(ind)%neighbour(right-1)=-1
      cell_glo(ind)%neighbour(rup-1)=-1
      cell_glo(ind)%neighbour(rdown-1)=-1
    ENDDO

 !! assignation des delta 
    DO nf=1,nb_face
      DO j=-2,jjm_glo
        DO i=1,iim_glo+3 
          ind=vertex_glo(i,j,nf)%ind
          IF (ind/=-1) THEN
            nf2=cell_glo(ind)%assign_face
            vertex_glo(i,j,nf)%delta=tab_index(nf,nf2,0)
          ENDIF
        ENDDO
      ENDDO
    ENDDO    

! 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
    
  END SUBROUTINE set_cell_area

 
  SUBROUTINE compute_extended_face_area_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+2
        CALL set_neighbour(i,jjm_glo,nf,lup-1)
        CALL set_neighbour(i,jjm_glo,nf,rup-1)
      ENDDO    
     
      DO j=-1,jjm_glo
        
        CALL set_neighbour(1,j,nf,left-1)
        CALL set_neighbour(1,j,nf,lup-1)

      ENDDO
      
    ENDDO
    
    DO nf=1,nb_face
      DO j=-2,jjm_glo+1
        DO i=0,iim_glo+3
          ind=vertex_glo(i,j,nf)%ind
          IF (ind/=-1) THEN
            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
          ENDIF
        ENDDO
      ENDDO
   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_area_bis



  SUBROUTINE set_cell_edge_area
  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)
          IF (ind2==-1) THEN
            cell_glo(ind)%edge(k)=-1
          ELSE
            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
        ENDIF
      ENDDO
    ENDDO
    
  END SUBROUTINE  set_cell_edge_area    
        
  SUBROUTINE set_cell_vertex_area
  IMPLICIT NONE
    INTEGER :: k,k2
    INTEGER :: ind,ind1,ind2
    INTEGER :: ng1,ng2
    INTEGER :: ne     
  
    DO ind=1,ncell_glo
      cell_glo(ind)%vertex(:)=-1
    ENDDO
    
    ne=0
    DO ind=1,ncell_glo
      DO k=0,5
        IF (cell_glo(ind)%vertex(k)==-1) THEN
          ind1=cell_glo(ind)%neighbour(k)
          IF (ind1/=-1) THEN
            DO ng1=0,5
              ind2=cell_glo(ind1)%neighbour(ng1)
              IF (ind2/=-1) THEN
                DO ng2=0,5
                  IF (cell_glo(ind2)%neighbour(ng2)==ind) THEN
                    DO k2=0,5
                      IF (cell_glo(ind)%neighbour(k2)==ind2) THEN
                        IF (k2==k+1 .OR. k2==k+2 .OR. k2+6==k+1 .AND. k2+6==k+2) THEN
                          IF (cell_glo(ind1)%vertex(ng1)==-1 .AND. cell_glo(ind2)%vertex(ng2)==-1) THEN
                            ne=ne+1
                            cell_glo(ind)%vertex(k)=ne 
                            cell_glo(ind1)%vertex(ng1)=ne
                            cell_glo(ind2)%vertex(ng2)=ne
                          ELSE IF (cell_glo(ind1)%vertex(ng1)==0) THEN
                            cell_glo(ind)%vertex(k)=cell_glo(ind2)%vertex(ng2)
                            cell_glo(ind1)%vertex(ng1)=cell_glo(ind2)%vertex(ng2)
                          ELSE IF (cell_glo(ind2)%vertex(ng2)==0) THEN
                            cell_glo(ind)%vertex(k)=cell_glo(ind1)%vertex(ng1)
                            cell_glo(ind2)%vertex(ng2)=cell_glo(ind1)%vertex(ng1)
                          ENDIF
                        ENDIF
                      ENDIF
                    ENDDO
                  ENDIF
                ENDDO
              ENDIF
            ENDDO
          ENDIF
        ENDIF
      ENDDO
    ENDDO                      
                     
 
  END SUBROUTINE set_cell_vertex_area

        
  SUBROUTINE set_vertex_edge_area
  IMPLICIT NONE
    INTEGER :: nf
    INTEGER :: i,j
    INTEGER :: ind,ind2
    INTEGER :: ne
    INTEGER :: k,l     
  
    DO nf=1,nb_face
      DO j=-1,jjm_glo+1
        DO i=0,iim_glo+2
          ind=vertex_glo(i,j,nf)%ind
          IF ( ind/=-1) THEN
            IF (ind/=-1) THEN
              DO k=0,5
                ind2=vertex_glo(i,j,nf)%neighbour(k)
                IF (ind2/=-1) THEN
                  DO l=0,5
                    ne=cell_glo(ind)%edge(l)
                    IF (ne/=-1) THEN
                      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
                    ENDIF
                  ENDDO
                ENDIF
              ENDDO
            ENDIF
          ENDIF
        ENDDO
      ENDDO
    ENDDO


  END SUBROUTINE set_vertex_edge_area
  
  SUBROUTINE check_area
    USE spherical_geom_mod
  IMPLICIT NONE
    INTEGER :: nf,i,j,ind
    REAL :: lon,lat
    
    DO nf=1,nb_face
      DO j=-2,jjm_glo
        DO i=1,iim_glo+3
          CALL xyz2lonlat(vertex_glo(i,j,nf)%xyz,lon,lat)
          PRINT*,"vertex_glo%xyz : ",nf, i,j,lon*180./Pi,lat*180/Pi
        ENDDO
      ENDDO
    ENDDO
    
    DO ind=1,ncell_glo
      CALL xyz2lonlat(cell_glo(ind)%xyz,lon,lat)
      PRINT*,"cell_glo%xyz : ",ind,lon*180./Pi,lat*180/Pi
    ENDDO  
  END  SUBROUTINE check_area


  SUBROUTINE compute_metric
  IMPLICIT NONE
  
    CALL allocate_metric
    CALL compute_face_area
    CALL set_index_area
    CALL set_cell_area
    CALL compute_extended_face_area_bis
    CALL set_cell_edge_area
    CALL set_cell_vertex_area
    CALL set_vertex_edge_area
!    CALL check_area
  END SUBROUTINE compute_metric

END MODULE metric_icosa_area_old_mod
