MODULE metric_icosa_area_mod
  USE genmod
  USE grid_param, ONLY : iim_glo, jjm_glo
  USE domain_param
  USE metric_var_mod  
  
  TYPE t_global_cells
    REAL(rstd) :: xyz(3)
    INTEGER    :: assign_face
    INTEGER    :: ind
  END TYPE t_global_cells
  
  TYPE t_island
    TYPE(t_global_cells),ALLOCATABLE ::  global_cells(:,:)
  END TYPE t_island

  TYPE(t_island),ALLOCATABLE,TARGET :: islands(:)

  INTEGER,ALLOCATABLE :: tile_island(:) ! the island number for each tile

  INTEGER,ALLOCATABLE    :: delta_glo(:)

  TYPE tile_t
    INTEGER :: vertex_i(4)
    INTEGER :: vertex_j(4)
    LOGICAL :: initialized
  END TYPE

  TYPE edge_t
    INTEGER :: tile
    INTEGER :: vertex(2)
  END TYPE edge_t


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 init_metric(connect)
   USE getin_mod
   USE mpipara, ONLY: is_mpi_root
   IMPLICIT NONE
    TYPE(edge_t),ALLOCATABLE,INTENT(INOUT) :: connect(:,:)
    CHARACTER(LEN=255), ALLOCATABLE :: str_connect(:,:)
    CHARACTER(LEN=255) :: area_metric

    area_metric="hexagon"
    PRINT*,"area_metric ",area_metric
    CALL getin("area_metric", area_metric)
    IF (TRIM(area_metric)=='hexagon') THEN
      CALL init_metric_hexa(str_connect)

    ELSE IF (TRIM(area_metric)=='extended_hexagon0') THEN
      CALL init_metric_extended_hexagon(str_connect,0)
    ELSE IF (TRIM(area_metric)=='extended_hexagon1') THEN
      CALL init_metric_extended_hexagon(str_connect,1)
    ELSE IF (TRIM(area_metric)=='extended_hexagon2') THEN
      CALL init_metric_extended_hexagon(str_connect,2)
    ELSE IF (TRIM(area_metric)=='extended_hexagon3') THEN
      CALL init_metric_extended_hexagon(str_connect,3)
    ELSE IF (TRIM(area_metric)=='extended_hexagon4') THEN
      CALL init_metric_extended_hexagon(str_connect,4)
    ELSE IF (TRIM(area_metric)=='extended_hexagon5') THEN
      CALL init_metric_extended_hexagon(str_connect,5)
    ELSE IF (TRIM(area_metric)=='extended_hexagon6') THEN
      CALL init_metric_extended_hexagon(str_connect,6)

    ELSE IF (TRIM(area_metric)=='rhombus') THEN
      CALL init_metric_rhombus(str_connect)

    ELSE IF (TRIM(area_metric)=='extended_rhombus0') THEN
      CALL init_metric_extended_rhombus(str_connect,0)
    ELSE IF (TRIM(area_metric)=='extended_rhombus1') THEN
      CALL init_metric_extended_rhombus(str_connect,1)
    ELSE IF (TRIM(area_metric)=='extended_rhombus2') THEN
      CALL init_metric_extended_rhombus(str_connect,2)
    ELSE IF (TRIM(area_metric)=='extended_rhombus3') THEN
      CALL init_metric_extended_rhombus(str_connect,3)
    ELSE IF (TRIM(area_metric)=='extended_rhombus4') THEN
      CALL init_metric_extended_rhombus(str_connect,4)
    ELSE IF (TRIM(area_metric)=='extendedrhombus5') THEN
      CALL init_metric_extended_rhombus(str_connect,5)
    ELSE IF (TRIM(area_metric)=='extended_rhombus6') THEN
      CALL init_metric_extended_rhombus(str_connect,6)

    ELSE IF (TRIM(area_metric)=='star') THEN
      CALL init_metric_star(str_connect)


    ELSE IF (TRIM(area_metric)=='user_defined') THEN
      CALL init_metric_user_defined(str_connect)
    ELSE
      IF (is_mpi_root) THEN 
        PRINT*,"Undefined parameter for area_metric : ",area_metric
        PRINT*,"Choices are  : ","<hexagon>, <extended_hexagonX>, <rhombus>, <extended_rhombusX>, <star> or <user_defined>"
      ENDIF
      STOP
    ENDIF
    
    CALL init_metric_generic(str_connect, connect)
  END SUBROUTINE init_metric
  
  SUBROUTINE init_metric_hexa(str_connect)
   USE getin_mod
   IMPLICIT NONE
    CHARACTER(LEN=255), ALLOCATABLE, INTENT(INOUT) :: str_connect(:,:)
    
    nb_face=3
    
    ALLOCATE(str_connect(4,2))
    str_connect(:,1) = (/ "1 ", "AB", "2 ", "AD" /)
    str_connect(:,2) = (/ "2 ", "CD", "3 ", "AD" /)

  END SUBROUTINE init_metric_hexa

  SUBROUTINE init_metric_rhombus(str_connect)
   USE getin_mod
   IMPLICIT NONE
    CHARACTER(LEN=255), ALLOCATABLE, INTENT(INOUT) :: str_connect(:,:)
    
    nb_face=1
    
    ALLOCATE(str_connect(4,0))

  END SUBROUTINE init_metric_rhombus

  SUBROUTINE init_metric_extended_rhombus(str_connect, rhombus_size)
   USE getin_mod
   IMPLICIT NONE
    CHARACTER(LEN=255), ALLOCATABLE, INTENT(INOUT) :: str_connect(:,:)
    INTEGER, INTENT(IN) :: rhombus_size
    CHARACTER(LEN=16) :: str1,str2
    INTEGER :: sz

    nb_face=1+rhombus_size
    
    ALLOCATE(str_connect(4,rhombus_size))

    DO sz=1, rhombus_size
      WRITE(str1,*) sz
      WRITE(str2,*) sz+1
      str1=ADJUSTL(str1)
      str2=ADJUSTL(str2)
      str_connect(:,sz) = (/ TRIM(str1), "BC", TRIM(str2), "AD" /)
    ENDDO

  END SUBROUTINE init_metric_extended_rhombus

SUBROUTINE init_metric_extended_hexagon(str_connect, hexa_size)
   USE getin_mod
   IMPLICIT NONE
    CHARACTER(LEN=255), ALLOCATABLE, INTENT(INOUT) :: str_connect(:,:)
    INTEGER, INTENT(IN) :: hexa_size
    CHARACTER(LEN=16) :: str1,str2
    INTEGER :: sz

    nb_face=3+2*hexa_size
    
    ALLOCATE(str_connect(4,2+2*hexa_size))
    
    str_connect(:,1) = (/ "1 ", "AB", "2 ", "AD" /)
    str_connect(:,2) = (/ "2 ", "CD", "3 ", "AD" /)

    DO sz=1, hexa_size
      WRITE(str1,*) sz*2
      WRITE(str2,*) sz*2+2
      str1=ADJUSTL(str1)
      str2=ADJUSTL(str2)
      str_connect(:,2+2*sz-1) = (/ TRIM(str1), "BC", TRIM(str2), "AD" /)

      WRITE(str1,*) sz*2+2
      WRITE(str2,*) sz*2+3
      str1=ADJUSTL(str1)
      str2=ADJUSTL(str2)
      str_connect(:,2+2*sz) = (/ TRIM(str1), "CD", TRIM(str2), "AD" /)
    ENDDO

  END SUBROUTINE init_metric_extended_hexagon

  SUBROUTINE init_metric_star(str_connect)
   USE getin_mod
   IMPLICIT NONE
    CHARACTER(LEN=255), ALLOCATABLE, INTENT(INOUT) :: str_connect(:,:)
    
    nb_face=6
    
    ALLOCATE(str_connect(4,5))
    str_connect(:,1) = (/ "1 ", "AB", "2 ", "AD" /)
    str_connect(:,2) = (/ "2 ", "AB", "3 ", "AD" /)
    str_connect(:,3) = (/ "3 ", "AB", "4 ", "AD" /)
    str_connect(:,4) = (/ "4 ", "AB", "5 ", "AD" /)
    str_connect(:,5) = (/ "5 ", "AB", "6 ", "AD" /)

  END SUBROUTINE init_metric_star

SUBROUTINE init_metric_user_defined(str_connect)
   USE getin_mod
   IMPLICIT NONE
    CHARACTER(LEN=255), ALLOCATABLE, INTENT(INOUT) :: str_connect(:,:)
    CHARACTER(LEN=255) :: str_nb_connects(4)
    INTEGER :: nb_connects
    CHARACTER(LEN=16) :: str_number
    INTEGER :: i

    nb_face=1
    CALL getin("nb_tiles",nb_face)
    
    nb_connects=0
    DO
      WRITE(str_number,*) nb_connects+1
      str_number=ADJUSTL(str_number)
      str_nb_connects(:)=""
      CALL getin("tiles_connection"//TRIM(str_number), str_nb_connects)
      IF (TRIM(str_nb_connects(1))=="") THEN
        EXIT
      ELSE
        nb_connects=nb_connects+1
      ENDIF
    ENDDO 

    !CALL getin("tiles_nb_connections", nb_connects)
    ALLOCATE(str_connect(4, nb_connects))

    DO i=1,nb_connects
      WRITE(str_number,*) i
      str_number=ADJUSTL(str_number)
      CALL getin("tiles_connection"//TRIM(str_number), str_connect(:,i))
    ENDDO

  END SUBROUTINE init_metric_user_defined

  SUBROUTINE init_metric_generic(str_connect, connect)
   USE getin_mod
   IMPLICIT NONE
    CHARACTER(LEN=255),INTENT(IN) :: str_connect(:,:)
    TYPE(edge_t),ALLOCATABLE,INTENT(INOUT) :: connect(:,:)
    INTEGER :: nb_connects
    INTEGER :: stat
    INTEGER :: tile1, tile2
    CHARACTER(LEN=2) :: edge1, edge2
    INTEGER :: i,j,n
    INTEGER :: max_island_num, island_num, min_val
    
    ALLOCATE(tile_island(nb_face))
    nb_connects = size(str_connect,2)
    ALLOCATE(connect(4,nb_connects))
    
    ! find different tile island (connected tiles)
    tile_island(:) = -1
    island_num=0
    DO i=1,nb_connects
      READ(str_connect(1,i),*,iostat=stat) tile1
      READ(str_connect(3,i),*,iostat=stat) tile2
      IF (tile_island(tile1)==-1 .AND. tile_island(tile2)==-1) THEN
        island_num=island_num+1
        tile_island(tile1) = island_num
        tile_island(tile2) = island_num
      ELSE IF (tile_island(tile1)==-1 .AND. tile_island(tile2)/=-1) THEN
        tile_island(tile1) = tile_island(tile2)
      ELSE IF (tile_island(tile1)/=-1 .AND. tile_island(tile2)==-1) THEN
        tile_island(tile2) = tile_island(tile1)
      ELSE IF (tile_island(tile1)/=-1 .AND. tile_island(tile2)/=-1) THEN
        DO j=1,nb_face
          IF (tile_island(j)==tile_island(tile2)) tile_island(j)=tile_island(tile1)
        ENDDO
      ENDIF
    ENDDO
    
    ! unconnected tile are standalone island
    DO j=1,nb_face
      IF (tile_island(j)==-1) THEN
        island_num=island_num+1
        tile_island(j) = island_num
      ENDIF
    ENDDO
    
    ! compress island_num to get the true island number
    max_island_num=MAXVAL(tile_island(:))
    island_num=0
    DO n=1,max_island_num
      min_val=max_island_num+1
      DO j=1,nb_face
        IF (tile_island(j) > island_num .AND. tile_island(j)<min_val) min_val=tile_island(j)
      ENDDO
      IF (min_val>=island_num+1) THEN
        island_num = island_num+1
        DO j=1,nb_face
          IF (tile_island(j)==min_val) tile_island(j)=island_num
        ENDDO
      ENDIF
    ENDDO

    


    DO i=1,nb_connects
      READ(str_connect(1,i),*,iostat=stat) tile1
      READ(str_connect(3,i),*,iostat=stat) tile2
      edge1=str_connect(2,i)
      edge2=str_connect(4,i)
      
      connect(1,i)%tile=tile1

      IF      (edge1=="AB" .OR. edge1=="ab" .OR. edge1=="Ab".OR. edge1=="aB" ) THEN
        connect(1,i)%vertex(1)=1 ;  connect(1,i)%vertex(2)=2 
      ELSE IF (edge1=="BA" .OR. edge1=="ba" .OR. edge1=="bA".OR. edge1=="Ba" ) THEN
        connect(1,i)%vertex(1)=2 ;  connect(1,i)%vertex(2)=1 
      ELSE IF (edge1=="BC" .OR. edge1=="bc" .OR. edge1=="Bc".OR. edge1=="bC" ) THEN
        connect(1,i)%vertex(1)=2 ;  connect(1,i)%vertex(2)=3 
      ELSE IF (edge1=="CB" .OR. edge1=="cb" .OR. edge1=="cB".OR. edge1=="Cb" ) THEN
        connect(1,i)%vertex(1)=3 ;  connect(1,i)%vertex(2)=2 
      ELSE IF (edge1=="CD" .OR. edge1=="cd" .OR. edge1=="Cd".OR. edge1=="cD" ) THEN
        connect(1,i)%vertex(1)=3 ;  connect(1,i)%vertex(2)=4 
      ELSE IF (edge1=="DC" .OR. edge1=="dc" .OR. edge1=="dC".OR. edge1=="Dc" ) THEN
        connect(1,i)%vertex(1)=4 ;  connect(1,i)%vertex(2)=3 
      ELSE IF (edge1=="DA" .OR. edge1=="da" .OR. edge1=="dA".OR. edge1=="Da" ) THEN
        connect(1,i)%vertex(1)=4 ;  connect(1,i)%vertex(2)=1 
      ELSE IF (edge1=="AD" .OR. edge1=="ad" .OR. edge1=="Ad".OR. edge1=="aD" ) THEN
        connect(1,i)%vertex(1)=1 ;  connect(1,i)%vertex(2)=4 
      ELSE
        PRINT*,"Unreconize edge type in init_metric", edge1," must be either AB, BA, BC, CB, CD, DC, DA, AD"
        STOP
      ENDIF

      connect(2,i)%tile=tile2

      IF      (edge2=="AB" .OR. edge2=="ab" .OR. edge2=="Ab".OR. edge2=="aB" ) THEN
        connect(2,i)%vertex(1)=1 ;  connect(2,i)%vertex(2)=2 
      ELSE IF (edge2=="BA" .OR. edge2=="ba" .OR. edge2=="bA".OR. edge2=="Ba" ) THEN
        connect(2,i)%vertex(1)=2 ;  connect(2,i)%vertex(2)=1 
      ELSE IF (edge2=="BC" .OR. edge2=="bc" .OR. edge2=="Bc".OR. edge2=="bC" ) THEN
        connect(2,i)%vertex(1)=2 ;  connect(2,i)%vertex(2)=3 
      ELSE IF (edge2=="CB" .OR. edge2=="cb" .OR. edge2=="cB".OR. edge2=="Cb" ) THEN
        connect(2,i)%vertex(1)=3 ;  connect(2,i)%vertex(2)=2 
      ELSE IF (edge2=="CD" .OR. edge2=="cd" .OR. edge2=="Cd".OR. edge2=="cD" ) THEN
        connect(2,i)%vertex(1)=3 ;  connect(2,i)%vertex(2)=4 
      ELSE IF (edge2=="DC" .OR. edge2=="dc" .OR. edge2=="dC".OR. edge2=="Dc" ) THEN
        connect(2,i)%vertex(1)=4 ;  connect(2,i)%vertex(2)=3 
      ELSE IF (edge2=="DA" .OR. edge2=="da" .OR. edge2=="dA".OR. edge2=="Da" ) THEN
        connect(2,i)%vertex(1)=4 ;  connect(2,i)%vertex(2)=1 
      ELSE IF (edge2=="AD" .OR. edge2=="ad" .OR. edge2=="Ad".OR. edge2=="aD" ) THEN
        connect(2,i)%vertex(1)=1 ;  connect(2,i)%vertex(2)=4 
      ELSE
        PRINT*,"Unreconize edge type in init_metric", edge2," must be either AB, BA, BC, CB, CD, DC, DA, AD"
        STOP
      ENDIF

    ENDDO

  END SUBROUTINE init_metric_generic

  

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

  SUBROUTINE compute_face_area(connect)
  USE spherical_geom_mod
  USE vector
  USE getin_mod
  use earth_const
  IMPLICIT NONE
    TYPE(edge_t),INTENT(IN) :: connect(:,:)
    
    REAL(rstd) :: lonA,lonB,lonC,lonD,lonE,lonF,lonO, testLonO, testLatO
    REAL(rstd) :: latA,latB,latC,latD,latE,latF,latO
    REAL(rstd) :: rotation
    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    :: offset_i, offset_j, ei_i,ei_j,ej_i,ej_j
    INTEGER    :: max_i,max_j,min_i,min_j,istart, jstart, iend, jend
    CHARACTER(LEN=16) :: str_number

    INTEGER :: i,j, ii, jj, nf, is
    
    TYPE(tile_t) :: tiles(nb_face)
    INTEGER :: tile1,tile2
    CHARACTER(LEN=2) :: edge1, edge2
    INTEGER :: stat
    INTEGER :: num_island
    INTEGER    :: iim_gc, jjm_gc
    TYPE(t_global_cells), POINTER :: global_cells(:,:)


        
    tiles(:)%initialized=.FALSE.

    CALL connect_tiles(tiles, connect)
    
    ! find number of island
    num_island=0
    DO i=1,nb_face
      IF (tile_island(i)>num_island) num_island=tile_island(i)
    ENDDO
    
    ALLOCATE(islands(num_island))

    DO is=1,num_island
    
      ! A point of the rectangle area
      lonO = 0.
      latO = 0.
      rotation = 0.
      length = 1000
      CALL getin("area_center_lon", lonO)
      CALL getin("area_center_lat", latO)
      CALL getin("area_rotation", rotation)
      CALL getin("area_radius", length)
     
      WRITE(str_number,*) is
      str_number=ADJUSTL(str_number)
      CALL getin("area_center_lon"//TRIM(str_number), lonO)
      CALL getin("area_center_lat"//TRIM(str_number), latO)
      CALL getin("area_rotation"//TRIM(str_number), rotation)
      CALL getin("area_radius"//TRIM(str_number), length)

      min_i = HUGE(min_i)
      max_i = -(HUGE(max_i)/2)
      min_j = HUGE(min_j)
      max_j = -(HUGE(max_j)/2)

      DO i=1,nb_face
        IF (tile_island(i)/=is) CYCLE

        DO j=1,4
          IF (tiles(i)%vertex_i(j)<min_i) min_i=tiles(i)%vertex_i(j)
          IF (tiles(i)%vertex_i(j)>max_i) max_i=tiles(i)%vertex_i(j)
          IF (tiles(i)%vertex_j(j)<min_j) min_j=tiles(i)%vertex_j(j)
          IF (tiles(i)%vertex_j(j)>max_j) max_j=tiles(i)%vertex_j(j)
        ENDDO
      ENDDO
      
      iim_gc=MAX(max_i-min_i+1, max_j-min_j+1)

      iim_gc=iim_gc + 12
      jjm_gc=iim_gc

      istart = min_i-(iim_gc-(max_i-min_i))/2
      jstart = min_j-(iim_gc-(max_i-min_i))/2
      iend=istart+iim_gc-1
      jend=jstart+jjm_gc-1

      ALLOCATE(islands(is)%global_cells(istart:iend,jstart:jend))
      global_cells => islands(is)%global_cells
      length=(length*iim_gc)/iim_glo
      


      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,-(length/2*1000)/radius,A_xyz1)
      CALL rotate_Oy(B_xyz2,-(length/2*1000)/radius,B_xyz1)
      CALL rotate_Oy(C_xyz2,-(length/2*1000)/radius,C_xyz1)
      CALL rotate_Oy(D_xyz2,-(length/2*1000)/radius,D_xyz1)
      CALL rotate_Oy(E_xyz2,-(length/2*1000)/radius,E_xyz1)
      CALL rotate_Oy(F_xyz2,-(length/2*1000)/radius,F_xyz1)
      CALL rotate_Oy(O_xyz2,-(length/2*1000)/radius,O_xyz1)             

      CALL rotate_Oz(A_xyz1, rotation*Pi/180.,A_xyz2)
      CALL rotate_Oz(B_xyz1, rotation*Pi/180., B_xyz2)
      CALL rotate_Oz(C_xyz1, rotation*Pi/180., C_xyz2)
      CALL rotate_Oz(D_xyz1, rotation*Pi/180., D_xyz2)
      CALL rotate_Oz(E_xyz1, rotation*Pi/180., E_xyz2)
      CALL rotate_Oz(F_xyz1, rotation*Pi/180., F_xyz2)
      CALL rotate_Oz(O_xyz1, rotation*Pi/180., 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)
    
      



      global_cells(:,:)%ind = 0
      !face 1
      global_cells(istart,jstart)%xyz             =  F_xyz
      global_cells(iend,jstart)%xyz               =  A_xyz
      global_cells(iend,jend)%xyz                 =  B_xyz
      global_cells(istart,jend)%xyz               =  O_xyz

      !ABCO
      DO i=istart+1,iend-1
         !EF
        CALL div_arc_bis(global_cells(istart,jstart)%xyz, global_cells(iend,jstart)%xyz, (i-istart)/(iim_gc-1.), global_cells(i,jstart)%xyz)
         !OA
        CALL div_arc_bis(global_cells(istart,jend)%xyz, global_cells(iend,jend)%xyz, (i-istart)/(iim_gc-1.), global_cells(i,jend)%xyz)
      ENDDO
      
      DO j=jstart+1,jend-1
        DO i=istart,iend
          CALL div_arc_bis(global_cells(i,jstart)%xyz, global_cells(i,jend)%xyz, (j-jstart)/(jjm_gc-1.), global_cells(i,j)%xyz)
        ENDDO
      ENDDO    
     
      DO nf=1,nb_face
        IF (tile_island(nf)/=is) CYCLE

        ei_i = (tiles(nf)%vertex_i(2)-tiles(nf)%vertex_i(1))/(iim_glo-1)
        ei_j = (tiles(nf)%vertex_j(2)-tiles(nf)%vertex_j(1))/(iim_glo-1)
        ej_i = (tiles(nf)%vertex_i(4)-tiles(nf)%vertex_i(1))/(jjm_glo-1)
        ej_j = (tiles(nf)%vertex_j(4)-tiles(nf)%vertex_j(1))/(jjm_glo-1)
        
        offset_i = tiles(nf)%vertex_i(1)
        offset_j = tiles(nf)%vertex_j(1)
        
        DO j=-2,jjm_glo+3
          DO i=-2,iim_glo+3
            vertex_glo(i,j,nf)%assign_i = offset_i + (i-1)*ei_i + (j-1)*ej_i 
            vertex_glo(i,j,nf)%assign_j = offset_j + (i-1)*ei_j + (j-1)*ej_j
            vertex_glo(i,j,nf)%xyz = global_cells(vertex_glo(i,j,nf)%assign_i,vertex_glo(i,j,nf)%assign_j)%xyz
          ENDDO
        ENDDO
        IF (ei_i==1 .AND. ei_j==0) delta_glo(nf)=0
        IF (ei_i==0 .AND. ei_j==1) delta_glo(nf)=1
        IF (ei_i==-1 .AND. ei_j==1) delta_glo(nf)=2
        IF (ei_i==-1 .AND. ei_j==0) delta_glo(nf)=3
        IF (ei_i==0 .AND. ei_j==-1) delta_glo(nf)=4
        IF (ei_i==1 .AND. ei_j==-1) delta_glo(nf)=5
      ENDDO
    ENDDO
  END SUBROUTINE compute_face_area

  SUBROUTINE connect_tiles(tiles, connect)
  IMPLICIT NONE  
    TYPE(tile_t),INTENT(INOUT),TARGET :: tiles(:)  
    TYPE(edge_t),INTENT(IN) :: connect(:,:)

    INTEGER :: nb_tiles
    INTEGER :: nb_connects
        
    INTEGER :: t1,t2,tswap
    INTEGER :: v1,v2, v3
    INTEGER :: c,k
    
    TYPE(tile_t), POINTER :: my_tiles(:)
    
    nb_tiles=size(tiles)
    nb_connects=size(connect,2)
    
    my_tiles=>tiles
    
    DO c=1,nb_connects
      t1 = connect(1,c)%tile
      t2 = connect(2,c)%tile
      IF (tiles(t1)%initialized .AND. tiles(t2)%initialized) THEN
        PRINT*, "WARNING : overhead connection"
        CYCLE
      ELSE IF (tiles(t1)%initialized .AND. .NOT. tiles(t2)%initialized) THEN
        ! good order
      ELSE IF (.NOT. tiles(t1)%initialized .AND. tiles(t2)%initialized) THEN
        tswap = t2
        t2=t1
        t1=tswap
      ELSE IF (.NOT. tiles(t1)%initialized .AND. .NOT. tiles(t2)%initialized) THEN
        tiles(t1)%initialized=.TRUE.
        tiles(t1)%vertex_i = [1, iim_glo, iim_glo,   1]    
        tiles(t1)%vertex_j = [1,   1,     jjm_glo, jjm_glo]
      ENDIF
      
      v1 = connect(2,c)%vertex(1)
      tiles(t2)%vertex_i(v1) = tiles(t1)%vertex_i(connect(1,c)%vertex(1))
      tiles(t2)%vertex_j(v1) = tiles(t1)%vertex_j(connect(1,c)%vertex(1))

      v2 = connect(2,c)%vertex(2)
      tiles(t2)%vertex_i(v2) = tiles(t1)%vertex_i(connect(1,c)%vertex(2))
      tiles(t2)%vertex_j(v2) = tiles(t1)%vertex_j(connect(1,c)%vertex(2))
      
      DO k=1,3
        CALL next_vertex(t2,v1,v2,v3)
        v1=v2
        v2=v3
      ENDDO
      tiles(t2)%initialized=.TRUE.
      
    ENDDO
    ! check for unconnected tiles
    DO c=1,nb_tiles
      IF (.NOT. tiles(c)%initialized) THEN
        tiles(c)%initialized=.TRUE.
        tiles(c)%vertex_i = [1, iim_glo, iim_glo,   1]    
        tiles(c)%vertex_j = [1,   1,     jjm_glo, jjm_glo]
      ENDIF
    ENDDO

  CONTAINS
    SUBROUTINE next_vertex(tile,v1,v2,v3)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: tile
      INTEGER, INTENT(IN) :: v1,v2
      INTEGER, INTENT(OUT) :: v3

      INTEGER :: ei_i,ei_j,ej_i,ej_j

      IF (v1==1 .AND. v2==2) THEN
        v3=MOD(v2+4,4)+1
        ei_i = (tiles(tile)%vertex_i(v2) - tiles(tile)%vertex_i(v1))/(iim_glo-1)
        ei_j = (tiles(tile)%vertex_j(v2) - tiles(tile)%vertex_j(v1))/(iim_glo-1)
        CALL get_eij(ei_i,ei_j,ej_i,ej_j)
        tiles(tile)%vertex_i(v3)=tiles(tile)%vertex_i(v2) + ej_i*(jjm_glo-1)
        tiles(tile)%vertex_j(v3)=tiles(tile)%vertex_j(v2) + ej_j*(jjm_glo-1)
      ELSE IF (v1==2 .AND. v2==3) THEN
        v3=MOD(v2+4,4)+1
        ei_i = (tiles(tile)%vertex_i(v2) - tiles(tile)%vertex_i(v1))/(jjm_glo-1)
        ei_j = (tiles(tile)%vertex_j(v2) - tiles(tile)%vertex_j(v1))/(jjm_glo-1)
        CALL get_eij(ei_i,ei_j,ej_i,ej_j)
        tiles(tile)%vertex_i(v3)=tiles(tile)%vertex_i(v2) + (-ei_i + ej_i)*(iim_glo-1)
        tiles(tile)%vertex_j(v3)=tiles(tile)%vertex_j(v2) + (-ei_j + ej_j)*(iim_glo-1)
      ELSE IF (v1==3 .AND. v2==4) THEN
        v3=MOD(v2+4,4)+1
        ei_i = (tiles(tile)%vertex_i(v2) - tiles(tile)%vertex_i(v1))/(iim_glo-1)
        ei_j = (tiles(tile)%vertex_j(v2) - tiles(tile)%vertex_j(v1))/(iim_glo-1)
        CALL get_eij(ei_i,ei_j,ej_i,ej_j)
        tiles(tile)%vertex_i(v3)=tiles(tile)%vertex_i(v2) + ej_i*(jjm_glo-1)
        tiles(tile)%vertex_j(v3)=tiles(tile)%vertex_j(v2) + ej_j*(jjm_glo-1)
      ELSE IF (v1==4 .AND. v2==1) THEN
        v3=MOD(v2+4,4)+1
        ei_i = (tiles(tile)%vertex_i(v2) - tiles(tile)%vertex_i(v1))/(jjm_glo-1)
        ei_j = (tiles(tile)%vertex_j(v2) - tiles(tile)%vertex_j(v1))/(jjm_glo-1)
        CALL get_eij(ei_i,ei_j,ej_i,ej_j)
        tiles(tile)%vertex_i(v3)=tiles(tile)%vertex_i(v2) + (-ei_i + ej_i)*(iim_glo-1)
        tiles(tile)%vertex_j(v3)=tiles(tile)%vertex_j(v2) + (-ei_j + ej_j)*(iim_glo-1)
      
      ELSE IF (v2==1 .AND. v1==2) THEN
        v3=MOD(v2-2+4,4)+1
        ei_i = (tiles(tile)%vertex_i(v2) - tiles(tile)%vertex_i(v1))/(iim_glo-1)
        ei_j = (tiles(tile)%vertex_j(v2) - tiles(tile)%vertex_j(v1))/(iim_glo-1)
        CALL get_eij(ei_i,ei_j,ej_i,ej_j)
        tiles(tile)%vertex_i(v3)=tiles(tile)%vertex_i(v2) + (-ej_i  )*(jjm_glo-1)
        tiles(tile)%vertex_j(v3)=tiles(tile)%vertex_j(v2) + (-ej_j  )*(jjm_glo-1)
      ELSE IF (v2==2 .AND. v1==3) THEN
        v3=MOD(v2-2+4,4)+1
        ei_i = (tiles(tile)%vertex_i(v2) - tiles(tile)%vertex_i(v1))/(jjm_glo-1)
        ei_j = (tiles(tile)%vertex_j(v2) - tiles(tile)%vertex_j(v1))/(jjm_glo-1)
        CALL get_eij(ei_i,ei_j,ej_i,ej_j)
        tiles(tile)%vertex_i(v3)=tiles(tile)%vertex_i(v2) + (ei_i - ej_i)*(iim_glo-1)
        tiles(tile)%vertex_j(v3)=tiles(tile)%vertex_j(v2) + (ei_j - ej_j)*(iim_glo-1)
      ELSE IF (v2==3 .AND. v1==4) THEN
        v3=MOD(v2-2+4,4)+1
        ei_i = (tiles(tile)%vertex_i(v2) - tiles(tile)%vertex_i(v1))/(iim_glo-1)
        ei_j = (tiles(tile)%vertex_j(v2) - tiles(tile)%vertex_j(v1))/(iim_glo-1)
        CALL get_eij(ei_i,ei_j,ej_i,ej_j)
        tiles(tile)%vertex_i(v3)=tiles(tile)%vertex_i(v2) + (-ej_i  )*(jjm_glo-1)
        tiles(tile)%vertex_j(v3)=tiles(tile)%vertex_j(v2) + (-ej_j  )*(jjm_glo-1)
      ELSE IF (v2==4 .AND. v1==1) THEN
        v3=MOD(v2-2,4)+1
        ei_i = (tiles(tile)%vertex_i(v2) - tiles(tile)%vertex_i(v1))/(jjm_glo-1)
        ei_j = (tiles(tile)%vertex_j(v2) - tiles(tile)%vertex_j(v1))/(jjm_glo-1)
        CALL get_eij(ei_i,ei_j,ej_i,ej_j)
        tiles(tile)%vertex_i(v3)=tiles(tile)%vertex_i(v2) + (ei_i - ej_i)*(iim_glo-1)
        tiles(tile)%vertex_j(v3)=tiles(tile)%vertex_j(v2) + (ei_j - ej_j)*(iim_glo-1)
      ENDIF

    END SUBROUTINE next_vertex

    SUBROUTINE get_eij(ei_i,ei_j,ej_i,ej_j)
    IMPLICIT NONE
      INTEGER, INTENT(IN) :: ei_i,ei_j
      INTEGER, INTENT(OUT) :: ej_i,ej_j

      IF (ei_i==1 .AND. ei_j==0) THEN
        ej_i=0
        ej_j=1
      ELSE IF (ei_i==0 .AND. ei_j==1) THEN
        ej_i=-1
        ej_j=1
      ELSE IF (ei_i==-1 .AND. ei_j==1) THEN
        ej_i=-1
        ej_j=0
      ELSE IF (ei_i==-1 .AND. ei_j==0) THEN
        ej_i=0
        ej_j=-1
      ELSE IF (ei_i==0 .AND. ei_j==-1) THEN
        ej_i=1
        ej_j=-1
      ELSE IF (ei_i==1 .AND. ei_j==-1) THEN 
        ej_i=1
        ej_j=0
      ENDIF
    END SUBROUTINE get_eij      

  END SUBROUTINE connect_tiles

  SUBROUTINE Set_index_area
  IMPLICIT NONE
    INTEGER :: delta
    INTEGER :: nf1,nf2,i
   
    DO nf1=1,nb_face
      DO nf2=1,nb_face
        delta=delta_glo(nf1)-delta_glo(nf2)
        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,ind1, ind2
    INTEGER :: nf,nf1,nf2,nfm1,nfp1
    INTEGER :: i,j, halo
    
    ind=0

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

! interieur des faces
    DO nf=1,nb_face     
      DO j=2,jjm_glo-1
        DO i=2,iim_glo-1
          CALL associate_vertex(nf,i,j)
        ENDDO
      ENDDO
    ENDDO

    DO halo=0,3
      DO nf=1,nb_face     
        DO j=1-halo,jjm_glo+halo
          CALL associate_vertex(nf,1-halo,j)
          CALL associate_vertex(nf,iim_glo+halo,j)
        ENDDO
      
        DO i=1-halo,iim_glo+halo
          CALL associate_vertex(nf,i,1-halo)
          CALL associate_vertex(nf,i,jjm_glo+halo)
        ENDDO
      ENDDO
    ENDDO

    ncell_glo=ind

!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Association des voisins !
!!!!!!!!!!!!!!!!!!!!!!!!!!!
    DO nf=1,nb_face     
      DO j=-1,jjm_glo+2
        DO i=-1,iim_glo+2
          ind1=vertex_glo(i,j,nf)%ind
          nf1=cell_glo(ind1)%assign_face

          ! right          
          ind2=vertex_glo(i+1,j,nf)%ind   
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind1)%neighbour(tab_index(nf,nf1,right-1))=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,left-1))=ind1
          vertex_glo(i,j,nf)%neighbour(right-1) = ind2

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

          ! left          
          ind2=vertex_glo(i-1,j,nf)%ind   
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind1)%neighbour(tab_index(nf,nf1,left-1))=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,right-1))=ind1
          vertex_glo(i,j,nf)%neighbour(left-1) = ind2

          ! ldown          
          ind2=vertex_glo(i,j-1,nf)%ind   
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind1)%neighbour(tab_index(nf,nf1,ldown-1))=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,rup-1))=ind1
          vertex_glo(i,j,nf)%neighbour(ldown-1) = ind2

          ! rdown          
          ind2=vertex_glo(i+1,j-1,nf)%ind   
          nf2=cell_glo(ind2)%assign_face
          cell_glo(ind1)%neighbour(tab_index(nf,nf1,rdown-1))=ind2
          cell_glo(ind2)%neighbour(tab_index(nf,nf2,lup-1))=ind1
          vertex_glo(i,j,nf)%neighbour(rdown-1) = ind2
        ENDDO
      ENDDO
    ENDDO

 !! assignation des delta 
    DO nf=1,nb_face
      DO j=-2,jjm_glo+3
        DO i=-2,iim_glo+3 
          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    

! 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

  CONTAINS
    
    SUBROUTINE associate_vertex(nf,i,j)
      INTEGER,INTENT(IN) :: i,j,nf
      INTEGER :: index
      TYPE(t_global_cells),POINTER :: global_cells(:,:)

      global_cells=>islands(tile_island(nf))%global_cells

      index = global_cells(vertex_glo(i,j,nf)%assign_i,vertex_glo(i,j,nf)%assign_j)%ind 
      IF (index==0) THEN
        ind=ind+1
        index=ind
        cell_glo(index)%assign_face=nf
        cell_glo(index)%assign_i=i
        cell_glo(index)%assign_j=j
        global_cells(vertex_glo(i,j,nf)%assign_i,vertex_glo(i,j,nf)%assign_j)%ind=index
      ENDIF
      vertex_glo(i,j,nf)%ind = index
    END SUBROUTINE associate_vertex

  END SUBROUTINE set_cell_area



  SUBROUTINE set_cell_edge_area
  IMPLICIT NONE
    INTEGER :: ind,k,ind2,ng,ne,i,j,nf
    
    DO ind=1,ncell_glo
      cell_glo(ind)%edge(:)=-1
    ENDDO

    ne=0
    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
            IF (cell_glo(ind)%edge(k)==-1) THEN
              ne=ne+1
              ind2=cell_glo(ind)%neighbour(k)
              DO ng=0,5
                IF (cell_glo(ind2)%neighbour(ng)==ind) THEN
                  cell_glo(ind)%edge(k)=ne
                  edge_glo(ne)%e1=ind
                  edge_glo(ne)%e2=ind2
                  cell_glo(ind2)%edge(ng)=ne
                ENDIF
              ENDDO
            ENDIF
          ENDDO
        ENDDO
      ENDDO
    ENDDO

  END SUBROUTINE set_cell_edge_area    
      

  SUBROUTINE set_cell_vertex_area
  IMPLICIT NONE
    INTEGER :: nf,i,j
    INTEGER :: k,k2
    INTEGER :: ind,ind1,ind2
    INTEGER :: ng
    INTEGER :: ne     
  
    DO ind=1,ncell_glo
      cell_glo(ind)%vertex(:)=-1
    ENDDO
    
    ne=0
    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
            IF (cell_glo(ind)%vertex(k)==-1) THEN
              ne=ne+1
              cell_glo(ind)%vertex(k)=ne
              ind1=cell_glo(ind)%neighbour(k)
              ind2=cell_glo(ind)%neighbour(MOD(k+1,6))
              
              DO ng=0,5
                IF (      (cell_glo(ind1)%neighbour(ng)==ind  .AND. cell_glo(ind1)%neighbour(MOD(ng+1,6))==ind2)      &
                  &  .OR. (cell_glo(ind1)%neighbour(ng)==ind2 .AND. cell_glo(ind1)%neighbour(MOD(ng+1,6))==ind )) THEN
                  cell_glo(ind1)%vertex(ng)=ne
                  EXIT
                ENDIF
              ENDDO

              DO ng=0,5
                IF (      (cell_glo(ind2)%neighbour(ng)==ind  .AND. cell_glo(ind2)%neighbour(MOD(ng+1,6))==ind1)  &
                   & .OR. (cell_glo(ind2)%neighbour(ng)==ind1 .AND. cell_glo(ind2)%neighbour(MOD(ng+1,6))==ind )) THEN
                  cell_glo(ind2)%vertex(ng)=ne
                  EXIT
                ENDIF
              ENDDO
              
            ENDIF
          ENDDO
        ENDDO
      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=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_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=-1,jjm_glo+2
        DO i=-1,iim_glo+2
          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
    TYPE(edge_t),ALLOCATABLE :: connect(:,:)

    CALL init_metric(connect)
    CALL allocate_metric
    CALL compute_face_area(connect)
    CALL set_index_area
    CALL set_cell_area
    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_mod
