!----------------------------------------------------------------------
! NEMO system team, System and Interface for oceanic RElocable Nesting
!----------------------------------------------------------------------
!
! MODULE: grid
!
! DESCRIPTION:
!> @brief grid manager <br/>
!>
!> @details
!> 
!> @author
!> J.Paul
! REVISION HISTORY:
!> @date Nov, 2013 - Initial Version
!
!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
!> @todo
!----------------------------------------------------------------------
MODULE grid
   USE netcdf
   USE kind                            ! F90 kind parameter
   USE fct                             ! basic usefull function
   USE global                          ! global parameter
   USE phycst                          ! physical constant
   USE logger                             ! log file manager
   USE file                            ! file manager
   USE var                             ! variable manager
   USE dim                             ! dimension manager
   USE dom                             ! domain manager
   USE iom                             ! I/O manager
   USE mpp                             ! MPP manager
   USE iom_mpp                         ! MPP I/O manager
   IMPLICIT NONE
   PRIVATE
   ! NOTE_avoid_public_variables_if_possible

   ! type and variable

   ! function and subroutine
   PUBLIC :: grid_check_dom !< check domain validity 
   PUBLIC :: grid_get_coarse_index !< get closest coarse grid indices of fine grid domain.
   PUBLIC :: grid_is_global !< check if grid is global or not
   PUBLIC :: grid_get_closest !< return closest coarse grid point from another point
   PUBLIC :: grid_distance !< compute grid distance to a point
   PUBLIC :: grid_get_fine_offset !< get fine grid offset
   PUBLIC :: grid_check_coincidence !< check fine and coarse grid coincidence
   PUBLIC :: grid_get_perio !< return NEMO periodicity index
   PUBLIC :: grid_get_pivot !< return NEMO pivot point index
   PUBLIC :: grid_add_ghost !< add ghost cell at boundaries.
   PUBLIC :: grid_del_ghost !< delete ghost cell at boundaries.
   PUBLIC :: grid_get_ghost !< return ghost cell factor
   PUBLIC :: grid_split_domain !< 
   PUBLIC :: grid_fill_small_dom !< 

   PRIVATE :: grid_get_coarse_index_ff
   PRIVATE :: grid_get_coarse_index_cf
   PRIVATE :: grid_get_coarse_index_fc
   PRIVATE :: grid_get_coarse_index_cc
   PRIVATE :: grid__get_ghost_f  
   PRIVATE :: grid__get_ghost_ll
   PRIVATE :: grid__check_corner
   
   INTERFACE  grid_get_ghost
      MODULE PROCEDURE grid__get_ghost_ll
      MODULE PROCEDURE grid__get_ghost_f
   END INTERFACE  grid_get_ghost

   INTERFACE  grid_get_coarse_index
      MODULE PROCEDURE grid_get_coarse_index_ff
      MODULE PROCEDURE grid_get_coarse_index_cf
      MODULE PROCEDURE grid_get_coarse_index_fc
      MODULE PROCEDURE grid_get_coarse_index_cc
   END INTERFACE grid_get_coarse_index

CONTAINS
   !-------------------------------------------------------------------
   !> @brief 
   !> This funtion return NEMO pivot point index of the input variable.
   !> - F-point : 0
   !> - T-point : 1
   !>
   !> @warning 
   !> - variable must be nav_lon or nav_lat
   !> - do not work with ORCA2 grid (T-point)
   !>
   !> @author J.Paul
   !> - Nov, 2013- Subroutine written
   !
   !> @todo 
   !> - improve check between T or F pivot.
   !
   !> @param[in] td_file : file structure
   !> @param[in] cd_varname : variable name
   !> @return NEMO pivot point index
   !-------------------------------------------------------------------
   !> @code
   INTEGER(i4) FUNCTION grid_get_pivot(td_file)
      IMPLICIT NONE
      ! Argument      
      TYPE(TFILE),       INTENT(IN) :: td_file

      ! local variable
      TYPE(TVAR)                        :: tl_var
      INTEGER(i4)                       :: il_varid
      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim

      ! loop indices
      INTEGER(i4) :: ji

      INTEGER(i4) :: it1
      INTEGER(i4) :: it2
      INTEGER(i4) :: jt1
      INTEGER(i4) :: jt2

      INTEGER(i4) :: if1
      INTEGER(i4) :: if2
      INTEGER(i4) :: jf1
      INTEGER(i4) :: jf2
      !----------------------------------------------------------------
      ! initialise
      grid_get_pivot=-1

      ! look for suitable variable
      il_varid=0
      DO ji=1,td_file%i_nvar
         IF( .NOT. ALL(td_file%t_var(ji)%t_dim(1:2)%l_use) ) CYCLE
         SELECT CASE(TRIM(fct_lower(td_file%t_var(ji)%c_stdname)) )
            CASE('longitude','latitude')
            CASE DEFAULT
               il_varid=ji
               EXIT
         END SELECT
      ENDDO

      IF( il_varid/=0 )THEN
         IF( ASSOCIATED(td_file%t_var(il_varid)%d_value) )THEN
            CALL logger_debug("GRID GET PIVOT: ASSOCIATED")
            tl_var=td_file%t_var(il_varid)
         ELSE
            ! read variable
            il_dim(:)=td_file%t_var(il_varid)%t_dim(:)%i_len

            CALL logger_debug("GRID GET PIVOT: read variable")
            tl_var=iom_read_var(td_file, td_file%t_var(il_varid)%c_name,     &
            &                   id_start=(/1,il_dim(2)-3,1,1/), &
            &                   id_count=(/3,4,1,1/) )
         ENDIF

         CALL logger_debug("GRID GET PIVOT: use variable "//TRIM(tl_var%c_name))

         IF( ASSOCIATED(tl_var%d_value) )THEN

            CALL logger_debug("GRID GET PIVOT: point "//TRIM(tl_var%c_point))
            ! T-point pivot !case of ORCA2, ORCA025, ORCA12 grid
            it1=1 ; jt1=4 
            it2=3 ; jt2=2 

            ! F-point pivot !case of ORCA05 grid
            if1=1 ; jf1=4 
            if2=2 ; jf2=3 

            SELECT CASE(TRIM(tl_var%c_point))
            CASE('T')
               IF( ABS(tl_var%d_value(it1,jt1,1,1)) == &
               &   ABS(tl_var%d_value(it2,jt2,1,1)) )THEN
                  CALL logger_info("GRID GET PIVOT: T-pivot")
                  grid_get_pivot=1
               ELSEIF( ABS(tl_var%d_value(if1,jf1,1,1)) == &
               &       ABS(tl_var%d_value(if2,jf2,1,1)) )THEN
                  CALL logger_info("GRID GET PIVOT: F-pivot")
                  grid_get_pivot=0
               ELSE
                  CALL logger_error("GRID GET PIVOT: something wrong when "//&
                  &  "computing pivot point")
               ENDIF
            CASE('U')
               IF( ABS(tl_var%d_value(it1  ,jt1,1,1)) == &
               &   ABS(tl_var%d_value(it2-1,jt2,1,1)) )THEN
                  CALL logger_info("GRID GET PIVOT: T-pivot")
                  grid_get_pivot=1
               ELSEIF( ABS(tl_var%d_value(if1  ,jf1,1,1)) == &
               &       ABS(tl_var%d_value(if2-1,jf2,1,1)) )THEN
                  CALL logger_info("GRID GET PIVOT: F-pivot")
                  grid_get_pivot=0
               ELSE
                  CALL logger_error("GRID GET PIVOT: something wrong when "//&
                  &  "computing pivot point")
               ENDIF
            CASE('V')
               IF( ABS(tl_var%d_value(it1,jt1  ,1,1)) == &
               &   ABS(tl_var%d_value(it2,jt2-1,1,1)) )THEN
                  CALL logger_info("GRID GET PIVOT: T-pivot")
                  grid_get_pivot=1
               ELSEIF( ABS(tl_var%d_value(if1,jf1  ,1,1)) == &
               &       ABS(tl_var%d_value(if2,jf2-1,1,1)) )THEN
                  CALL logger_info("GRID GET PIVOT: F-pivot")
                  grid_get_pivot=0
               ELSE
                  CALL logger_error("GRID GET PIVOT: something wrong when "//&
                  &  "computing pivot point")
               ENDIF
            CASE('F')
               IF( ABS(tl_var%d_value(it1  ,jt1  ,1,1)) == &
               &   ABS(tl_var%d_value(it2-1,jt2-1,1,1)) )THEN
                  CALL logger_info("GRID GET PIVOT: T-pivot")
                  grid_get_pivot=1
               ELSEIF( ABS(tl_var%d_value(if1  ,jf1  ,1,1)) == &
               &       ABS(tl_var%d_value(if2-1,jf2-1,1,1)) )THEN
                  CALL logger_info("GRID GET PIVOT: F-pivot")
                  grid_get_pivot=0
               ELSE
                  CALL logger_error("GRID GET PIVOT: something wrong when "//&
                  &  "computing pivot point")
               ENDIF
            END SELECT
         ELSE
            CALL logger_error("GRID GET PIVOT: can't compute pivot point. "//&
            &  "no value associated to variable "//TRIM(tl_var%c_name) )
         ENDIF

      ELSE
         CALL logger_error("GRID GET PIVOT: no suitable variable to compute "//&
         &              "pivot point in file "//TRIM(td_file%c_name))
      ENDIF

   END FUNCTION grid_get_pivot
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief 
   !> This funtion return NEMO periodicity index of the input file.
   !> The variable used must be on T point.
   !>
   !> @note the NEMO periodicity index can't be compute from coordinates file,
   !> neither with mpp files.
   !>
   !> 0: closed boundaries
   !> 1: cyclic east-west boundary
   !> 2: symmetric boundary condition across the equator
   !> 3: North fold boundary (with a F-point pivot)
   !> 4: North fold boundary (with a F-point pivot) and cyclic east-west boundary
   !> 5: North fold boundary (with a T-point pivot)
   !> 6: North fold boundary (with a T-point pivot) and cyclic east-west boundary
   !>
   !> @author J.Paul
   !> - Nov, 2013- Subroutine written
   !
   !> @todo 
   !> - improve check between T or F pivot.
   !> - manage mpp case (read only border files)
   !
   !> @param[in] td_file : file structure
   !> @param[in] id_pivot : pivot point
   !> @return NEMO periodicity index
   !-------------------------------------------------------------------
   !> @code
   INTEGER(i4) FUNCTION grid_get_perio(td_file, id_pivot)
      IMPLICIT NONE

      ! Argument      
      TYPE(TFILE), INTENT(IN) :: td_file
      INTEGER(i4), INTENT(IN) :: id_pivot

      ! local variable
      TYPE(TVAR)                        :: tl_var
      INTEGER(i4)                       :: il_varid
      INTEGER(i4), DIMENSION(ip_maxdim) :: il_dim

      ! loop indices
      INTEGER(i4) :: ji
      !----------------------------------------------------------------

      ! initialise
      grid_get_perio=-1

      IF( id_pivot < 0 .OR. id_pivot > 1 )THEN
         CALL logger_error("GRID GET PERIO: invalid pivot point index. "//&
         &  "you should use grid_get_pivot to compute it")
      ENDIF

      ! look for suitable variable
      il_varid=0
      DO ji=1,td_file%i_nvar
         IF( .NOT. ALL(td_file%t_var(ji)%t_dim(1:2)%l_use) ) CYCLE
         SELECT CASE(TRIM(fct_lower(td_file%t_var(ji)%c_stdname)) )
            CASE('longitude','latitude')
            CASE DEFAULT
               il_varid=ji
               EXIT
         END SELECT
      ENDDO

      IF( il_varid==0 )THEN

         CALL logger_error("GRID GET PERIO: no suitable variable to compute "//&
         &              " periodicity in file "//TRIM(td_file%c_name))
      ELSE
         il_dim(:)=td_file%t_var(il_varid)%t_dim(:)%i_len

         IF( ASSOCIATED(td_file%t_var(il_varid)%d_value) )THEN
            tl_var=td_file%t_var(il_varid)
         ELSE
            ! read variable
            tl_var=iom_read_var(td_file, td_file%t_var(il_varid)%c_name, &
            &                   id_start=(/1,1,1,1/), &
            &                   id_count=(/il_dim(1),il_dim(2),1,1/) )
         ENDIF

         IF(ALL(tl_var%d_value(    1    ,    :    ,1,1)/=tl_var%d_fill).AND.&
         &  ALL(tl_var%d_value(il_dim(1),    :    ,1,1)/=tl_var%d_fill).AND.&
         &  ALL(tl_var%d_value(    :    ,    1    ,1,1)/=tl_var%d_fill).AND.&
         &  ALL(tl_var%d_value(    :    ,il_dim(2),1,1)/=tl_var%d_fill))THEN
         ! no boundary closed
            CALL logger_warn("GRID GET PERIO: can't determined periodicity. "//&
            &             "there is no boundary closed for variable "//&
            &              TRIM(tl_var%c_name)//" in file "//&
            &              TRIM(td_file%c_name) )
         ELSE
            ! check periodicity
            IF(ANY(tl_var%d_value(   1     ,:,1,1)/=tl_var%d_fill).OR.&
            &  ANY(tl_var%d_value(il_dim(1),:,1,1)/=tl_var%d_fill))THEN
            ! East-West cyclic (1,4,6)

               IF( ANY(tl_var%d_value(:, 1, 1, 1) /= tl_var%d_fill) )THEN
               ! South boundary not closed 

                  CALL logger_error("GRID GET PERIO: should have been an "//&
                  &              "impossible case")
                  CALL logger_debug("GRID GET PERIO: East_West cyclic")
                  CALL logger_debug("GRID GET PERIO: South boundary not closed")

               ELSE
               ! South boundary closed (1,4,6)
                  CALL logger_info("GRID GET PERIO: South boundary closed")

                  IF(ANY(tl_var%d_value(:,il_dim(2),1,1)/=tl_var%d_fill) )THEN
                  ! North boundary not closed (4,6)
                     CALL logger_info("GRID GET PERIO: North boundary not closed")
                     ! check pivot
                     SELECT CASE(id_pivot)
                        CASE(0)
                           ! F pivot
                           grid_get_perio=4
                        CASE(1)
                           ! T pivot
                           grid_get_perio=6
                        CASE DEFAULT
                           CALL logger_error("GRID GET PERIO: invalid pivot ")
                     END SELECT
                  ELSE
                  ! North boundary closed
                     CALL logger_info("GRID GET PERIO: North boundary closed")
                     grid_get_perio=1 ! North and South boundaries closed
                  ENDIF

               ENDIF

            ELSE
            ! East-West boundaries closed (0,2,3,5)
               CALL logger_info("GRID GET PERIO: East West boundaries closed")

               IF( ANY(tl_var%d_value(:, 1, 1, 1) /= tl_var%d_fill) )THEN
               ! South boundary not closed (2)
                  CALL logger_info("GRID GET PERIO: South boundary not closed")

                  IF(ANY(tl_var%d_value(:,il_dim(2),1,1)/=tl_var%d_fill))THEN
                  ! North boundary not closed
                     CALL logger_error("GRID GET PERIO: should have been "//&
                     &              "an impossible case")
                     CALL logger_debug("GRID GET PERIO: East West boundaries "//&
                     &              "closed")
                     CALL logger_debug("GRID GET PERIO: South boundary not closed")
                     CALL logger_debug("GRID GET PERIO: North boundary not closed")
                  ELSE
                  ! North boundary closed
                     grid_get_perio=2   ! East-West and North boundaries closed 
                  ENDIF

               ELSE
               ! South boundary closed (0,3,5)
                  CALL logger_info("GRID GET PERIO: South boundary closed")

                  IF(ANY(tl_var%d_value(:,il_dim(2),1,1)/=tl_var%d_fill))THEN
                  ! North boundary not closed (3,5)
                     CALL logger_info("GRID GET PERIO: North boundary not closed")
                     ! check pivot
                     SELECT CASE(id_pivot)
                        CASE(0)
                           ! F pivot
                           grid_get_perio=3
                        CASE(1)
                           ! T pivot
                           grid_get_perio=5
                        CASE DEFAULT
                           CALL logger_error("GRID GET PERIO: invalid pivot")
                     END SELECT
                  ELSE
                  ! North boundary closed   
                     CALL logger_info("GRID GET PERIO: North boundary closed")
                     grid_get_perio=0   ! all boundary closed
                  ENDIF

               ENDIF

            ENDIF

         ENDIF
      ENDIF

   END FUNCTION grid_get_perio
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This subroutine check domain validity.
   !
   !> @details
   !> If maximum latitude greater than 88°N, program will stop. 
   !> It is not able to manage north fold boundary for now.
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] cd_coord : coordinate file 
   !> @param[in] id_imin : i-direction lower left  point indice  
   !> @param[in] id_imax : i-direction upper right point indice 
   !> @param[in] id_jmin : j-direction lower left  point indice 
   !> @param[in] id_jmax : j-direction upper right point indice 
   !>
   !> @todo
   !> - use domain instead of start count
   !-------------------------------------------------------------------
   !> @code
   SUBROUTINE grid_check_dom(td_coord, id_imin, id_imax, id_jmin, id_jmax)
      IMPLICIT NONE
      ! Argument      
      TYPE(TFILE), INTENT(IN) :: td_coord
      INTEGER(i4), INTENT(IN) :: id_imin
      INTEGER(i4), INTENT(IN) :: id_imax
      INTEGER(i4), INTENT(IN) :: id_jmin
      INTEGER(i4), INTENT(IN) :: id_jmax

      ! local variable
      TYPE(TVAR)                        :: tl_var

      TYPE(TFILE)                       :: tl_coord

      TYPE(TMPP)                        :: tl_mppcoord

      TYPE(TDOM)                        :: tl_dom

      ! loop indices
      !----------------------------------------------------------------

      IF( id_jmin >= id_jmax )THEN

         CALL logger_fatal("GRID CHECK DOM: invalid domain. "//&
         &  "can not create configuration with north pole.")

      ELSE

         IF( td_coord%i_id == 0 )THEN
            CALL logger_error("GRID CHECK DOM: can not check domain. "//&
            &  " file "//TRIM(td_coord%c_name)//" not opened." )
         ELSE

            IF( id_imin == id_imax .AND. td_coord%i_ew < 0 )THEN
               CALL logger_fatal("GRID CHECK DOM: invalid domain."//&
               &  " can not create east-west cyclic fine grid"//&
               &  " inside closed coarse grid")
            ENDIF

            !1- read domain
            tl_coord=td_coord
            CALL iom_open(tl_coord)

            !1-1 compute domain
            tl_dom=dom_init( tl_coord,        &
            &                 id_imin, id_imax,&
            &                 id_jmin, id_jmax )
            
            !1-2 close file
            CALL iom_close(tl_coord)

            !1-3 read variables on domain (ugly way to do it, have to work on it)
            !1-3-1 init mpp structure
            tl_mppcoord=mpp_init(tl_coord)   

            CALL file_clean(tl_coord)

            !1-3-2 get processor to be used
            CALL mpp_get_use( tl_mppcoord, tl_dom )

            !1-3-3 open mpp files
            CALL iom_mpp_open(tl_mppcoord)

            !1-3-4 read variable value on domain
            tl_var=iom_mpp_read_var(tl_mppcoord,'latitude',td_dom=tl_dom)

            !1-3-5 close mpp files
            CALL iom_mpp_close(tl_mppcoord)

            !1-3-6 clean structure
            CALL mpp_clean(tl_mppcoord)

            IF( MAXVAL(tl_var%d_value(:,:,:,:), &
            &          tl_var%d_value(:,:,:,:)/= tl_var%d_fill) >= 88.0 )THEN
               
               CALL logger_debug("GRID CHECK DOM: max latitude "//&
               &  TRIM(fct_str(MAXVAL(tl_var%d_value(:,:,:,:)))) )
               CALL logger_fatal("GRID CHECK DOM: invalid domain. "//&
               &  "can not create configuration too close from north pole.")

            ENDIF

            ! clean
            CALL var_clean(tl_var)

         ENDIF


      ENDIF

   END SUBROUTINE grid_check_dom
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This function get closest coarse grid indices of fine grid domain.
   !
   !> @details
   !>
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] td_coord0 : coarse grid coordinate structure
   !> @param[in] td_coord1 : fine grid coordinate structure
   !> @return coarse grid indices (/ (/imin0, imax0/), (/jmin0, jmax0/) /)
   !> @todo
   !> - use domain instead of start count
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid_get_coarse_index_ff( td_coord0, td_coord1, &
   &                                  id_rho )
      IMPLICIT NONE
      ! Argument
      TYPE(TFILE), INTENT(IN) :: td_coord0
      TYPE(TFILE), INTENT(IN) :: td_coord1
      INTEGER(i4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho

      ! function
      INTEGER(i4), DIMENSION(2,2,2) :: grid_get_coarse_index_ff

      ! local variable
      TYPE(TVAR)  :: tl_lon0
      TYPE(TVAR)  :: tl_lat0
      TYPE(TVAR)  :: tl_lon1
      TYPE(TVAR)  :: tl_lat1

      INTEGER(i4) , DIMENSION(:), ALLOCATABLE :: il_rho

      INTEGER(i4), DIMENSION(ip_maxdim) :: il_start
      INTEGER(i4), DIMENSION(ip_maxdim) :: il_count
      INTEGER(i4), DIMENSION(2)         :: il_xghost0
      INTEGER(i4), DIMENSION(2)         :: il_xghost1

      INTEGER(i4) :: il_imin0
      INTEGER(i4) :: il_imax0
      INTEGER(i4) :: il_jmin0
      INTEGER(i4) :: il_jmax0

      INTEGER(i4) :: il_imin1
      INTEGER(i4) :: il_imax1
      INTEGER(i4) :: il_jmin1
      INTEGER(i4) :: il_jmax1

      ! loop indices
      !----------------------------------------------------------------

      ! init
      grid_get_coarse_index_ff(:,:,:)=0

      ALLOCATE(il_rho(ig_ndim))
      il_rho(:)=1
      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)

      IF( td_coord0%i_id == 0 .OR. td_coord1%i_id == 0 )THEN
         CALL logger_error("GRID GET COARSE INDEX: can not get corase "//&
         &  "grid indices. file "//TRIM(td_coord0%c_name)//" and/or "//&
         &   TRIM(td_coord1%c_name)//" not opened." )
      ELSE
         !1- Coarse grid
         ! read coarse longitue and latitude
         tl_lon0=iom_read_var(td_coord0,'longitude')
         tl_lat0=iom_read_var(td_coord0,'latitude')

         ! get ghost cell factor on coarse grid
         il_xghost0(:)=grid_get_ghost( tl_lon0, tl_lat0 )

         il_imin0=1+il_xghost0(1)*ig_ghost
         il_jmin0=1+il_xghost0(2)*ig_ghost

         il_imax0=tl_lon0%t_dim(1)%i_len-il_xghost0(1)*ig_ghost
         il_jmax0=tl_lon0%t_dim(2)%i_len-il_xghost0(2)*ig_ghost

         CALL var_clean(tl_lon0)
         CALL var_clean(tl_lat0)

         ! read coarse longitue and latitude without ghost cell
         il_start(:)=(/il_imin0,il_jmin0,1,1/)
         il_count(:)=(/il_imax0-il_imin0+1, &
         &             il_jmax0-il_jmin0+1, &
         &             tl_lon0%t_dim(3)%i_len, &
         &             tl_lon0%t_dim(4)%i_len /)

         tl_lon0=iom_read_var(td_coord0,'longitude',il_start(:), il_count(:))
         tl_lat0=iom_read_var(td_coord0,'latitude' ,il_start(:), il_count(:))
         
         !2- Fine grid
         ! read fine longitue and latitude
         tl_lon1=iom_read_var(td_coord1,'longitude')
         tl_lat1=iom_read_var(td_coord1,'latitude')

         ! get ghost cell factor on fine grid
         il_xghost1(:)=grid_get_ghost( tl_lon1, tl_lat1 )

         il_imin1=1+il_xghost1(1)*ig_ghost
         il_jmin1=1+il_xghost1(2)*ig_ghost

         il_imax1=tl_lon1%t_dim(1)%i_len-il_xghost1(1)*ig_ghost
         il_jmax1=tl_lon1%t_dim(2)%i_len-il_xghost1(2)*ig_ghost

         CALL var_clean(tl_lon1)
         CALL var_clean(tl_lat1)

         ! read fine longitue and latitude without ghost cell
         il_start(:)=(/il_imin1,il_jmin1,1,1/)
         il_count(:)=(/il_imax1-il_imin1+1, &
         &             il_jmax1-il_jmin1+1, &
         &             tl_lon1%t_dim(3)%i_len, &
         &             tl_lon1%t_dim(4)%i_len /)

         tl_lon1=iom_read_var(td_coord1,'longitude',il_start(:), il_count(:))

         tl_lat1=iom_read_var(td_coord1,'latitude' ,il_start(:), il_count(:))
 
         !3- compute
         
         grid_get_coarse_index_ff(:,:,:)=grid_get_coarse_index(tl_lon0,tl_lat0,&
         &                                                     tl_lon1,tl_lat1,&
         &                                                     il_rho(:) )

         il_imin0=grid_get_coarse_index_ff(1,1,1)-il_xghost0(1)*ig_ghost
         il_imax0=grid_get_coarse_index_ff(1,2,1)+il_xghost0(1)*ig_ghost
         il_jmin0=grid_get_coarse_index_ff(2,1,1)-il_xghost0(2)*ig_ghost
         il_jmax0=grid_get_coarse_index_ff(2,2,1)+il_xghost0(2)*ig_ghost

         grid_get_coarse_index_ff(1,1,1)=il_imin0
         grid_get_coarse_index_ff(1,2,1)=il_imax0
         grid_get_coarse_index_ff(2,1,1)=il_jmin0
         grid_get_coarse_index_ff(2,2,1)=il_jmax0

         CALL var_clean(tl_lon0)
         CALL var_clean(tl_lat0)         
         CALL var_clean(tl_lon1)
         CALL var_clean(tl_lat1)         

      ENDIF

   END FUNCTION grid_get_coarse_index_ff
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This function get closest coarse grid indices of fine grid domain.
   !
   !> @details
   !>
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] td_longitude0 : coarse grid longitude
   !> @param[in] td_latitude0  : coarse grid latitude
   !> @param[in] td_coord1 : fine grid coordinate structure
   !> @return coarse grid indices (/ (/imin0, imax0/), (/jmin0, jmax0/) /)
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid_get_coarse_index_cf( td_lon0, td_lat0, td_coord1, &
   &                                  id_rho )
      IMPLICIT NONE
      ! Argument
      TYPE(TVAR ), INTENT(IN) :: td_lon0
      TYPE(TVAR ), INTENT(IN) :: td_lat0
      TYPE(TFILE), INTENT(IN) :: td_coord1
      INTEGER(i4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho

      ! function
      INTEGER(i4), DIMENSION(2,2,2) :: grid_get_coarse_index_cf

      ! local variable
      TYPE(TVAR)  :: tl_lon1
      TYPE(TVAR)  :: tl_lat1

      INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho

      INTEGER(i4), DIMENSION(ip_maxdim) :: il_start
      INTEGER(i4), DIMENSION(ip_maxdim) :: il_count
      INTEGER(i4), DIMENSION(2)         :: il_xghost

      INTEGER(i4) :: il_imin1
      INTEGER(i4) :: il_imax1
      INTEGER(i4) :: il_jmin1
      INTEGER(i4) :: il_jmax1

      ! loop indices
      !----------------------------------------------------------------

      ! init
      grid_get_coarse_index_cf(:,:,:)=0

      ALLOCATE(il_rho(ig_ndim) )
      il_rho(:)=1
      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)

      IF( td_coord1%i_id == 0 )THEN
         CALL logger_error("GRID GET COARSE INDEX: file "//&
         &   TRIM(td_coord1%c_name)//" not opened." )

      ELSE IF( .NOT. ASSOCIATED(td_lon0%d_value) .OR. &
      &        .NOT. ASSOCIATED(td_lat0%d_value) )THEN

         CALL logger_error("GRID GET COARSE INDEX: some coarse grid"//&
         &                 " coordinate value are not associated.")         

      ELSE
         
         !1- Fine grid
         ! read fine longitue and latitude
         tl_lon1=iom_read_var(td_coord1,'longitude')
         tl_lat1=iom_read_var(td_coord1,'latitude')

         ! get ghost cell factor on fine grid
         il_xghost(:)=grid_get_ghost( tl_lon1, tl_lat1 )

         il_imin1=1+il_xghost(1)*ig_ghost
         il_jmin1=1+il_xghost(2)*ig_ghost

         il_imax1=tl_lon1%t_dim(1)%i_len-il_xghost(1)*ig_ghost
         il_jmax1=tl_lon1%t_dim(2)%i_len-il_xghost(2)*ig_ghost

         CALL var_clean(tl_lon1)
         CALL var_clean(tl_lat1)

         ! read fine longitue and latitude without ghost cell
         il_start(:)=(/il_imin1,il_jmin1,1,1/)
         il_count(:)=(/il_imax1-il_imin1+1, &
         &             il_jmax1-il_jmin1+1, &
         &             tl_lon1%t_dim(3)%i_len, &
         &             tl_lon1%t_dim(4)%i_len /)

         tl_lon1=iom_read_var(td_coord1,'longitude',il_start(:), il_count(:))
         tl_lat1=iom_read_var(td_coord1,'latitude' ,il_start(:), il_count(:))
         
         !3- compute
         grid_get_coarse_index_cf(:,:,:)=grid_get_coarse_index(td_lon0,td_lat0,&
         &                                                     tl_lon1,tl_lat1,&
         &                                                     il_rho(:) )

         CALL var_clean(tl_lon1)
         CALL var_clean(tl_lat1)         

      ENDIF

   END FUNCTION grid_get_coarse_index_cf
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This function get closest coarse grid indices of fine grid domain.
   !
   !> @details
   !>
   !> @warning use ghost cell so can not be used on extracted domain without
   !> ghost cell
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] td_coord0 : coarse grid coordinate structure
   !> @param[in] td_lon1 : fine grid longitude
   !> @param[in] td_lat1 : fine grid latitude
   !> @return coarse grid indices (/ (/imin0, imax0/), (/jmin0, jmax0/) /)
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid_get_coarse_index_fc( td_coord0, td_lon1, td_lat1, &
   &                                  id_rho )
      IMPLICIT NONE
      ! Argument
      TYPE(TFILE), INTENT(IN) :: td_coord0
      TYPE(TVAR ), INTENT(IN) :: td_lon1
      TYPE(TVAR ), INTENT(IN) :: td_lat1
      INTEGER(i4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho

      ! function
      INTEGER(i4), DIMENSION(2,2,2) :: grid_get_coarse_index_fc

      ! local variable
      TYPE(TVAR)  :: tl_lon0
      TYPE(TVAR)  :: tl_lat0

      INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho

      INTEGER(i4), DIMENSION(ip_maxdim) :: il_start
      INTEGER(i4), DIMENSION(ip_maxdim) :: il_count
      INTEGER(i4), DIMENSION(2)         :: il_xghost

      INTEGER(i4) :: il_imin0
      INTEGER(i4) :: il_imax0
      INTEGER(i4) :: il_jmin0
      INTEGER(i4) :: il_jmax0


      ! loop indices
      !----------------------------------------------------------------

      ! init
      grid_get_coarse_index_fc(:,:,:)=0

      ALLOCATE(il_rho(ig_ndim))
      il_rho(:)=1
      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)

      IF( td_coord0%i_id == 0 )THEN
         CALL logger_error("GRID GET COARSE INDEX: file "//&
         &                 TRIM(td_coord0%c_name)//" not opened." )

      ELSE IF( .NOT. ASSOCIATED(td_lon1%d_value) .OR. &
      &        .NOT. ASSOCIATED(td_lat1%d_value) )THEN

         CALL logger_error("GRID GET COARSE INDEX: some fine grid"//&
         &                 " coordinate value are not associated.")

      ELSE
         ! read coarse longitue and latitude
         tl_lon0=iom_read_var(td_coord0,'longitude')
         tl_lat0=iom_read_var(td_coord0,'latitude')

         ! get ghost cell factor on coarse grid
         il_xghost(:)=grid_get_ghost( tl_lon0, tl_lat0 )

         il_imin0=1+il_xghost(1)*ig_ghost
         il_jmin0=1+il_xghost(2)*ig_ghost

         il_imax0=tl_lon0%t_dim(1)%i_len-il_xghost(1)*ig_ghost
         il_jmax0=tl_lon0%t_dim(2)%i_len-il_xghost(2)*ig_ghost

         CALL var_clean(tl_lon0)
         CALL var_clean(tl_lat0)

         ! read coarse longitue and latitude without ghost cell
         il_start(:)=(/il_imin0,il_jmin0,1,1/)
         il_count(:)=(/il_imax0-il_imin0+1, &
         &             il_jmax0-il_jmin0+1, &
         &             tl_lon0%t_dim(3)%i_len, &
         &             tl_lon0%t_dim(4)%i_len /)

         tl_lon0=iom_read_var(td_coord0,'longitude',il_start(:), il_count(:))
         tl_lat0=iom_read_var(td_coord0,'latitude' ,il_start(:), il_count(:))

         grid_get_coarse_index_fc(:,:,:)=grid_get_coarse_index(tl_lon0,tl_lat0,&
         &                                                     td_lon1,td_lat1,&
         &                                                     il_rho(:) )

         ! remove ghost cell
         il_imin0=grid_get_coarse_index_fc(1,1,1)+il_xghost(1)*ig_ghost
         il_imax0=grid_get_coarse_index_fc(1,2,1)+il_xghost(1)*ig_ghost
         il_jmin0=grid_get_coarse_index_fc(2,1,1)+il_xghost(2)*ig_ghost
         il_jmax0=grid_get_coarse_index_fc(2,2,1)+il_xghost(2)*ig_ghost

         grid_get_coarse_index_fc(1,1,1)=il_imin0
         grid_get_coarse_index_fc(1,2,1)=il_imax0
         grid_get_coarse_index_fc(2,1,1)=il_jmin0
         grid_get_coarse_index_fc(2,2,1)=il_jmax0

         CALL var_clean(tl_lon0)
         CALL var_clean(tl_lat0)

      ENDIF

   END FUNCTION grid_get_coarse_index_fc
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This function get closest coarse grid indices of fine grid domain.
   !
   !> @details
   !>
   !> @warning use ghost cell so can not be used on extracted domain without
   !> ghost cell
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] td_lon0 : coarse grid longitude
   !> @param[in] td_lat0 : coarse grid latitude
   !> @param[in] td_lon1 : fine grid longitude
   !> @param[in] td_lat1 : fine grid latitude
   !> @return coarse grid indices (/ (/imin0, imax0/), (/jmin0, jmax0/) /)
   !>
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid_get_coarse_index_cc( td_lon0, td_lat0, td_lon1, td_lat1, &
   &                                  id_rho )
      IMPLICIT NONE
      ! Argument
      TYPE(TVAR) , INTENT(IN) :: td_lon0
      TYPE(TVAR) , INTENT(IN) :: td_lat0
      TYPE(TVAR) , INTENT(IN) :: td_lon1
      TYPE(TVAR) , INTENT(IN) :: td_lat1
      INTEGER(i4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho

      ! function
      INTEGER(i4), DIMENSION(2,2,2) :: grid_get_coarse_index_cc

      ! local variable
      REAL(dp)    :: dl_lon1_ll
      REAL(dp)    :: dl_lon1_ul
      REAL(dp)    :: dl_lon1_lr
      REAL(dp)    :: dl_lon1_ur

      REAL(dp)    :: dl_lat1_ll
      REAL(dp)    :: dl_lat1_ul
      REAL(dp)    :: dl_lat1_lr
      REAL(dp)    :: dl_lat1_ur

      REAL(dp)    :: dl_dlon
      REAL(dp)    :: dl_dlat

      INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho

      INTEGER(i4), DIMENSION(2) :: il_ill
      INTEGER(i4), DIMENSION(2) :: il_ilr
      INTEGER(i4), DIMENSION(2) :: il_iul
      INTEGER(i4), DIMENSION(2) :: il_iur

      INTEGER(i4) :: il_ew0 
      INTEGER(i4) :: il_imin0
      INTEGER(i4) :: il_imax0
      INTEGER(i4) :: il_jmin0
      INTEGER(i4) :: il_jmax0

      INTEGER(i4) :: il_ew1 
      INTEGER(i4) :: il_imin1
      INTEGER(i4) :: il_imax1
      INTEGER(i4) :: il_jmin1
      INTEGER(i4) :: il_jmax1

      INTEGER(i4) :: il_imin
      INTEGER(i4) :: il_imax
      INTEGER(i4) :: il_jmin
      INTEGER(i4) :: il_jmax      

      INTEGER(i4), DIMENSION(2,2) :: il_offset

      ! loop indices
      INTEGER(i4) :: ji
      INTEGER(i4) :: jj
      !----------------------------------------------------------------

      ! init
      grid_get_coarse_index_cc(:,:,:)=0

      ALLOCATE( il_rho(ig_ndim) )
      il_rho(:)=1
      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)

      IF( .NOT. ASSOCIATED(td_lon0%d_value) .OR. &
      &   .NOT. ASSOCIATED(td_lat0%d_value) .OR. &
      &   .NOT. ASSOCIATED(td_lon1%d_value) .OR. &
      &   .NOT. ASSOCIATED(td_lat1%d_value) )THEN
         CALL logger_error("GRID GET COARSE INDEX: some fine or coarse grid"//&
         &                 " coordinate value not associated.")
      ELSE

         IF( grid_is_global(td_lon1, td_lat1) )THEN

            IF( grid_is_global(td_lon0, td_lat0) )THEN
               CALL logger_trace("GRID GET COARSE INDEX: fine grid is global ")
               grid_get_coarse_index_cc(:,:,1) = 1
               grid_get_coarse_index_cc(:,:,2) = 0
            ELSE
               CALL logger_error("GRID GET COARSE INDEX: fine grid is "//&
               &                 "global, coarse grid not.")
            ENDIF

         ELSE

            ! "global" coarse grid indice
            il_imin0=1
            il_jmin0=1

            il_imax0=td_lon0%t_dim(1)%i_len
            il_jmax0=td_lon0%t_dim(2)%i_len

            ! get east west overlap for coarse grid
            il_ew0=dom_get_ew_overlap(td_lon0)
            IF( il_ew0 >= 0 )THEN
               ! last point before overlap
               il_imax0=il_imax0-il_ew0
            ENDIF

            ! "global" fine grid indice
            il_imin1=1
            il_jmin1=1

            il_imax1=td_lon1%t_dim(1)%i_len
            il_jmax1=td_lon1%t_dim(2)%i_len

            ! get east west overlap for coarse grid
            il_ew1=dom_get_ew_overlap(td_lon1)
            IF( il_ew1 >= 0 )THEN
               ! last point before overlap
               il_imax1=il_imax1-il_ew1
            ENDIF

            ! get indices for each corner
            !1- search lower left corner indices
            dl_lon1_ll=td_lon1%d_value( il_imin1, il_jmin1, 1, 1 )
            dl_lat1_ll=td_lat1%d_value( il_imin1, il_jmin1, 1, 1 )

            dl_dlon=ABS(td_lon1%d_value(il_imin1+1,il_jmin1  ,1,1)-dl_lon1_ll)
            dl_dlat=ABS(td_lat1%d_value(il_imin1  ,il_jmin1+1,1,1)-dl_lat1_ll)

!            CALL logger_debug("GRID GET COARSE INDEX: lon1 ll "//&
!            &  TRIM(fct_str(dl_lon1_ll)) )
!            CALL logger_debug("GRID GET COARSE INDEX: lat1 ll "//&
!            &  TRIM(fct_str(dl_lat1_ll)) )
!
!            CALL logger_debug("GRID GET COARSE INDEX: lon0 min "//&
!            &  TRIM(fct_str(minval(td_lon0%d_value(2:,2:,:,:)))) )
!            CALL logger_debug("GRID GET COARSE INDEX: lon0 max "//&
!            &  TRIM(fct_str(maxval(td_lon0%d_value(2:,2:,:,:)))) )
!
!            CALL logger_debug("GRID GET COARSE INDEX: lat0 min "//&
!            &  TRIM(fct_str(minval(td_lat0%d_value(2:,2:,:,:)))) )
!            CALL logger_debug("GRID GET COARSE INDEX: lat0 max "//&
!            &  TRIM(fct_str(maxval(td_lat0%d_value(2:,2:,:,:)))) )

            ! look for closest point on coarse grid
            il_ill(:)= grid_get_closest(td_lon0%d_value(il_imin0:il_imax0, &
            &                                           il_jmin0:il_jmax0, &
            &                                           1,1), &
            &                           td_lat0%d_value(il_imin0:il_imax0, &
            &                                           il_jmin0:il_jmax0, &
            &                                           1,1), &
            &                           dl_lon1_ll, dl_lat1_ll   )

            ! coarse grid point should be south west of fine grid domain
            ji = il_ill(1)
            jj = il_ill(2)

            IF( ABS(td_lon0%d_value(ji,jj,1,1)-dl_lon1_ll) > dl_dlon*1.e-3 )THEN
               IF(td_lon0%d_value(ji,jj,1,1) > dl_lon1_ll ) il_ill(1)=il_ill(1)-1
            ENDIF
            IF( ABS(td_lat0%d_value(ji,jj,1,1)-dl_lat1_ll) > dl_dlat*1.e-3 )THEN
               IF(td_lat0%d_value(ji,jj,1,1) > dl_lat1_ll ) il_ill(2)=il_ill(2)-1
            ENDIF

            !2- search upper left corner indices
            dl_lon1_ul=td_lon1%d_value( il_imin1, il_jmax1, 1, 1 )
            dl_lat1_ul=td_lat1%d_value( il_imin1, il_jmax1, 1, 1 )

            dl_dlon=ABS(td_lon1%d_value(il_imin1+1,il_jmax1  ,1,1)-dl_lon1_ll)
            dl_dlat=ABS(td_lat1%d_value(il_imin1  ,il_jmax1-1,1,1)-dl_lat1_ll)
            
            ! look for closest point on coarse grid
            il_iul(:)= grid_get_closest(td_lon0%d_value(il_imin0:il_imax0, &
            &                                           il_jmin0:il_jmax0, &
            &                                           1,1), &
            &                           td_lat0%d_value(il_imin0:il_imax0, &
            &                                           il_jmin0:il_jmax0, &
            &                                           1,1), &
            &                           dl_lon1_ul, dl_lat1_ul   )

            ! coarse grid point should be north west of fine grid domain
            ji = il_iul(1)
            jj = il_iul(2)

            IF( ABS(td_lon0%d_value(ji,jj,1,1)-dl_lon1_ul) > dl_dlon*1.e-3 )THEN
               IF(td_lon0%d_value(ji,jj,1,1) > dl_lon1_ul ) il_iul(1)=il_iul(1)-1
            ENDIF
            IF( ABS(td_lat0%d_value(ji,jj,1,1)-dl_lat1_ul) > dl_dlat*1.e-3 )THEN
               IF(td_lat0%d_value(ji,jj,1,1) < dl_lat1_ul ) il_iul(2)=il_iul(2)+1
            ENDIF

            !3- search lower right corner indices
            dl_lon1_lr=td_lon1%d_value( il_imax1, il_jmin1, 1, 1 )
            dl_lat1_lr=td_lat1%d_value( il_imax1, il_jmin1, 1, 1 )

            dl_dlon=ABS(td_lon1%d_value(il_imax1-1,il_jmin1  ,1,1)-dl_lon1_ll)
            dl_dlat=ABS(td_lat1%d_value(il_imax1  ,il_jmin1+1,1,1)-dl_lat1_ll)
            
            ! look for closest point on coarse grid
            il_ilr(:)= grid_get_closest(td_lon0%d_value(il_imin0:il_imax0, &
            &                                           il_jmin0:il_jmax0, &
            &                                           1,1), &
            &                           td_lat0%d_value(il_imin0:il_imax0, &
            &                                           il_jmin0:il_jmax0, &
            &                                           1,1), &
            &                           dl_lon1_lr, dl_lat1_lr   )

            ! coarse grid point should be south east of fine grid domain
            ji = il_ilr(1)
            jj = il_ilr(2)
            IF( ABS(td_lon0%d_value(ji,jj,1,1)-dl_lon1_lr) > dl_dlon*1.e-3 )THEN
               IF( td_lon0%d_value(ji,jj,1,1) < dl_lon1_lr ) il_ilr(1)=il_ilr(1)+1
            ENDIF
            IF( ABS(td_lat0%d_value(ji,jj,1,1)-dl_lat1_lr) > dl_dlat*1.e-3 )THEN
               IF( td_lat0%d_value(ji,jj,1,1) > dl_lat1_lr ) il_ilr(2)=il_ilr(2)-1
            ENDIF

            !4- search upper right corner indices
            dl_lon1_ur=td_lon1%d_value( il_imax1, il_jmax1, 1, 1 )
            dl_lat1_ur=td_lat1%d_value( il_imax1, il_jmax1, 1, 1 )

            dl_dlon=ABS(td_lon1%d_value(il_imax1-1,il_jmax1  ,1,1)-dl_lon1_ll)
            dl_dlat=ABS(td_lat1%d_value(il_imax1  ,il_jmax1-1,1,1)-dl_lat1_ll)
            
            ! look for closest point on coarse grid
            il_iur(:)= grid_get_closest(td_lon0%d_value(il_imin0:il_imax0, &
            &                                           il_jmin0:il_jmax0, &
            &                                           1,1), &
            &                           td_lat0%d_value(il_imin0:il_imax0, &
            &                                           il_jmin0:il_jmax0, &
            &                                           1,1), &
            &                           dl_lon1_ur, dl_lat1_ur   )

            ! coarse grid point should be north east fine grid domain
            ji = il_iur(1)
            jj = il_iur(2)
            IF( ABS(td_lon0%d_value(ji,jj,1,1)-dl_lon1_ur) > dl_dlon*1.e-3 )THEN
               IF( td_lon0%d_value(ji,jj,1,1) < dl_lon1_ur ) il_iur(1)=il_iur(1)+1
            ENDIF
            IF( ABS(td_lat0%d_value(ji,jj,1,1)-dl_lat1_ur) > dl_dlat*1.e-3 )THEN
               IF( td_lat0%d_value(ji,jj,1,1) < dl_lat1_ur ) il_iur(2)=il_iur(2)+1
            ENDIF

            ! coarse grid indices
            il_imin = il_imin0-1+MIN(il_ill(1), il_iul(1))
            il_imax = il_imin0-1+MAX(il_ilr(1), il_iur(1))

            IF( il_imax <= il_ew0 )THEN
               il_imax = td_lon0%t_dim(1)%i_len - il_ew0 + il_imax 
            ENDIF

            il_jmin = il_jmin0-1+MIN(il_ill(2), il_ilr(2))
            il_jmax = il_jmin0-1+MAX(il_iul(2), il_iur(2))

            il_offset(:,:)= grid_get_fine_offset( td_lon0%d_value( :,:,1,1 ), &
            &                                     td_lat0%d_value( :,:,1,1 ), &
            &                                     il_imin, il_jmin, &
            &                                     il_imax, il_jmax, &
            &                                     td_lon1%d_value( :,:,1,1 ), &
            &                                     td_lat1%d_value( :,:,1,1 ), &
            &                                     il_rho(:) )

            grid_get_coarse_index_cc(1,1,2) = il_offset(1,1)
            grid_get_coarse_index_cc(1,2,2) = il_offset(1,2)

            grid_get_coarse_index_cc(2,1,2) = il_offset(2,1)
            grid_get_coarse_index_cc(2,2,2) = il_offset(2,2) 

            ! special case if east west overlap
            IF( il_ew1 >= 0 )THEN
               CALL logger_debug("GRID GET COARSE INDEX: East-West overlap "//&
               &                 "found for fine grid " )

               il_imin = 1
               il_imax = 1

               grid_get_coarse_index_cc(1,1,2) = 0
               grid_get_coarse_index_cc(1,2,2) = 0
            ENDIF

         ENDIF

         IF( il_imin == il_imax ) il_imax=td_lon0%t_dim(1)%i_len
         IF( il_jmin == il_jmax ) il_jmax=td_lon0%t_dim(2)%i_len

         grid_get_coarse_index_cc(1,1,1) = il_imin
         grid_get_coarse_index_cc(1,2,1) = il_imax

         grid_get_coarse_index_cc(2,1,1) = il_jmin
         grid_get_coarse_index_cc(2,2,1) = il_jmax
 
      ENDIF

   END FUNCTION grid_get_coarse_index_cc
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This function check if grid is global or not
   !
   !> @details
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] td_lon : longitude structure 
   !> @param[in] td_lat : latitude structure 
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid_is_global(td_lon, td_lat)
      IMPLICIT NONE
      ! Argument      
      TYPE(TVAR), INTENT(IN) :: td_lon
      TYPE(TVAR), INTENT(IN) :: td_lat

      ! function
      LOGICAL :: grid_is_global
      
      ! local variable
      INTEGER(i4)               :: il_ew
      INTEGER(i4)               :: il_south
      INTEGER(i4)               :: il_north

      REAL(dp)                  :: dl_lat_min
      REAL(dp)                  :: dl_lat_max

      ! loop indices
      !----------------------------------------------------------------

      ! init
      grid_is_global=.FALSE.

      IF( ANY( td_lon%t_dim(:)%i_len /= td_lat%t_dim(:)%i_len )  )THEN
         CALL logger_fatal("GRID IS GLOBAL: dimension of longitude and "//&
         &              " latitude differ")
      ENDIF

      IF( .NOT. ASSOCIATED(td_lon%d_value) .OR. &
      &   .NOT. ASSOCIATED(td_lat%d_value) )THEN
         CALL logger_error("GRID IS GLOBAL: na value associated to "//&
         &              " longitude or latitude strucutre")      
      ELSE

         il_south=1
         il_north=td_lon%t_dim(2)%i_len

         dl_lat_min=MINVAL(td_lat%d_value(:,il_south,1,1))
         dl_lat_max=MAXVAL(td_lat%d_value(:,il_north,1,1))

         IF( dl_lat_min < -77.0 .AND. dl_lat_max > 89.0 )THEN

            il_ew=td_lon%i_ew
            IF( il_ew >= 0 )THEN

               grid_is_global=.TRUE.

            ENDIF

         ENDIF
      ENDIF

   END FUNCTION grid_is_global
   !> @endcode

   !-------------------------------------------------------------------
   !> @brief This function return coarse grid indices of the closest point
   !> from fine grid point (lon1,lat1) 
   !> 
   !
   !> @details
   !
   !> @note overlap band should have been already removed from coarse grid table 
   !> of longitude and latitude, before running this function
   !>
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] dd_lon0 : coarse grid table of longitude
   !> @param[in] dd_lat0 : coarse grid table of latitude
   !> @param[in] dd_lon1 : fine   grid longitude
   !> @param[in] dd_lat1 : fine   grid latitude
   !> @return coarse grid indices of closest point of fine grid point
   !>
   !> @todo 
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid_get_closest( dd_lon0, dd_lat0, dd_lon1, dd_lat1  )
      IMPLICIT NONE
      ! Argument
      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon0
      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat0
      REAL(dp),                 INTENT(IN) :: dd_lon1
      REAL(dp),                 INTENT(IN) :: dd_lat1

      ! function
      INTEGER(i4), DIMENSION(2) :: grid_get_closest

      ! local variable
      INTEGER(i4)                              :: il_iinf
      INTEGER(i4)                              :: il_imid
      INTEGER(i4)                              :: il_isup
      INTEGER(i4)                              :: il_jinf
      INTEGER(i4)                              :: il_jmid
      INTEGER(i4)                              :: il_jsup
      INTEGER(i4), DIMENSION(2)                :: il_shape
      INTEGER(i4), DIMENSION(1)                :: il_ind
    
      LOGICAL                                  :: ll_north
      LOGICAL                                  :: ll_continue

      REAL(dp)                                 :: dl_lon1
      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_dist
      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon0

      ! loop indices
      !----------------------------------------------------------------

      IF( ANY( SHAPE(dd_lon0(:,:)) /= SHAPE(dd_lat0(:,:)) ) )THEN
         CALL logger_fatal("GRID GET CLOSEST: dimension of longitude and "//&
         &              " latitude differ")
      ENDIF

      il_shape(:)=SHAPE(dd_lon0(:,:))
      
      ALLOCATE( dl_lon0(il_shape(1),il_shape(2)) ) 
      
      dl_lon0(:,:) = dd_lon0(:,:)
      WHERE(dd_lon0(:,:) < 0 ) dl_lon0(:,:) = dd_lon0(:,:) + 360.

      dl_lon1 = dd_lon1
      IF( dd_lon1 < 0 ) dl_lon1 = dd_lon1 + 360.

      !1- first, use dichotomy to reduce domain
      il_iinf = 1              ; il_jinf = 1
      il_isup = il_shape(1)    ; il_jsup = il_shape(2)

      il_shape(1)= il_isup - il_iinf + 1
      il_shape(2)= il_jsup - il_jinf + 1

      ll_north=.FALSE.
      ll_continue=.TRUE.

      !1-1 look for meridian 0°/360°
      il_jmid = il_jinf + INT(il_shape(2)/2)
      il_ind(:) = MAXLOC( dl_lon0(:,il_jmid), dl_lon0(:,il_jmid) <= 360._dp )

      il_imid=il_ind(1)

      IF( dl_lon1 == dl_lon0(il_imid,il_jmid) .AND. &
      &   dd_lat1 == dd_lat0(il_imid,il_jmid) )THEN

         il_iinf = il_imid ;  il_isup = il_imid
         il_jinf = il_jmid ;  il_jsup = il_jmid

         ll_continue=.FALSE.

      ELSE
         IF( dl_lon1 < dl_lon0(il_isup,il_jmid) .AND. &
         &   il_imid /= il_isup )THEN

            ! point east
            il_iinf = il_imid
      
         ELSE IF( dl_lon1 > dl_lon0(il_iinf,il_jmid) .AND. &
         &        il_imid /= il_iinf )THEN

            ! point west
            il_isup = il_imid

         ENDIF

         il_shape(1)= il_isup - il_iinf + 1
         il_shape(2)= il_jsup - il_jinf + 1

         il_imid = il_iinf + INT(il_shape(1)/2) 
         il_jmid = il_jinf + INT(il_shape(2)/2)

         ! exit if too close from north fold (safer)
         IF( dd_lat0(il_imid,il_jmid) > 50.0 ) ll_north=.TRUE.

         ! exit when close enough of point
         IF( ANY(il_shape(:) < 10 ) ) ll_continue=.FALSE.
      ENDIF

      !1-2
      DO WHILE( ll_continue .AND. .NOT. ll_north )

         IF( dl_lon1 == dl_lon0(il_imid,il_jmid) .AND. &
         &   dd_lat1 == dd_lat0(il_imid,il_jmid) )THEN

            il_iinf = il_imid ;  il_isup = il_imid
            il_jinf = il_jmid ;  il_jsup = il_jmid

            ll_continue=.FALSE.

         ELSE
            IF( dl_lon1 > dl_lon0(il_imid,il_jmid) )THEN

               ! point east
               il_iinf = il_imid
        
            ELSE IF(dl_lon1 < dl_lon0(il_imid,il_jmid) )THEN

               ! point west
               il_isup = il_imid

            ENDIF


            IF( dd_lat1 > dd_lat0(il_imid,il_jmid) )THEN
               
               ! point north
               il_jinf = il_jmid

            ELSE IF(dd_lat1 < dd_lat0(il_imid,il_jmid) )THEN

               ! point south
               il_jsup = il_jmid
            
            ENDIF

            il_shape(1)= il_isup - il_iinf + 1
            il_shape(2)= il_jsup - il_jinf + 1

            il_imid = il_iinf + INT(il_shape(1)/2) 
            il_jmid = il_jinf + INT(il_shape(2)/2)

            ! exit if too close from north fold (safer)
            IF( dd_lat0(il_imid,il_jmid) > 50.0 ) ll_north=.TRUE.

            ! exit when close enough of point
            IF( ANY(il_shape(:) < 10 ) ) ll_continue=.FALSE.
         ENDIF
         
      ENDDO

      !2- then find closest point by computing distances
      il_shape(1)= il_isup - il_iinf + 1
      il_shape(2)= il_jsup - il_jinf + 1

      ALLOCATE( dl_dist(il_shape(1), il_shape(2)) )

      dl_dist(:,:)=grid_distance(dl_lon0(il_iinf:il_isup,il_jinf:il_jsup), &
      &                          dd_lat0(il_iinf:il_isup,il_jinf:il_jsup), &
      &                          dl_lon1, dd_lat1 )

      grid_get_closest(:)=MINLOC(dl_dist(:,:),dl_dist(:,:)/=NF90_FILL_DOUBLE)

      grid_get_closest(1)=grid_get_closest(1)+il_iinf-1
      grid_get_closest(2)=grid_get_closest(2)+il_jinf-1

      DEALLOCATE( dl_dist )
      DEALLOCATE( dl_lon0 )

   END FUNCTION grid_get_closest
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This function compute the distance between a point A and 
   !> points of a grid  
   !
   !> @details
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] dd_lon : grid longitude table
   !> @param[in] dd_lat : grid latitude  table
   !> @param[in] dd_lonA : longitude of point A
   !> @param[in] dd_latA : latitude  of point A
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid_distance(dd_lon, dd_lat, dd_lonA, dd_latA)
      IMPLICIT NONE
      ! Argument      
      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon
      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat
      REAL(dp),                 INTENT(IN) :: dd_lonA
      REAL(dp),                 INTENT(IN) :: dd_latA

      ! function
      REAL(dp), DIMENSION(SIZE(dd_lon(:,:),DIM=1),&
      &                   SIZE(dd_lon(:,:),DIM=2)) :: grid_distance

      ! local variable
      INTEGER(i4), DIMENSION(2) :: il_shape

      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon
      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lat
      REAL(dp)                                 :: dl_lonA
      REAL(dp)                                 :: dl_latA

      REAL(dp)                                 :: dl_tmp

      ! loop indices
      INTEGER(i4) :: ji
      INTEGER(i4) :: jj
      !----------------------------------------------------------------

      IF( ANY( SHAPE(dd_lon(:,:)) /= SHAPE(dd_lat(:,:)) ) )THEN
         CALL logger_fatal("GRID DISTANCE: dimension of longitude and "//&
         &              " latitude differ")
      ENDIF
      il_shape(:)=SHAPE(dd_lon(:,:))
      
      ALLOCATE(dl_lon(il_shape(1),il_shape(2)))
      ALLOCATE(dl_lat(il_shape(1),il_shape(2)))

      dl_lon(:,:) = dd_lon(:,:)
      dl_lonA     = dd_lonA

      WHERE(dd_lon(:,:) < 0 ) dl_lon(:,:) = dd_lon(:,:) + 360.
      IF(   dd_lonA     < 0 ) dl_lonA     = dd_lonA     + 360.
      
      dl_lonA = dd_lonA * dg_deg2rad
      dl_latA = dd_latA * dg_deg2rad

      dl_lon(:,:) = dl_lon(:,:) * dg_deg2rad
      dl_lat(:,:) = dd_lat(:,:) * dg_deg2rad

      grid_distance(:,:)=NF90_FILL_DOUBLE

      DO jj=1,il_shape(2)
         DO ji=1,il_shape(1)
            IF( dl_lon(ji,jj) == dl_lonA .AND. &
            &   dl_lat(ji,jj) == dl_lATA )THEN
               grid_distance(ji,jj)=0.0
            ELSE
               dl_tmp= SIN(dl_latA)*SIN(dl_lat(ji,jj)) + &
               &       COS(dl_latA)*COS(dl_lat(ji,jj))*COS(dl_lon(ji,jj)-dl_lonA)
               ! check to avoid mistake with ACOS
               IF( dl_tmp < -1.0 ) dl_tmp = -1.0
               IF( dl_tmp >  1.0 ) dl_tmp =  1.0
               grid_distance(ji,jj)=ACOS(dl_tmp)*dg_rearth
            ENDIF
         ENDDO
      ENDDO

      DEALLOCATE(dl_lon)
      DEALLOCATE(dl_lat)

   END FUNCTION grid_distance
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This function get fine grid offset.
   !
   !> @details
   !> offset value could be 0,1,..,rho-1
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] dd_lon0 : coarse grid longitude table 
   !> @param[in] dd_lat0 : coarse grid latitude  table
   !> @param[in] dd_lon1 : fine   grid longitude table 
   !> @param[in] dd_lat1 : fine   grid latitude  table
   !> @param[in] id_imin0 : coarse grid lower left corner i-indice of fine grid domain
   !> @param[in] id_jmin0 : coarse grid lower left corner j-indice of fine grid domain
   !> @param[in] id_imax0 : coarse grid upper right corner i-indice of fine grid domain
   !> @param[in] id_jmax0 : coarse grid upper right corner j-indice of fine grid domain
   !> @param[in] id_rhoi : i-direction refinement factor
   !> @param[in] id_rhoj : j-direction refinement factor
   !> @return offset table (/ (/i_offset_left,i_offset_right!/),(/j_offset_lower,j_offset_upper/) /)
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid_get_fine_offset( dd_lon0, dd_lat0, &
   &                              id_imin0, id_jmin0, id_imax0, id_jmax0, &
   &                              dd_lon1, dd_lat1, id_rho )
      IMPLICIT NONE
      ! Argument
      REAL(dp)   , DIMENSION(:,:), INTENT(IN) :: dd_lon0
      REAL(dp)   , DIMENSION(:,:), INTENT(IN) :: dd_lat0
      REAL(dp)   , DIMENSION(:,:), INTENT(IN) :: dd_lon1
      REAL(dp)   , DIMENSION(:,:), INTENT(IN) :: dd_lat1

      INTEGER(i4),                 INTENT(IN) :: id_imin0
      INTEGER(i4),                 INTENT(IN) :: id_jmin0
      INTEGER(i4),                 INTENT(IN) :: id_imax0
      INTEGER(i4),                 INTENT(IN) :: id_jmax0

      INTEGER(i4), DIMENSION(:)  , INTENT(IN) :: id_rho

      ! function
      INTEGER(i4), DIMENSION(2,2) :: grid_get_fine_offset

      ! local variable
      INTEGER(i4), DIMENSION(2) :: il_shape0
      INTEGER(i4), DIMENSION(2) :: il_shape1
      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon0
      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon1
      
      REAL(dp) :: dl_dlon
      REAL(dp) :: dl_dlat

      ! loop indices
      INTEGER(i4) :: ji
      INTEGER(i4) :: jj

      INTEGER(i4) :: ii
      INTEGER(i4) :: ij
      !----------------------------------------------------------------
      IF( ANY( SHAPE(dd_lon0(:,:)) /= SHAPE(dd_lat0(:,:)) ) )THEN
         CALL logger_fatal("GRID GET FINE OFFSET: dimension of coarse "//&
         &              "longitude and latitude differ")
      ENDIF

      IF( ANY( SHAPE(dd_lon1(:,:)) /= SHAPE(dd_lat1(:,:)) ) )THEN
         CALL logger_fatal("GRID GET FINE OFFSET: dimension of fine "//&
         &              "longitude and latitude differ")
      ENDIF      

      il_shape0(:)=SHAPE(dd_lon0(:,:))
      ALLOCATE( dl_lon0(il_shape0(1),il_shape0(2)) )

      dl_lon0(:,:)=dd_lon0(:,:)
      WHERE( dd_lon0(:,:) < 0 ) dl_lon0(:,:)=dd_lon0(:,:)+360.

      il_shape1(:)=SHAPE(dd_lon1(:,:))
      ALLOCATE( dl_lon1(il_shape1(1),il_shape1(2)) )

      dl_lon1(:,:)=dd_lon1(:,:)
      WHERE( dd_lon1(:,:) < 0 ) dl_lon1(:,:)=dd_lon1(:,:)+360.         

      grid_get_fine_offset(:,:)=-1

      ! look for i-direction left offset 
      IF( dl_lon1(1,1) < dl_lon0(id_imin0+1,id_jmin0) )THEN
         DO ji=1,id_rho(jp_I)+2
            dl_dlon=ABS(dl_lon1(ji+1,1)-dl_lon1(ji,1))*1.e-3
            IF( dl_lon1(ji,1) > dl_lon0(id_imin0+1,id_jmin0) + dl_dlon )THEN
               grid_get_fine_offset(1,1)=(id_rho(jp_I)+1)-ji+MOD(id_rho(jp_I),2)
               EXIT
            ENDIF
         ENDDO
      ELSE
         CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
         &                 " not match fine grid lower left corner.")
      ENDIF

      ! look for i-direction right offset
      IF( dl_lon1(il_shape1(1),1) > dl_lon0(id_imax0-1,id_jmin0) )THEN
         DO ji=1,id_rho(jp_I)+2
            ii=il_shape1(1)-ji+1
            dl_dlon=ABS(dl_lon1(ii,1)-dl_lon1(ii-1,1))*1.e-3
            IF( dl_lon1(ii,1) < dl_lon0(id_imax0-1,id_jmin0) - dl_dlon )THEN
               grid_get_fine_offset(1,2)=(id_rho(jp_I)+1)-ji+MOD(id_rho(jp_I),2)
               EXIT
            ENDIF
         ENDDO
      ELSE
         CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
         &                 " not match fine grid lower right corner.")
      ENDIF

      ! look for j-direction lower offset 
      IF( dd_lat1(1,1) < dd_lat0(id_imin0,id_jmin0+1) )THEN
         DO jj=1,id_rho(jp_J)+2
            dl_dlat=ABS(dd_lat1(1,jj+1)-dd_lat1(1,jj))*1.e-3
            IF( dd_lat1(1,jj) > dd_lat0(id_imin0,id_jmin0+1) + dl_dlat )THEN
               grid_get_fine_offset(2,1)=(id_rho(jp_J)+1)-jj+MOD(id_rho(jp_J),2)
               EXIT
            ENDIF
         ENDDO
      ELSE
         CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
         &                 " not match fine grid upper left corner.")
      ENDIF

      ! look for j-direction upper offset 
      IF( dd_lat1(1,il_shape1(2)) > dd_lat0(id_imin0,id_jmax0-1) )THEN
         DO jj=1,id_rho(jp_J)+2
            ij=il_shape1(2)-jj+1
            dl_dlat=ABS(dd_lat1(1,ij)-dd_lat1(1,ij-1))*1.e-3
            IF( dd_lat1(1,ij) < dd_lat0(id_imin0,id_jmax0-1) - dl_dlat )THEN
               grid_get_fine_offset(2,2)=(id_rho(jp_J)+1)-jj+MOD(id_rho(jp_J),2)
               EXIT
            ENDIF
         ENDDO
      ELSE
         CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//&
         &                 " not match fine grid upper right corner.")
      ENDIF

      DEALLOCATE( dl_lon0 )
      DEALLOCATE( dl_lon1 )

   END FUNCTION grid_get_fine_offset
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This function check if ghost cell are used or not, and return ghost
   !> cell factor (0,1) in i- and j-direction.
   !
   !> @details
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] td_lon  : grid longitude sturcture 
   !> @param[in] td_lat  : grid latitude  structure
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid__get_ghost_ll( td_lon, td_lat )
      IMPLICIT NONE
      ! Argument
      TYPE(TVAR), INTENT(IN) :: td_lon
      TYPE(TVAR), INTENT(IN) :: td_lat

      ! function
      INTEGER(i4), DIMENSION(2) :: grid__get_ghost_ll

      ! local variable
      INTEGER(i4) :: il_ew
      ! loop indices
      !----------------------------------------------------------------
      ! init
      grid__get_ghost_ll(:)=0

      IF( grid_is_global(td_lon, td_lat) )THEN
         grid__get_ghost_ll(:)=0
      ELSE
         grid__get_ghost_ll(2)=1

         il_ew=td_lon%i_ew
         IF( il_ew < 0 )THEN
            grid__get_ghost_ll(1)=1
         ELSE
            grid__get_ghost_ll(1)=0
         ENDIF
      ENDIF

   END FUNCTION grid__get_ghost_ll
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This function check if ghost cell are used or not, and return ghost
   !> cell factor (0,1) in i- and j-direction.
   !
   !> @details
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] td_file : file sturcture 
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid__get_ghost_f( td_file )
      IMPLICIT NONE
      ! Argument
      TYPE(TFILE), INTENT(IN) :: td_file

      ! function
      INTEGER(i4), DIMENSION(2) :: grid__get_ghost_f

      ! local variable
      TYPE(TVAR)  :: tl_lon
      TYPE(TVAR)  :: tl_lat

      INTEGER(i4) :: il_lonid
      INTEGER(i4) :: il_latid
      ! loop indices
      INTEGER(i4)  :: ji
      !----------------------------------------------------------------
      ! init
      grid__get_ghost_f(:)=0

      IF( td_file%i_id == 0 )THEN
         CALL logger_error("GRID GET GHOST: file "//&
         &                 TRIM(td_file%c_name)//" not opened." )

      ELSE

         IF( ASSOCIATED(td_file%t_var) )THEN
            ! read coarse longitue and latitude
            il_lonid=var_get_id(td_file%t_var(:),'longitude')
            il_latid=var_get_id(td_file%t_var(:),'latitude')

            print *,'file ',trim(td_file%c_name),td_file%i_ew
            DO ji=1,td_file%i_nvar
               print *,ji,trim(td_file%t_var(ji)%c_name),': ',td_file%t_var(ji)%i_ew
            ENDDO
            print *,'lonid ',il_lonid
            print *,'latid ',il_latid
            IF( il_lonid /=0 .AND. il_latid /= 0 )THEN
               tl_lon=iom_read_var(td_file,il_lonid)
               print *,'lon ',tl_lon%i_ew
               tl_lat=iom_read_var(td_file,il_latid)
               print *,'lat ',tl_lat%i_ew
               ! get ghost cell factor on coarse grid
               grid__get_ghost_f(:)=grid_get_ghost( tl_lon, tl_lat )
            ELSE
               CALL logger_error("GRID GET GHOST: can not find "//&
               &           "longitude or latitude "//&
               &           "in file "//TRIM(td_file%c_name))
            ENDIF
         ELSE
           CALL logger_error("GRID GET GHOST: no variable "//&
           &           "associated to file "//TRIM(td_file%c_name))
         ENDIF

      ENDIF

   END FUNCTION grid__get_ghost_f
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This subroutine check fine and coarse grid coincidence
   !
   !> @details
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] td_coord0 : coarse grid coordinate file structure 
   !> @param[in] td_coord1 : fine   grid coordinate file structure 
   !> @param[in] id_imin0 : coarse grid lower left  corner i-indice of fine grid domain 
   !> @param[in] id_imax0 : coarse grid upper right corner i-indice of fine grid domain
   !> @param[in] id_jmin0 : coarse grid lower left  corner j-indice of fine grid domain 
   !> @param[in] id_jmax0 : coarse grid upper right corner j-indice of fine grid domain  
   !> @param[in] id_rho   : table of refinement factor 
   !-------------------------------------------------------------------
   !> @code
   SUBROUTINE grid_check_coincidence( td_coord0, td_coord1, &
   &                                  id_imin0, id_imax0, &
   &                                  id_jmin0, id_jmax0, &
   &                                  id_rho )
      IMPLICIT NONE
      
      ! Argument      
      TYPE(TFILE), INTENT(IN) :: td_coord0
      TYPE(TFILE), INTENT(IN) :: td_coord1
      INTEGER(i4), INTENT(IN) :: id_imin0
      INTEGER(i4), INTENT(IN) :: id_imax0
      INTEGER(i4), INTENT(IN) :: id_jmin0
      INTEGER(i4), INTENT(IN) :: id_jmax0
      INTEGER(i4), DIMENSION(:), INTENT(IN) :: id_rho

      ! local variable
      INTEGER(i4) :: il_imid1
      INTEGER(i4) :: il_jmid1
      
      INTEGER(i4) :: il_ew0
      INTEGER(i4) :: il_ew1

      INTEGER(i4) :: il_imin1
      INTEGER(i4) :: il_imax1
      INTEGER(i4) :: il_jmin1
      INTEGER(i4) :: il_jmax1

      INTEGER(i4), DIMENSION(2) :: il_indC
      INTEGER(i4), DIMENSION(2) :: il_indF
      INTEGER(i4), DIMENSION(2) :: il_iind
      INTEGER(i4), DIMENSION(2) :: il_jind

      REAL(dp)    :: dl_lon0
      REAL(dp)    :: dl_lat0
      REAL(dp)    :: dl_lon1
      REAL(dp)    :: dl_lat1

      REAL(dp)    :: dl_lon1p
      REAL(dp)    :: dl_lat1p

      REAL(dp)    :: dl_dlon
      REAL(dp)    :: dl_dlat

      LOGICAL     :: ll_coincidence

      TYPE(TVAR)  :: tl_lon0
      TYPE(TVAR)  :: tl_lat0
      TYPE(TVAR)  :: tl_lon1
      TYPE(TVAR)  :: tl_lat1

      TYPE(TFILE) :: tl_coord0

      TYPE(TMPP)  :: tl_mppcoord0

      TYPE(TDOM)  :: tl_dom0

      ! loop indices
      INTEGER(i4) :: ji
      INTEGER(i4) :: jj
      !----------------------------------------------------------------

      ll_coincidence=.TRUE.

      ! read coarse longitue and latitude on domain
      tl_coord0=td_coord0
      CALL iom_open(tl_coord0)

      !2-1 compute domain
      tl_dom0=dom_init( tl_coord0,         &
      &                 id_imin0, id_imax0,&
      &                 id_jmin0, id_jmax0 )

      !2-2 close file
      CALL iom_close(tl_coord0)

      !2-3 read variables on domain (ugly way to do it, have to work on it)
      !2-3-1 init mpp structure
      tl_mppcoord0=mpp_init(tl_coord0)

      CALL file_clean(tl_coord0)

      !2-3-2 get processor to be used
      CALL mpp_get_use( tl_mppcoord0, tl_dom0 )

      !2-3-3 open mpp files
      CALL iom_mpp_open(tl_mppcoord0)

      !2-3-4 read variable value on domain
      tl_lon0=iom_mpp_read_var(tl_mppcoord0,'longitude',td_dom=tl_dom0)
      tl_lat0=iom_mpp_read_var(tl_mppcoord0,'latitude' ,td_dom=tl_dom0)

      !2-3-5 close mpp files
      CALL iom_mpp_close(tl_mppcoord0)

      !2-3-6 clean structure
      CALL mpp_clean(tl_mppcoord0)

      ! read fine longitue and latitude
      tl_lon1=iom_read_var(td_coord1,'longitude')
      tl_lat1=iom_read_var(td_coord1,'latitude')
      
      CALL logger_debug("GRID CHECK COINCIDENCE:"//&
      &        " fine   grid "//TRIM(td_coord1%c_name) )
      CALL logger_debug("GRID CHECK COINCIDENCE:"//&
      &        " coarse grid "//TRIM(td_coord0%c_name) )

      !1- check domain
      !1-1 check global grid
      IF( .NOT. grid_is_global(tl_lon0, tl_lat0) )THEN
         IF( grid_is_global(tl_lon1, tl_lat1) )THEN

            ll_coincidence=.FALSE.
            CALL logger_fatal("GRID CHECK COINCIDENCE:"//&
            &        " fine   grid is global,"//&
            &        " coarse grid is not ")

         ELSE
            !1-2 ew overlap
            il_ew1=tl_lon1%i_ew
            IF( il_ew1 >= 0 )THEN

               il_ew0=tl_lon0%i_ew
               IF( il_ew0 < 0 )THEN
                  CALL logger_fatal("GRID CHECK COINCIDENCE: "//&
                  &        "fine grid has east west overlap,"//&
                  &        " coarse grid not ")
               ENDIF

               il_jmin1=1+ig_ghost
               il_jmax1=tl_lon1%t_dim(2)%i_len-ig_ghost

               ll_coincidence=grid__check_lat(&
               &                     tl_lat0%d_value(1,:,1,1),&
               &                     tl_lat1%d_value(1,il_jmin1:il_jmax1,1,1),&
               &                     id_rho(jp_J) )

            ELSE
               !1-3 other case
               il_imin1=1+ig_ghost
               il_jmin1=1+ig_ghost

               il_imax1=tl_lon1%t_dim(1)%i_len-ig_ghost
               il_jmax1=tl_lon1%t_dim(2)%i_len-ig_ghost

               ll_coincidence=grid__check_corner(&
               &                      tl_lon0%d_value(:,:,1,1),&
               &                      tl_lat0%d_value(:,:,1,1),&
               &                      tl_lon1%d_value(il_imin1:il_imax1, &
               &                                      il_jmin1:il_jmax1, &
               &                                      1,1),&
               &                      tl_lat1%d_value(il_imin1:il_imax1, &
               &                                      il_jmin1:il_jmax1, &
               &                                      1,1) )

            ENDIF
            
         ENDIF

         IF( .NOT. ll_coincidence )THEN
            CALL logger_fatal("GRID CHECK COINCIDENCE: no coincidence "//&
            &              "between fine grid and coarse grid. invalid domain" )
         ENDIF

      ENDIF
 
      !2- check refinement factor
      ! select point in middle of fine grid
      il_imid1=INT(tl_lon1%t_dim(1)%i_len*0.5)
      il_jmid1=INT(tl_lon1%t_dim(2)%i_len*0.5)
      
      dl_lon1=tl_lon1%d_value(il_imid1, il_jmid1,1,1)
      dl_lat1=tl_lat1%d_value(il_imid1, il_jmid1,1,1)

      ! select closest point on coarse grid
      il_indC(:)=grid_get_closest(tl_lon0%d_value(:,:,1,1),&
      &                           tl_lat0%d_value(:,:,1,1),&
      &                           dl_lon1, dl_lat1   )

      IF( ANY(il_indC(:)==0) )THEN
         CALL logger_fatal("GRID CHECK COINCIDENCE: can not find valid "//&
         &              "coarse grid indices. invalid domain" )
      ENDIF

      dl_lon0=tl_lon0%d_value(il_indC(1),il_indC(2),1,1)
      dl_lat0=tl_lat0%d_value(il_indC(1),il_indC(2),1,1)

      ! look for closest fine grid point from selected coarse grid point
      il_iind(:)=MAXLOC( tl_lon1%d_value(:,:,1,1), &
      &                  tl_lon1%d_value(:,:,1,1) <= dl_lon0)

      il_jind(:)=MAXLOC( tl_lat1%d_value(:,:,1,1), &
      &                  tl_lat1%d_value(:,:,1,1) <= dl_lat0 )

      il_indF(1)=il_iind(1)
      il_indF(2)=il_jind(2)

      IF( ANY(il_indF(:)==0) )THEN
         CALL logger_fatal("GRID CHECK COINCIDENCE: can not find valid "//&
         &              "fine grid indices. invalid domain" )
      ENDIF

      dl_lon1=tl_lon1%d_value(il_indF(1),il_indF(2),1,1)
      dl_lat1=tl_lat1%d_value(il_indF(1),il_indF(2),1,1)

      !2-1 check i-direction refinement factor
      DO ji=1,MIN(3,il_imid1)

         IF( il_indF(1)+ji*id_rho(jp_I)+1 > tl_lon1%t_dim(1)%i_len )THEN
            CALL logger_debug("GRID CHECK COINCIDENCE: tl_lon1%t_dim(1)%i_len "//TRIM(fct_str(tl_lon1%t_dim(1)%i_len)))
            CALL logger_debug("GRID CHECK COINCIDENCE: il_indF(1)+ji*id_rhoi+1 "//TRIM(fct_str(il_indF(1)+ji*id_rho(jp_I)+1)))
            CALL logger_debug("GRID CHECK COINCIDENCE: il_indF(1) "//TRIM(fct_str(il_indF(1))))
            CALL logger_debug("GRID CHECK COINCIDENCE: id_rhoi "//TRIM(fct_str(id_rho(jp_I))))
            CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//&
            &  " to check i-direction refinement factor ")
            EXIT
         ELSE
            dl_lon1=tl_lon1%d_value(il_indF(1)+ji*id_rho(jp_I),il_indF(2),1,1)
            dl_lon0=tl_lon0%d_value(il_indC(1)+ji,il_indC(2),1,1)

            dl_lon1p=tl_lon1%d_value(il_indF(1)+ji*id_rho(jp_I)+1,il_indF(2),1,1)

            dl_dlon=ABS(dl_lon1p-dl_lon1)*1.e-3

            SELECT CASE(MOD(id_rho(jp_I),2))

            CASE(0)

               IF( dl_lon1 >= dl_lon0 .OR. dl_lon0 >= dl_lon1p )THEN
                  ll_coincidence=.FALSE.
                  CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
                  &  "i-direction refinement factor ("//&
                  &   TRIM(fct_str(id_rho(jp_I)))//&
                  &   ") between fine grid and coarse grid ")
               ENDIF

            CASE DEFAULT         
            
               IF( ABS(dl_lon1 - dl_lon0) > dl_dlon )THEN
                  ll_coincidence=.FALSE.
                  CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
                  &  "i-direction refinement factor ("//&
                  &   TRIM(fct_str(id_rho(jp_I)))//&
                  &  ") between fine grid and coarse grid ")
               ENDIF
            
            END SELECT
         ENDIF

      ENDDO

      !2-2 check j-direction refinement factor
      DO jj=1,MIN(3,il_jmid1)

         IF( il_indF(2)+jj*id_rho(jp_J)+1 > tl_lat1%t_dim(2)%i_len )THEN
            CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//&
            &  " to check j-direction refinement factor ")
            EXIT
         ELSE      
            dl_lat1=tl_lat1%d_value(il_indF(1),il_indF(2)+jj*id_rho(jp_J),1,1)
            dl_lat0=tl_lat0%d_value(il_indC(1),il_indC(2)+jj,1,1)

            dl_lat1p=tl_lat1%d_value(il_indF(1),il_indF(2)+jj*id_rho(jp_J)+1,1,1)

            dl_dlat=ABS(dl_lat1p-dl_lat1)*1.e-3

            SELECT CASE(MOD(id_rho(jp_J),2))

            CASE(0)
               
               IF( dl_lat1 >= dl_lat0 .OR. dl_lat0 >= dl_lat1p )THEN
                  ll_coincidence=.FALSE.
                  CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
                  &  "j-direction refinement factor ("//&
                  &   TRIM(fct_str(id_rho(jp_J)))//&
                  &  ") between fine grid and coarse grid ")
               ENDIF

            CASE DEFAULT

               IF( ABS(dl_lat1-dl_lat0) > dl_dlat )THEN
                  ll_coincidence=.FALSE.
                  CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//&
                  &  "j-direction refinement factor ("//&
                  &   TRIM(fct_str(id_rho(jp_J)))//&
                  &  ") between fine grid and coarse grid ")
               ENDIF

            END SELECT
         ENDIF

      ENDDO

      IF( .NOT. ll_coincidence )THEN
         CALL logger_fatal("GRID CHECK COINCIDENCE: no coincidence "//&
         &              "between fine and coarse grid: "//&
         &              "invalid refinement factor" )
      ENDIF

   END SUBROUTINE grid_check_coincidence
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This function check that fine grid is 
   !> inside coarse grid 
   !
   !> @details
   !> 
   !> @note deltalon and delatlat are used only to avoid issue due to 
   !> cubic interpolation approximation on the firsts grid points  
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] dd_lon0 : table of coarse grid longitude
   !> @param[in] dd_lat0 : table of coarse grid latitude
   !> @param[in] dd_lon1 : table of fine   grid longitude 
   !> @param[in] dd_lat1 : table of fine   grid latitude
   !> @return logical, fine grid is inside coarse grid
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid__check_corner(dd_lon0, dd_lat0, &
   &                           dd_lon1, dd_lat1 )
   IMPLICIT NONE
      ! Argument      
      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon0
      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat0
      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lon1
      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_lat1

      ! function
      LOGICAL :: grid__check_corner

      ! local variable
      INTEGER(i4), DIMENSION(2) :: il_shape0
      INTEGER(i4), DIMENSION(2) :: il_shape1

      INTEGER(i4) :: il_imin0
      INTEGER(i4) :: il_jmin0
      INTEGER(i4) :: il_imax0
      INTEGER(i4) :: il_jmax0

      INTEGER(i4) :: il_imin1
      INTEGER(i4) :: il_jmin1
      INTEGER(i4) :: il_imax1
      INTEGER(i4) :: il_jmax1

      REAL(dp)    :: dl_lon0
      REAL(dp)    :: dl_lat0

      REAL(dp)    :: dl_lon1
      REAL(dp)    :: dl_lat1

      REAL(dp)    :: dl_dlon
      REAL(dp)    :: dl_dlat
      ! loop indices
      !----------------------------------------------------------------

      ! init
      grid__check_corner=.TRUE.

      il_shape0=SHAPE(dd_lon0(:,:))
      il_shape1=SHAPE(dd_lon1(:,:))

      !1- check if fine grid inside coarse grid domain
      il_imin0=1 ; il_imax0=il_shape0(1)
      il_jmin0=1 ; il_jmax0=il_shape0(2)

      il_imin1=1 ; il_imax1=il_shape1(1)
      il_jmin1=1 ; il_jmax1=il_shape1(2)

      ! check lower left corner
      dl_lon0 = dd_lon0(il_imin0, il_jmin0  )
      dl_lat0 = dd_lat0(il_imin0, il_jmin0  )

      dl_lon1 = dd_lon1(il_imin1, il_jmin1)
      dl_lat1 = dd_lat1(il_imin1, il_jmin1)

      dl_dlon=ABS(dd_lon1(il_imin1+1,il_jmin1  )-dl_lon1)*1.e-3
      dl_dlat=ABS(dd_lat1(il_imin1  ,il_jmin1+1)-dl_lat1)*1.e-3

      IF( (ABS(dl_lon1-dl_lon0)>dl_dlon) .AND. (dl_lon1 < dl_lon0 ) .OR. & 
      &   (ABS(dl_lat1-dl_lat0)>dl_dlat) .AND. (dl_lat1 < dl_lat0 ) )THEN

         CALL logger_error("GRID CHECK COINCIDENCE: fine grid lower left "//&
         &     "corner  not north east of coarse grid (imin,jmin) ")
         CALL logger_debug(" fine   grid lower left ( "//&
         &              TRIM(fct_str(dl_lon1))//","//&
         &              TRIM(fct_str(dl_lat1))//")" )
         CALL logger_debug(" coarse grid lower left ( "//&
         &              TRIM(fct_str(dl_lon0))//","//&
         &              TRIM(fct_str(dl_lat0))//")" )
         grid__check_corner=.FALSE.

      ENDIF

      ! check upper left corner
      dl_lon0 = dd_lon0(il_imin0, il_jmax0  )
      dl_lat0 = dd_lat0(il_imin0, il_jmax0  )

      dl_lon1 = dd_lon1(il_imin1, il_jmax1)
      dl_lat1 = dd_lat1(il_imin1, il_jmax1)

      dl_dlon=ABS(dd_lon1(il_imin1+1,il_jmax1  )-dl_lon1)*1.e-3
      dl_dlat=ABS(dd_lat1(il_imin1  ,il_jmax1-1)-dl_lat1)*1.e-3

      IF( (ABS(dl_lon1-dl_lon0)>dl_dlon) .AND. (dl_lon1 < dl_lon0) .OR. &
      &   (ABS(dl_lat1-dl_lat0)>dl_dlat) .AND. (dl_lat1 > dl_lat0) )THEN

         CALL logger_error("GRID CHECK COINCIDENCE: fine grid upper left "//&
         &     "corner not south east of coarse grid (imin,jmax) ")
         CALL logger_debug(" fine   grid upper left ("//&
         &              TRIM(fct_str(dl_lon1))//","//&
         &              TRIM(fct_str(dl_lat1))//")")
         CALL logger_debug(" coasre grid upper left ("//&
         &              TRIM(fct_str(dl_lon0))//","//&
         &              TRIM(fct_str(dl_lat0))//")")
         grid__check_corner=.FALSE.

      ENDIF

      ! check lower right corner
      dl_lon0 = dd_lon0(il_imax0, il_jmin0  )
      dl_lat0 = dd_lat0(il_imax0, il_jmin0  )

      dl_lon1 = dd_lon1(il_imax1, il_jmin1)
      dl_lat1 = dd_lat1(il_imax1, il_jmin1)

      dl_dlon=ABS(dd_lon1(il_imax1-1,il_jmin1  )-dl_lon1)*1.e-3
      dl_dlat=ABS(dd_lat1(il_imax1  ,il_jmin1+1)-dl_lat1)*1.e-3

      IF( (ABS(dl_lon1-dl_lon0)>dl_dlon) .AND. (dl_lon1 > dl_lon0) .OR. &
      &   (ABS(dl_lat1-dl_lat0)>dl_dlat) .AND. (dl_lat1 < dl_lat0) )THEN

         CALL logger_error("GRID CHECK COINCIDENCE: fine grid lower right "//&
         &     "corner not north west west of coarse grid (imax,jmin) ")
         CALL logger_debug(" fine   grid lower right ( "//&
         &              TRIM(fct_str(dl_lon1))//","//&
         &              TRIM(fct_str(dl_lat1))//")" )
         CALL logger_debug(" coarse grid lower right ( "//&
         &              TRIM(fct_str(dl_lon0))//","//&
         &              TRIM(fct_str(dl_lat0))//")" )   
         grid__check_corner=.FALSE.

      ENDIF

      ! check upper right corner
      dl_lon0 = dd_lon0(il_imax0, il_jmax0  )
      dl_lat0 = dd_lat0(il_imax0, il_jmax0  )

      dl_lon1 = dd_lon1(il_imax1, il_jmax1)
      dl_lat1 = dd_lat1(il_imax1, il_jmax1)

      dl_dlon=ABS(dd_lon1(il_imax1-1,il_jmax1  )-dl_lon1)*1.e-3
      dl_dlat=ABS(dd_lat1(il_imax1  ,il_jmax1-1)-dl_lat1)*1.e-3

      IF( (ABS(dl_lon1-dl_lon0)>dl_dlon) .AND. (dl_lon1 > dl_lon0) .OR. &
      &   (ABS(dl_lat1-dl_lat0)>dl_dlat) .AND. (dl_lat1 > dl_lat0) )THEN

         CALL logger_error("GRID CHECK COINCIDENCE: fine grid upper right "//&
         &     "corner not south west of coarse grid (imax,jmax) ")
         CALL logger_debug(" fine   grid upper right ( "//&
         &              TRIM(fct_str(dl_lon1))//","//&
         &              TRIM(fct_str(dl_lat1))//")" )
         CALL logger_debug(" fine   imax1 jmax1 ( "//&
         &              TRIM(fct_str(il_imax1))//","//&
         &              TRIM(fct_str(il_jmax1))//")" )
         CALL logger_debug(" coarse grid upper right ( "//&
         &              TRIM(fct_str(dl_lon0))//","//&
         &              TRIM(fct_str(dl_lat0))//")" )    
         CALL logger_debug(" fine   imax0 jmax0 ( "//&
         &              TRIM(fct_str(il_imax0))//","//&
         &              TRIM(fct_str(il_jmax0))//")" )
         grid__check_corner=.FALSE.

      ENDIF      

   END FUNCTION grid__check_corner
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This function check that fine grid latitude are 
   !> inside coarse grid latitude
   !
   !> @details
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] dd_lat0 : table of coarse grid latitude 
   !> @param[in] dd_lat1 : table of fine grid latitude 
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid__check_lat(dd_lat0, dd_lat1, id_rhoj)
      IMPLICIT NONE
      ! Argument 
      REAL(dp), DIMENSION(:), INTENT(IN) :: dd_lat0
      REAL(dp), DIMENSION(:), INTENT(IN) :: dd_lat1
      INTEGER(i4)           , INTENT(IN) :: id_rhoj

      ! function
      LOGICAL :: grid__check_lat

      ! local variable
      INTEGER(i4), DIMENSION(1) :: il_shape0
      INTEGER(i4), DIMENSION(1) :: il_shape1

      INTEGER(i4) :: il_jmin0
      INTEGER(i4) :: il_jmax0

      INTEGER(i4) :: il_jmin1
      INTEGER(i4) :: il_jmax1

      REAL(dp)    :: dl_dlat
      ! loop indices
      !----------------------------------------------------------------

      ! init
      grid__check_lat=.TRUE.

      il_shape0(:)=SHAPE(dd_lat0(:))
      il_shape1(:)=SHAPE(dd_lat1(:))

      !1- check if fine grid inside coarse grid domain
      il_jmin0=1+1 ; il_jmax0=il_shape0(1)-1

      il_jmin1=1+id_rhoj ; il_jmax1=il_shape1(1)-id_rhoj

      dl_dlat=ABS(dd_lat1(il_jmin1+1)-dd_lat1(il_jmin1))*1.e-3

      ! check lower left fine grid
      IF( ABS(dd_lat1(il_jmin1)-dd_lat0(il_jmin0)) > dl_dlat .AND. &
      &   dd_lat1(il_jmin1) < dd_lat0(il_jmin0) )THEN

         CALL logger_error("GRID CHECK COINCIDENCE: fine grid lower point"//&
         &     " not north of coarse grid (jmin) ")
         CALL logger_debug(" fine grid lower point ( "//&
         &              TRIM(fct_str(dd_lat1(il_jmin1)))//")" )
         CALL logger_debug(" coarse grid lower point ( "//&
         &              TRIM(fct_str(dd_lat0(il_jmin0)))//")" )
         grid__check_lat=.FALSE.

      ENDIF

      dl_dlat=ABS(dd_lat1(il_jmax1-1)-dd_lat1(il_jmax1))*1.e-3

      ! check upper left fine grid
      IF( ABS(dd_lat1(il_jmax1)-dd_lat0(il_jmax0)) > dl_dlat .AND. &
      &   dd_lat1(il_jmax1) > dd_lat0(il_jmax0) )THEN

         CALL logger_error("GRID CHECK COINCIDENCE: fine grid upper point"//&
         &     " not south of coarse grid (jmax) ")
         CALL logger_debug(" fine grid upper point ("//&
         &              TRIM(fct_str(dd_lat1(il_jmax1)))//")")
         CALL logger_debug(" coasre grid upper point ("//&
         &              TRIM(fct_str(dd_lat0(il_jmax0)))//")")
         grid__check_lat=.FALSE.

      ENDIF
      
   END FUNCTION grid__check_lat
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief
   !> This subroutine add ghost cell at boundaries.
   !> 
   !> @author J.Paul
   !> - Nov, 2013-Initial version
   !
   !> @param[inout] td_var : table of variable structure 
   !> @param[in] id_ighost : i-direction ghost cell factor 
   !> @param[in] id_jghost : j-direction ghost cell factor 
   !-------------------------------------------------------------------
   !> @code
   SUBROUTINE grid_add_ghost(td_var, id_ighost, id_jghost)
      IMPLICIT NONE
      ! Argument
      TYPE(TVAR) , INTENT(INOUT) :: td_var
      INTEGER(i4), INTENT(IN   ) :: id_ighost
      INTEGER(i4), INTENT(IN   ) :: id_jghost

      ! local variable
      INTEGER(i4) :: il_imin
      INTEGER(i4) :: il_jmin
      INTEGER(i4) :: il_imax
      INTEGER(i4) :: il_jmax

      REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
      
      TYPE(TVAR) :: tl_var

      ! loop indices
      !----------------------------------------------------------------

      IF( ALL(td_var%t_dim(1:2)%l_use) )THEN

         CALL logger_warn( "ADD GHOST: dimension change in variable "//&
         &              TRIM(td_var%c_name) )

         ! copy variable
         tl_var=td_var

         CALL var_del_value(td_var)

         ! compute indice to fill center
         il_imin=1+id_ighost*ig_ghost
         il_jmin=1+id_jghost*ig_ghost

         il_imax=il_imin+tl_var%t_dim(1)%i_len-1
         il_jmax=il_jmin+tl_var%t_dim(2)%i_len-1

         ! compute new dimension
         td_var%t_dim(1)%i_len = tl_var%t_dim(1)%i_len + 2*id_ighost*ig_ghost
         td_var%t_dim(2)%i_len = tl_var%t_dim(2)%i_len + 2*id_jghost*ig_ghost

         ALLOCATE(dl_value(td_var%t_dim(1)%i_len, &
         &                 td_var%t_dim(2)%i_len, &
         &                 td_var%t_dim(3)%i_len, &
         &                 td_var%t_dim(4)%i_len) )

         dl_value(:,:,:,:)=tl_var%d_fill

         dl_value(il_imin:il_imax, &
         &        il_jmin:il_jmax, &
         &              :,:)  =  tl_var%d_value(:,:,:,:)

         ! add variable value
         CALL var_add_value(td_var,dl_value(:,:,:,:))

         ! save variable type
         td_var%i_type=tl_var%i_type
         
         DEALLOCATE( dl_value )

         CALL var_clean(tl_var)

      ENDIF

   END SUBROUTINE grid_add_ghost
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief
   !> This subroutine delete ghost cell at boundaries.
   !> 
   !> @author J.Paul
   !> - Nov, 2013-Initial version
   !
   !> @param[inout] td_var : table of variable structure 
   !> @param[in] id_ighost : i-direction ghost cell factor 
   !> @param[in] id_jghost : j-direction ghost cell factor 
   !-------------------------------------------------------------------
   !> @code
   SUBROUTINE grid_del_ghost(td_var, id_ighost, id_jghost)
      IMPLICIT NONE
      ! Argument
      TYPE(TVAR) , INTENT(INOUT) :: td_var
      INTEGER(i4), INTENT(IN   ) :: id_ighost
      INTEGER(i4), INTENT(IN   ) :: id_jghost

      ! local variable
      INTEGER(i4) :: il_imin
      INTEGER(i4) :: il_jmin
      INTEGER(i4) :: il_imax
      INTEGER(i4) :: il_jmax

      REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
      
      TYPE(TVAR) :: tl_var

      ! loop indices
      !----------------------------------------------------------------

      IF( ALL(td_var%t_dim(1:2)%l_use) )THEN

         CALL logger_warn( "DEL GHOST: dimension change in variable "//&
         &              TRIM(td_var%c_name) )

         ! copy variable
         tl_var=td_var

         CALL var_del_value(td_var)

         ! compute indice to get center
         il_imin=1+id_ighost*ig_ghost
         il_jmin=1+id_jghost*ig_ghost

         il_imax=tl_var%t_dim(1)%i_len-id_ighost*ig_ghost
         il_jmax=tl_var%t_dim(2)%i_len-id_jghost*ig_ghost

         ! compute new dimension
         td_var%t_dim(1)%i_len = tl_var%t_dim(1)%i_len - 2*id_ighost*ig_ghost
         td_var%t_dim(2)%i_len = tl_var%t_dim(2)%i_len - 2*id_jghost*ig_ghost

         ALLOCATE(dl_value(td_var%t_dim(1)%i_len, &
         &                 td_var%t_dim(2)%i_len, &
         &                 td_var%t_dim(3)%i_len, &
         &                 td_var%t_dim(4)%i_len) )

         dl_value(:,:,:,:)=tl_var%d_fill

         dl_value(:,:,:,:)  =  tl_var%d_value(il_imin:il_imax, &
         &                                    il_jmin:il_jmax, &
         &                                    :,:)

         ! add variable value
         CALL var_add_value(td_var,dl_value(:,:,:,:))

         ! save variable type
         td_var%i_type=tl_var%i_type
         
         DEALLOCATE( dl_value )

         CALL var_clean(tl_var)

      ENDIF

   END SUBROUTINE grid_del_ghost
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This subroutine fill small closed sea with fill value. 
   !
   !> @details
   !> the minimum size (nbumber of point) of closed sea to be kept could be
   !> sepcify with id_minsize.
   !> By default only the biggest sea is preserve.
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[inout] td_var : variable structure
   !> @param[in] id_mask : domain mask (from grid_split_domain)
   !> @param[in] id_minsize : minimum size of sea to be kept
   !-------------------------------------------------------------------
   !> @code
   SUBROUTINE grid_fill_small_dom(td_var, id_mask, id_minsize)
      IMPLICIT NONE
      ! Argument      
      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var
      INTEGER(i4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_mask
      INTEGER(i4),                 INTENT(IN   ), OPTIONAL :: id_minsize

      ! local variable
      INTEGER(i4)                              :: il_ndom
      INTEGER(i4)                              :: il_minsize
      INTEGER(i4), DIMENSION(2)                :: il_shape
      INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_tmp

      ! loop indices
      INTEGER(i4) :: ji
      INTEGER(i4) :: jk
      INTEGER(i4) :: jl
      !----------------------------------------------------------------

      il_shape(:)=SHAPE(id_mask(:,:))
      IF( ANY(il_shape(:) /= td_var%t_dim(1:2)%i_len) )THEN
         CALL logger_error("GRID FILL SMALL DOM: variable and mask "//&
         &              "dimension differ")
      ELSE

         il_ndom=MINVAL(id_mask(:,:))

         ALLOCATE( il_tmp(il_shape(1),il_shape(2)) )
         il_tmp(:,:)=0
         DO ji=-1,il_ndom,-1
            WHERE( id_mask(:,:)==ji ) 
               il_tmp(:,:)=SUM(id_mask(:,:),id_mask(:,:)==ji)/ji
            END WHERE
         ENDDO

         il_minsize=MAXVAL(il_tmp(:,:))
         IF( PRESENT(id_minsize) ) il_minsize=id_minsize

         DO jl=1,td_var%t_dim(4)%i_len
            DO jk=1,td_var%t_dim(3)%i_len
               WHERE( il_tmp(:,:) < il_minsize ) 
                  td_var%d_value(:,:,jk,jl)=td_var%d_fill
               END WHERE
            ENDDO
         ENDDO

         DEALLOCATE( il_tmp )

      ENDIF

   END SUBROUTINE grid_fill_small_dom
   !> @endcode
   !-------------------------------------------------------------------
   !> @brief This subroutine compute closed sea domain.
   !
   !> @details
   !> to each domain is associated a negative value id (from -1 to ...)
   !
   !> @author J.Paul
   !> - Nov, 2013- Initial Version
   !
   !> @param[in] td_var : variable strucutre 
   !> @param[in] id_level : level
   !> @return domain mask  
   !-------------------------------------------------------------------
   !> @code
   FUNCTION grid_split_domain(td_var, id_level)
      IMPLICIT NONE
      ! Argument
      TYPE(TVAR) , INTENT(IN) :: td_var
      INTEGER(i4), INTENT(IN), OPTIONAL :: id_level

      ! function
      INTEGER(i4), DIMENSION(td_var%t_dim(1)%i_len, &
      &                      td_var%t_dim(2)%i_len ) :: grid_split_domain

      ! local variable
      INTEGER(i4)                              :: il_domid
      INTEGER(i4)                              :: il_level
      INTEGER(i4), DIMENSION(2)                :: il_shape
      INTEGER(i4), DIMENSION(2)                :: il_ind
      INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_mask
      INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_tmp

      LOGICAL                                  :: ll_full

      ! loop indices
      INTEGER(i4) :: ji
      INTEGER(i4) :: jim
      INTEGER(i4) :: jip
      INTEGER(i4) :: jj
      INTEGER(i4) :: jjm
      INTEGER(i4) :: jjp
      !----------------------------------------------------------------
      il_level=1
      IF( PRESENT(id_level) ) il_level=id_level

      ! init
      il_domid=-1

      il_shape(:)=td_var%t_dim(1:2)%i_len
      ALLOCATE( il_mask(il_shape(1),il_shape(2)) )
      il_mask(:,:)=0
      WHERE( td_var%d_value(:,:,il_level,1)/=td_var%d_fill ) il_mask(:,:)=1

      il_ind(:)=MAXLOC( il_mask(:,:) )
      DO WHILE( il_mask(il_ind(1), il_ind(2)) == 1 )

         il_mask(il_ind(1),il_ind(2))=il_domid
         ll_full=.FALSE.

         ALLOCATE( il_tmp(il_shape(1),il_shape(2)) )

         DO WHILE( .NOT. ll_full )
            il_tmp(:,:)=0

            ll_full=.TRUE.
            DO jj=1,il_shape(2)
               DO ji=1,il_shape(1)
                  IF( il_mask(ji,jj)==il_domid )THEN
                     jim=MAX(1,ji-1)   ;  jip=MIN(il_shape(1),ji+1)
                     jjm=MAX(1,jj-1)   ;  jjp=MIN(il_shape(2),jj+1)
                     
                     WHERE( il_mask(jim:jip,jjm:jjp)==1 )
                        il_mask(jim:jip,jjm:jjp)=il_domid
                        il_tmp(jim:jip,jjm:jjp)=1
                     END WHERE
                  ENDIF
               ENDDO
            ENDDO
            IF( SUM(il_tmp(:,:))/=0 ) ll_full=.FALSE.

         ENDDO
         DEALLOCATE( il_tmp )

         il_ind(:)=MAXLOC( il_mask(:,:) )
         il_domid=il_domid-1

      ENDDO

      ! save result
      grid_split_domain(:,:)=il_mask(:,:)

      DEALLOCATE( il_mask )

      CALL logger_info("GRID SPLIT DOMAIN: "//TRIM( fct_str(ABS(il_domid+1)) )//&
      &             " domain found" ) 

   END FUNCTION grid_split_domain
   !> @endcode
!   !-------------------------------------------------------------------
!   !> @brief This function 
!   !
!   !> @details
!   !
!   !> @author J.Paul
!   !> - Nov, 2013- Initial Version
!   !
!   !> @param[in] 
!   !-------------------------------------------------------------------
!   !> @code
!   FUNCTION grid_()
!      IMPLICIT NONE
!      ! Argument      
!      ! function
!      ! local variable
!      ! loop indices
!      !----------------------------------------------------------------
!
!   END FUNCTION grid_
!   !> @endcode
!   !-------------------------------------------------------------------
!   !> @brief This subroutine 
!   !
!   !> @details
!   !
!   !> @author J.Paul
!   !> - Nov, 2013- Initial Version
!   !
!   !> @param[in] 
!   !-------------------------------------------------------------------
!   !> @code
!   SUBROUTINE grid_()
!      IMPLICIT NONE
!      ! Argument      
!      ! local variable
!      ! loop indices
!      !----------------------------------------------------------------
!
!   END SUBROUTINE grid_
!   !> @endcode
END MODULE grid