Index: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/TOOLS/SIREN/src/extrap.f90
===================================================================
--- branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/TOOLS/SIREN/src/extrap.f90	(revision 4763)
+++ branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/TOOLS/SIREN/src/extrap.f90	(revision 5214)
@@ -7,7 +7,86 @@
 ! DESCRIPTION:
 !> @brief 
-!> This module 
+!> This module manage extrapolation.
 !>
 !> @details
+!>    Extrapolation method to be used is specify inside variable 
+!>    strcuture, as array of string character.<br/> 
+!>    - td_var\%c_extrap(1) string character is the interpolation name choose between:
+!>       - 'dist_weight'
+!>       - 'min_error'
+!>
+!>    @note Extrapolation method could be specify for each variable in namelist _namvar_,
+!>    defining string character _cn\_varinfo_. By default _dist_weight_.<br/>
+!>    Example:
+!>       - cn_varinfo='varname1:dist_weight', 'varname2:min_error'
+!>
+!>    to detect point to be extrapolated:<br/>
+!> @code
+!>    il_detect(:,:,:)=extrap_detect(td_var, [td_level], [id_offset,] [id_rho,] [id_ext]) 
+!> @endcode
+!>       - il_detect(:,:,:) is 3D array of point to be extrapolated
+!>       - td_var  is coarse grid variable to be extrapolated
+!>       - td_level is fine grid array of level (see vgrid_get_level) [optional]
+!>       - id_offset is array of offset between fine and coarse grid [optional]
+!>       - id_rho    is array of refinment factor [optional]
+!>       - id_ext    is array of number of points to be extrapolated [optional]
+!>
+!>    to extrapolate variable:<br/>
+!> @code
+!>    CALL extrap_fill_value( td_var, [td_level], [id_offset], [id_rho], [id_iext], [id_jext], [id_kext], [id_radius], [id_maxiter])
+!> @endcode
+!>       - td_var  is coarse grid variable to be extrapolated
+!>       - td_level is fine grid array of level (see vgrid_get_level) [optional]
+!>       - id_offset is array of offset between fine and coarse grid [optional]
+!>       - id_rho    is array of refinment factor [optional]
+!>       - id_iext   is number of points to be extrapolated in i-direction [optional]
+!>       - id_jext   is number of points to be extrapolated in j-direction [optional]
+!>       - id_kext   is number of points to be extrapolated in k-direction [optional]
+!>       - id_radius is radius of the halo used to compute extrapolation [optional]
+!>       - id_maxiter is maximum number of iteration [optional]
+!>
+!>    to add extraband to the variable (to be extrapolated):<br/>
+!> @code
+!>    CALL extrap_add_extrabands(td_var, [id_isize,] [id_jsize] )
+!> @endcode
+!>       - td_var is variable structure
+!>       - id_isize : i-direction size of extra bands [optional]
+!>       - id_jsize : j-direction size of extra bands [optional]
+!>
+!>    to delete extraband of a variable:<br/>
+!> @code
+!>    CALL extrap_del_extrabands(td_var, [id_isize,] [id_jsize] )
+!> @endcode
+!>       - td_var is variable structure
+!>       - id_isize : i-direction size of extra bands [optional]
+!>       - id_jsize : j-direction size of extra bands [optional]
+!>
+!>    to compute first derivative of 1D array:<br/>
+!> @code
+!>    dl_value(:)=extrap_deriv_1D( dd_value(:), dd_fill, [ld_discont] )
+!> @endcode
+!>       - dd_value is 1D array of variable
+!>       - dd_fill is FillValue of variable
+!>       - ld_discont is logical to take into account longitudinal East-West discontinuity [optional]
+!>
+!>    to compute first derivative of 2D array:<br/>
+!> @code
+!>    dl_value(:,:)=extrap_deriv_2D( dd_value(:,:), dd_fill, cd_dim, [ld_discont] )
+!> @endcode
+!>       - dd_value is 2D array of variable
+!>       - dd_fill is FillValue of variable
+!>       - cd_dim is character to compute derivative on first (I) or second (J) dimension
+!>       - ld_discont is logical to take into account longitudinal East-West discontinuity [optional]
+!>
+!>    to compute first derivative of 3D array:<br/>
+!> @code
+!>    dl_value(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, cd_dim, [ld_discont] )
+!> @endcode
+!>       - dd_value is 3D array of variable
+!>       - dd_fill is FillValue of variable
+!>       - cd_dim is character to compute derivative on first (I), second (J), or third (K) dimension
+!>       - ld_discont is logical to take into account longitudinal East-West discontinuity [optional]
+!>
+!> @warning _FillValue must not be zero (use var_chg_FillValue())
 !>
 !> @author
@@ -15,47 +94,50 @@
 ! REVISION HISTORY:
 !> @date Nov, 2013 - Initial Version
+!> @date September, 2014
+!> - add header
 !>
-!> @note WARNING: FillValue must not be zero (use var_chg_FillValue)
+!> @todo
+!> - create module for each extrapolation method
 !>
 !> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-!> @todo
-!> - something wrong when computing point to be extralopated
-!> - take care of ew value in variable structure
 !----------------------------------------------------------------------
 MODULE extrap
    USE netcdf                          ! nf90 library
-   USE kind
-   USE phycst
-   USE global
-   USE fct
-   USE logger
-   USE dim
-   USE att
-   USE var
+   USE kind                            ! F90 kind parameter
+   USE phycst                          ! physical constant
+   USE global                          ! global variable
+   USE fct                             ! basic useful function
+   USE date                            ! date manager
+   USE logger                          ! log file manager
+   USE att                             ! attribute manager
+   USE dim                             ! dimension manager
+   USE var                             ! variable manager
 
    IMPLICIT NONE
-   PRIVATE
    ! NOTE_avoid_public_variables_if_possible
 
    ! type and variable
+   PRIVATE :: im_maxiter   !< default maximum number of iteration 
+   PRIVATE :: im_minext    !< default minumum number of point to extrapolate
+   PRIVATE :: im_mincubic  !< default minumum number of point to extrapolate for cubic interpolation
 
    ! function and subroutine
    PUBLIC :: extrap_detect         !< detected point to be extrapolated
    PUBLIC :: extrap_fill_value     !< extrapolate value over detected point 
-   PUBLIC :: extrap_add_extrabands !< 
-   PUBLIC :: extrap_del_extrabands !< 
-   PUBLIC :: extrap_deriv_1D !< 
-   PUBLIC :: extrap_deriv_2D !< 
-
-   PRIVATE :: extrap__detect_wrapper !< detected point to be extrapolated
-   PRIVATE :: extrap__detect         !< detected point to be extrapolated
-   PRIVATE :: extrap__fill_value_wrapper     !< extrapolate value over detected point 
-   PRIVATE :: extrap__fill_value     !< extrapolate value over detected point 
-   PRIVATE :: extrap__3D
-   PRIVATE :: extrap_deriv_3D
-   PRIVATE :: extrap__3D_min_error_coef
-   PRIVATE :: extrap__3D_min_error_fill
-   PRIVATE :: extrap__3D_dist_weight_coef
-   PRIVATE :: extrap__3D_dist_weight_fill
+   PUBLIC :: extrap_add_extrabands !< add extraband to the variable (to be extrapolated) 
+   PUBLIC :: extrap_del_extrabands !< delete extraband of the variable 
+   PUBLIC :: extrap_deriv_1D       !< compute first derivative of 1D array 
+   PUBLIC :: extrap_deriv_2D       !< compute first derivative of 2D array 
+   PUBLIC :: extrap_deriv_3D       !< compute first derivative of 3D array
+
+   PRIVATE :: extrap__detect_wrapper      ! detected point to be extrapolated wrapper
+   PRIVATE :: extrap__detect              ! detected point to be extrapolated
+   PRIVATE :: extrap__fill_value_wrapper  ! extrapolate value over detected point wrapper 
+   PRIVATE :: extrap__fill_value          ! extrapolate value over detected point 
+   PRIVATE :: extrap__3D                  !
+   PRIVATE :: extrap__3D_min_error_coef   !
+   PRIVATE :: extrap__3D_min_error_fill   !
+   PRIVATE :: extrap__3D_dist_weight_coef !
+   PRIVATE :: extrap__3D_dist_weight_fill ! 
 
    INTEGER(i4), PARAMETER :: im_maxiter = 10 !< default maximum number of iteration
@@ -64,5 +146,5 @@
 
    INTERFACE extrap_detect
-      MODULE PROCEDURE extrap__detect_wrapper !< detected point to be extrapolated
+      MODULE PROCEDURE extrap__detect_wrapper     !< detected point to be extrapolated
    END INTERFACE extrap_detect
 
@@ -74,23 +156,33 @@
    !-------------------------------------------------------------------
    !> @brief
-   !> This function detected point to be extrapolated.
+   !> This function detected point to be extrapolated, given variable structure.
    !> 
    !> @details 
+   !> optionaly, you could sepcify fine grid level, refinment factor (default 1), 
+   !> offset between fine and coarse grid (default compute from refinment factor
+   !> as offset=(rho-1)/2), number of point to be extrapolated in each direction
+   !> (default im_minext).<br/>
    !>
+   !> First coarsening fine grid level, if need be, then select point near 
+   !> grid point already inform.
+   !>
+   !> @note point to be extrapolated are selected using FillValue, 
+   !> so to avoid mistake FillValue should not be zero (use var_chg_FillValue)
+   !> 
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
+   !> - November, 2013- Initial Version
    !
-   !> @param[in] td_var : variable to extrapolate
-   !> @param[in] id_iext : number of points to be extrapolated in i-direction
-   !> @param[in] id_jext : number of points to be extrapolated in j-direction
-   !> @param[in] id_kext : number of points to be extrapolated in k-direction
-   !> @return
-   !-------------------------------------------------------------------
-   !> @code
+   !> @param[in] td_var0   coarse grid variable to extrapolate
+   !> @param[in] td_level1 fine grid array of level
+   !> @param[in] id_offset array of offset between fine and coarse grid 
+   !> @param[in] id_rho    array of refinment factor 
+   !> @param[in] id_ext    array of number of points to be extrapolated
+   !> @return array of point to be extrapolated
+   !-------------------------------------------------------------------
    FUNCTION extrap__detect( td_var0, td_level1, &
    &                        id_offset, id_rho, id_ext )
       IMPLICIT NONE
       ! Argument
-      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var0
+      TYPE(TVAR) ,                 INTENT(IN   ) :: td_var0
       TYPE(TVAR) , DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: td_level1
       INTEGER(i4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset
@@ -106,6 +198,7 @@
       CHARACTER(LEN=lc)                                :: cl_level
 
-      INTEGER(i4)                                      :: il_varid
+      INTEGER(i4)                                      :: il_ind
       INTEGER(i4)      , DIMENSION(:,:,:), ALLOCATABLE :: il_detect
+      INTEGER(i4)      , DIMENSION(:,:,:), ALLOCATABLE :: il_tmp
       INTEGER(i4)      , DIMENSION(:,:)  , ALLOCATABLE :: il_offset
       INTEGER(i4)      , DIMENSION(:,:)  , ALLOCATABLE :: il_level1
@@ -143,9 +236,10 @@
       ALLOCATE( il_offset(ip_maxdim,2) )
       il_offset(:,:)=0
-      il_offset(jp_I,:)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5)
-      il_offset(jp_J,:)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5)
       IF( PRESENT(id_offset) )THEN
          il_offset(1:SIZE(id_offset(:,:),DIM=1),&
          &         1:SIZE(id_offset(:,:),DIM=2) )= id_offset(:,:)
+      ELSE
+         il_offset(jp_I,:)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5)
+         il_offset(jp_J,:)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5)
       ENDIF
 
@@ -160,19 +254,26 @@
 
       ! select point already inform
-      WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=1
+      DO jk0=1,td_var0%t_dim(3)%i_len
+         DO jj0=1,td_var0%t_dim(2)%i_len
+            DO ji0=1,td_var0%t_dim(1)%i_len
+               IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill ) il_detect(ji0,jj0,jk0)=1
+            ENDDO
+         ENDDO
+      ENDDO
  
       IF( PRESENT(td_level1) )THEN
          SELECT CASE(TRIM(td_var0%c_point))
-         CASE DEFAULT !'T'
-            cl_level='tlevel'
-         CASE('U')
-            cl_level='ulevel'
-         CASE('V')
-            cl_level='vlevel'
-         CASE('F')
-            cl_level='flevel'
+            CASE DEFAULT !'T'
+               cl_level='tlevel'
+            CASE('U')
+               cl_level='ulevel'
+            CASE('V')
+               cl_level='vlevel'
+            CASE('F')
+               cl_level='flevel'
          END SELECT
-         il_varid=var_get_id(td_level1(:),TRIM(cl_level))
-         IF( il_varid == 0 )THEN
+
+         il_ind=var_get_index(td_level1(:),TRIM(cl_level))
+         IF( il_ind == 0 )THEN
             CALL logger_error("EXTRAP DETECT: can not compute point to be "//&
             &     "extrapolated for variable "//TRIM(td_var0%c_name)//&
@@ -180,84 +281,66 @@
             &     "level for variable point "//TRIM(TRIM(td_var0%c_point)))
          ELSE
-            print *,'read ',trim(cl_level)
-            tl_var1=td_level1(il_varid)
+            tl_var1=var_copy(td_level1(il_ind))
+
+            ALLOCATE( il_level1_G0( il_dim0(1), il_dim0(2)) )
+            IF( ALL(tl_var1%t_dim(1:2)%i_len == il_dim0(1:2)) )THEN
+
+               ! variable to be extrapolated use same resolution than level
+               il_level1_G0(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)
+               
+            ELSE
+               ! variable to be extrapolated do not use same resolution than level
+               ALLOCATE( il_level1(tl_var1%t_dim(1)%i_len, &
+               &                   tl_var1%t_dim(2)%i_len) )
+               ! match fine grid vertical level with coarse grid
+               il_level1(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)/il_rho(jp_K)
+
+               ALLOCATE( il_extra(ip_maxdim,2) )
+               ! coarsening fine grid level
+               il_extra(jp_I,1)=CEILING(REAL(il_rho(jp_I)-1,dp)*0.5_dp)
+               il_extra(jp_I,2)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5_dp)
+
+               il_extra(jp_J,1)=CEILING(REAL(il_rho(jp_J)-1,dp)*0.5_dp)
+               il_extra(jp_J,2)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5_dp)
+
+               DO jj0=1,td_var0%t_dim(2)%i_len
+                  
+                  jj1=(jj0-1)*il_rho(jp_J)+1-il_offset(jp_J,1)
+
+                  jj1m=MAX( jj1-il_extra(jp_J,1), 1 )
+                  jj1p=MIN( jj1+il_extra(jp_J,2), &
+                  &         tl_var1%t_dim(2)%i_len-il_offset(jp_J,2) )
+                  
+                  DO ji0=1,td_var0%t_dim(1)%i_len
+
+                     ji1=(ji0-1)*il_rho(jp_I)+1-id_offset(jp_I,1)
+
+                     ji1m=MAX( ji1-il_extra(jp_I,1), 1 )
+                     ji1p=MIN( ji1+il_extra(jp_I,2), &
+                     &         tl_var1%t_dim(1)%i_len-id_offset(jp_I,2) )
+               
+                     il_level1_G0(ji0,jj0)=MAXVAL(il_level1(ji1m:ji1p,jj1m:jj1p))
+
+                  ENDDO
+               ENDDO
+
+               ! clean
+               DEALLOCATE( il_extra )
+               DEALLOCATE( il_level1 )
+
+            ENDIF
+
+            ! look for sea point
+            DO jk0=1,td_var0%t_dim(3)%i_len
+               WHERE( il_level1_G0(:,:) >= jk0)
+                  il_detect(:,:,jk0)=1
+               END WHERE
+            ENDDO
+
+            ! clean
+            DEALLOCATE( il_level1_G0 )
+            CALL var_clean(tl_var1)
+
          ENDIF
-
-         ALLOCATE( il_level1_G0( il_dim0(1), il_dim0(2)) )
-         IF( ALL(tl_var1%t_dim(1:2)%i_len == il_dim0(1:2)) )THEN
-
-            ! variable to be extrapolated use same resolution than level
-            il_level1_G0(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)
-            
-         ELSE
-            ! variable to be extrapolated do not use same resolution than level
-            ALLOCATE( il_level1(tl_var1%t_dim(1)%i_len, &
-            &                   tl_var1%t_dim(2)%i_len) )
-            ! match fine grid vertical level with coarse grid
-            il_level1(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)/il_rho(jp_K)
-
-            ALLOCATE( il_extra(ig_ndim,2) )
-            ! coarsening fine grid level
-            il_extra(jp_I,1)=CEILING(REAL(il_rho(jp_I)-1,dp)/2._dp)
-            il_extra(jp_I,2)=FLOOR(REAL(il_rho(jp_I)-1,dp)/2._dp)
-
-            il_extra(jp_J,1)=CEILING(REAL(il_rho(jp_J)-1,dp)/2._dp)
-            il_extra(jp_J,2)=FLOOR(REAL(il_rho(jp_J)-1,dp)/2._dp)
-
-            DO jj0=1,td_var0%t_dim(2)%i_len
-               
-               jj1=(jj0-1)*il_rho(jp_J)+1-il_offset(jp_J,1)
-
-               jj1m=MAX( jj1-il_extra(jp_J,1), 1 )
-               jj1p=MIN( jj1+il_extra(jp_J,2), &
-               &         tl_var1%t_dim(2)%i_len-il_offset(jp_J,2) )
-               
-               DO ji0=1,td_var0%t_dim(1)%i_len
-
-                  ji1=(ji0-1)*il_rho(jp_I)+1-id_offset(jp_I,1)
-
-                  ji1m=MAX( ji1-il_extra(jp_I,1), 1 )
-                  ji1p=MIN( ji1+il_extra(jp_I,2), &
-                  &         tl_var1%t_dim(1)%i_len-id_offset(jp_I,2) )
-            
-                  il_level1_G0(ji0,jj0)=MAXVAL(il_level1(ji1m:ji1p,jj1m:jj1p))
-               ENDDO
-            ENDDO
-
-            !il_level1_G0(:,:)=0
-            !DO jj1=1,tl_var1%t_dim(2)%i_len
-      
-            !   jj0=INT(REAL((jj1+il_offset(jp_J,1)-1)-1,dp)/REAL(il_rho(jp_J),dp)) +1
-
-            !   DO ji1=1,tl_var1%t_dim(1)%i_len
-
-            !      ji0=INT(REAL((ji1+il_offset(jp_I,1)-1)-1,dp)/REAL(il_rho(jp_I),dp)) +1
-
-            !      il_level1_G0(ji0,jj0)=MAX(il_level1_G0(ji0,jj0),il_level1(ji1,jj1))
-            !      
-            !   ENDDO
-            !ENDDO
-
-            ! clean
-            DEALLOCATE( il_extra )
-            DEALLOCATE( il_level1 )
-
-         ENDIF
-
-         ! look for sea point
-         !il_detect(:,:,1)=0
-         DO jk0=1,td_var0%t_dim(3)%i_len
-            WHERE( il_level1_G0(:,:) >= jk0)
-               il_detect(:,:,jk0)=1
-            END WHERE
-            !il_detect(:,:,jk0)=il_level1_G0(:,:)
-            !WHERE( td_var0%d_value(:,:,jk0,1) /= td_var0%d_fill )
-            !   il_detect(:,:,1)=jk0-1
-            !END WHERE
-         ENDDO
-
-         ! clean
-         DEALLOCATE( il_level1_G0 )
-
       ENDIF
 
@@ -265,19 +348,37 @@
       DEALLOCATE( il_offset )
  
-      !! select extra point depending on interpolation method
-      !! compute point near grid point already inform
-      !FORALL( ji0=1:il_dim0(1), &
-      !&       jj0=1:il_dim0(2), &
-      !&       jk0=1:il_dim0(3), &
-      !&       il_detect(ji0,jj0,jk0) == 1 )
-
-      !   il_detect(MAX(1,ji0-il_ext(jp_I)):MIN(ji0+il_ext(jp_I),il_dim0(1)),&
-      !   &         MAX(1,jj0-il_ext(jp_J)):MIN(jj0+il_ext(jp_J),il_dim0(2)),&
-      !   &         MAX(1,jk0-il_ext(jp_K)):MIN(jk0+il_ext(jp_K),il_dim0(3)) )=1 
-
-      !END FORALL
+      ALLOCATE( il_tmp(il_dim0(1),&
+      &                il_dim0(2),&
+      &                il_dim0(3)) )
+      il_tmp(:,:,:)=il_detect(:,:,:)
+      ! select extra point depending on interpolation method
+      ! compute point near grid point already inform
+      DO jk0=1,il_dim0(3)
+         DO jj0=1,il_dim0(2)
+            DO ji0=1,il_dim0(1)
+
+               IF( il_tmp(ji0,jj0,jk0) == 1 )THEN
+                  il_detect( &
+                  &  MAX(1,ji0-il_ext(jp_I)):MIN(ji0+il_ext(jp_I),il_dim0(1)),&
+                  &  MAX(1,jj0-il_ext(jp_J)):MIN(jj0+il_ext(jp_J),il_dim0(2)),&
+                  &  MAX(1,jk0-il_ext(jp_K)):MIN(jk0+il_ext(jp_K),il_dim0(3)) &
+                  &  ) = 1 
+               ENDIF
+
+            ENDDO
+         ENDDO
+      ENDDO
       
+      ! clean
+      DEALLOCATE( il_tmp )
+
       ! do not compute grid point already inform
-      WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=0
+      DO jk0=1,td_var0%t_dim(3)%i_len
+         DO jj0=1,td_var0%t_dim(2)%i_len
+            DO ji0=1,td_var0%t_dim(1)%i_len
+               IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill ) il_detect(ji0,jj0,jk0)=0
+            ENDDO
+         ENDDO
+      ENDDO
 
       ! save result
@@ -288,23 +389,22 @@
       DEALLOCATE( il_ext )
       DEALLOCATE( il_detect )
+      DEALLOCATE( il_rho )
 
    END FUNCTION extrap__detect
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This function detected point to be extrapolated.
+   !> This function sort variable to be extrapolated, depending on number of
+   !> dimentsion, then detected point to be extrapolated.
    !> 
-   !> @details 
+   !> @author J.Paul
+   !> - November, 2013- Initial Version
    !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in] td_var : variable to extrapolate
-   !> @param[in] id_iext : number of points to be extrapolated in i-direction
-   !> @param[in] id_jext : number of points to be extrapolated in j-direction
-   !> @param[in] id_kext : number of points to be extrapolated in k-direction
-   !> @return
-   !-------------------------------------------------------------------
-   !> @code
+   !> @param[in] td_var    coarse grid variable to extrapolate
+   !> @param[in] td_level  fine grid array of level
+   !> @param[in] id_offset array of offset between fine and coarse grid 
+   !> @param[in] id_rho    array of refinment factor 
+   !> @param[in] id_ext    array of number of points to be extrapolated
+   !> @return 3D array of point to be extrapolated
+   !-------------------------------------------------------------------
    FUNCTION extrap__detect_wrapper( td_var, td_level, &
    &                                id_offset, id_rho, id_ext )
@@ -312,6 +412,6 @@
       IMPLICIT NONE
       ! Argument
-      TYPE(TVAR) ,                         INTENT(INOUT) :: td_var
-      TYPE(TVAR) , DIMENSION(:)          , INTENT(IN   ), OPTIONAL :: td_level
+      TYPE(TVAR) ,                 INTENT(IN   ) :: td_var
+      TYPE(TVAR) , DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: td_level
       INTEGER(i4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset
       INTEGER(i4), DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_rho
@@ -335,5 +435,5 @@
       ELSE IF( ALL(td_var%t_dim(1:3)%l_use) )THEN
          
-         ! detect point to be interpolated on I-J-K
+         ! detect point to be extrapolated on I-J-K
          CALL logger_debug(" EXTRAP DETECT: detect point "//&
          &              " for variable "//TRIM(td_var%c_name) )
@@ -346,5 +446,5 @@
       ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
 
-         ! detect point to be interpolated on I-J
+         ! detect point to be extrapolated on I-J
          CALL logger_debug(" EXTRAP DETECT: detect horizontal point "//&
          &              " for variable "//TRIM(td_var%c_name) )
@@ -357,5 +457,5 @@
       ELSE IF( td_var%t_dim(3)%l_use )THEN
          
-         ! detect point to be interpolated on K
+         ! detect point to be extrapolated on K
          CALL logger_debug(" EXTRAP DETECT: detect vertical point "//&
          &              " for variable "//TRIM(td_var%c_name) )
@@ -373,98 +473,32 @@
       
    END FUNCTION extrap__detect_wrapper
-   !> @endcode
-!   !-------------------------------------------------------------------
-!   !> @brief
-!   !> This function detected point to be extrapolated.
-!   !> 
-!   !> @details 
-!   !>
-!   !> @author J.Paul
-!   !> - Nov, 2013- Initial Version
-!   !
-!   !> @param[in] td_var : variable to extrapolate
-!   !> @param[in] id_iext : number of points to be extrapolated in i-direction
-!   !> @param[in] id_jext : number of points to be extrapolated in j-direction
-!   !> @param[in] id_kext : number of points to be extrapolated in k-direction
-!   !> @return 
-!   !
-!   !> @todo
-!   !-------------------------------------------------------------------
-!   !> @code
-!   FUNCTION extrap__detect(td_var, id_iext, id_jext, id_kext)
-!      IMPLICIT NONE
-!      ! Argument
-!      TYPE(TVAR) ,                   INTENT(INOUT) :: td_var
-!      !INTEGER(i4), DIMENSION(:,:,:), INTENT(OUT  ) :: id_detect
-!      INTEGER(i4),                   INTENT(IN   ), OPTIONAL :: id_iext
-!      INTEGER(i4),                   INTENT(IN   ), OPTIONAL :: id_jext
-!      INTEGER(i4),                   INTENT(IN   ), OPTIONAL :: id_kext
-!
-!      ! function
-!      INTEGER(i4), DIMENSION(td_var%t_dim(1)%i_len, &
-!      &                      td_var%t_dim(2)%i_len, &
-!      &                      td_var%t_dim(3)%i_len) :: extrap__detect
-!
-!      ! local variable
-!      INTEGER(i4) :: il_iext
-!      INTEGER(i4) :: il_jext
-!      INTEGER(i4) :: il_kext
-!
-!      INTEGER(i4), DIMENSION(3) :: il_dim
-!
-!      ! loop indices
-!      INTEGER(i4) :: ji
-!      INTEGER(i4) :: jj
-!      INTEGER(i4) :: jk
-!      !----------------------------------------------------------------
-!
-!      ! optional argument
-!      il_iext=im_minext
-!      IF( PRESENT(id_iext) ) il_iext=id_iext
-!      il_jext=im_minext
-!      IF( PRESENT(id_jext) ) il_jext=id_jext
-!      il_kext=im_minext
-!      IF( PRESENT(id_kext) ) il_kext=id_kext
-!
-!      ! init
-!      extrap__detect(:,:,:)=0
-!
-!      il_dim(:)=td_var%t_dim(1:3)%i_len
-!
-!      ! compute point near grid point already inform
-!      FORALL( ji=1:il_dim(1), &
-!      &       jj=1:il_dim(2), &
-!      &       jk=1:il_dim(3), &
-!      &       td_var%d_value(ji,jj,jk,1) /= td_var%d_fill )
-!
-!         extrap__detect(MAX(1,ji-il_iext):MIN(ji+il_iext,il_dim(1)),&
-!         &              MAX(1,jj-il_jext):MIN(jj+il_jext,il_dim(2)),&
-!         &              MAX(1,jk-il_kext):MIN(jk+il_kext,il_dim(3)) )=1 
-!
-!
-!      END FORALL
-!      
-!      ! do not compute grid point already inform
-!      WHERE( td_var%d_value(:,:,:,1) /= td_var%d_fill ) extrap__detect(:,:,:)=0
-!
-!   END FUNCTION extrap__detect
-!   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This subroutine
+   !> This subroutine select method to be used for extrapolation.
+   !> If need be, increase number of points to be extrapolated.
+   !> Finally launch extrap__fill_value.
    !> 
-   !> @details 
+   !> @details
+   !> optionaly, you could specify :<br/>
+   !>  - refinment factor (default 1)
+   !>  - offset between fine and coarse grid (default compute from refinment factor
+   !> as offset=(rho-1)/2) 
+   !>  - number of point to be extrapolated in each direction (default im_minext)
+   !>  - radius of the halo used to compute extrapolation
+   !>  - maximum number of iteration
    !>
    !> @author J.Paul
    !> - Nov, 2013- Initial Version
    !
-   !> @param[inout] td_var : variable structure
-   !> @param[in] id_iext : number of points to be extrapolated in i-direction
-   !> @param[in] id_jext : number of points to be extrapolated in j-direction
-   !> @param[in] id_kext : number of points to be extrapolated in k-direction
-   !> @param[in] id_extend :  radius of the box used to compute extrapolation 
-   !> @param[in] id_maxiter : maximum nuber of iteration
-   !-------------------------------------------------------------------
-   !> @code
+   !> @param[inout] td_var    variable structure
+   !> @param[in] td_level     fine grid array of level
+   !> @param[in] id_offset    array of offset between fine and coarse grid 
+   !> @param[in] id_rho       array of refinment factor 
+   !> @param[in] id_iext      number of points to be extrapolated in i-direction
+   !> @param[in] id_jext      number of points to be extrapolated in j-direction
+   !> @param[in] id_kext      number of points to be extrapolated in k-direction
+   !> @param[in] id_radius    radius of the halo used to compute extrapolation 
+   !> @param[in] id_maxiter   maximum number of iteration
+   !-------------------------------------------------------------------
    SUBROUTINE extrap__fill_value_wrapper( td_var, td_level, &
    &                                      id_offset,        &
@@ -496,5 +530,5 @@
       !----------------------------------------------------------------
       IF( .NOT. ASSOCIATED(td_var%d_value) )THEN
-         CALL logger_error("EXTRAP FILL VALUE: no table of value "//&
+         CALL logger_error("EXTRAP FILL VALUE: no value "//&
          &  "associted to variable "//TRIM(td_var%c_name) )
       ELSE
@@ -542,6 +576,5 @@
          IF( (il_iext /= 0 .AND. td_var%t_dim(1)%l_use) .OR. &
          &   (il_jext /= 0 .AND. td_var%t_dim(2)%l_use) .OR. &
-         &   (il_kext /= 0 .AND. td_var%t_dim(3)%l_use) .OR. &
-         &   PRESENT(td_level) )THEN
+         &   (il_kext /= 0 .AND. td_var%t_dim(3)%l_use) )THEN
 
             ! number of point use to compute box
@@ -577,23 +610,28 @@
 
    END SUBROUTINE extrap__fill_value_wrapper
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This subroutine
+   !> This subroutine compute point to be extrapolated, then extrapolate point.
    !> 
    !> @details 
+   !> optionaly, you could specify :<br/>
+   !>  - refinment factor (default 1)
+   !>  - offset between fine and coarse grid (default compute from refinment factor
+   !> as offset=(rho-1)/2) 
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
+   !> - November, 2013- Initial Version
    !
-   !> @param[inout] td_var : variable structure
-   !> @param[in] cd_method : extrapolation method
-   !> @param[in] id_iext : number of points to be extrapolated in i-direction
-   !> @param[in] id_jext : number of points to be extrapolated in j-direction
-   !> @param[in] id_kext : number of points to be extrapolated in k-direction
-   !> @param[in] id_radius : radius of the halo used to compute extrapolation
-   !> @param[in] id_maxiter : maximum nuber of iteration
-   !-------------------------------------------------------------------
-   !> @code
+   !> @param[inout] td_var    variable structure
+   !> @param[in] cd_method    extrapolation method
+   !> @param[in] id_iext      number of points to be extrapolated in i-direction
+   !> @param[in] id_jext      number of points to be extrapolated in j-direction
+   !> @param[in] id_kext      number of points to be extrapolated in k-direction
+   !> @param[in] id_radius    radius of the halo used to compute extrapolation
+   !> @param[in] id_maxiter   maximum number of iteration
+   !> @param[in] td_level     fine grid array of level
+   !> @param[in] id_offset    array of offset between fine and coarse grid 
+   !> @param[in] id_rho       array of refinment factor
+   !-------------------------------------------------------------------
    SUBROUTINE extrap__fill_value( td_var, cd_method, &
    &                              id_iext, id_jext, id_kext, &
@@ -619,11 +657,8 @@
 
       INTEGER(i4), DIMENSION(:,:,:)  , ALLOCATABLE :: il_detect
-      INTEGER(i4)                                  :: il_radius
-      INTEGER(i4)                                  :: il_iter
 
       TYPE(TATT)                                   :: tl_att
 
       ! loop indices
-      INTEGER(i4) :: jl
       !----------------------------------------------------------------
 
@@ -637,5 +672,4 @@
       &                                 id_rho,           &
       &                                 id_ext=(/id_iext, id_jext, id_kext/) )
-
       !2- add attribute to variable
       cl_extrap=fct_concat(td_var%c_extrap(:))
@@ -643,65 +677,50 @@
       CALL var_move_att(td_var, tl_att)
 
-      CALL logger_warn(" EXTRAP FILL: "//&
+      CALL att_clean(tl_att)
+
+      CALL logger_info(" EXTRAP FILL: "//&
          &              TRIM(fct_str(SUM(il_detect(:,:,:))))//&
          &              " point(s) to extrapolate " )
 
       !3- extrapolate
-      DO jl=1,td_var%t_dim(4)%i_len
-        
-!            td_var%d_value(:,:,:,jl)=il_detect(:,:,:)
-         il_iter=1
-         DO WHILE( ANY(il_detect(:,:,:)==1) )
-            ! change extend value to minimize number of iteration
-            il_radius=id_radius+(il_iter/id_maxiter)
-            
-            CALL logger_debug(" EXTRAP FILL VALUE: "//&
-            &              TRIM(fct_str(SUM(il_detect(:,:,:))))//&
-            &              " points to extrapolate " )
-            
-            CALL extrap__3D(td_var%d_value(:,:,:,jl), td_var%d_fill,    &
-            &               il_detect(:,:,:),                           &
-            &               cd_method, il_radius )
-               
-            il_iter=il_iter+1
-         ENDDO
-
-      ENDDO
- 
-      IF( SUM(il_detect(:,:,:)) /= 0 )THEN
-         CALL logger_warn(" EXTRAP FILL: still "//&
-         &              TRIM(fct_str(SUM(il_detect(:,:,:))))//&
-         &              " point(s) to extrapolate " )
-      ENDIF
+      CALL extrap__3D(td_var%d_value(:,:,:,:), td_var%d_fill,    &
+      &               il_detect(:,:,:),                           &
+      &               cd_method, id_radius, id_maxiter  )
 
       DEALLOCATE(il_detect)
 
    END SUBROUTINE extrap__fill_value
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This subroutine
+   !> This subroutine compute point to be extrapolated in 3D array.
    !> 
    !> @details 
+   !> in case of 'min_error' method:<br/>
+   !>    - compute derivative in i-, j- and k- direction
+   !>    - compute minimum error coefficient (distance to center of halo)
+   !>    - compute extrapolatd values by calculated minimum error using taylor expansion 
+   !> in case of 'dist_weight' method:<br/>
+   !>    - compute distance weight coefficient (inverse of distance to center of halo) 
+   !>    - compute extrapolatd values using Inverse Distance Weighting
    !>
    !> @author J.Paul
    !> - Nov, 2013- Initial Version
    !
-   !> @param[inout] dd_value : 3D table of variable to be extrapolated
-   !> @param[in] dd_fill : FillValue of variable
-   !> @param[inout] id_detect : table of point to extrapolate
-   !> @param[in] id_ext : number of point use to compute box 
-   !> @param[in] cd_method : extrapolation method
-   !-------------------------------------------------------------------
-   !> @code
+   !> @param[inout] dd_value  3D array of variable to be extrapolated
+   !> @param[in] dd_fill      FillValue of variable
+   !> @param[inout] id_detect array of point to extrapolate
+   !> @param[in] cd_method    extrapolation method
+   !> @param[in] id_radius    radius of the halo used to compute extrapolation
+   !-------------------------------------------------------------------
    SUBROUTINE extrap__3D( dd_value, dd_fill, id_detect,&
-   &                      cd_method, id_ext )
+   &                      cd_method, id_radius, id_maxiter )
       IMPLICIT NONE
       ! Argument
-      REAL(dp)   , DIMENSION(:,:,:), INTENT(INOUT) :: dd_value
+      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_ext
+      INTEGER(i4),                   INTENT(IN   ) :: id_radius
+      INTEGER(i4),                   INTENT(IN   ) :: id_maxiter
 
       ! local variable
@@ -712,7 +731,11 @@
       INTEGER(i4) :: il_kmin
       INTEGER(i4) :: il_kmax
-
-      INTEGER(i4), DIMENSION(3) :: il_shape
+      INTEGER(i4) :: il_iter
+      INTEGER(i4) :: il_radius
+
+      INTEGER(i4), DIMENSION(4) :: il_shape
       INTEGER(i4), DIMENSION(3) :: il_dim
+
+      INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect
 
       REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx
@@ -725,187 +748,217 @@
       INTEGER(i4) :: jj
       INTEGER(i4) :: jk
+      INTEGER(i4) :: jl
       !----------------------------------------------------------------
 
       il_shape(:)=SHAPE(dd_value)
 
+      ALLOCATE( il_detect( il_shape(1), il_shape(2), il_shape(3)) )
+
       SELECT CASE(TRIM(cd_method))
-         CASE('min_error')
-
-            ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) ) 
-            ALLOCATE( dl_dfdy(il_shape(1), il_shape(2), il_shape(3)) ) 
-            ALLOCATE( dl_dfdz(il_shape(1), il_shape(2), il_shape(3)) ) 
-
-
-            ! compute derivative in i-direction
-            dl_dfdx(:,:,:)=dd_fill
-            IF( il_shape(1) > 1 )THEN
-               dl_dfdx(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'I' )
-            ENDIF
-
-            ! compute derivative in i-direction
-            dl_dfdy(:,:,:)=dd_fill
-            IF( il_shape(2) > 1 )THEN
-               dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'J' )
-            ENDIF
-               
-            ! compute derivative in i-direction
-            dl_dfdz(:,:,:)=dd_fill
-            IF( il_shape(3) > 1 )THEN
-               dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'K' )
-            ENDIF
-      
-            il_dim(1)=2*id_ext+1
-            IF( il_shape(1) < 2*id_ext+1 ) il_dim(1)=1
-            il_dim(2)=2*id_ext+1
-            IF( il_shape(2) < 2*id_ext+1 ) il_dim(2)=1
-            il_dim(3)=2*id_ext+1
-            IF( il_shape(3) < 2*id_ext+1 ) il_dim(3)=1
-
-            ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 
-
-            dl_coef(:,:,:)=extrap__3D_min_error_coef(dd_value( 1:il_dim(1), &
-            &                                                  1:il_dim(2), &
-            &                                                  1:il_dim(3)))
-
-            DO jk=1,il_shape(3)
-               IF( ALL(id_detect(:,:,jk) == 0) ) CYCLE
-               DO jj=1,il_shape(2)
-                  IF( ALL(id_detect(:,jj,jk) == 0) ) CYCLE
-                  DO ji=1,il_shape(1)
-
-                     IF( id_detect(ji,jj,jk) == 1 )THEN
-                       
-                        il_imin=MAX(ji-id_ext,1)
-                        il_imax=MIN(ji+id_ext,il_shape(1))
-                        IF( il_dim(1) == 1 )THEN
-                           il_imin=ji
-                           il_imax=ji
+      CASE('min_error')
+         DO jl=1,il_shape(4)
+
+            ! intitialise table of poitn to be extrapolated
+            il_detect(:,:,:)=id_detect(:,:,:)
+
+            il_iter=1
+            DO WHILE( ANY(il_detect(:,:,:)==1) )
+               ! change extend value to minimize number of iteration
+               il_radius=id_radius+(il_iter/id_maxiter)
+
+               ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) ) 
+               ALLOCATE( dl_dfdy(il_shape(1), il_shape(2), il_shape(3)) ) 
+               ALLOCATE( dl_dfdz(il_shape(1), il_shape(2), il_shape(3)) ) 
+
+               ! compute derivative in i-direction
+               dl_dfdx(:,:,:)=dd_fill
+               IF( il_shape(1) > 1 )THEN
+                  dl_dfdx(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'I' )
+               ENDIF
+
+               ! compute derivative in j-direction
+               dl_dfdy(:,:,:)=dd_fill
+               IF( il_shape(2) > 1 )THEN
+                  dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'J' )
+               ENDIF
+ 
+               ! compute derivative in k-direction
+               dl_dfdz(:,:,:)=dd_fill
+               IF( il_shape(3) > 1 )THEN
+                  dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'K' )
+               ENDIF
+ 
+               il_dim(1)=2*il_radius+1
+               IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1
+               il_dim(2)=2*il_radius+1
+               IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1
+               il_dim(3)=2*il_radius+1
+               IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1
+
+               ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 
+
+               dl_coef(:,:,:)=extrap__3D_min_error_coef(dd_value( 1:il_dim(1), &
+               &                                                  1:il_dim(2), &
+               &                                                  1:il_dim(3), &
+               &                                                  jl ))
+
+               DO jk=1,il_shape(3)
+                  IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE
+                  DO jj=1,il_shape(2)
+                     IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE
+                     DO ji=1,il_shape(1)
+
+                        IF( il_detect(ji,jj,jk) == 1 )THEN
+                          
+                           il_imin=MAX(ji-il_radius,1)
+                           il_imax=MIN(ji+il_radius,il_shape(1))
+                           IF( il_dim(1) == 1 )THEN
+                              il_imin=ji
+                              il_imax=ji
+                           ENDIF
+
+                           il_jmin=MAX(jj-il_radius,1)
+                           il_jmax=MIN(jj+il_radius,il_shape(2))
+                           IF( il_dim(2) == 1 )THEN
+                              il_jmin=jj
+                              il_jmax=jj
+                           ENDIF
+
+                           il_kmin=MAX(jk-il_radius,1)
+                           il_kmax=MIN(jk+il_radius,il_shape(3))
+                           IF( il_dim(3) == 1 )THEN
+                              il_kmin=jk
+                              il_kmax=jk
+                           ENDIF
+
+                           dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( &
+                           &  dd_value( il_imin:il_imax, &
+                           &            il_jmin:il_jmax, &
+                           &            il_kmin:il_kmax,jl ), dd_fill, il_radius, &
+                           &  dl_dfdx(  il_imin:il_imax, &
+                           &            il_jmin:il_jmax, &
+                           &            il_kmin:il_kmax ), &
+                           &  dl_dfdy(  il_imin:il_imax, &
+                           &            il_jmin:il_jmax, &
+                           &            il_kmin:il_kmax ), &
+                           &  dl_dfdz(  il_imin:il_imax, &
+                           &            il_jmin:il_jmax, &
+                           &            il_kmin:il_kmax ), &
+                           &  dl_coef(:,:,:) )
+
+                           IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN
+                              il_detect(ji,jj,jk)= 0
+                           ENDIF
+
                         ENDIF
 
-                        il_jmin=MAX(jj-id_ext,1)
-                        il_jmax=MIN(jj+id_ext,il_shape(2))
-                        IF( il_dim(2) == 1 )THEN
-                           il_jmin=jj
-                           il_jmax=jj
-                        ENDIF
-
-                        il_kmin=MAX(jk-id_ext,1)
-                        il_kmax=MIN(jk+id_ext,il_shape(3))
-                        IF( il_dim(3) == 1 )THEN
-                           il_kmin=jk
-                           il_kmax=jk
-                        ENDIF
-
-                        dd_value(ji,jj,jk)=extrap__3D_min_error_fill( &
-                        &  dd_value( il_imin:il_imax, &
-                        &            il_jmin:il_jmax, &
-                        &            il_kmin:il_kmax ), dd_fill, id_ext, &
-                        &  dl_dfdx(  il_imin:il_imax, &
-                        &            il_jmin:il_jmax, &
-                        &            il_kmin:il_kmax ), &
-                        &  dl_dfdy(  il_imin:il_imax, &
-                        &            il_jmin:il_jmax, &
-                        &            il_kmin:il_kmax ), &
-                        &  dl_dfdz(  il_imin:il_imax, &
-                        &            il_jmin:il_jmax, &
-                        &            il_kmin:il_kmax ), &
-                        &  dl_coef(:,:,:) )
-
-                        IF( dd_value(ji,jj,jk) /= dd_fill )THEN
-                           id_detect(ji,jj,jk)= 0
-                        ENDIF
-
-                     ENDIF
-
+                     ENDDO
                   ENDDO
                ENDDO
+
+               DEALLOCATE( dl_dfdx )
+               DEALLOCATE( dl_dfdy )
+               DEALLOCATE( dl_dfdz )
+               DEALLOCATE( dl_coef )
+
+               il_iter=il_iter+1
             ENDDO
-
-            IF( ALLOCATED(dl_dfdx) ) DEALLOCATE( dl_dfdx )
-            IF( ALLOCATED(dl_dfdy) ) DEALLOCATE( dl_dfdy )
-            IF( ALLOCATED(dl_dfdz) ) DEALLOCATE( dl_dfdz )
-            IF( ALLOCATED(dl_coef) ) DEALLOCATE( dl_coef )
-
-         CASE DEFAULT ! 'dist_weight'
-
-            il_dim(1)=2*id_ext+1
-            IF( il_shape(1) < 2*id_ext+1 ) il_dim(1)=1
-            il_dim(2)=2*id_ext+1
-            IF( il_shape(2) < 2*id_ext+1 ) il_dim(2)=1
-            il_dim(3)=2*id_ext+1
-            IF( il_shape(3) < 2*id_ext+1 ) il_dim(3)=1
-            
-            ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) )
-
-            dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), &
-            &                                                   1:il_dim(2), &
-            &                                                   1:il_dim(3)) )
-
-            DO jk=1,il_shape(3)
-               IF( ALL(id_detect(:,:,jk) == 0) ) CYCLE
-               DO jj=1,il_shape(2)
-                  IF( ALL(id_detect(:,jj,jk) == 0) ) CYCLE
-                  DO ji=1,il_shape(1)
-
-                     IF( id_detect(ji,jj,jk) == 1 )THEN
-                        
-                        il_imin=MAX(ji-id_ext,1)
-                        il_imax=MIN(ji+id_ext,il_shape(1))
-                        IF( il_dim(1) == 1 )THEN
-                           il_imin=ji
-                           il_imax=ji
+         ENDDO
+
+      CASE DEFAULT ! 'dist_weight'
+         DO jl=1,il_shape(4)
+
+            ! intitialise table of poitn to be extrapolated
+            il_detect(:,:,:)=id_detect(:,:,:)
+
+            il_iter=1
+            DO WHILE( ANY(il_detect(:,:,:)==1) )
+               ! change extend value to minimize number of iteration
+               il_radius=id_radius+(il_iter/id_maxiter)
+
+               il_dim(1)=2*il_radius+1
+               IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1
+               il_dim(2)=2*il_radius+1
+               IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1
+               il_dim(3)=2*il_radius+1
+               IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1
+               
+               ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) )
+
+               dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), &
+               &                                                   1:il_dim(2), &
+               &                                                   1:il_dim(3), &
+               &                                                   jl ) )
+
+               DO jk=1,il_shape(3)
+                  IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE
+                  DO jj=1,il_shape(2)
+                     IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE
+                     DO ji=1,il_shape(1)
+
+                        IF( il_detect(ji,jj,jk) == 1 )THEN
+                           
+                           il_imin=MAX(ji-il_radius,1)
+                           il_imax=MIN(ji+il_radius,il_shape(1))
+                           IF( il_dim(1) == 1 )THEN
+                              il_imin=ji
+                              il_imax=ji
+                           ENDIF
+
+                           il_jmin=MAX(jj-il_radius,1)
+                           il_jmax=MIN(jj+il_radius,il_shape(2))
+                           IF( il_dim(2) == 1 )THEN
+                              il_jmin=jj
+                              il_jmax=jj
+                           ENDIF
+
+                           il_kmin=MAX(jk-il_radius,1)
+                           il_kmax=MIN(jk+il_radius,il_shape(3))
+                           IF( il_dim(3) == 1 )THEN
+                              il_kmin=jk
+                              il_kmax=jk
+                           ENDIF
+
+                           dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( &
+                           &  dd_value( il_imin:il_imax, &
+                           &            il_jmin:il_jmax, &
+                           &            il_kmin:il_kmax, &
+                           &            jl), dd_fill, il_radius, &
+                           &  dl_coef(:,:,:) )
+
+                           IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN
+                              il_detect(ji,jj,jk)= 0
+                           ENDIF
+
                         ENDIF
 
-                        il_jmin=MAX(jj-id_ext,1)
-                        il_jmax=MIN(jj+id_ext,il_shape(2))
-                        IF( il_dim(2) == 1 )THEN
-                           il_jmin=jj
-                           il_jmax=jj
-                        ENDIF
-
-                        il_kmin=MAX(jk-id_ext,1)
-                        il_kmax=MIN(jk+id_ext,il_shape(3))
-                        IF( il_dim(3) == 1 )THEN
-                           il_kmin=jk
-                           il_kmax=jk
-                        ENDIF
-
-                        dd_value(ji,jj,jk)=extrap__3D_dist_weight_fill( &
-                        &  dd_value( il_imin:il_imax, &
-                        &            il_jmin:il_jmax, &
-                        &            il_kmin:il_kmax ), dd_fill, id_ext, &
-                        &  dl_coef(:,:,:) )
-
-                        IF( dd_value(ji,jj,jk) /= dd_fill )THEN
-                           id_detect(ji,jj,jk)= 0
-                        ENDIF
-
-                     ENDIF
-
+                     ENDDO
                   ENDDO
                ENDDO
+
+               DEALLOCATE( dl_coef )
+               il_iter=il_iter+1
             ENDDO
-
-            IF( ALLOCATED(dl_coef) ) DEALLOCATE( dl_coef )
-            
+         ENDDO            
       END SELECT
 
+      DEALLOCATE( il_detect )
+
    END SUBROUTINE extrap__3D
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This function compute derivative of a function in i- and
-   !> j-direction
+   !> This function compute derivative of 1D array.
    !> 
    !> @details 
+   !> optionaly you could specify to take into account east west discontinuity
+   !> (-180° 180° or 0° 360° for longitude variable)
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
+   !> - November, 2013- Initial Version
    !
-   !> @param[in] dd_value : 1D table of variable to be extrapolated
-   !> @param[in] dd_fill : FillValue of variable
-   !> @param[in] cd_dim : compute derivative on first (I) or second (J) dimension 
-   !-------------------------------------------------------------------
-   !> @code
+   !> @param[in] dd_value     1D array of variable to be extrapolated
+   !> @param[in] dd_fill      FillValue of variable
+   !> @param[in] ld_discont   logical to take into account east west discontinuity 
+   !-------------------------------------------------------------------
    PURE FUNCTION extrap_deriv_1D( dd_value, dd_fill, ld_discont )
 
@@ -1003,20 +1056,22 @@
 
    END FUNCTION extrap_deriv_1D
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This function compute derivative of a function in i- and
-   !> j-direction
-   !> 
+   !> This function compute derivative of 2D array.
+   !> you have to specify in which direction derivative have to be computed:
+   !> first (I) or second (J) dimension. 
+   !>
    !> @details 
+   !> optionaly you could specify to take into account east west discontinuity
+   !> (-180° 180° or 0° 360° for longitude variable)
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
+   !> - November, 2013- Initial Version
    !
-   !> @param[in] dd_value : 2D table of variable to be extrapolated
-   !> @param[in] dd_fill : FillValue of variable
-   !> @param[in] cd_dim : compute derivative on first (I) or second (J) dimension 
-   !-------------------------------------------------------------------
-   !> @code
+   !> @param[in] dd_value     2D array of variable to be extrapolated
+   !> @param[in] dd_fill      FillValue of variable
+   !> @param[in] cd_dim       compute derivative on first (I) or second (J) dimension 
+   !> @param[in] ld_discont   logical to take into account east west discontinuity 
+   !-------------------------------------------------------------------
    FUNCTION extrap_deriv_2D( dd_value, dd_fill, cd_dim, ld_discont )
 
@@ -1123,5 +1178,5 @@
             END WHERE
 
-      ENDDO
+         ENDDO
 
       CASE('J')
@@ -1187,20 +1242,22 @@
 
    END FUNCTION extrap_deriv_2D
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This function compute derivative of a function in i- and
-   !> j-direction
+   !> This function compute derivative of 3D array.
+   !> you have to specify in which direction derivative have to be computed:
+   !> first (I), second (J) or third (K) dimension.
    !> 
    !> @details 
+   !> optionaly you could specify to take into account east west discontinuity
+   !> (-180° 180° or 0° 360° for longitude variable)
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
+   !> - November, 2013- Initial Version
    !
-   !> @param[inout] dd_value : 3D table of variable to be extrapolated
-   !> @param[in] dd_fill : FillValue of variable
-   !> @param[in] cd_dim : compute derivative on first (I) second (J) or third (K) dimension   
-   !-------------------------------------------------------------------
-   !> @code
+   !> @param[inout] dd_value  3D array of variable to be extrapolated
+   !> @param[in] dd_fill      FillValue of variable
+   !> @param[in] cd_dim       compute derivative on first (I) second (J) or third (K) dimension   
+   !> @param[in] ld_discont   logical to take into account east west discontinuity
+   !-------------------------------------------------------------------
    PURE FUNCTION extrap_deriv_3D( dd_value, dd_fill, cd_dim, ld_discont )
 
@@ -1431,26 +1488,18 @@
 
    END FUNCTION extrap_deriv_3D
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This function compute extrapolatd values by calculated minimum error using
-   !> taylor expansion
+   !> This function compute coefficient for min_error extrapolation.
    !> 
    !> @details 
+   !> coefficients are  "grid distance" to the center of the box choosed to compute
+   !> extrapolation.
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
+   !> - November, 2013- Initial Version
    !
-   !> @param[in] dd_value : 3D table of variable to be extrapolated
-   !> @param[in] dd_fill : FillValue of variable
-   !> @param[in] dd_ideriv : derivative of function in i-direction
-   !> @param[in] dd_jderiv : derivative of function in j-direction
-   !> @param[in] dd_kderiv : derivative of function in k-direction
-   !> @param[in] id_ji  : i-direction indices table
-   !> @param[in] id_jj  : j-direction indices table 
-   !> @param[in] id_ii  : i-direction indices of the point to extrapolate 
-   !> @param[in] id_ij  : j-direction indices of the point to extrapolate
-   !-------------------------------------------------------------------
-   !> @code
+   !> @param[in] dd_value  3D array of variable to be extrapolated
+   !> @return 3D array of coefficient for minimum error extrapolation
+   !-------------------------------------------------------------------
    PURE FUNCTION extrap__3D_min_error_coef( dd_value )
 
@@ -1514,24 +1563,22 @@
 
    END FUNCTION extrap__3D_min_error_coef
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This function compute extrapolatd values by calculated minimum error using
+   !> This function compute extrapolatd value by calculated minimum error using
    !> taylor expansion
    !> 
-   !> @details 
+   !> @author J.Paul
+   !> - November, 2013- Initial Version
    !>
-   !> @author J.Paul
-   !> - Nov, 2013- Initial Version
-   !
-   !> @param[in] dd_value : 3D table of variable to be extrapolated
-   !> @param[in] dd_fill : FillValue of variable
-   !> @param[in] dd_dfdx : derivative of function in i-direction
-   !> @param[in] dd_dfdy : derivative of function in j-direction
-   !> @param[in] dd_dfdz : derivative of function in k-direction
-   !> @param[in] dd_coef : 
-   !-------------------------------------------------------------------
-   !> @code
-   PURE FUNCTION extrap__3D_min_error_fill( dd_value, dd_fill, id_ext, &
+   !> @param[in] dd_value  3D array of variable to be extrapolated
+   !> @param[in] dd_fill   FillValue of variable
+   !> @param[in] id_radius radius of the halo used to compute extrapolation
+   !> @param[in] dd_dfdx   derivative of function in i-direction
+   !> @param[in] dd_dfdy   derivative of function in j-direction
+   !> @param[in] dd_dfdz   derivative of function in k-direction
+   !> @param[in] dd_coef   array of coefficient for min_error extrapolation 
+   !> @return extrapolatd value
+   !-------------------------------------------------------------------
+   PURE FUNCTION extrap__3D_min_error_fill( dd_value, dd_fill, id_radius, &
    &                                        dd_dfdx, dd_dfdy, dd_dfdz, &
    &                                        dd_coef )
@@ -1540,5 +1587,5 @@
       REAL(dp)   , DIMENSION(:,:,:), INTENT(IN) :: dd_value
       REAL(dp)   ,                   INTENT(IN) :: dd_fill
-      INTEGER(i4),                   INTENT(IN) :: id_ext
+      INTEGER(i4),                   INTENT(IN) :: id_radius
       REAL(dp)   , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdx
       REAL(dp)   , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdy
@@ -1564,5 +1611,5 @@
       extrap__3D_min_error_fill=dd_fill
 
-      il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_ext*2))
+      il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2))
 
       IF( COUNT(dd_value(:,:,:) /= dd_fill) >= il_min )THEN
@@ -1602,18 +1649,18 @@
 
    END FUNCTION extrap__3D_min_error_fill
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This function compute extrapolatd values by calculated minimum error using
-   !> taylor expansion
+   !> This function compute coefficient for inverse distance weighted method 
    !> 
    !> @details 
+   !> coefficients are inverse "grid distance" to the center of the box choosed to compute
+   !> extrapolation.
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
+   !> - November, 2013- Initial Version
    !
-   !> @param[in] dd_value : 3D table of variable to be extrapolated
-   !-------------------------------------------------------------------
-   !> @code
+   !> @param[in] dd_value  3D array of variable to be extrapolated
+   !> @return 3D array of coefficient for inverse distance weighted extrapolation
+   !-------------------------------------------------------------------
    PURE FUNCTION extrap__3D_dist_weight_coef( dd_value )
 
@@ -1677,27 +1724,27 @@
 
    END FUNCTION extrap__3D_dist_weight_coef
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
-   !> This function compute extrapolatd values by calculated minimum error using
-   !> taylor expansion
+   !> This function compute extrapolatd value using inverse distance weighted
+   !> method
    !> 
    !> @details 
    !>
    !> @author J.Paul
-   !> - Nov, 2013- Initial Version
+   !> - November, 2013- Initial Version
    !
-   !> @param[in] dd_value : 3D table of variable to be extrapolated
-   !> @param[in] dd_fill : FillValue of variable
-   !> @param[in] dd_coef :  
-   !-------------------------------------------------------------------
-   !> @code
-   FUNCTION extrap__3D_dist_weight_fill( dd_value, dd_fill, id_ext, &
-   &                                          dd_coef )
+   !> @param[in] dd_value  3D array of variable to be extrapolated
+   !> @param[in] dd_fill   FillValue of variable
+   !> @param[in] id_radius radius of the halo used to compute extrapolation
+   !> @param[in] dd_coef   3D array of coefficient for inverse distance weighted extrapolation 
+   !> @return extrapolatd value
+   !-------------------------------------------------------------------
+   FUNCTION extrap__3D_dist_weight_fill( dd_value, dd_fill, id_radius, &
+   &                                     dd_coef )
       IMPLICIT NONE
       ! Argument
       REAL(dp)   , DIMENSION(:,:,:), INTENT(IN) :: dd_value
       REAL(dp)   ,                   INTENT(IN) :: dd_fill
-      INTEGER(i4),                   INTENT(IN) :: id_ext
+      INTEGER(i4),                   INTENT(IN) :: id_radius
       REAL(dp)   , DIMENSION(:,:,:), INTENT(IN) :: dd_coef
 
@@ -1722,5 +1769,5 @@
       extrap__3D_dist_weight_fill=dd_fill
 
-      il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_ext*2))
+      il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2))
 
       IF( COUNT(dd_value(:,:,:)/= dd_fill) >= il_min )THEN
@@ -1733,18 +1780,22 @@
          dl_coef(:,:,:)=0
 
-         FORALL( ji=1:il_shape(1), &
-         &       jj=1:il_shape(2), &
-         &       jk=1:il_shape(3), &
-         &       dd_value(ji,jj,jk) /= dd_fill )
-
-            ! compute factor
-            dl_value(ji,jj,jk)=dd_coef(ji,jj,jk)*dd_value(ji,jj,jk)
-            dl_coef(ji,jj,jk)=dd_coef(ji,jj,jk)
-
-         END FORALL
+         DO jk=1,il_shape(3)
+            DO jj=1,il_shape(2)
+               DO ji=1,il_shape(1)
+
+                  IF( dd_value(ji,jj,jk) /= dd_fill )THEN
+                     ! compute factor
+                     dl_value(ji,jj,jk)=dd_coef(ji,jj,jk)*dd_value(ji,jj,jk)
+                     dl_coef(ji,jj,jk)=dd_coef(ji,jj,jk)
+                  ENDIF
+
+               ENDDO
+            ENDDO
+         ENDDO
 
          ! return value
          IF( SUM( dl_coef(:,:,:) ) /= 0 )THEN
-            extrap__3D_dist_weight_fill=SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) )
+            extrap__3D_dist_weight_fill = &
+            &  SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) )
          ENDIF
 
@@ -1755,5 +1806,4 @@
 
    END FUNCTION extrap__3D_dist_weight_fill
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
@@ -1761,14 +1811,16 @@
    !> extraband of N points at north,south,east and west boundaries.
    !> 
+   !> @details
+   !> optionaly you could specify size of extra bands in i- and j-direction
+   !>
    !> @author J.Paul
-   !> - 2013-Initial version
+   !> - November, 2013-Initial version
    !
-   !> @param[inout] td_var : variable 
-   !> @param[in] id_isize : i-direction size of extra bands (default=im_minext)
-   !> @param[in] id_jsize : j-direction size of extra bands (default=im_minext)
+   !> @param[inout] td_var variable 
+   !> @param[in] id_isize  i-direction size of extra bands (default=im_minext)
+   !> @param[in] id_jsize  j-direction size of extra bands (default=im_minext)
    !> @todo
    !> - invalid special case for grid with north fold
    !-------------------------------------------------------------------
-   !> @code
    SUBROUTINE extrap_add_extrabands(td_var, id_isize, id_jsize )
       IMPLICIT NONE
@@ -1821,4 +1873,5 @@
       dl_value(:,:,:,:)=td_var%d_value(:,:,:,:)
 
+
       td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len + 2*il_isize
       td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len + 2*il_jsize
@@ -1855,5 +1908,4 @@
 
    END SUBROUTINE extrap_add_extrabands
-   !> @endcode
    !-------------------------------------------------------------------
    !> @brief
@@ -1861,16 +1913,18 @@
    !> of N points at north,south,east and west boundaries.
    !> 
+   !> @details
+   !> optionaly you could specify size of extra bands in i- and j-direction
+   !>
    !> @author J.Paul
-   !> - 2013-Initial version
-   !
-   !> @param[inout] td_var : variable 
-   !> @param[in] id_isize : i-direction size of extra bands (default=im_minext)
-   !> @param[in] id_jsize : j-direction size of extra bands (default=im_minext)
-   !-------------------------------------------------------------------
-   !> @code
+   !> - November, 2013-Initial version
+   !>
+   !> @param[inout] td_var variable 
+   !> @param[in] id_isize  i-direction size of extra bands (default=im_minext)
+   !> @param[in] id_jsize  j-direction size of extra bands (default=im_minext)
+   !-------------------------------------------------------------------
    SUBROUTINE extrap_del_extrabands(td_var, id_isize, id_jsize )
       IMPLICIT NONE
       ! Argument
-      TYPE(TVAR) , INTENT(INOUT)  :: td_var
+      TYPE(TVAR) , INTENT(INOUT) :: td_var
       INTEGER(i4), INTENT(IN   ), OPTIONAL :: id_isize
       INTEGER(i4), INTENT(IN   ), OPTIONAL :: id_jsize
@@ -1881,5 +1935,5 @@
       INTEGER(i4) :: il_isize
       INTEGER(i4) :: il_jsize
-      
+ 
       INTEGER(i4) :: il_imin
       INTEGER(i4) :: il_imax
@@ -1935,50 +1989,3 @@
 
    END SUBROUTINE extrap_del_extrabands
-   !> @endcode
-!   !-------------------------------------------------------------------
-!   !> @brief
-!   !> This function
-!   !> 
-!   !> @details 
-!   !>
-!   !> @author J.Paul
-!   !> - Nov, 2013- Initial Version
-!   !
-!   !> @param[in]
-!   !> @param[out] 
-!   !-------------------------------------------------------------------
-!   !> @code
-!   FUNCTION extrap_( )
-!      IMPLICIT NONE
-!      ! Argument
-!
-!      ! local variable
-!
-!      ! loop indices
-!      !----------------------------------------------------------------
-!   END FUNCTION extrap_
-!   !> @endcode
-!   !-------------------------------------------------------------------
-!   !> @brief
-!   !> This subroutine
-!   !> 
-!   !> @details 
-!   !>
-!   !> @author J.Paul
-!   !> - Nov, 2013- Initial Version
-!   !
-!   !> @param[in]
-!   !> @param[out] 
-!   !-------------------------------------------------------------------
-!   !> @code
-!   SUBROUTINE extrap_( )
-!      IMPLICIT NONE
-!      ! Argument
-!
-!      ! local variable
-!
-!      ! loop indices
-!      !----------------------------------------------------------------
-!   END SUBROUTINE extrap_
-!   !> @endcode
 END MODULE extrap
