Index: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/TOOLS/SIREN/src/interp.f90
===================================================================
--- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/TOOLS/SIREN/src/interp.f90	(revision 5006)
+++ branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/TOOLS/SIREN/src/interp.f90	(revision 5602)
@@ -8,19 +8,73 @@
 !> @brief 
 !> This module manage interpolation on regular grid.
-!> @note It is used to work on ORCA grid, as we work only with grid indices.
+!>
+!> @details Interpolation method to be used is specify inside variable 
+!>    strcuture, as array of string character.<br/> 
+!>    - td_var\%c_interp(1) string character is the interpolation name choose between:
+!>       - 'nearest'
+!>       - 'cubic  '
+!>       - 'linear '
+!>    - td_var\%c_interp(2) string character is an operation to be used 
+!>    on interpolated value.<br/>
+!>          operation have to be mulitplication '*' or division '/'.<br/>
+!>          coefficient have to be refinement factor following i-direction 'rhoi', 
+!>          j-direction 'rhoj', or k-direction 'rhok'.<br/>
+!>
+!>          Examples: '*rhoi', '/rhoj'.
+!> 
+!>    @note Those informations are read from namelist or variable configuration file (default).<br/>
+!>    Interplation method could be specify for each variable in namelist _namvar_,
+!>    defining string character _cn\_varinfo_.<br/>
+!>    Example:
+!>       - cn_varinfo='varname1:cubic/rhoi', 'varname2:linear' 
+!>
+!>    to create mixed grid (with coarse grid point needed to compute
+!> interpolation):<br/>
+!> @code
+!>    CALL interp_create_mixed_grid( td_var, td_mix [,id_rho] )
+!> @endcode
+!>       - td_var is coarse grid variable (should be extrapolated)
+!>       - td_mix is mixed grid variable structure [output]
+!>       - id_rho is array of refinment factor [optional]
+!>
+!>    to detected point to be interpolated:<br/>
+!> @code
+!>    il_detect(:,:,:)=interp_detect( td_mix [,id_rho] )
+!> @endcode
+!>       - il_detect(:,:,:) is 3D array of detected point to be interpolated
+!>       - td_mix is mixed grid variable
+!>       - id_rho is array of refinement factor [optional]
+!>
+!>    to interpolate variable value:<br/>
+!> @code
+!>    CALL  interp_fill_value( td_var [,id_rho] [,id_offset] )
+!> @endcode
+!>       - td_var is variable structure
+!>       - id_rho is array of refinement factor [optional]
+!>       - id_offset is array of offset between fine and coarse grid [optional]
+!>
+!>    to clean mixed grid (remove points added on mixed grid to compute interpolation):<br/>
+!> @code
+!>    CALL interp_clean_mixed_grid( td_mix, td_var, id_rho )
+!> @endcode
+!>       - td_mix is mixed grid variable structure
+!>       - td_var is variable structure [output]
+!>       - id_rho is array of refinement factor [optional]
+!>       - id_offset is array of offset between fine and coarse grid [optional]
+!>
+!> @note It use to work on ORCA grid, as we work only with grid indices.
 !>
 !> @warning due to the use of second derivative when using cubic interpolation
-!> you should add 2 extrabands
+!> you should add at least 2 extrabands.
 !>
-!> @details
 !> @author
 !> J.Paul
 ! REVISION HISTORY:
-!> @date Nov, 2013 - Initial Version
-!
+!> @date November, 2013 - Initial Version
+!> @date September, 2014
+!> - add header
+!> - use interpolation method modules
+!>
 !> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-!> @todo
-!> - interp 3D
-!> - see issue when fill value is zero for check cubic_fill..
 !----------------------------------------------------------------------
 MODULE interp
@@ -29,22 +83,17 @@
    USE global                          ! global variable
    USE kind                            ! F90 kind parameter
-   USE logger                             ! log file manager
+   USE logger                          ! log file manager
    USE fct                             ! basic useful function
-!   USE date                            ! date manager
+   USE date                            ! date manager
    USE att                             ! attribute manager
    USE dim                             ! dimension manager
    USE var                             ! variable manager
-!   USE file                            ! file manager
-!   USE iom                             ! I/O manager
-!   USE dom                             ! domain manager
    USE grid                            ! grid manager
    USE extrap                          ! extrapolation manager
-!   USE interp                          ! interpolation manager
-!   USE filter                          ! filter manager
-!   USE mpp                             ! MPP manager
-!   USE iom_mpp                         ! MPP I/O manager
+   USE interp_cubic                    !   cubic interpolation manager
+   USE interp_linear                   !  linear interpolation manager
+   USE interp_nearest                  ! nearest interpolation manager
 
    IMPLICIT NONE
-   PRIVATE
    ! NOTE_avoid_public_variables_if_possible
 
@@ -53,34 +102,16 @@
    ! function and subroutine
    PUBLIC :: interp_detect             !< detected point to be interpolated
-   PUBLIC :: interp_fill_value         !< interpolate value over detectected point
+   PUBLIC :: interp_fill_value         !< interpolate value
    PUBLIC :: interp_create_mixed_grid  !< create mixed grid
    PUBLIC :: interp_clean_mixed_grid   !< clean mixed grid
 
-   PRIVATE :: interp__detect           !< detected point to be interpolated
-   PRIVATE :: interp__detect_wrapper   !< detected point to be interpolated
-   PRIVATE :: interp__fill_value_wrapper !< interpolate value over detectected point
-   PRIVATE :: interp__fill_value       !< interpolate value over detectected point
-   PRIVATE :: interp__clean_even_grid  !< clean even mixed grid
-   PRIVATE :: interp__del_offset       !< remove offset from interpolated grid
-   PRIVATE :: interp__check_method     !< check if interpolation method available
-!   PRIVATE :: interp__3D               !< interpolate 3D grid
-   PRIVATE :: interp__2D               !< interpolate 2D grid
-   PRIVATE :: interp__1D               !< interpolate 1D grid 
-   PRIVATE :: interp__2D_cubic_coef    !< compute coefficient for bicubic interpolation 
-   PRIVATE :: interp__2D_cubic_fill    !< compute bicubic interpolation 
-   PRIVATE :: interp__2D_linear_coef   !< compute coefficient for bilinear interpolation
-   PRIVATE :: interp__2D_linear_fill   !< compute bilinear interpolation 
-   PRIVATE :: interp__2D_nearest_fill  !< compute nearest interpolation 
-   PRIVATE :: interp__1D_cubic_coef    !< compute coefficient for cubic interpolation 
-   PRIVATE :: interp__1D_cubic_fill    !< compute cubic interpolation 
-   PRIVATE :: interp__1D_linear_coef   !< compute coefficient for linear interpolation
-   PRIVATE :: interp__1D_linear_fill   !< compute linear interpolation 
-   PRIVATE :: interp__1D_nearest_fill  !< compute nearest interpolation 
-!   PRIVATE :: interp__longitude
+   PRIVATE :: interp__detect              ! detected point to be interpolated
+   PRIVATE :: interp__detect_wrapper      ! detected point to be interpolated
+   PRIVATE :: interp__fill_value_wrapper  ! interpolate value over detectected point
+   PRIVATE :: interp__fill_value          ! interpolate value over detectected point
+   PRIVATE :: interp__clean_even_grid     ! clean even mixed grid
+   PRIVATE :: interp__check_method        ! check if interpolation method available
    
    TYPE TINTERP
-      !CHARACTER(LEN=lc) :: c_name = 'unknown'   !< interpolation method name
-      !CHARACTER(LEN=lc) :: c_factor = 'unknown' !< interpolation factor
-      !CHARACTER(LEN=lc) :: c_divisor = 'unknown'   !< interpolation divisor
       CHARACTER(LEN=lc) :: c_name   = '' !< interpolation method name
       CHARACTER(LEN=lc) :: c_factor = '' !< interpolation factor
@@ -89,9 +120,9 @@
 
    INTERFACE interp_detect
-      MODULE PROCEDURE interp__detect_wrapper !< detected point to be interpolated
+      MODULE PROCEDURE interp__detect_wrapper 
    END INTERFACE interp_detect
 
    INTERFACE interp_fill_value
-      MODULE PROCEDURE interp__fill_value_wrapper !< detected point to be interpolated
+      MODULE PROCEDURE interp__fill_value_wrapper 
    END INTERFACE interp_fill_value
 
@@ -99,16 +130,16 @@
    !-------------------------------------------------------------------
    !> @brief
-   !> This function check if interpolation method available.
+   !> This function check if interpolation method is available.
    !> 
    !> @details 
+   !> check if name of interpolation method is present in global list of string
+   !> character cp_interp_list (see global.f90).
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
+   !> - November, 2013- Initial Version
    !
-   !> @param[in] cd_method : interpolation method
+   !> @param[in] cd_method interpolation method
    !> @return 
-   !> @todo see extrap_detect
-   !-------------------------------------------------------------------
-   !> @code
+   !-------------------------------------------------------------------
    FUNCTION interp__check_method( cd_method )
       IMPLICIT NONE
@@ -130,6 +161,6 @@
 
       interp__check_method=.FALSE.
-      DO ji=1,ig_ninterp
-         cl_interp=fct_lower(cg_interp_list(ji))
+      DO ji=1,ip_ninterp
+         cl_interp=fct_lower(cp_interp_list(ji))
          IF( TRIM(cl_interp) == TRIM(cl_method) )THEN
             interp__check_method=.TRUE.
@@ -139,5 +170,4 @@
 
    END FUNCTION interp__check_method
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
@@ -149,11 +179,10 @@
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in] td_mix : mixed grid variable (to interpolate)
-   !> @param[in] id_rho : table of refinement factor
-   !> @return table of detected point to be interpolated 
-   !-------------------------------------------------------------------
-   !> @code
+   !> - November, 2013- Initial Version
+   !>
+   !> @param[in] td_mix mixed grid variable (to interpolate)
+   !> @param[in] id_rho array of refinement factor
+   !> @return 3D array of detected point to be interpolated 
+   !-------------------------------------------------------------------
    FUNCTION interp__detect_wrapper( td_mix, id_rho )
       IMPLICIT NONE
@@ -168,5 +197,4 @@
 
       ! local variable
-
       ! loop indices
       !----------------------------------------------------------------
@@ -208,5 +236,4 @@
 
    END FUNCTION interp__detect_wrapper
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
@@ -217,11 +244,10 @@
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
+   !> - November, 2013- Initial Version
    !
-   !> @param[in] td_mix : mixed grid variable (to interpolate)
-   !> @param[in] id_rho : table of refinement factor
-   !> @return table of detected point to be interpolated 
-   !-------------------------------------------------------------------
-   !> @code
+   !> @param[in] td_mix mixed grid variable (to interpolate)
+   !> @param[in] id_rho array of refinement factor
+   !> @return 3D array of detected point to be interpolated 
+   !-------------------------------------------------------------------
    FUNCTION interp__detect( td_mix, id_rho )
       IMPLICIT NONE
@@ -253,5 +279,5 @@
       ALLOCATE( il_rho(ip_maxdim) )
       il_rho(:)=1
-      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
+      IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
 
       ! special case for even refinement on ARAKAWA-C grid
@@ -289,4 +315,5 @@
       il_dim(:)=td_mix%t_dim(1:3)%i_len
 
+      ! init
       interp__detect(:,:,:)=1
 
@@ -297,26 +324,32 @@
 
       ! do not compute point near fill value
-      FORALL( ji=1:il_dim(1):il_rho(jp_I), &
-      &       jj=1:il_dim(2):il_rho(jp_J), &
-      &       jk=1:il_dim(3):il_rho(jp_K), &
-      &       td_mix%d_value(ji,jj,jk,1) == td_mix%d_fill )
-
-         ! i-direction
-         interp__detect(MAX(1,ji-il_xextra):MIN(ji+il_xextra,il_dim(1)),&
-         &              MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
-         &              MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0 
-         ! j-direction
-         interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
-         &              MAX(1,jj-il_yextra):MIN(jj+il_yextra,il_dim(2)),&
-         &              MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0
-         ! k-direction
-         interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
-         &              MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
-         &              MAX(1,jk-il_zextra):MIN(jk+il_zextra,il_dim(3)) )=0         
-
-      END FORALL
+      DO jk=1,il_dim(3),il_rho(jp_K)
+         DO jj=1,il_dim(2),il_rho(jp_J)
+            DO ji=1,il_dim(1),il_rho(jp_I)
+
+               IF( td_mix%d_value(ji,jj,jk,1) == td_mix%d_fill )THEN
+
+                  ! i-direction
+                  interp__detect(MAX(1,ji-il_xextra):MIN(ji+il_xextra,il_dim(1)),&
+                  &              MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
+                  &              MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0 
+                  ! j-direction
+                  interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
+                  &              MAX(1,jj-il_yextra):MIN(jj+il_yextra,il_dim(2)),&
+                  &              MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0
+                  ! k-direction
+                  interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
+                  &              MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
+                  &              MAX(1,jk-il_zextra):MIN(jk+il_zextra,il_dim(3)) )=0         
+
+               ENDIF
+
+            ENDDO
+         ENDDO
+      ENDDO
+
+      DEALLOCATE( il_rho )
 
    END FUNCTION interp__detect
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
@@ -330,13 +363,10 @@
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in] td_var : coarse grid variable (should be extrapolated)
-   !> @param[out] td_mix : mixed grid variable
-   !> @param[in] id_rho  : table of refinment factor 
-   !>
-   !> @todo
-   !-------------------------------------------------------------------
-   !> @code
+   !> - November, 2013- Initial Version
+   !>
+   !> @param[in] td_var    coarse grid variable (should be extrapolated)
+   !> @param[out] td_mix   mixed grid variable
+   !> @param[in] id_rho    array of refinment factor (default 1) 
+   !-------------------------------------------------------------------
    SUBROUTINE interp_create_mixed_grid( td_var, td_mix, id_rho )
       IMPLICIT NONE
@@ -359,5 +389,5 @@
       ALLOCATE(il_rho(ip_maxdim))
       il_rho(:)=1
-      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
+      IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
 
       ! special case for even refinement on ARAKAWA-C grid
@@ -381,5 +411,5 @@
 
       ! copy variable
-      td_mix=td_var
+      td_mix=var_copy(td_var)
 
       ! compute new dimension length
@@ -408,6 +438,7 @@
       &     td_var%d_value(:,:,:,:)
 
+      DEALLOCATE(il_rho)
+
    END SUBROUTINE interp_create_mixed_grid
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
@@ -418,16 +449,13 @@
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]  td_mix : mixed grid variable
-   !> @param[in] id_rho  : table of refinment factor 
-   !>
-   !> @todo
-   !-------------------------------------------------------------------
-   !> @code
+   !> - November, 2013- Initial Version
+   !>
+   !> @param[inout] td_mix mixed grid variable
+   !> @param[in] id_rho    array of refinment factor 
+   !-------------------------------------------------------------------
    SUBROUTINE interp__clean_even_grid( td_mix, id_rho )
       IMPLICIT NONE
       ! Argument
-      TYPE(TVAR) , INTENT(INOUT) :: td_mix 
+      TYPE(TVAR) ,               INTENT(INOUT) :: td_mix 
       INTEGER(I4), DIMENSION(:), INTENT(IN   ), OPTIONAL :: id_rho
 
@@ -468,9 +496,9 @@
       END SELECT
 
-      ! copy variable
-      tl_mix=td_mix
-
       ! remove some point only if refinement in some direction is even
       IF( ANY(ll_even(:)) )THEN
+
+         ! copy variable
+         tl_mix=var_copy(td_mix)
 
          ALLOCATE( ll_mask( tl_mix%t_dim(1)%i_len, &
@@ -554,5 +582,8 @@
 
             td_mix%d_value(:,:,:,:)=RESHAPE( dl_vect(:), &
-            &                                SHAPE=td_mix%t_dim(:)%i_len )
+            &                                SHAPE=(/td_mix%t_dim(1)%i_len, &
+            &                                        td_mix%t_dim(2)%i_len, &
+            &                                        td_mix%t_dim(3)%i_len, &
+            &                                        td_mix%t_dim(4)%i_len/) )
 
             DEALLOCATE( dl_vect )
@@ -562,35 +593,38 @@
          DEALLOCATE( ll_mask )
 
+         ! clean
+         CALL var_clean(tl_mix)
+
       ENDIF
 
-      CALL var_clean(tl_mix)
+      ! clean
+      DEALLOCATE(il_rho)
 
    END SUBROUTINE interp__clean_even_grid
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This subroutine save interpolated value over domain. 
-   !> And so remove points added on mixed grid 
-   !> to compute interpolation
+   !> This subroutine remove points added on mixed grid 
+   !> to compute interpolation. And save interpolated value over domain. 
    !> 
    !> @details 
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in] td_var : table of mixed grid variable (to interpolate)
-   !> @param[in] id_rho : table of refinement factor
-   !> @return 
-   !-------------------------------------------------------------------
-   !> @code
+   !> - November, 2013- Initial Version
+   !> @date September, 2014
+   !> - use offset to save useful domain
+   !>
+   !> @param[in] td_mix    mixed grid variable structure
+   !> @param[out] td_var   variable structure
+   !> @param[in] id_rho    array of refinement factor (default 1)
+   !> @param[in] id_offset 2D array of offset between fine and coarse grid
+   !-------------------------------------------------------------------
    SUBROUTINE interp_clean_mixed_grid( td_mix, td_var, &
-   &                                   id_rho, &
-   &                                   id_offset )
+   &                                   id_rho, id_offset )
       IMPLICIT NONE
       ! Argument
       TYPE(TVAR)                 , INTENT(IN   ) :: td_mix
       TYPE(TVAR)                 , INTENT(  OUT) :: td_var
-      INTEGER(I4), DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_rho
-      INTEGER(I4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset
+      INTEGER(I4), DIMENSION(:)  , INTENT(IN   ) :: id_rho
+      INTEGER(I4), DIMENSION(2,2), INTENT(IN   ) :: id_offset
 
       ! local variable
@@ -605,6 +639,4 @@
       INTEGER(i4) :: il_jmax1
 
-      INTEGER(i4), DIMENSION(2,2) :: il_offset
-
       REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
 
@@ -613,15 +645,6 @@
       ! loop indices
       !----------------------------------------------------------------
-      il_offset(:,:)=0 
-      IF( PRESENT(id_offset) )THEN
-         IF( ANY( SHAPE(id_offset(:,:)) /= SHAPE(il_offset(:,:)) ) )THEN
-            CALL logger_error("INTERP CLEAN MIXED GRID: invalid dimension of"//&
-            &                 " offset table")
-         ELSE
-            il_offset(:,:)=id_offset(:,:)
-         ENDIF
-      ENDIF
       ! copy mixed variable in temporary structure
-      tl_mix=td_mix
+      tl_mix=var_copy(td_mix)
 
       ! remove unusefull points over mixed grid for even refinement
@@ -629,131 +652,62 @@
 
       ! copy cleaned mixed variable
-      td_var=tl_mix
-      IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
-
-         ! delete table of value
-         CALL var_del_value(td_var)
-
-         ! compute domain indices
-         il_imin0=1  ; il_imax0=td_var%t_dim(1)%i_len
-         il_jmin0=1  ; il_jmax0=td_var%t_dim(2)%i_len
-   
-         il_imin1=il_imin0+il_offset(1,1)
-         il_jmin1=il_jmin0+il_offset(2,1)
-   
-         il_imax1=il_imax0-il_offset(1,2)
-         il_jmax1=il_jmax0-il_offset(2,2)
-
-         SELECT CASE(TRIM(td_var%c_point))
-         CASE('U')
-            il_imin1=il_imin0+(il_offset(1,1)-1)
-            il_imax1=il_imax0-(il_offset(1,2)+1)
-         CASE('V')
-            il_jmin1=il_jmin0+(il_offset(2,1)-1)
-            il_jmax1=il_jmax0-(il_offset(2,2)+1)
-         CASE('F')
-            il_imin1=il_imin0+(il_offset(1,1)-1)
-            il_imax1=il_imax0-(il_offset(1,2)+1)
-
-            il_jmin1=il_jmin0+(il_offset(2,1)-1)
-            il_jmax1=il_jmax0-(il_offset(2,2)+1)
-         END SELECT
-
-         ! compute new dimension
-         td_var%t_dim(1)%i_len=il_imax1-il_imin1+1
-         td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1
+      td_var=var_copy(tl_mix)
+
+      ! delete array of value
+      CALL var_del_value(td_var)
+
+      ! compute domain indices in i-direction
+      il_imin0=1  ; il_imax0=td_var%t_dim(1)%i_len
  
-         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( 1:td_var%t_dim(1)%i_len, &
-         &         :,:,:) = tl_mix%d_value( il_imin1:il_imax1, &
-         &                                  il_jmin1:il_jmax1, &
-         &                                      :, : )
-
-         ! add variable value
-         CALL var_add_value(td_var,dl_value(:,:,:,:))
-
-         ! save variable type
-         td_var%i_type=tl_mix%i_type
-
-         DEALLOCATE(dl_value)
+      IF( td_var%t_dim(1)%l_use )THEN
+         il_imin1=il_imin0+id_offset(Jp_I,1)
+         il_imax1=il_imax0-id_offset(Jp_I,2)
+      ELSE
+
+         il_imin1=il_imin0
+         il_imax1=il_imax0
 
       ENDIF
-      
+
+      ! compute domain indices in j-direction
+      il_jmin0=1  ; il_jmax0=td_var%t_dim(2)%i_len
+ 
+      IF( td_var%t_dim(2)%l_use )THEN
+         il_jmin1=il_jmin0+id_offset(Jp_J,1)
+         il_jmax1=il_jmax0-id_offset(Jp_J,2)
+      ELSE
+
+         il_jmin1=il_jmin0
+         il_jmax1=il_jmax0
+ 
+      ENDIF
+
+      ! compute new dimension
+      td_var%t_dim(1)%i_len=il_imax1-il_imin1+1
+      td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1
+ 
+      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( 1:td_var%t_dim(1)%i_len, &
+      &         1:td_var%t_dim(2)%i_len, &
+      &         :,:) = tl_mix%d_value( il_imin1:il_imax1, &
+      &                                il_jmin1:il_jmax1, &
+      &                                     :, : )
+
+      ! add variable value
+      CALL var_add_value(td_var,dl_value(:,:,:,:))
+
+      DEALLOCATE(dl_value)
+ 
+      ! clean
       CALL var_clean(tl_mix)
 
    END SUBROUTINE interp_clean_mixed_grid
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This subroutine save interpolated value over domain. 
-   !> And so remove points added on mixed grid 
-   !> to compute interpolation
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in] td_var : table of mixed grid variable (to interpolate)
-   !> @param[in] id_offset : table of offset
-   !> @return 
-   !-------------------------------------------------------------------
-   !> @code
-   SUBROUTINE interp__del_offset( td_var, id_offset )
-      IMPLICIT NONE
-      ! Argument
-      TYPE(TVAR)                 , INTENT(INOUT) :: td_var
-      INTEGER(i4), DIMENSION(:,:), INTENT(IN   ) :: id_offset
-
-      ! local variable
-      INTEGER(i4) :: il_imin1
-      INTEGER(i4) :: il_jmin1
-      INTEGER(i4) :: il_imax1
-      INTEGER(i4) :: il_jmax1
-
-      REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
-
-      ! loop indices
-      !----------------------------------------------------------------
-   
-      IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
-
-         il_imin1=1+id_offset(1,1)
-         il_jmin1=1+id_offset(2,1)
-
-         il_imax1=td_var%t_dim(1)%i_len-id_offset(2,1)
-         il_jmax1=td_var%t_dim(2)%i_len-id_offset(2,2)
-
-         ! compute new dimension
-         td_var%t_dim(1)%i_len=il_imax1-il_imin1+1
-         td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1
-         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(:,:,:,:)=td_var%d_value( il_imin1:il_imax1, &
-         &                                 il_jmin1:il_jmax1, &
-         &                                 :,: )
-
-         ! delete table of value
-         CALL var_del_value(td_var)
-
-         ! add variable value
-         CALL var_add_value(td_var,dl_value(:,:,:,:))
-
-         DEALLOCATE(dl_value)
-
-      ENDIF
-
-   END SUBROUTINE interp__del_offset
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine interpolate detected point.
+   !> This subroutine interpolate variable value.
    !> 
    !> @details 
@@ -762,10 +716,10 @@
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in] td_var : variable to be interpolated
-   !> @param[in] id_rho : table of refinement factor
-   !-------------------------------------------------------------------
-   !> @code
+   !> - November, 2013- Initial Version
+   !>
+   !> @param[inout] td_var variable structure
+   !> @param[in] id_rho    array of refinement factor
+   !> @param[in] id_offset 2D array of offset between fine and coarse grid
+   !-------------------------------------------------------------------
    SUBROUTINE interp__fill_value_wrapper( td_var, & 
    &                                      id_rho, &
@@ -774,30 +728,18 @@
       ! Argument
       TYPE(TVAR) ,                 INTENT(INOUT) :: td_var
-      INTEGER(I4), DIMENSION(:),   INTENT(IN   ), OPTIONAL :: id_rho
+      INTEGER(I4), DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_rho
       INTEGER(I4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset
 
       ! local variable
-      CHARACTER(LEN=lc)                 :: cl_method
-      INTEGER(i4)      , DIMENSION(:), ALLOCATABLE  :: il_rho
-      INTEGER(i4)      , DIMENSION(2,2) :: il_offset
+      CHARACTER(LEN=lc)                            :: cl_method
+      INTEGER(i4)      , DIMENSION(:), ALLOCATABLE :: il_rho
+      INTEGER(i4)      , DIMENSION(2,2)            :: il_offset
 
       ! loop indices
       !----------------------------------------------------------------
 
-      SELECT CASE(TRIM(td_var%c_interp(1)))
-         CASE('cubic','linear','nearest')
-            cl_method=TRIM(td_var%c_interp(1))
-         CASE DEFAULT
-            CALL logger_warn("INTERP FILL: interpolation method unknown."//&
-            &  " use linear interpolation")
-            cl_method='linear'
-
-            ! update variable structure value
-            td_var%c_interp(1)='linear'
-      END SELECT
-
       ALLOCATE( il_rho(ip_maxdim) )
       il_rho(:)=1
-      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
+      IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
       IF( ANY(il_rho(:) < 0) )THEN
          CALL logger_error("INTERP FILL VALUE: invalid "//&
@@ -808,5 +750,5 @@
       IF( PRESENT(id_offset) )THEN
          IF( ANY(SHAPE(id_offset(:,:)) /= (/2,2/)) )THEN
-            CALL logger_error("INTERP FILL VALUE: invalid table of offset")
+            CALL logger_error("INTERP FILL VALUE: invalid array of offset")
          ELSE
             il_offset(:,:)=id_offset(:,:)
@@ -817,4 +759,15 @@
       &   (il_rho(jp_J) /= 1 .AND. td_var%t_dim(2)%l_use) .OR. &
       &   (il_rho(jp_K) /= 1 .AND. td_var%t_dim(3)%l_use) )THEN
+
+         SELECT CASE(TRIM(td_var%c_interp(1)))
+            CASE('cubic','linear','nearest')
+               cl_method=TRIM(td_var%c_interp(1))
+            CASE DEFAULT
+               CALL logger_warn("INTERP FILL VALUE: interpolation method unknown."//&
+               &  " use linear interpolation")
+               cl_method='linear'
+               ! update variable structure value
+               td_var%c_interp(1)='linear'
+         END SELECT
 
          CALL logger_info("INTERP FILL: interpolate "//TRIM(td_var%c_name)//&
@@ -826,6 +779,5 @@
  
          CALL interp__fill_value( td_var, cl_method, & 
-         &                        il_rho(:), &
-         &                        il_offset(:,:) )
+         &                        il_rho(:), il_offset(:,:) )
 
          SELECT CASE(TRIM(td_var%c_interp(2)))
@@ -864,26 +816,27 @@
          END SELECT
 
+      ELSE
+         td_var%c_interp(:)=''
       ENDIF
 
+      DEALLOCATE(il_rho)
+
    END SUBROUTINE interp__fill_value_wrapper
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
    !> This subroutine interpolate value over mixed grid.
    !> 
-   !> @details 
-   !> 
-   !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[inout] td_var : variable
-   !> @param[in] cd_method : interpolation method 
-   !> @param[in] id_rho : table of refinment factor 
-   !-------------------------------------------------------------------
-   !> @code
+   !> - November, 2013- Initial Version
+   !> @date September, 2014
+   !> - use interpolation method modules
+   !>
+   !> @param[inout] td_var variable structure
+   !> @param[in] cd_method interpolation method 
+   !> @param[in] id_rho    array of refinment factor 
+   !> @param[in] id_offset 2D array of offset between fine and coarse grid
+   !-------------------------------------------------------------------
    SUBROUTINE interp__fill_value( td_var, cd_method, &
-   &                              id_rho, &
-   &                              id_offset )
+   &                              id_rho, id_offset )
       IMPLICIT NONE
       ! Argument
@@ -891,5 +844,5 @@
       CHARACTER(LEN=*)                , INTENT(IN   ) :: cd_method
       INTEGER(I4)     , DIMENSION(:)  , INTENT(IN   ) :: id_rho
-      INTEGER(I4)     , DIMENSION(:,:), INTENT(IN   ) :: id_offset
+      INTEGER(I4)     , DIMENSION(2,2), INTENT(IN   ) :: id_offset
 
       ! local variable
@@ -908,10 +861,6 @@
 
       TYPE(TATT)                                       :: tl_att
-
+ 
       ! loop indices
-      INTEGER(i4) :: ji
-      INTEGER(i4) :: jj
-      INTEGER(i4) :: jk
-      INTEGER(i4) :: jl
       !----------------------------------------------------------------
 
@@ -935,4 +884,7 @@
       tl_att=att_init('interpolation',cl_interp)
       CALL var_move_att(tl_mix, tl_att)
+
+      ! clean
+      CALL att_clean(tl_att)
 
       ! special case for even refinement on ARAKAWA-C grid
@@ -972,79 +924,20 @@
 
       !3- interpolate
-      DO jl=1,tl_mix%t_dim(4)%i_len
-         IF( il_rho(jp_K) /= 1 )THEN
-!            CALL interp__3D(tl_mix%d_value(:,:,:,jl), tl_mix%d_fill, &
-!              &             il_detect(:,:,:), cd_method,             &
-!              &             il_rhoi, il_rhoj, il_rhok,               &
-!              &             ll_even(:), ll_discont )
-             CALL logger_error("INTERP FILL: can not interpolate "//&
-             &                "vertically for now ")
-         ENDIF
-
-         IF( ANY(il_detect(:,:,:)==1) )THEN
-            ! I-J plan
-            DO jk=1,tl_mix%t_dim(3)%i_len
-              CALL interp__2D(tl_mix%d_value(:,:,jk,jl), tl_mix%d_fill, &
-              &               il_detect(:,:,jk), cd_method,             &
-              &               il_rho(jp_I), il_rho(jp_J),               &
-              &               ll_even(1:2), ll_discont)
-            ENDDO
-            IF( ALL(il_detect(:,:,:)==0) ) CYCLE
-            IF( il_rho(jp_K) /= 1 )THEN
-               ! I-K plan
-               DO jj=1,tl_mix%t_dim(2)%i_len
-                  CALL interp__2D(tl_mix%d_value(:,jj,:,jl), tl_mix%d_fill, &
-                  &               il_detect(:,jj,:), cd_method,             &
-                  &               il_rho(jp_J), il_rho(jp_K),               &
-                  &               ll_even(1:3:2), ll_discont  )
-               ENDDO
-               IF( ALL(il_detect(:,:,:)==0) ) CYCLE
-               ! J-K plan
-               DO ji=1,tl_mix%t_dim(1)%i_len
-                  CALL interp__2D(tl_mix%d_value(ji,:,:,jl), tl_mix%d_fill, &
-                  &               il_detect(ji,:,:), cd_method,             &
-                  &               il_rho(jp_J), il_rho(jp_K),               &
-                  &               ll_even(2:3), ll_discont )
-               ENDDO
-
-            ENDIF
-            IF( ANY(il_detect(:,:,:)==1) )THEN
-               ! I direction
-               DO jk=1,tl_mix%t_dim(3)%i_len
-                  DO jj=1,tl_mix%t_dim(2)%i_len
-                     CALL interp__1D( tl_mix%d_value(:,jj,jk,jl),    &
-                     &                tl_mix%d_fill,                 &
-                     &                il_detect(:,jj,jk), cd_method, &
-                     &                il_rho(jp_I),                  &
-                     &                ll_even(1), ll_discont )
-                  ENDDO
-               ENDDO
-               IF( ALL(il_detect(:,:,:)==0) ) CYCLE
-               ! J direction
-               DO jk=1,tl_mix%t_dim(3)%i_len
-                  DO ji=1,tl_mix%t_dim(1)%i_len
-                     CALL interp__1D( tl_mix%d_value(ji,:,jk,jl),    &
-                     &                tl_mix%d_fill,                 &
-                     &                il_detect(ji,:,jk), cd_method, &
-                     &                il_rho(jp_J),                  &
-                     &                ll_even(2), ll_discont )
-                  ENDDO
-               ENDDO
-               IF( il_rho(jp_K) /= 1 )THEN
-                  IF( ALL(il_detect(:,:,:)==0) ) CYCLE
-                  ! K direction
-                  DO jj=1,tl_mix%t_dim(2)%i_len
-                     DO ji=1,tl_mix%t_dim(1)%i_len
-                        CALL interp__1D( tl_mix%d_value(ji,jj,:,jl),    &
-                        &                tl_mix%d_fill,                 & 
-                        &                il_detect(ji,jj,:), cd_method, &
-                        &                il_rho(jp_K),                  &
-                        &                ll_even(3), ll_discont  )
-                     ENDDO
-                  ENDDO
-               ENDIF
-            ENDIF
-         ENDIF
-      ENDDO
+      CALL logger_debug("INTERP 2D: interpolation method "//TRIM(cd_method)//&
+      &  " discont "//TRIM(fct_str(ll_discont)) )
+      SELECT CASE(TRIM(cd_method))
+      CASE('cubic')
+         CALL interp_cubic_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, &
+           &                    il_detect(:,:,:),                        &
+           &                    il_rho(:), ll_even(:), ll_discont )
+      CASE('nearest')
+         CALL interp_nearest_fill(tl_mix%d_value(:,:,:,:), &
+              &                   il_detect(:,:,:),         &
+              &                   il_rho(:) )
+      CASE DEFAULT ! linear
+         CALL interp_linear_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, &
+              &                  il_detect(:,:,:),                        &
+              &                  il_rho(:), ll_even(:), ll_discont )
+      END SELECT         
 
       IF( ANY(il_detect(:,:,:)==1) )THEN
@@ -1056,1267 +949,10 @@
       !4- save useful domain (remove offset)
       CALL interp_clean_mixed_grid( tl_mix, td_var, &
-      &                             id_rho(:),      &
-      &                             id_offset(:,:) )      
+      &                             id_rho(:), id_offset(:,:)  )  
 
       ! clean variable structure
+      DEALLOCATE(il_rho)
       CALL var_clean(tl_mix)
 
    END SUBROUTINE interp__fill_value
-   !> @endcode
-!   !-------------------------------------------------------------------
-!   !> @brief
-!   !> This subroutine 
-!   !> 
-!   !> @details 
-!   !>
-!   !> @author J.Paul
-!   !> - Nov, 2013- Initial Version
-!   !
-!   !> @param[in]
-!   !> @param[out] 
-!   !-------------------------------------------------------------------
-!   !> @code
-!   SUBROUTINE interp__3D( dd_value,  dd_fill,   &
-!      &                   id_detect, cd_method, &
-!      &                   id_rhoi, id_rhoj, id_rhok, &
-!      &                   ld_even,  ld_discont )
-!      IMPLICIT NONE
-!      ! Argument
-!      REAL(dp)        , DIMENSION(:,:,:), INTENT(INOUT) :: dd_value 
-!      REAL(dp)                          , INTENT(IN   ) :: dd_fill 
-!      INTEGER(I4)     , DIMENSION(:,:,:), INTENT(INOUT) :: id_detect
-!      CHARACTER(LEN=*)                  , INTENT(IN   ) :: cd_method
-!      INTEGER(I4)                       , INTENT(IN   ) :: id_rhoi
-!      INTEGER(I4)                       , INTENT(IN   ) :: id_rhoj
-!      INTEGER(I4)                       , INTENT(IN   ) :: id_rhok
-!      LOGICAL         , DIMENSION(:)    , INTENT(IN   ) :: ld_even
-!      LOGICAL                           , INTENT(IN   ), OPTIONAL :: ld_discont
-!
-!      ! local variable
-!      INTEGER(i4), DIMENSION(3)                  :: il_shape
-!      INTEGER(i4), DIMENSION(3)                  :: il_dim
-!
-!      REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_coarse
-!      REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx
-!      REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdy
-!      REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_d2fdxy
-!
-!      LOGICAL                                    :: ll_discont
-!
-!      ! loop indices
-!!      INTEGER(i4) :: ji
-!!      INTEGER(i4) :: jj
-!!      INTEGER(i4) :: ii
-!!      INTEGER(i4) :: ij
-!      !----------------------------------------------------------------
-!      ll_discont=.FALSE.
-!      IF( PRESENT(ld_discont) ) ll_discont=ld_discont
-!
-!      il_shape(:)=SHAPE(dd_value)
-!      il_dim(:)=(il_shape(:)-1)/2
-!
-!      ALLOCATE( dl_coarse(il_dim(1),il_dim(2),il_dim(3)) )
-!      ! value on coarse grid
-!      dl_coarse(:,:,:)=dd_value( 1:il_shape(1):id_rhoi, &
-!      &                        1:il_shape(2):id_rhoj, &
-!      &                        1:il_shape(3):id_rhok )
-!      SELECT CASE(TRIM(cd_method))
-!
-!         CASE('cubic')
-!            
-!            ALLOCATE( dl_dfdx(  il_dim(1),il_dim(2),il_dim(3)), &
-!            &         dl_dfdy(  il_dim(1),il_dim(2),il_dim(3)), &
-!            &         dl_d2fdxy(il_dim(1),il_dim(2),il_dim(3)) )
-!
-!!            ! compute derivative on coarse grid
-!!            dl_dfdx(:,:)=extrap_deriv_2D(dl_coarse(:,:,:), dd_fill, 'I')
-!!            dl_dfdy(:,:)=extrap_deriv_2D(dl_coarse(:,:,:), dd_fill, 'J')
-!!
-!!            ! compute cross derivative on coarse grid
-!!            dl_d2fdxy(:,:)=extrap_deriv_2D(dl_dfdx(:,:,:), dd_fill, 'J')
-!!
-!!            DO jj=1,il_shape(2)-1,id_rhoj
-!!               ij=((jj-1)/id_rhoj)+1
-!!               DO ji=1,il_shape(1)-1,id_rhoi
-!!                  ii=((ji-1)/id_rhoi)+1
-!!            
-!!                  IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE
-!!                  ! compute bicubic coefficient
-!!                  dl_coef(:)=interp__2D_cubic_coef(dl_coarse(ii:ii+1,ij:ij+1),&
-!!                  &                                dl_dfdx(  ii:ii+1,ij:ij+1),&
-!!                  &                                dl_dfdy(  ii:ii+1,ij:ij+1),&
-!!                  &                                dl_d2fdxy(ii:ii+1,ij:ij+1) )
-!!
-!!                  ! compute value on detetected point
-!!                  CALL interp__2D_cubic_fill(dl_coef(:), &
-!!                  &                          dd_value( ji:ji+id_rhoi,   &
-!!                  &                                    jj:jj+id_rhoj ), &
-!!                  &                          id_detect(ji:ji+id_rhoi,   &
-!!                  &                                    jj:jj+id_rhoj )  )
-!!
-!!               ENDDO
-!!            ENDDO
-!
-!            DEALLOCATE( dl_dfdx, &
-!            &           dl_dfdy, &
-!            &           dl_d2fdxy )
-!
-!         CASE('nearest')
-!
-!         CASE DEFAULT ! linear 
-!
-!!            DO jj=1,il_shape(2)-1,id_rhoj
-!!               ij=((jj-1)/id_rhoj)+1
-!!               DO ji=1,il_shape(1)-1,id_rhoi
-!!                  ii=((ji-1)/id_rhoi)+1            
-!!
-!!                  IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE
-!!                  ! compute bilinear coefficient
-!!                  dl_coef(:)=interp__2D_linear_coef(dl_coarse(ii:ii+1,ij:ij+1))
-!!
-!!                  ! compute value on detetected point
-!!                  CALL interp__2D_linear_fill(dl_coef(:), &
-!!                  &                          dd_value( ji:ji+id_rhoi,   &
-!!                  &                                    jj:jj+id_rhoj ), &
-!!                  &                          id_detect(ji:ji+id_rhoi,   &
-!!                  &                                    jj:jj+id_rhoj )  )
-!!
-!!               ENDDO
-!!            ENDDO
-!
-!      END SELECT
-!
-!      DEALLOCATE( dl_coarse )
-!
-!   END SUBROUTINE interp__3D
-!   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   SUBROUTINE interp__2D( dd_value,  dd_fill,   &
-      &                   id_detect, cd_method, &
-      &                   id_rhoi, id_rhoj,     &
-      &                   ld_even,  ld_discont )
-
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value 
-      REAL(dp)                        , INTENT(IN   ) :: dd_fill 
-      INTEGER(I4)     , DIMENSION(:,:), INTENT(INOUT) :: id_detect
-      CHARACTER(LEN=*)                , INTENT(IN   ) :: cd_method
-      INTEGER(I4)                     , INTENT(IN   ) :: id_rhoi
-      INTEGER(I4)                     , INTENT(IN   ) :: id_rhoj
-      LOGICAL         , DIMENSION(:)  , INTENT(IN   ) :: ld_even
-      LOGICAL                         , INTENT(IN   ), OPTIONAL :: ld_discont
-
-      ! local variable
-      INTEGER(I4)                              :: il_xextra
-      INTEGER(I4)                              :: il_yextra
-      INTEGER(i4), DIMENSION(2)                :: il_shape
-      INTEGER(i4), DIMENSION(2)                :: il_dim
-
-      REAL(dp)                                 :: dl_min 
-      REAL(dp)                                 :: dl_max 
-      REAL(dp)   , DIMENSION(:)  , ALLOCATABLE :: dl_coef
-      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_coarse
-      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_tmp
-      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_dfdx
-      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_dfdy
-      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_d2fdxy
-
-      LOGICAL                                  :: ll_discont
-      ! loop indices
-      INTEGER(i4) :: ji
-      INTEGER(i4) :: jj
-      INTEGER(i4) :: ii
-      INTEGER(i4) :: ij
-
-      !----------------------------------------------------------------
-      ll_discont=.FALSE.
-      IF( PRESENT(ld_discont) ) ll_discont=ld_discont
-
-      CALL logger_debug("INTERP 2D: interpolation method "//TRIM(cd_method)//&
-      &  " discont "//TRIM(fct_str(ll_discont)) )
-
-      il_shape(:)=SHAPE(dd_value)
-
-      ! compute coarse grid dimension
-      il_xextra=id_rhoi-1
-      il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi
-
-      il_yextra=id_rhoj-1
-      il_dim(2)=(il_shape(2)+il_yextra)/id_rhoj
-
-      ALLOCATE( dl_coarse(il_dim(1),il_dim(2)) )
-
-      ! value on coarse grid
-      dl_coarse(:,:)=dd_value( 1:il_shape(1):id_rhoi, &
-      &                        1:il_shape(2):id_rhoj )
-
-      SELECT CASE(TRIM(cd_method))
-
-         CASE('cubic')
-            
-            ALLOCATE( dl_dfdx(  il_dim(1),il_dim(2)), &
-            &         dl_dfdy(  il_dim(1),il_dim(2)), &
-            &         dl_d2fdxy(il_dim(1),il_dim(2)) )
-
-            ! compute derivative on coarse grid
-            dl_dfdx(:,:)=extrap_deriv_2D(dl_coarse(:,:), dd_fill, 'I', ll_discont)
-            dl_dfdy(:,:)=extrap_deriv_2D(dl_coarse(:,:), dd_fill, 'J', ll_discont)
-
-            ! compute cross derivative on coarse grid
-            dl_d2fdxy(:,:)=extrap_deriv_2D(dl_dfdx(:,:), dd_fill, 'J', ll_discont)
-
-            ALLOCATE( dl_tmp(2,2) )
-            ALLOCATE( dl_coef(16) )
-
-            DO jj=1,il_shape(2)-1,id_rhoj
-               ij=((jj-1)/id_rhoj)+1
-               DO ji=1,il_shape(1)-1,id_rhoi
-                  ii=((ji-1)/id_rhoi)+1
-            
-                  IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill) .OR. &
-                  &   ANY(  dl_dfdx(ii:ii+1,ij:ij+1)==dd_fill) .OR. &
-                  &   ANY(  dl_dfdy(ii:ii+1,ij:ij+1)==dd_fill) .OR. &
-                  &   ANY(dl_d2fdxy(ii:ii+1,ij:ij+1)==dd_fill) ) CYCLE
-
-                  dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1)
-                  IF( ll_discont )THEN
-
-                     dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
-                     dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
-                     IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
-                        WHERE( dl_tmp(:,:) < 0_dp ) 
-                           dl_tmp(:,:) = dl_tmp(:,:)+360._dp
-                        END WHERE
-                     ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
-                        WHERE( dl_tmp(:,:) > 180_dp ) 
-                           dl_tmp(:,:) = dl_tmp(:,:)-180._dp
-                        END WHERE
-                     ENDIF
-                  ENDIF
-
-                  ! compute bicubic coefficient
-                  dl_coef(:)=interp__2D_cubic_coef(dl_tmp(:,:),&
-                  &                                dl_dfdx(   ii:ii+1,ij:ij+1),&
-                  &                                dl_dfdy(   ii:ii+1,ij:ij+1),&
-                  &                                dl_d2fdxy( ii:ii+1,ij:ij+1),&
-                  &                                dd_fill )
-
-                  ! compute value on detetected point
-                  CALL interp__2D_cubic_fill(dd_value( ji:ji+id_rhoi,   &
-                  &                                    jj:jj+id_rhoj ), &
-                  &                          id_detect(ji:ji+id_rhoi,   &
-                  &                                    jj:jj+id_rhoj ), &
-                  &                          dl_coef(:), dd_fill,       &
-                  &                          ld_even(:), id_rhoi, id_rhoj )
-
-                  IF( ll_discont )THEN
-                     WHERE( dd_value( ji:ji+id_rhoi, &
-                        &             jj:jj+id_rhoj ) >= 180._dp .AND. &
-                        &   dd_value( ji:ji+id_rhoi, &
-                        &             jj:jj+id_rhoj ) /= dd_fill )
-                        dd_value( ji:ji+id_rhoi, &
-                        &         jj:jj+id_rhoj ) = &
-                        &           dd_value( ji:ji+id_rhoi, &
-                        &                     jj:jj+id_rhoj ) - 360._dp
-                     END WHERE
-                  ENDIF
-
-               ENDDO
-            ENDDO
-
-            DEALLOCATE(dl_coef)
-            DEALLOCATE(dl_tmp )
-
-            DEALLOCATE(dl_dfdx, &
-            &          dl_dfdy, &
-            &          dl_d2fdxy )
-
-         CASE('nearest') 
-
-            DO jj=1,il_shape(2)-1,id_rhoj
-               ij=((jj-1)/id_rhoj)+1
-               DO ji=1,il_shape(1)-1,id_rhoi
-                  ii=((ji-1)/id_rhoi)+1            
-
-                  ! compute value on detetected point
-                  CALL interp__2D_nearest_fill(dd_value( ji:ji+id_rhoi,   &
-                  &                                      jj:jj+id_rhoj ), &
-                  &                            id_detect(ji:ji+id_rhoi,   &
-                  &                                      jj:jj+id_rhoj )  )
-
-               ENDDO
-            ENDDO
-
-         CASE DEFAULT ! linear 
-
-            ALLOCATE( dl_coef(4)  )
-            ALLOCATE( dl_tmp(2,2) )
-
-            DO jj=1,il_shape(2)-1,id_rhoj
-               ij=((jj-1)/id_rhoj)+1
-               DO ji=1,il_shape(1)-1,id_rhoi
-                  ii=((ji-1)/id_rhoi)+1            
-
-                  IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE
-
-                  dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1)
-                  IF( ll_discont )THEN
-
-                     dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
-                     dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
-                     IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
-                        WHERE( dl_tmp(:,:) < 0_dp ) 
-                           dl_tmp(:,:) = dl_tmp(:,:)+360._dp
-                        END WHERE
-                     ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
-                        WHERE( dl_tmp(:,:) > 180_dp ) 
-                           dl_tmp(:,:) = dl_tmp(:,:)-180._dp
-                        END WHERE
-                     ENDIF
-                  ENDIF
-
-                  ! compute bilinear coefficient
-                  dl_coef(:)=interp__2D_linear_coef(dl_tmp(:,:), dd_fill)
-
-                  ! compute value on detetected point
-                  CALL interp__2D_linear_fill(dd_value( ji:ji+id_rhoi,   &
-                  &                                     jj:jj+id_rhoj ), &
-                  &                           id_detect(ji:ji+id_rhoi,   &
-                  &                                     jj:jj+id_rhoj ), &
-                  &                           dl_coef(:), dd_fill,       &
-                  &                           ld_even(:), id_rhoi, id_rhoj )
-
-                  IF( ll_discont )THEN
-                     WHERE( dd_value( ji:ji+id_rhoi, &
-                        &             jj:jj+id_rhoj ) >= 180._dp .AND. &
-                        &   dd_value( ji:ji+id_rhoi, &
-                        &             jj:jj+id_rhoj ) /= dd_fill )
-                        dd_value( ji:ji+id_rhoi, &
-                        &         jj:jj+id_rhoj ) = &
-                        &           dd_value( ji:ji+id_rhoi, &
-                        &                     jj:jj+id_rhoj ) - 360._dp
-                     END WHERE
-                  ENDIF
-
-               ENDDO
-            ENDDO
-
-            DEALLOCATE(dl_coef)
-            DEALLOCATE(dl_tmp )
-
-      END SELECT
-
-      DEALLOCATE( dl_coarse )
-
-   END SUBROUTINE interp__2D
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   SUBROUTINE interp__1D( dd_value,  dd_fill,   &
-      &                   id_detect, cd_method, &
-      &                   id_rhoi,              &
-      &                   ld_even, ld_discont )
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp)        , DIMENSION(:), INTENT(INOUT) :: dd_value 
-      REAL(dp)                      , INTENT(IN   ) :: dd_fill 
-      INTEGER(I4)     , DIMENSION(:), INTENT(INOUT) :: id_detect
-      CHARACTER(LEN=*)              , INTENT(IN   ) :: cd_method
-      INTEGER(I4)                   , INTENT(IN   ) :: id_rhoi
-      LOGICAL                       , INTENT(IN   ) :: ld_even
-      LOGICAL                       , INTENT(IN   ), OPTIONAL :: ld_discont
-
-      ! local variable
-      INTEGER(i4), DIMENSION(1)              :: il_shape
-      INTEGER(i4), DIMENSION(1)              :: il_dim
-      INTEGER(I4)                            :: il_xextra
-
-      REAL(dp)                               :: dl_min 
-      REAL(dp)                               :: dl_max 
-      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_coarse
-      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_tmp
-      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_dfdx
-
-      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_coef
-
-      LOGICAL                                :: ll_discont
-
-      ! loop indices
-      INTEGER(i4) :: ji
-      INTEGER(i4) :: ii
-      !----------------------------------------------------------------
-      ll_discont=.FALSE.
-      IF( PRESENT(ld_discont) ) ll_discont=ld_discont
-
-      CALL logger_debug("INTERP 1D: interpolation method "//TRIM(cd_method)//&
-      &  " discont "//TRIM(fct_str(ll_discont)) )
-
-      il_shape(:)=SHAPE(dd_value)
-
-      il_xextra=id_rhoi-1
-      il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi
-
-      ALLOCATE( dl_coarse(il_dim(1)) )
-      ! value on coarse grid
-      dl_coarse(:)=dd_value( 1:il_shape(1):id_rhoi )
-
-      SELECT CASE(TRIM(cd_method))
-
-         CASE('cubic')
-            
-            ALLOCATE( dl_dfdx( il_dim(1)) )
-
-            ! compute derivative on coarse grid
-            dl_dfdx(:)=extrap_deriv_1D(dl_coarse(:), dd_fill, ll_discont)
-
-            ALLOCATE( dl_coef(4))
-            ALLOCATE( dl_tmp(2) )
-
-            DO ji=1,il_shape(1)-1,id_rhoi
-               ii=((ji-1)/id_rhoi)+1
-            
-               IF( ANY( dl_tmp(:)==dd_fill ) .OR. &
-               &   ANY(dl_dfdx(:)==dd_fill ) ) CYCLE
-
-               dl_tmp(:)=dl_coarse(ii:ii+1)
-               IF( ll_discont )THEN
-
-                  dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
-                  dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
-                  IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
-                     WHERE( dl_tmp(:) < 0_dp ) 
-                        dl_tmp(:) = dl_tmp(:)+360._dp
-                     END WHERE
-                  ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
-                     WHERE( dl_tmp(:) > 180_dp ) 
-                        dl_tmp(:) = dl_tmp(:)-180._dp
-                     END WHERE
-                  ENDIF
-               ENDIF
-
-               ! compute bicubic coefficient
-               dl_coef(:)=interp__1D_cubic_coef(dl_tmp(:),&
-               &                                dl_dfdx(   ii:ii+1), &
-               &                                dd_fill )
-
-               ! compute value on detetected point
-               CALL interp__1D_cubic_fill(dd_value( ji:ji+id_rhoi),&
-               &                          id_detect(ji:ji+id_rhoi),&
-               &                          dl_coef(:), dd_fill,     &
-               &                          ld_even, id_rhoi )
-
-               IF( ll_discont )THEN
-                  WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. &
-                     &   dd_value( ji:ji+id_rhoi ) /= dd_fill )
-                     dd_value( ji:ji+id_rhoi ) = &
-                     &           dd_value( ji:ji+id_rhoi ) - 360._dp
-                  END WHERE
-               ENDIF
-
-            ENDDO
-
-            DEALLOCATE(dl_coef)
-            DEALLOCATE(dl_tmp )
-
-         CASE('nearest') 
-
-            DO ji=1,il_shape(1)-1,id_rhoi
-               ii=((ji-1)/id_rhoi)+1            
-
-               ! compute value on detetected point
-               CALL interp__1D_nearest_fill(dd_value( ji:ji+id_rhoi), &
-               &                            id_detect(ji:ji+id_rhoi)  )
-
-            ENDDO
-
-         CASE DEFAULT ! linear 
-
-            ALLOCATE(dl_coef(2))
-            ALLOCATE( dl_tmp(2) )
-
-            DO ji=1,il_shape(1)-1,id_rhoi
-               ii=((ji-1)/id_rhoi)+1            
-
-               IF( ANY(dl_coarse(ii:ii+1)==dd_fill ) ) CYCLE
-
-               dl_tmp(:)=dl_coarse(ii:ii+1)
-               IF( ll_discont )THEN
-
-                  dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
-                  dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
-                  IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
-                     WHERE( dl_tmp(:) < 0_dp ) 
-                        dl_tmp(:) = dl_tmp(:)+360._dp
-                     END WHERE
-                  ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
-                     WHERE( dl_tmp(:) > 180_dp ) 
-                        dl_tmp(:) = dl_tmp(:)-180._dp
-                     END WHERE
-                  ENDIF
-               ENDIF
-
-               ! compute bilinear coefficient
-               dl_coef(:)=interp__1D_linear_coef(dl_tmp(:), dd_fill)
-
-               ! compute value on detetected point
-               CALL interp__1D_linear_fill( dd_value( ji:ji+id_rhoi),&
-               &                            id_detect(ji:ji+id_rhoi),&
-               &                            dl_coef(:), dd_fill,     &
-               &                            ld_even, id_rhoi  )
-
-               IF( ll_discont )THEN
-                  WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. &
-                     &   dd_value( ji:ji+id_rhoi ) /= dd_fill )
-                     dd_value( ji:ji+id_rhoi ) = &
-                     &           dd_value( ji:ji+id_rhoi ) - 360._dp
-                  END WHERE
-               ENDIF
-
-            ENDDO
-
-            DEALLOCATE(dl_coef)
-
-      END SELECT
-
-      DEALLOCATE( dl_coarse )
-
-   END SUBROUTINE interp__1D
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   FUNCTION interp__2D_cubic_coef( dd_value, &
-      &                            dd_dfdx,  &
-      &                            dd_dfdy,  &
-      &                            dd_d2fdxy,&
-      &                            dd_fill )
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_value 
-      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdx  
-      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdy  
-      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_d2fdxy
-      REAL(dp)                , INTENT(IN) :: dd_fill
-
-      ! function
-      REAL(dp), DIMENSION(16) :: interp__2D_cubic_coef
-
-      ! local variable
-      REAL(dp), DIMENSION(16,16), PARAMETER :: dp_matrix = RESHAPE( &
-      & (/ 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,& 
-           0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,& 
-          -3 , 3 , 0 , 0 ,-2 ,-1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,& 
-           2 ,-2 , 0 , 0 , 1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,& 
-           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,& 
-           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 ,& 
-           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,-3 , 3 , 0 , 0 ,-2 ,-1 , 0 , 0 ,& 
-           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 ,-2 , 0 , 0 , 1 , 1 , 0 , 0 ,& 
-          -3 , 0 , 3 , 0 , 0 , 0 , 0 , 0 ,-2 , 0 ,-1 , 0 , 0 , 0 , 0 , 0 ,& 
-           0 , 0 , 0 , 0 ,-3 , 0 , 3 , 0 , 0 , 0 , 0 , 0 ,-2 , 0 ,-1 , 0 ,& 
-           9 ,-9 ,-9 , 9 , 6 , 3 ,-6 ,-3 , 6 ,-6 , 3 ,-3 , 4 , 2 , 2 , 1 ,& 
-          -6 , 6 , 6 ,-6 ,-3 ,-3 , 3 , 3 ,-4 , 4 ,-2 , 2 ,-2 ,-2 ,-1 ,-1 ,& 
-           2 , 0 ,-2 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 ,& 
-           0 , 0 , 0 , 0 , 2 , 0 ,-2 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 ,& 
-          -6 , 6 , 6 ,-6 ,-4 ,-2 , 4 , 2 ,-3 , 3 ,-3 , 3 ,-2 ,-1 ,-2 ,-1 ,& 
-           4 ,-4 ,-4 , 4 , 2 , 2 ,-2 ,-2 , 2 ,-2 , 2 ,-2 , 1 , 1 , 1 , 1 /), &
-      & (/ 16, 16 /) )
-
-      REAL(dp), DIMENSION(16) :: dl_vect
-
-      !----------------------------------------------------------------
-      ! init
-      interp__2D_cubic_coef(:)=dd_fill
-
-      IF( ANY(SHAPE( dd_value(:,:))/= 2) .OR. &
-      &   ANY(SHAPE(  dd_dfdx(:,:))/= 2) .OR. &
-      &   ANY(SHAPE(  dd_dfdy(:,:))/= 2) .OR. &
-      &   ANY(SHAPE(dd_d2fdxy(:,:))/= 2) )THEN
-   
-         CALL logger_error("INTERP CUBIC COEF: invalid dimension of "//&
-         &              "input tables. shape should be (/2,2/)")
-
-      ELSEIF( ANY( dd_value(:,:) == dd_fill) .OR. &
-      &       ANY(  dd_dfdx(:,:) == dd_fill) .OR. &
-      &       ANY(  dd_dfdy(:,:) == dd_fill) .OR. &
-      &       ANY(dd_d2fdxy(:,:) == dd_fill) )THEN
-      
-         CALL logger_warn("INTERP CUBIC COEF: fill value detected. "//&
-         &  "can not compute coefficient ")
-
-      ELSE
-
-         dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. )
-         dl_vect( 5: 8)=PACK(dd_dfdx(:,:),.TRUE. )
-         dl_vect( 9:12)=PACK(dd_dfdy(:,:),.TRUE. )
-         dl_vect(13:16)=PACK(dd_d2fdxy(:,:),.TRUE. )
-
-         interp__2D_cubic_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))
-
-      ENDIF
-   END FUNCTION interp__2D_cubic_coef
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   SUBROUTINE interp__2D_cubic_fill( dd_value, id_detect, dd_coef, dd_fill, &
-   &                                 ld_even, id_rhoi, id_rhoj )
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value 
-      INTEGER(i4)     , DIMENSION(:,:), INTENT(INOUT) :: id_detect
-      REAL(dp)        , DIMENSION(:)  , INTENT(IN   ) :: dd_coef 
-      REAL(dp)                        , INTENT(IN   ) :: dd_fill 
-      LOGICAL         , DIMENSION(:),   INTENT(IN   ) :: ld_even
-      INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoi
-      INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoj
-
-      ! local variable
-      INTEGER(i4), DIMENSION(2)  :: il_shape
-
-      REAL(dp)   , DIMENSION(16) :: dl_vect
-      REAL(dp)                   :: dl_dx
-      REAL(dp)                   :: dl_dy
-      REAL(dp)                   :: dl_x
-      REAL(dp)                   :: dl_x2
-      REAL(dp)                   :: dl_x3
-      REAL(dp)                   :: dl_y
-      REAL(dp)                   :: dl_y2
-      REAL(dp)                   :: dl_y3
-
-      ! loop indices
-      INTEGER(i4) :: ji
-      INTEGER(i4) :: jj
-      !----------------------------------------------------------------
-
-      IF( SIZE(dd_coef(:)) /= 16 )THEN
-         CALL logger_error("INTERP CUBIC FILL: invalid dimension of "//&
-         &              "coef table. shape should be (/16/)")
-      ELSEIF( ANY(   dd_coef(:)==dd_fill ) )THEN
-         CALL logger_error("INTERP CUBIC FILL: fill value detected in coef . "//&
-         &              "can not compute interpolation.")
-      ELSE
-
-         il_shape(:)=SHAPE(dd_value(:,:))
-
-         dl_dx=1./REAL(id_rhoi)
-         dl_dy=1./REAL(id_rhoj)
-
-         DO jj=1,il_shape(2)
-
-            IF( ld_even(2) )THEN
-               dl_y=(jj-1)*dl_dy - dl_dy*0.5
-            ELSE ! odd refinement
-               dl_y=(jj-1)*dl_dy
-            ENDIF
-            dl_y2=dl_y*dl_y
-            dl_y3=dl_y2*dl_y
-
-            DO ji=1,il_shape(1)
-               
-               IF(id_detect(ji,jj)==1)THEN
-                  
-                  IF( ld_even(1) )THEN
-                     dl_x=(ji-1)*dl_dx - dl_dx*0.5 
-                  ELSE ! odd refinement
-                     dl_x=(ji-1)*dl_dx 
-                  ENDIF
-                  dl_x2=dl_x*dl_x
-                  dl_x3=dl_x2*dl_x
-
-                  dl_vect(:)=(/1._dp, dl_x      , dl_x2      , dl_x3      , &
-                  &            dl_y , dl_x*dl_y , dl_x2*dl_y , dl_x3*dl_y , &
-                  &            dl_y2, dl_x*dl_y2, dl_x2*dl_y2, dl_x3*dl_y2, &
-                  &            dl_y3, dl_x*dl_y3, dl_x2*dl_y3, dl_x3*dl_y3 /)
-
-                  dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dl_vect(:))
-                  id_detect(ji,jj)=0
-
-               ENDIF
-
-            ENDDO
-         ENDDO
-
-      ENDIF
-
-   END SUBROUTINE interp__2D_cubic_fill
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   FUNCTION interp__2D_linear_coef( dd_value, dd_fill )
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp), DIMENSION(:,:)  , INTENT(IN) :: dd_value 
-      REAL(dp)                  , INTENT(IN) :: dd_fill 
-
-      ! function
-      REAL(dp), DIMENSION(4) :: interp__2D_linear_coef
-
-      ! local variable
-
-      REAL(dp), DIMENSION(4,4), PARAMETER :: dp_matrix = RESHAPE( &
-      & (/  1 , 0 , 0 , 0 ,& 
-           -1 , 1 , 0 , 0 ,& 
-           -1 , 0 , 1 , 0 ,& 
-            1 ,-1 ,-1 , 1 /), &
-      & (/ 4, 4 /) )
-
-      REAL(dp), DIMENSION(4) :: dl_vect
-
-      !----------------------------------------------------------------
-
-      IF( ANY(SHAPE(dd_value(:,:))/= 2) )THEN
-         CALL logger_error("INTERP LINEAR COEF: invalid dimension of "//&
-         &              "input tables. shape should be (/2,2/)")
-      ELSEIF( ANY(dd_value(:,:)==dd_fill) )THEN
-         CALL logger_error("INTERP LINEAR COEF: fill value detected. "//&
-         &              "can not compute coefficient.")
-      ELSE
-
-         dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. )
-
-         interp__2D_linear_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))
-
-      ENDIF
-   END FUNCTION interp__2D_linear_coef
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   SUBROUTINE interp__2D_linear_fill( dd_value, id_detect, dd_coef, dd_fill, &
-   &                                  ld_even, id_rhoi, id_rhoj )
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value 
-      INTEGER(i4)     , DIMENSION(:,:), INTENT(INOUT) :: id_detect
-      REAL(dp)        , DIMENSION(:)  , INTENT(IN   ) :: dd_coef 
-      REAL(dp)                        , INTENT(IN   ) :: dd_fill 
-      LOGICAL         , DIMENSION(:),   INTENT(IN   ) :: ld_even
-      INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoi
-      INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoj
-
-      ! local variable
-      INTEGER(i4), DIMENSION(2) :: il_shape
-
-      REAL(dp)   , DIMENSION(4) :: dl_vect
-      REAL(dp)                  :: dl_dx
-      REAL(dp)                  :: dl_dy
-      REAL(dp)                  :: dl_x
-      REAL(dp)                  :: dl_y
-
-      ! loop indices
-      INTEGER(i4) :: ji
-      INTEGER(i4) :: jj
-      !----------------------------------------------------------------
-
-      IF( SIZE(dd_coef(:)) /= 4 )THEN
-         CALL logger_error("INTERP LINEAR FILL: invalid dimension of "//&
-         &              "coef table. shape should be (/4/)")
-      ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN
-         CALL logger_error("INTERP LINEAR FILL: fill value detected in coef. "//&
-         &              "can not compute interpolation.")
-      ELSE
-
-         il_shape(:)=SHAPE(dd_value(:,:))
-
-         dl_dx=1./REAL(id_rhoi)
-         dl_dy=1./REAL(id_rhoj)
-
-         DO jj=1,il_shape(2)
-
-            IF( ld_even(2) )THEN
-               dl_y=(jj-1)*dl_dy - dl_dy*0.5
-            ELSE ! odd refinement
-               dl_y=(jj-1)*dl_dy
-            ENDIF
-
-            DO ji=1,il_shape(1)
-               
-               IF(id_detect(ji,jj)==1)THEN
-
-                  IF( ld_even(1) )THEN
-                     dl_x=(ji-1)*dl_dx - dl_dx*0.5 
-                  ELSE ! odd refinement
-                     dl_x=(ji-1)*dl_dx 
-                  ENDIF
-
-                  dl_vect(:)=(/1._dp, dl_x, dl_y, dl_x*dl_y /)
-
-                  dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dl_vect(:))
-                  id_detect(ji,jj)=0
-
-               ENDIF
-
-            ENDDO
-         ENDDO
-
-      ENDIF
-
-   END SUBROUTINE interp__2D_linear_fill
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   SUBROUTINE interp__2D_nearest_fill( dd_value, id_detect )
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp)   , DIMENSION(:,:)  , INTENT(INOUT) :: dd_value 
-      INTEGER(i4), DIMENSION(:,:)  , INTENT(INOUT) :: id_detect
-
-      ! local variable
-      INTEGER(i4), DIMENSION(2) :: il_shape
-
-      INTEGER(i4) :: il_i1
-      INTEGER(i4) :: il_i2
-      INTEGER(i4) :: il_j1
-      INTEGER(i4) :: il_j2
-
-      INTEGER(i4) :: il_half1
-      INTEGER(i4) :: il_half2
-
-      ! loop indices
-      INTEGER(i4) :: ji
-      INTEGER(i4) :: jj
-      !----------------------------------------------------------------
-
-      il_shape(:)=SHAPE(dd_value(:,:))
-
-      il_i1=1
-      il_i2=il_shape(1)
-
-      il_j1=1
-      il_j2=il_shape(2)
-
-      il_half1=CEILING(il_shape(1)*0.5)
-      il_half2=CEILING(il_shape(2)*0.5)
-
-      DO jj=1,il_half2
-
-         DO ji=1,il_half1
-            
-            ! lower left point
-            IF(id_detect(ji,jj)==1)THEN
-
-               dd_value( ji,jj)=dd_value(il_i1,il_j1)
-               id_detect(ji,jj)=0
-
-            ENDIF
-
-            ! lower right point
-            IF(id_detect(il_shape(1)-ji+1,jj)==1)THEN
-
-               dd_value( il_shape(1)-ji+1,jj)=dd_value(il_i2,il_j1)
-               id_detect(il_shape(1)-ji+1,jj)=0
-
-            ENDIF
-
-            ! upper left point
-            IF(id_detect(ji,il_shape(2)-jj+1)==1)THEN
-
-               dd_value( ji,il_shape(2)-jj+1)=dd_value(il_i1,il_j2)
-               id_detect(ji,il_shape(2)-jj+1)=0
-
-            ENDIF            
-
-            ! upper right point
-            IF(id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)==1)THEN
-
-               dd_value( il_shape(1)-ji+1,il_shape(2)-jj+1)=dd_value(il_i2,il_j2)
-               id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)=0
-
-            ENDIF            
-
-         ENDDO
-
-      ENDDO
-
-   END SUBROUTINE interp__2D_nearest_fill
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   FUNCTION interp__1D_cubic_coef( dd_value, &
-      &                            dd_dfdx,  &
-      &                            dd_fill )
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp), DIMENSION(:)  , INTENT(IN) :: dd_value 
-      REAL(dp), DIMENSION(:)  , INTENT(IN) :: dd_dfdx  
-      REAL(dp)                , INTENT(IN) :: dd_fill  
-
-      ! function
-      REAL(dp), DIMENSION(4) :: interp__1D_cubic_coef
-
-      ! local variable
-
-      REAL(dp), DIMENSION(4,4), PARAMETER :: dp_matrix = RESHAPE( &
-      & (/  1 , 0 , 0 , 0 ,& 
-           -1 , 1 , 0 , 0 ,& 
-           -3 , 3 ,-2 ,-1 ,& 
-            2 ,-2 , 1 , 1  /), &
-      & (/ 4, 4 /) )
-
-      REAL(dp), DIMENSION(4) :: dl_vect
-
-      !----------------------------------------------------------------
-      IF( SIZE(dd_value(:))/= 2 .OR. &
-      &   SIZE(dd_dfdx(:) )/= 2  )THEN
-   
-         CALL logger_error("INTERP CUBIC COEF: invalid dimension of "//&
-         &              "input tables. shape should be (/2,2/)")
-
-      ELSEIF( ANY(dd_value(:)==dd_fill) .OR. &
-      &       ANY( dd_dfdx(:)==dd_fill) )THEN
-         CALL logger_error("INTERP CUBIC COEF: fill value detected. "//&
-         &              "can not compute coefficient.")
-      ELSE
-
-         dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. )
-         dl_vect( 3: 4)=PACK(dd_dfdx(:),.TRUE. )
-
-         interp__1D_cubic_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))
-
-      ENDIF
-
-   END FUNCTION interp__1D_cubic_coef
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   SUBROUTINE interp__1D_cubic_fill( dd_value, id_detect, dd_coef, dd_fill, &
-   &                                 ld_even, id_rhoi )
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp)        , DIMENSION(:), INTENT(INOUT) :: dd_value 
-      INTEGER(i4)     , DIMENSION(:), INTENT(INOUT) :: id_detect
-      REAL(dp)        , DIMENSION(:), INTENT(IN   ) :: dd_coef 
-      REAL(dp)                      , INTENT(IN   ) :: dd_fill 
-      LOGICAL         ,               INTENT(IN   ) :: ld_even
-      INTEGER(I4)     ,               INTENT(IN   ) :: id_rhoi
-
-      ! local variable
-      INTEGER(i4), DIMENSION(1)  :: il_shape
-
-      REAL(dp)   , DIMENSION(4) :: dl_vect
-      REAL(dp)                   :: dl_dx
-      REAL(dp)                   :: dl_x
-      REAL(dp)                   :: dl_x2
-      REAL(dp)                   :: dl_x3
-
-      ! loop indices
-      INTEGER(i4) :: ji
-      !----------------------------------------------------------------
-
-      IF( SIZE(dd_coef(:)) /= 4 )THEN
-         CALL logger_error("INTERP CUBIC FILL: invalid dimension of "//&
-         &              "coef table. shape should be (/4/)")
-      !ELSEIF( ANY(dd_value(:)==dd_fill .AND. id_detect(:)==0 ) .OR. &
-      ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN
-         CALL logger_error("INTERP CUBIC FILL: fill value detected. "//&
-         &              "can not compute interpolation")
-      ELSE
-
-         il_shape(:)=SHAPE(dd_value(:))
-
-         dl_dx=1./REAL(id_rhoi)
-         
-         DO ji=1,il_shape(1)
-            
-            IF(id_detect(ji)==1)THEN
-
-               IF( ld_even )THEN
-                  dl_x=(ji-1)*dl_dx - dl_dx*0.5 
-               ELSE ! odd refinement
-                  dl_x=(ji-1)*dl_dx 
-               ENDIF
-               dl_x2=dl_x*dl_x
-               dl_x3=dl_x2*dl_x
-
-               dl_vect(:)=(/1._dp, dl_x, dl_x2, dl_x3 /)
-
-               dd_value(ji)=DOT_PRODUCT(dd_coef(:),dl_vect(:))
-               id_detect(ji)=0
-
-            ENDIF
-
-         ENDDO
-
-      ENDIF
-
-   END SUBROUTINE interp__1D_cubic_fill
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   FUNCTION interp__1D_linear_coef( dd_value, dd_fill )
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp), DIMENSION(:)  , INTENT(IN) :: dd_value 
-      REAL(dp)                , INTENT(IN) :: dd_fill  
-
-      ! function
-      REAL(dp), DIMENSION(2) :: interp__1D_linear_coef
-
-      ! local variable
-
-      REAL(dp), DIMENSION(2,2), PARAMETER :: dp_matrix = RESHAPE( &
-      & (/  1 , 0 ,& 
-           -1 , 1  /), &
-      &  (/ 2, 2 /) )
-
-      REAL(dp), DIMENSION(2) :: dl_vect
-
-      !----------------------------------------------------------------
-
-      IF( ANY(SHAPE(dd_value(:))/= 2) )THEN
-         CALL logger_error("INTERP LINEAR COEF: invalid dimension of "//&
-         &              "input tables. shape should be (/2/)")
-      ELSEIF( ANY(dd_value(:)==dd_fill) )THEN
-         CALL logger_error("INTERP LINEAR COEF: fill value detected. "//&
-         &              "can not compute coefficient.")
-      ELSE
-         
-         dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. )
-
-         interp__1D_linear_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))
-
-      ENDIF
-
-   END FUNCTION interp__1D_linear_coef
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   SUBROUTINE interp__1D_linear_fill( dd_value, id_detect, dd_coef, dd_fill, &
-   &                                  ld_even, id_rhoi )
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp)        , DIMENSION(:), INTENT(INOUT) :: dd_value 
-      INTEGER(i4)     , DIMENSION(:), INTENT(INOUT) :: id_detect
-      REAL(dp)        , DIMENSION(:), INTENT(IN   ) :: dd_coef 
-      REAL(dp)                      , INTENT(IN   ) :: dd_fill 
-      LOGICAL         ,               INTENT(IN   ) :: ld_even
-      INTEGER(I4)     ,               INTENT(IN   ) :: id_rhoi      
-
-      ! local variable
-      INTEGER(i4), DIMENSION(1) :: il_shape
-
-      REAL(dp)   , DIMENSION(2) :: dl_vect
-      REAL(dp)                  :: dl_dx
-      REAL(dp)                  :: dl_x
-
-      ! loop indices
-      INTEGER(i4) :: ji
-      !----------------------------------------------------------------
-
-      IF( SIZE(dd_coef(:)) /= 2 )THEN
-         CALL logger_error("INTERP LINEAR FILL: invalid dimension of "//&
-         &              "coef table. shape should be (/2/)")
-      !ELSEIF( ANY(dd_value(:)==dd_fill .AND. id_detect(:)==0 ) .OR. &
-      ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN
-         CALL logger_error("INTERP LINEAR FILL: fill value detected. "//&
-         &              "can not compute interpolation")
-      ELSE
-
-         il_shape(:)=SHAPE(dd_value)
-
-         dl_dx=1./REAL(id_rhoi)
-
-         DO ji=1,il_shape(1)
-            
-            IF(id_detect(ji)==1)THEN
-
-               IF( ld_even )THEN
-                  dl_x=(ji-1)*dl_dx - dl_dx*0.5 
-               ELSE ! odd refinement
-                  dl_x=(ji-1)*dl_dx 
-               ENDIF
-
-               dl_vect(:)=(/1._dp, dl_x /)
-
-               dd_value(ji)=DOT_PRODUCT(dd_coef(:),dl_vect(:))
-               id_detect(ji)=0
-
-            ENDIF
-
-         ENDDO
-
-      ENDIF
-
-   END SUBROUTINE interp__1D_linear_fill
-   !> @endcode
-   !-------------------------------------------------------------------
-   !> @brief
-   !> This subroutine 
-   !> 
-   !> @details 
-   !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in]
-   !> @param[out] 
-   !-------------------------------------------------------------------
-   !> @code
-   SUBROUTINE interp__1D_nearest_fill( dd_value, id_detect )
-      IMPLICIT NONE
-      ! Argument
-      REAL(dp)   , DIMENSION(:), INTENT(INOUT) :: dd_value 
-      INTEGER(i4), DIMENSION(:), INTENT(INOUT) :: id_detect
-
-      ! local variable
-      INTEGER(i4), DIMENSION(1) :: il_shape
-
-      INTEGER(i4) :: il_i1
-      INTEGER(i4) :: il_i2
-
-      INTEGER(i4) :: il_half1
-      
-      ! loop indices
-      INTEGER(i4) :: ji
-      !----------------------------------------------------------------
-
-      il_shape(:)=SHAPE(dd_value)
-
-      il_i1=1
-      il_i2=il_shape(1)
-
-      il_half1=CEILING(il_shape(1)*0.5)
-      
-      DO ji=1,il_half1
-         
-         ! lower left point
-         IF(id_detect(ji)==1)THEN
-
-            dd_value( ji)=dd_value(il_i1)
-            id_detect(ji)=0
-
-         ENDIF
-
-         ! lower right point
-         IF(id_detect(il_shape(1)-ji+1)==1)THEN
-
-            dd_value( il_shape(1)-ji+1)=dd_value(il_i2)
-            id_detect(il_shape(1)-ji+1)=0
-
-         ENDIF
-
-      ENDDO
-
-   END SUBROUTINE interp__1D_nearest_fill
-   !> @endcode
 END MODULE interp
