New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5240 for branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/TOOLS/SIREN/src/extrap.f90 – NEMO

Ignore:
Timestamp:
2015-04-29T12:17:12+02:00 (9 years ago)
Author:
davestorkey
Message:

Update UKMO nn_etau_revision branch with trunk changes to rev 5107.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/TOOLS/SIREN/src/extrap.f90

    r4213 r5240  
    77! DESCRIPTION: 
    88!> @brief  
    9 !> This module  
     9!> This module manage extrapolation. 
    1010!> 
    1111!> @details 
     12!>    Extrapolation method to be used is specify inside variable  
     13!>    strcuture, as array of string character.<br/>  
     14!>    - td_var\%c_extrap(1) string character is the interpolation name choose between: 
     15!>       - 'dist_weight' 
     16!>       - 'min_error' 
     17!> 
     18!>    @note Extrapolation method could be specify for each variable in namelist _namvar_, 
     19!>    defining string character _cn\_varinfo_. By default _dist_weight_.<br/> 
     20!>    Example: 
     21!>       - cn_varinfo='varname1:dist_weight', 'varname2:min_error' 
     22!> 
     23!>    to detect point to be extrapolated:<br/> 
     24!> @code 
     25!>    il_detect(:,:,:)=extrap_detect(td_var, [td_level], [id_offset,] [id_rho,] [id_ext])  
     26!> @endcode 
     27!>       - il_detect(:,:,:) is 3D array of point to be extrapolated 
     28!>       - td_var  is coarse grid variable to be extrapolated 
     29!>       - td_level is fine grid array of level (see vgrid_get_level) [optional] 
     30!>       - id_offset is array of offset between fine and coarse grid [optional] 
     31!>       - id_rho    is array of refinment factor [optional] 
     32!>       - id_ext    is array of number of points to be extrapolated [optional] 
     33!> 
     34!>    to extrapolate variable:<br/> 
     35!> @code 
     36!>    CALL extrap_fill_value( td_var, [td_level], [id_offset], [id_rho], [id_iext], [id_jext], [id_kext], [id_radius], [id_maxiter]) 
     37!> @endcode 
     38!>       - td_var  is coarse grid variable to be extrapolated 
     39!>       - td_level is fine grid array of level (see vgrid_get_level) [optional] 
     40!>       - id_offset is array of offset between fine and coarse grid [optional] 
     41!>       - id_rho    is array of refinment factor [optional] 
     42!>       - id_iext   is number of points to be extrapolated in i-direction [optional] 
     43!>       - id_jext   is number of points to be extrapolated in j-direction [optional] 
     44!>       - id_kext   is number of points to be extrapolated in k-direction [optional] 
     45!>       - id_radius is radius of the halo used to compute extrapolation [optional] 
     46!>       - id_maxiter is maximum number of iteration [optional] 
     47!> 
     48!>    to add extraband to the variable (to be extrapolated):<br/> 
     49!> @code 
     50!>    CALL extrap_add_extrabands(td_var, [id_isize,] [id_jsize] ) 
     51!> @endcode 
     52!>       - td_var is variable structure 
     53!>       - id_isize : i-direction size of extra bands [optional] 
     54!>       - id_jsize : j-direction size of extra bands [optional] 
     55!> 
     56!>    to delete extraband of a variable:<br/> 
     57!> @code 
     58!>    CALL extrap_del_extrabands(td_var, [id_isize,] [id_jsize] ) 
     59!> @endcode 
     60!>       - td_var is variable structure 
     61!>       - id_isize : i-direction size of extra bands [optional] 
     62!>       - id_jsize : j-direction size of extra bands [optional] 
     63!> 
     64!>    to compute first derivative of 1D array:<br/> 
     65!> @code 
     66!>    dl_value(:)=extrap_deriv_1D( dd_value(:), dd_fill, [ld_discont] ) 
     67!> @endcode 
     68!>       - dd_value is 1D array of variable 
     69!>       - dd_fill is FillValue of variable 
     70!>       - ld_discont is logical to take into account longitudinal East-West discontinuity [optional] 
     71!> 
     72!>    to compute first derivative of 2D array:<br/> 
     73!> @code 
     74!>    dl_value(:,:)=extrap_deriv_2D( dd_value(:,:), dd_fill, cd_dim, [ld_discont] ) 
     75!> @endcode 
     76!>       - dd_value is 2D array of variable 
     77!>       - dd_fill is FillValue of variable 
     78!>       - cd_dim is character to compute derivative on first (I) or second (J) dimension 
     79!>       - ld_discont is logical to take into account longitudinal East-West discontinuity [optional] 
     80!> 
     81!>    to compute first derivative of 3D array:<br/> 
     82!> @code 
     83!>    dl_value(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, cd_dim, [ld_discont] ) 
     84!> @endcode 
     85!>       - dd_value is 3D array of variable 
     86!>       - dd_fill is FillValue of variable 
     87!>       - cd_dim is character to compute derivative on first (I), second (J), or third (K) dimension 
     88!>       - ld_discont is logical to take into account longitudinal East-West discontinuity [optional] 
     89!> 
     90!> @warning _FillValue must not be zero (use var_chg_FillValue()) 
    1291!> 
    1392!> @author 
     
    1594! REVISION HISTORY: 
    1695!> @date Nov, 2013 - Initial Version 
     96!> @date September, 2014 
     97!> - add header 
    1798!> 
    18 !> @note WARNING: FillValue must not be zero (use var_chg_FillValue) 
     99!> @todo 
     100!> - create module for each extrapolation method 
    19101!> 
    20102!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    21 !> @todo 
    22 !> - something wrong when computing point to be extralopated 
    23 !> - take care of ew value in variable structure 
    24103!---------------------------------------------------------------------- 
    25104MODULE extrap 
    26105   USE netcdf                          ! nf90 library 
    27    USE kind 
    28    USE phycst 
    29    USE global 
    30    USE fct 
    31    USE logger 
    32    USE dim 
    33    USE att 
    34    USE var 
     106   USE kind                            ! F90 kind parameter 
     107   USE phycst                          ! physical constant 
     108   USE global                          ! global variable 
     109   USE fct                             ! basic useful function 
     110   USE date                            ! date manager 
     111   USE logger                          ! log file manager 
     112   USE att                             ! attribute manager 
     113   USE dim                             ! dimension manager 
     114   USE var                             ! variable manager 
    35115 
    36116   IMPLICIT NONE 
    37    PRIVATE 
    38117   ! NOTE_avoid_public_variables_if_possible 
    39118 
    40119   ! type and variable 
     120   PRIVATE :: im_maxiter   !< default maximum number of iteration  
     121   PRIVATE :: im_minext    !< default minumum number of point to extrapolate 
     122   PRIVATE :: im_mincubic  !< default minumum number of point to extrapolate for cubic interpolation 
    41123 
    42124   ! function and subroutine 
    43125   PUBLIC :: extrap_detect         !< detected point to be extrapolated 
    44126   PUBLIC :: extrap_fill_value     !< extrapolate value over detected point  
    45    PUBLIC :: extrap_add_extrabands !<  
    46    PUBLIC :: extrap_del_extrabands !<  
    47    PUBLIC :: extrap_deriv_1D !<  
    48    PUBLIC :: extrap_deriv_2D !<  
    49  
    50    PRIVATE :: extrap__detect_wrapper !< detected point to be extrapolated 
    51    PRIVATE :: extrap__detect         !< detected point to be extrapolated 
    52    PRIVATE :: extrap__fill_value_wrapper     !< extrapolate value over detected point  
    53    PRIVATE :: extrap__fill_value     !< extrapolate value over detected point  
    54    PRIVATE :: extrap__3D 
    55    PRIVATE :: extrap_deriv_3D 
    56    PRIVATE :: extrap__3D_min_error_coef 
    57    PRIVATE :: extrap__3D_min_error_fill 
    58    PRIVATE :: extrap__3D_dist_weight_coef 
    59    PRIVATE :: extrap__3D_dist_weight_fill 
     127   PUBLIC :: extrap_add_extrabands !< add extraband to the variable (to be extrapolated)  
     128   PUBLIC :: extrap_del_extrabands !< delete extraband of the variable  
     129   PUBLIC :: extrap_deriv_1D       !< compute first derivative of 1D array  
     130   PUBLIC :: extrap_deriv_2D       !< compute first derivative of 2D array  
     131   PUBLIC :: extrap_deriv_3D       !< compute first derivative of 3D array 
     132 
     133   PRIVATE :: extrap__detect_wrapper      ! detected point to be extrapolated wrapper 
     134   PRIVATE :: extrap__detect              ! detected point to be extrapolated 
     135   PRIVATE :: extrap__fill_value_wrapper  ! extrapolate value over detected point wrapper  
     136   PRIVATE :: extrap__fill_value          ! extrapolate value over detected point  
     137   PRIVATE :: extrap__3D                  ! 
     138   PRIVATE :: extrap__3D_min_error_coef   ! 
     139   PRIVATE :: extrap__3D_min_error_fill   ! 
     140   PRIVATE :: extrap__3D_dist_weight_coef ! 
     141   PRIVATE :: extrap__3D_dist_weight_fill !  
    60142 
    61143   INTEGER(i4), PARAMETER :: im_maxiter = 10 !< default maximum number of iteration 
     
    64146 
    65147   INTERFACE extrap_detect 
    66       MODULE PROCEDURE extrap__detect_wrapper !< detected point to be extrapolated 
     148      MODULE PROCEDURE extrap__detect_wrapper     !< detected point to be extrapolated 
    67149   END INTERFACE extrap_detect 
    68150 
     
    74156   !------------------------------------------------------------------- 
    75157   !> @brief 
    76    !> This function detected point to be extrapolated. 
     158   !> This function detected point to be extrapolated, given variable structure. 
    77159   !>  
    78160   !> @details  
     161   !> optionaly, you could sepcify fine grid level, refinment factor (default 1),  
     162   !> offset between fine and coarse grid (default compute from refinment factor 
     163   !> as offset=(rho-1)/2), number of point to be extrapolated in each direction 
     164   !> (default im_minext).<br/> 
    79165   !> 
     166   !> First coarsening fine grid level, if need be, then select point near  
     167   !> grid point already inform. 
     168   !> 
     169   !> @note point to be extrapolated are selected using FillValue,  
     170   !> so to avoid mistake FillValue should not be zero (use var_chg_FillValue) 
     171   !>  
    80172   !> @author J.Paul 
    81    !> - Nov, 2013- Initial Version 
     173   !> - November, 2013- Initial Version 
    82174   ! 
    83    !> @param[in] td_var : variable to extrapolate 
    84    !> @param[in] id_iext : number of points to be extrapolated in i-direction 
    85    !> @param[in] id_jext : number of points to be extrapolated in j-direction 
    86    !> @param[in] id_kext : number of points to be extrapolated in k-direction 
    87    !> @return 
    88    !------------------------------------------------------------------- 
    89    !> @code 
     175   !> @param[in] td_var0   coarse grid variable to extrapolate 
     176   !> @param[in] td_level1 fine grid array of level 
     177   !> @param[in] id_offset array of offset between fine and coarse grid  
     178   !> @param[in] id_rho    array of refinment factor  
     179   !> @param[in] id_ext    array of number of points to be extrapolated 
     180   !> @return array of point to be extrapolated 
     181   !------------------------------------------------------------------- 
    90182   FUNCTION extrap__detect( td_var0, td_level1, & 
    91183   &                        id_offset, id_rho, id_ext ) 
    92184      IMPLICIT NONE 
    93185      ! Argument 
    94       TYPE(TVAR) ,                 INTENT(INOUT) :: td_var0 
     186      TYPE(TVAR) ,                 INTENT(IN   ) :: td_var0 
    95187      TYPE(TVAR) , DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: td_level1 
    96188      INTEGER(i4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset 
     
    106198      CHARACTER(LEN=lc)                                :: cl_level 
    107199 
    108       INTEGER(i4)                                      :: il_varid 
     200      INTEGER(i4)                                      :: il_ind 
    109201      INTEGER(i4)      , DIMENSION(:,:,:), ALLOCATABLE :: il_detect 
     202      INTEGER(i4)      , DIMENSION(:,:,:), ALLOCATABLE :: il_tmp 
    110203      INTEGER(i4)      , DIMENSION(:,:)  , ALLOCATABLE :: il_offset 
    111204      INTEGER(i4)      , DIMENSION(:,:)  , ALLOCATABLE :: il_level1 
     
    143236      ALLOCATE( il_offset(ip_maxdim,2) ) 
    144237      il_offset(:,:)=0 
    145       il_offset(jp_I,:)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5) 
    146       il_offset(jp_J,:)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5) 
    147238      IF( PRESENT(id_offset) )THEN 
    148239         il_offset(1:SIZE(id_offset(:,:),DIM=1),& 
    149240         &         1:SIZE(id_offset(:,:),DIM=2) )= id_offset(:,:) 
     241      ELSE 
     242         il_offset(jp_I,:)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5) 
     243         il_offset(jp_J,:)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5) 
    150244      ENDIF 
    151245 
     
    160254 
    161255      ! select point already inform 
    162       WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=1 
     256      DO jk0=1,td_var0%t_dim(3)%i_len 
     257         DO jj0=1,td_var0%t_dim(2)%i_len 
     258            DO ji0=1,td_var0%t_dim(1)%i_len 
     259               IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill ) il_detect(ji0,jj0,jk0)=1 
     260            ENDDO 
     261         ENDDO 
     262      ENDDO 
    163263  
    164264      IF( PRESENT(td_level1) )THEN 
    165265         SELECT CASE(TRIM(td_var0%c_point)) 
    166          CASE DEFAULT !'T' 
    167             cl_level='tlevel' 
    168          CASE('U') 
    169             cl_level='ulevel' 
    170          CASE('V') 
    171             cl_level='vlevel' 
    172          CASE('F') 
    173             cl_level='flevel' 
     266            CASE DEFAULT !'T' 
     267               cl_level='tlevel' 
     268            CASE('U') 
     269               cl_level='ulevel' 
     270            CASE('V') 
     271               cl_level='vlevel' 
     272            CASE('F') 
     273               cl_level='flevel' 
    174274         END SELECT 
    175          il_varid=var_get_id(td_level1(:),TRIM(cl_level)) 
    176          IF( il_varid == 0 )THEN 
     275 
     276         il_ind=var_get_index(td_level1(:),TRIM(cl_level)) 
     277         IF( il_ind == 0 )THEN 
    177278            CALL logger_error("EXTRAP DETECT: can not compute point to be "//& 
    178279            &     "extrapolated for variable "//TRIM(td_var0%c_name)//& 
     
    180281            &     "level for variable point "//TRIM(TRIM(td_var0%c_point))) 
    181282         ELSE 
    182             print *,'read ',trim(cl_level) 
    183             tl_var1=td_level1(il_varid) 
     283            tl_var1=var_copy(td_level1(il_ind)) 
     284 
     285            ALLOCATE( il_level1_G0( il_dim0(1), il_dim0(2)) ) 
     286            IF( ALL(tl_var1%t_dim(1:2)%i_len == il_dim0(1:2)) )THEN 
     287 
     288               ! variable to be extrapolated use same resolution than level 
     289               il_level1_G0(:,:)=INT(tl_var1%d_value(:,:,1,1),i4) 
     290                
     291            ELSE 
     292               ! variable to be extrapolated do not use same resolution than level 
     293               ALLOCATE( il_level1(tl_var1%t_dim(1)%i_len, & 
     294               &                   tl_var1%t_dim(2)%i_len) ) 
     295               ! match fine grid vertical level with coarse grid 
     296               il_level1(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)/il_rho(jp_K) 
     297 
     298               ALLOCATE( il_extra(ip_maxdim,2) ) 
     299               ! coarsening fine grid level 
     300               il_extra(jp_I,1)=CEILING(REAL(il_rho(jp_I)-1,dp)*0.5_dp) 
     301               il_extra(jp_I,2)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5_dp) 
     302 
     303               il_extra(jp_J,1)=CEILING(REAL(il_rho(jp_J)-1,dp)*0.5_dp) 
     304               il_extra(jp_J,2)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5_dp) 
     305 
     306               DO jj0=1,td_var0%t_dim(2)%i_len 
     307                   
     308                  jj1=(jj0-1)*il_rho(jp_J)+1-il_offset(jp_J,1) 
     309 
     310                  jj1m=MAX( jj1-il_extra(jp_J,1), 1 ) 
     311                  jj1p=MIN( jj1+il_extra(jp_J,2), & 
     312                  &         tl_var1%t_dim(2)%i_len-il_offset(jp_J,2) ) 
     313                   
     314                  DO ji0=1,td_var0%t_dim(1)%i_len 
     315 
     316                     ji1=(ji0-1)*il_rho(jp_I)+1-id_offset(jp_I,1) 
     317 
     318                     ji1m=MAX( ji1-il_extra(jp_I,1), 1 ) 
     319                     ji1p=MIN( ji1+il_extra(jp_I,2), & 
     320                     &         tl_var1%t_dim(1)%i_len-id_offset(jp_I,2) ) 
     321                
     322                     il_level1_G0(ji0,jj0)=MAXVAL(il_level1(ji1m:ji1p,jj1m:jj1p)) 
     323 
     324                  ENDDO 
     325               ENDDO 
     326 
     327               ! clean 
     328               DEALLOCATE( il_extra ) 
     329               DEALLOCATE( il_level1 ) 
     330 
     331            ENDIF 
     332 
     333            ! look for sea point 
     334            DO jk0=1,td_var0%t_dim(3)%i_len 
     335               WHERE( il_level1_G0(:,:) >= jk0) 
     336                  il_detect(:,:,jk0)=1 
     337               END WHERE 
     338            ENDDO 
     339 
     340            ! clean 
     341            DEALLOCATE( il_level1_G0 ) 
     342            CALL var_clean(tl_var1) 
     343 
    184344         ENDIF 
    185  
    186          ALLOCATE( il_level1_G0( il_dim0(1), il_dim0(2)) ) 
    187          IF( ALL(tl_var1%t_dim(1:2)%i_len == il_dim0(1:2)) )THEN 
    188  
    189             ! variable to be extrapolated use same resolution than level 
    190             il_level1_G0(:,:)=INT(tl_var1%d_value(:,:,1,1),i4) 
    191              
    192          ELSE 
    193             ! variable to be extrapolated do not use same resolution than level 
    194             ALLOCATE( il_level1(tl_var1%t_dim(1)%i_len, & 
    195             &                   tl_var1%t_dim(2)%i_len) ) 
    196             ! match fine grid vertical level with coarse grid 
    197             il_level1(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)/il_rho(jp_K) 
    198  
    199             ALLOCATE( il_extra(ig_ndim,2) ) 
    200             ! coarsening fine grid level 
    201             il_extra(jp_I,1)=CEILING(REAL(il_rho(jp_I)-1,dp)/2._dp) 
    202             il_extra(jp_I,2)=FLOOR(REAL(il_rho(jp_I)-1,dp)/2._dp) 
    203  
    204             il_extra(jp_J,1)=CEILING(REAL(il_rho(jp_J)-1,dp)/2._dp) 
    205             il_extra(jp_J,2)=FLOOR(REAL(il_rho(jp_J)-1,dp)/2._dp) 
    206  
    207             DO jj0=1,td_var0%t_dim(2)%i_len 
    208                 
    209                jj1=(jj0-1)*il_rho(jp_J)+1-il_offset(jp_J,1) 
    210  
    211                jj1m=MAX( jj1-il_extra(jp_J,1), 1 ) 
    212                jj1p=MIN( jj1+il_extra(jp_J,2), & 
    213                &         tl_var1%t_dim(2)%i_len-il_offset(jp_J,2) ) 
    214                 
    215                DO ji0=1,td_var0%t_dim(1)%i_len 
    216  
    217                   ji1=(ji0-1)*il_rho(jp_I)+1-id_offset(jp_I,1) 
    218  
    219                   ji1m=MAX( ji1-il_extra(jp_I,1), 1 ) 
    220                   ji1p=MIN( ji1+il_extra(jp_I,2), & 
    221                   &         tl_var1%t_dim(1)%i_len-id_offset(jp_I,2) ) 
    222              
    223                   il_level1_G0(ji0,jj0)=MAXVAL(il_level1(ji1m:ji1p,jj1m:jj1p)) 
    224                ENDDO 
    225             ENDDO 
    226  
    227             !il_level1_G0(:,:)=0 
    228             !DO jj1=1,tl_var1%t_dim(2)%i_len 
    229        
    230             !   jj0=INT(REAL((jj1+il_offset(jp_J,1)-1)-1,dp)/REAL(il_rho(jp_J),dp)) +1 
    231  
    232             !   DO ji1=1,tl_var1%t_dim(1)%i_len 
    233  
    234             !      ji0=INT(REAL((ji1+il_offset(jp_I,1)-1)-1,dp)/REAL(il_rho(jp_I),dp)) +1 
    235  
    236             !      il_level1_G0(ji0,jj0)=MAX(il_level1_G0(ji0,jj0),il_level1(ji1,jj1)) 
    237             !       
    238             !   ENDDO 
    239             !ENDDO 
    240  
    241             ! clean 
    242             DEALLOCATE( il_extra ) 
    243             DEALLOCATE( il_level1 ) 
    244  
    245          ENDIF 
    246  
    247          ! look for sea point 
    248          !il_detect(:,:,1)=0 
    249          DO jk0=1,td_var0%t_dim(3)%i_len 
    250             WHERE( il_level1_G0(:,:) >= jk0) 
    251                il_detect(:,:,jk0)=1 
    252             END WHERE 
    253             !il_detect(:,:,jk0)=il_level1_G0(:,:) 
    254             !WHERE( td_var0%d_value(:,:,jk0,1) /= td_var0%d_fill ) 
    255             !   il_detect(:,:,1)=jk0-1 
    256             !END WHERE 
    257          ENDDO 
    258  
    259          ! clean 
    260          DEALLOCATE( il_level1_G0 ) 
    261  
    262345      ENDIF 
    263346 
     
    265348      DEALLOCATE( il_offset ) 
    266349  
    267       !! select extra point depending on interpolation method 
    268       !! compute point near grid point already inform 
    269       !FORALL( ji0=1:il_dim0(1), & 
    270       !&       jj0=1:il_dim0(2), & 
    271       !&       jk0=1:il_dim0(3), & 
    272       !&       il_detect(ji0,jj0,jk0) == 1 ) 
    273  
    274       !   il_detect(MAX(1,ji0-il_ext(jp_I)):MIN(ji0+il_ext(jp_I),il_dim0(1)),& 
    275       !   &         MAX(1,jj0-il_ext(jp_J)):MIN(jj0+il_ext(jp_J),il_dim0(2)),& 
    276       !   &         MAX(1,jk0-il_ext(jp_K)):MIN(jk0+il_ext(jp_K),il_dim0(3)) )=1  
    277  
    278       !END FORALL 
     350      ALLOCATE( il_tmp(il_dim0(1),& 
     351      &                il_dim0(2),& 
     352      &                il_dim0(3)) ) 
     353      il_tmp(:,:,:)=il_detect(:,:,:) 
     354      ! select extra point depending on interpolation method 
     355      ! compute point near grid point already inform 
     356      DO jk0=1,il_dim0(3) 
     357         DO jj0=1,il_dim0(2) 
     358            DO ji0=1,il_dim0(1) 
     359 
     360               IF( il_tmp(ji0,jj0,jk0) == 1 )THEN 
     361                  il_detect( & 
     362                  &  MAX(1,ji0-il_ext(jp_I)):MIN(ji0+il_ext(jp_I),il_dim0(1)),& 
     363                  &  MAX(1,jj0-il_ext(jp_J)):MIN(jj0+il_ext(jp_J),il_dim0(2)),& 
     364                  &  MAX(1,jk0-il_ext(jp_K)):MIN(jk0+il_ext(jp_K),il_dim0(3)) & 
     365                  &  ) = 1  
     366               ENDIF 
     367 
     368            ENDDO 
     369         ENDDO 
     370      ENDDO 
    279371       
     372      ! clean 
     373      DEALLOCATE( il_tmp ) 
     374 
    280375      ! do not compute grid point already inform 
    281       WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=0 
     376      DO jk0=1,td_var0%t_dim(3)%i_len 
     377         DO jj0=1,td_var0%t_dim(2)%i_len 
     378            DO ji0=1,td_var0%t_dim(1)%i_len 
     379               IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill ) il_detect(ji0,jj0,jk0)=0 
     380            ENDDO 
     381         ENDDO 
     382      ENDDO 
    282383 
    283384      ! save result 
     
    288389      DEALLOCATE( il_ext ) 
    289390      DEALLOCATE( il_detect ) 
     391      DEALLOCATE( il_rho ) 
    290392 
    291393   END FUNCTION extrap__detect 
    292    !> @endcode 
    293394   !------------------------------------------------------------------- 
    294395   !> @brief 
    295    !> This function detected point to be extrapolated. 
     396   !> This function sort variable to be extrapolated, depending on number of 
     397   !> dimentsion, then detected point to be extrapolated. 
    296398   !>  
    297    !> @details  
     399   !> @author J.Paul 
     400   !> - November, 2013- Initial Version 
    298401   !> 
    299    !> @author J.Paul 
    300    !> - Nov, 2013- Initial Version 
    301    ! 
    302    !> @param[in] td_var : variable to extrapolate 
    303    !> @param[in] id_iext : number of points to be extrapolated in i-direction 
    304    !> @param[in] id_jext : number of points to be extrapolated in j-direction 
    305    !> @param[in] id_kext : number of points to be extrapolated in k-direction 
    306    !> @return 
    307    !------------------------------------------------------------------- 
    308    !> @code 
     402   !> @param[in] td_var    coarse grid variable to extrapolate 
     403   !> @param[in] td_level  fine grid array of level 
     404   !> @param[in] id_offset array of offset between fine and coarse grid  
     405   !> @param[in] id_rho    array of refinment factor  
     406   !> @param[in] id_ext    array of number of points to be extrapolated 
     407   !> @return 3D array of point to be extrapolated 
     408   !------------------------------------------------------------------- 
    309409   FUNCTION extrap__detect_wrapper( td_var, td_level, & 
    310410   &                                id_offset, id_rho, id_ext ) 
     
    312412      IMPLICIT NONE 
    313413      ! Argument 
    314       TYPE(TVAR) ,                         INTENT(INOUT) :: td_var 
    315       TYPE(TVAR) , DIMENSION(:)          , INTENT(IN   ), OPTIONAL :: td_level 
     414      TYPE(TVAR) ,                 INTENT(IN   ) :: td_var 
     415      TYPE(TVAR) , DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: td_level 
    316416      INTEGER(i4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset 
    317417      INTEGER(i4), DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_rho 
     
    335435      ELSE IF( ALL(td_var%t_dim(1:3)%l_use) )THEN 
    336436          
    337          ! detect point to be interpolated on I-J-K 
     437         ! detect point to be extrapolated on I-J-K 
    338438         CALL logger_debug(" EXTRAP DETECT: detect point "//& 
    339439         &              " for variable "//TRIM(td_var%c_name) ) 
     
    346446      ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN 
    347447 
    348          ! detect point to be interpolated on I-J 
     448         ! detect point to be extrapolated on I-J 
    349449         CALL logger_debug(" EXTRAP DETECT: detect horizontal point "//& 
    350450         &              " for variable "//TRIM(td_var%c_name) ) 
     
    357457      ELSE IF( td_var%t_dim(3)%l_use )THEN 
    358458          
    359          ! detect point to be interpolated on K 
     459         ! detect point to be extrapolated on K 
    360460         CALL logger_debug(" EXTRAP DETECT: detect vertical point "//& 
    361461         &              " for variable "//TRIM(td_var%c_name) ) 
     
    373473       
    374474   END FUNCTION extrap__detect_wrapper 
    375    !> @endcode 
    376 !   !------------------------------------------------------------------- 
    377 !   !> @brief 
    378 !   !> This function detected point to be extrapolated. 
    379 !   !>  
    380 !   !> @details  
    381 !   !> 
    382 !   !> @author J.Paul 
    383 !   !> - Nov, 2013- Initial Version 
    384 !   ! 
    385 !   !> @param[in] td_var : variable to extrapolate 
    386 !   !> @param[in] id_iext : number of points to be extrapolated in i-direction 
    387 !   !> @param[in] id_jext : number of points to be extrapolated in j-direction 
    388 !   !> @param[in] id_kext : number of points to be extrapolated in k-direction 
    389 !   !> @return  
    390 !   ! 
    391 !   !> @todo 
    392 !   !------------------------------------------------------------------- 
    393 !   !> @code 
    394 !   FUNCTION extrap__detect(td_var, id_iext, id_jext, id_kext) 
    395 !      IMPLICIT NONE 
    396 !      ! Argument 
    397 !      TYPE(TVAR) ,                   INTENT(INOUT) :: td_var 
    398 !      !INTEGER(i4), DIMENSION(:,:,:), INTENT(OUT  ) :: id_detect 
    399 !      INTEGER(i4),                   INTENT(IN   ), OPTIONAL :: id_iext 
    400 !      INTEGER(i4),                   INTENT(IN   ), OPTIONAL :: id_jext 
    401 !      INTEGER(i4),                   INTENT(IN   ), OPTIONAL :: id_kext 
    402 ! 
    403 !      ! function 
    404 !      INTEGER(i4), DIMENSION(td_var%t_dim(1)%i_len, & 
    405 !      &                      td_var%t_dim(2)%i_len, & 
    406 !      &                      td_var%t_dim(3)%i_len) :: extrap__detect 
    407 ! 
    408 !      ! local variable 
    409 !      INTEGER(i4) :: il_iext 
    410 !      INTEGER(i4) :: il_jext 
    411 !      INTEGER(i4) :: il_kext 
    412 ! 
    413 !      INTEGER(i4), DIMENSION(3) :: il_dim 
    414 ! 
    415 !      ! loop indices 
    416 !      INTEGER(i4) :: ji 
    417 !      INTEGER(i4) :: jj 
    418 !      INTEGER(i4) :: jk 
    419 !      !---------------------------------------------------------------- 
    420 ! 
    421 !      ! optional argument 
    422 !      il_iext=im_minext 
    423 !      IF( PRESENT(id_iext) ) il_iext=id_iext 
    424 !      il_jext=im_minext 
    425 !      IF( PRESENT(id_jext) ) il_jext=id_jext 
    426 !      il_kext=im_minext 
    427 !      IF( PRESENT(id_kext) ) il_kext=id_kext 
    428 ! 
    429 !      ! init 
    430 !      extrap__detect(:,:,:)=0 
    431 ! 
    432 !      il_dim(:)=td_var%t_dim(1:3)%i_len 
    433 ! 
    434 !      ! compute point near grid point already inform 
    435 !      FORALL( ji=1:il_dim(1), & 
    436 !      &       jj=1:il_dim(2), & 
    437 !      &       jk=1:il_dim(3), & 
    438 !      &       td_var%d_value(ji,jj,jk,1) /= td_var%d_fill ) 
    439 ! 
    440 !         extrap__detect(MAX(1,ji-il_iext):MIN(ji+il_iext,il_dim(1)),& 
    441 !         &              MAX(1,jj-il_jext):MIN(jj+il_jext,il_dim(2)),& 
    442 !         &              MAX(1,jk-il_kext):MIN(jk+il_kext,il_dim(3)) )=1  
    443 ! 
    444 ! 
    445 !      END FORALL 
    446 !       
    447 !      ! do not compute grid point already inform 
    448 !      WHERE( td_var%d_value(:,:,:,1) /= td_var%d_fill ) extrap__detect(:,:,:)=0 
    449 ! 
    450 !   END FUNCTION extrap__detect 
    451 !   !> @endcode 
    452475   !------------------------------------------------------------------- 
    453476   !> @brief 
    454    !> This subroutine 
     477   !> This subroutine select method to be used for extrapolation. 
     478   !> If need be, increase number of points to be extrapolated. 
     479   !> Finally launch extrap__fill_value. 
    455480   !>  
    456    !> @details  
     481   !> @details 
     482   !> optionaly, you could specify :<br/> 
     483   !>  - refinment factor (default 1) 
     484   !>  - offset between fine and coarse grid (default compute from refinment factor 
     485   !> as offset=(rho-1)/2)  
     486   !>  - number of point to be extrapolated in each direction (default im_minext) 
     487   !>  - radius of the halo used to compute extrapolation 
     488   !>  - maximum number of iteration 
    457489   !> 
    458490   !> @author J.Paul 
    459491   !> - Nov, 2013- Initial Version 
    460492   ! 
    461    !> @param[inout] td_var : variable structure 
    462    !> @param[in] id_iext : number of points to be extrapolated in i-direction 
    463    !> @param[in] id_jext : number of points to be extrapolated in j-direction 
    464    !> @param[in] id_kext : number of points to be extrapolated in k-direction 
    465    !> @param[in] id_extend :  radius of the box used to compute extrapolation  
    466    !> @param[in] id_maxiter : maximum nuber of iteration 
    467    !------------------------------------------------------------------- 
    468    !> @code 
     493   !> @param[inout] td_var    variable structure 
     494   !> @param[in] td_level     fine grid array of level 
     495   !> @param[in] id_offset    array of offset between fine and coarse grid  
     496   !> @param[in] id_rho       array of refinment factor  
     497   !> @param[in] id_iext      number of points to be extrapolated in i-direction 
     498   !> @param[in] id_jext      number of points to be extrapolated in j-direction 
     499   !> @param[in] id_kext      number of points to be extrapolated in k-direction 
     500   !> @param[in] id_radius    radius of the halo used to compute extrapolation  
     501   !> @param[in] id_maxiter   maximum number of iteration 
     502   !------------------------------------------------------------------- 
    469503   SUBROUTINE extrap__fill_value_wrapper( td_var, td_level, & 
    470504   &                                      id_offset,        & 
     
    496530      !---------------------------------------------------------------- 
    497531      IF( .NOT. ASSOCIATED(td_var%d_value) )THEN 
    498          CALL logger_error("EXTRAP FILL VALUE: no table of value "//& 
     532         CALL logger_error("EXTRAP FILL VALUE: no value "//& 
    499533         &  "associted to variable "//TRIM(td_var%c_name) ) 
    500534      ELSE 
     
    542576         IF( (il_iext /= 0 .AND. td_var%t_dim(1)%l_use) .OR. & 
    543577         &   (il_jext /= 0 .AND. td_var%t_dim(2)%l_use) .OR. & 
    544          &   (il_kext /= 0 .AND. td_var%t_dim(3)%l_use) .OR. & 
    545          &   PRESENT(td_level) )THEN 
     578         &   (il_kext /= 0 .AND. td_var%t_dim(3)%l_use) )THEN 
    546579 
    547580            ! number of point use to compute box 
     
    577610 
    578611   END SUBROUTINE extrap__fill_value_wrapper 
    579    !> @endcode 
    580612   !------------------------------------------------------------------- 
    581613   !> @brief 
    582    !> This subroutine 
     614   !> This subroutine compute point to be extrapolated, then extrapolate point. 
    583615   !>  
    584616   !> @details  
     617   !> optionaly, you could specify :<br/> 
     618   !>  - refinment factor (default 1) 
     619   !>  - offset between fine and coarse grid (default compute from refinment factor 
     620   !> as offset=(rho-1)/2)  
    585621   !> 
    586622   !> @author J.Paul 
    587    !> - Nov, 2013- Initial Version 
     623   !> - November, 2013- Initial Version 
    588624   ! 
    589    !> @param[inout] td_var : variable structure 
    590    !> @param[in] cd_method : extrapolation method 
    591    !> @param[in] id_iext : number of points to be extrapolated in i-direction 
    592    !> @param[in] id_jext : number of points to be extrapolated in j-direction 
    593    !> @param[in] id_kext : number of points to be extrapolated in k-direction 
    594    !> @param[in] id_radius : radius of the halo used to compute extrapolation 
    595    !> @param[in] id_maxiter : maximum nuber of iteration 
    596    !------------------------------------------------------------------- 
    597    !> @code 
     625   !> @param[inout] td_var    variable structure 
     626   !> @param[in] cd_method    extrapolation method 
     627   !> @param[in] id_iext      number of points to be extrapolated in i-direction 
     628   !> @param[in] id_jext      number of points to be extrapolated in j-direction 
     629   !> @param[in] id_kext      number of points to be extrapolated in k-direction 
     630   !> @param[in] id_radius    radius of the halo used to compute extrapolation 
     631   !> @param[in] id_maxiter   maximum number of iteration 
     632   !> @param[in] td_level     fine grid array of level 
     633   !> @param[in] id_offset    array of offset between fine and coarse grid  
     634   !> @param[in] id_rho       array of refinment factor 
     635   !------------------------------------------------------------------- 
    598636   SUBROUTINE extrap__fill_value( td_var, cd_method, & 
    599637   &                              id_iext, id_jext, id_kext, & 
     
    619657 
    620658      INTEGER(i4), DIMENSION(:,:,:)  , ALLOCATABLE :: il_detect 
    621       INTEGER(i4)                                  :: il_radius 
    622       INTEGER(i4)                                  :: il_iter 
    623659 
    624660      TYPE(TATT)                                   :: tl_att 
    625661 
    626662      ! loop indices 
    627       INTEGER(i4) :: jl 
    628663      !---------------------------------------------------------------- 
    629664 
     
    637672      &                                 id_rho,           & 
    638673      &                                 id_ext=(/id_iext, id_jext, id_kext/) ) 
    639  
    640674      !2- add attribute to variable 
    641675      cl_extrap=fct_concat(td_var%c_extrap(:)) 
     
    643677      CALL var_move_att(td_var, tl_att) 
    644678 
    645       CALL logger_warn(" EXTRAP FILL: "//& 
     679      CALL att_clean(tl_att) 
     680 
     681      CALL logger_info(" EXTRAP FILL: "//& 
    646682         &              TRIM(fct_str(SUM(il_detect(:,:,:))))//& 
    647683         &              " point(s) to extrapolate " ) 
    648684 
    649685      !3- extrapolate 
    650       DO jl=1,td_var%t_dim(4)%i_len 
    651          
    652 !            td_var%d_value(:,:,:,jl)=il_detect(:,:,:) 
    653          il_iter=1 
    654          DO WHILE( ANY(il_detect(:,:,:)==1) ) 
    655             ! change extend value to minimize number of iteration 
    656             il_radius=id_radius+(il_iter/id_maxiter) 
    657              
    658             CALL logger_debug(" EXTRAP FILL VALUE: "//& 
    659             &              TRIM(fct_str(SUM(il_detect(:,:,:))))//& 
    660             &              " points to extrapolate " ) 
    661              
    662             CALL extrap__3D(td_var%d_value(:,:,:,jl), td_var%d_fill,    & 
    663             &               il_detect(:,:,:),                           & 
    664             &               cd_method, il_radius ) 
    665                 
    666             il_iter=il_iter+1 
    667          ENDDO 
    668  
    669       ENDDO 
    670   
    671       IF( SUM(il_detect(:,:,:)) /= 0 )THEN 
    672          CALL logger_warn(" EXTRAP FILL: still "//& 
    673          &              TRIM(fct_str(SUM(il_detect(:,:,:))))//& 
    674          &              " point(s) to extrapolate " ) 
    675       ENDIF 
     686      CALL extrap__3D(td_var%d_value(:,:,:,:), td_var%d_fill,    & 
     687      &               il_detect(:,:,:),                           & 
     688      &               cd_method, id_radius, id_maxiter  ) 
    676689 
    677690      DEALLOCATE(il_detect) 
    678691 
    679692   END SUBROUTINE extrap__fill_value 
    680    !> @endcode 
    681693   !------------------------------------------------------------------- 
    682694   !> @brief 
    683    !> This subroutine 
     695   !> This subroutine compute point to be extrapolated in 3D array. 
    684696   !>  
    685697   !> @details  
     698   !> in case of 'min_error' method:<br/> 
     699   !>    - compute derivative in i-, j- and k- direction 
     700   !>    - compute minimum error coefficient (distance to center of halo) 
     701   !>    - compute extrapolatd values by calculated minimum error using taylor expansion  
     702   !> in case of 'dist_weight' method:<br/> 
     703   !>    - compute distance weight coefficient (inverse of distance to center of halo)  
     704   !>    - compute extrapolatd values using Inverse Distance Weighting 
    686705   !> 
    687706   !> @author J.Paul 
    688707   !> - Nov, 2013- Initial Version 
    689708   ! 
    690    !> @param[inout] dd_value : 3D table of variable to be extrapolated 
    691    !> @param[in] dd_fill : FillValue of variable 
    692    !> @param[inout] id_detect : table of point to extrapolate 
    693    !> @param[in] id_ext : number of point use to compute box  
    694    !> @param[in] cd_method : extrapolation method 
    695    !------------------------------------------------------------------- 
    696    !> @code 
     709   !> @param[inout] dd_value  3D array of variable to be extrapolated 
     710   !> @param[in] dd_fill      FillValue of variable 
     711   !> @param[inout] id_detect array of point to extrapolate 
     712   !> @param[in] cd_method    extrapolation method 
     713   !> @param[in] id_radius    radius of the halo used to compute extrapolation 
     714   !------------------------------------------------------------------- 
    697715   SUBROUTINE extrap__3D( dd_value, dd_fill, id_detect,& 
    698    &                      cd_method, id_ext ) 
     716   &                      cd_method, id_radius, id_maxiter ) 
    699717      IMPLICIT NONE 
    700718      ! Argument 
    701       REAL(dp)   , DIMENSION(:,:,:), INTENT(INOUT) :: dd_value 
     719      REAL(dp)   , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_value 
    702720      REAL(dp)   ,                   INTENT(IN   ) :: dd_fill 
    703721      INTEGER(i4), DIMENSION(:,:,:), INTENT(INOUT) :: id_detect 
    704722      CHARACTER(LEN=*),              INTENT(IN   ) :: cd_method 
    705       INTEGER(i4),                   INTENT(IN   ) :: id_ext 
     723      INTEGER(i4),                   INTENT(IN   ) :: id_radius 
     724      INTEGER(i4),                   INTENT(IN   ) :: id_maxiter 
    706725 
    707726      ! local variable 
     
    712731      INTEGER(i4) :: il_kmin 
    713732      INTEGER(i4) :: il_kmax 
    714  
    715       INTEGER(i4), DIMENSION(3) :: il_shape 
     733      INTEGER(i4) :: il_iter 
     734      INTEGER(i4) :: il_radius 
     735 
     736      INTEGER(i4), DIMENSION(4) :: il_shape 
    716737      INTEGER(i4), DIMENSION(3) :: il_dim 
     738 
     739      INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect 
    717740 
    718741      REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx 
     
    725748      INTEGER(i4) :: jj 
    726749      INTEGER(i4) :: jk 
     750      INTEGER(i4) :: jl 
    727751      !---------------------------------------------------------------- 
    728752 
    729753      il_shape(:)=SHAPE(dd_value) 
    730754 
     755      ALLOCATE( il_detect( il_shape(1), il_shape(2), il_shape(3)) ) 
     756 
    731757      SELECT CASE(TRIM(cd_method)) 
    732          CASE('min_error') 
    733  
    734             ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) )  
    735             ALLOCATE( dl_dfdy(il_shape(1), il_shape(2), il_shape(3)) )  
    736             ALLOCATE( dl_dfdz(il_shape(1), il_shape(2), il_shape(3)) )  
    737  
    738  
    739             ! compute derivative in i-direction 
    740             dl_dfdx(:,:,:)=dd_fill 
    741             IF( il_shape(1) > 1 )THEN 
    742                dl_dfdx(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'I' ) 
    743             ENDIF 
    744  
    745             ! compute derivative in i-direction 
    746             dl_dfdy(:,:,:)=dd_fill 
    747             IF( il_shape(2) > 1 )THEN 
    748                dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'J' ) 
    749             ENDIF 
    750                 
    751             ! compute derivative in i-direction 
    752             dl_dfdz(:,:,:)=dd_fill 
    753             IF( il_shape(3) > 1 )THEN 
    754                dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'K' ) 
    755             ENDIF 
    756        
    757             il_dim(1)=2*id_ext+1 
    758             IF( il_shape(1) < 2*id_ext+1 ) il_dim(1)=1 
    759             il_dim(2)=2*id_ext+1 
    760             IF( il_shape(2) < 2*id_ext+1 ) il_dim(2)=1 
    761             il_dim(3)=2*id_ext+1 
    762             IF( il_shape(3) < 2*id_ext+1 ) il_dim(3)=1 
    763  
    764             ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) )  
    765  
    766             dl_coef(:,:,:)=extrap__3D_min_error_coef(dd_value( 1:il_dim(1), & 
    767             &                                                  1:il_dim(2), & 
    768             &                                                  1:il_dim(3))) 
    769  
    770             DO jk=1,il_shape(3) 
    771                IF( ALL(id_detect(:,:,jk) == 0) ) CYCLE 
    772                DO jj=1,il_shape(2) 
    773                   IF( ALL(id_detect(:,jj,jk) == 0) ) CYCLE 
    774                   DO ji=1,il_shape(1) 
    775  
    776                      IF( id_detect(ji,jj,jk) == 1 )THEN 
    777                         
    778                         il_imin=MAX(ji-id_ext,1) 
    779                         il_imax=MIN(ji+id_ext,il_shape(1)) 
    780                         IF( il_dim(1) == 1 )THEN 
    781                            il_imin=ji 
    782                            il_imax=ji 
     758      CASE('min_error') 
     759         DO jl=1,il_shape(4) 
     760 
     761            ! intitialise table of poitn to be extrapolated 
     762            il_detect(:,:,:)=id_detect(:,:,:) 
     763 
     764            il_iter=1 
     765            DO WHILE( ANY(il_detect(:,:,:)==1) ) 
     766               ! change extend value to minimize number of iteration 
     767               il_radius=id_radius+(il_iter/id_maxiter) 
     768 
     769               ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) )  
     770               ALLOCATE( dl_dfdy(il_shape(1), il_shape(2), il_shape(3)) )  
     771               ALLOCATE( dl_dfdz(il_shape(1), il_shape(2), il_shape(3)) )  
     772 
     773               ! compute derivative in i-direction 
     774               dl_dfdx(:,:,:)=dd_fill 
     775               IF( il_shape(1) > 1 )THEN 
     776                  dl_dfdx(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'I' ) 
     777               ENDIF 
     778 
     779               ! compute derivative in j-direction 
     780               dl_dfdy(:,:,:)=dd_fill 
     781               IF( il_shape(2) > 1 )THEN 
     782                  dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'J' ) 
     783               ENDIF 
     784  
     785               ! compute derivative in k-direction 
     786               dl_dfdz(:,:,:)=dd_fill 
     787               IF( il_shape(3) > 1 )THEN 
     788                  dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'K' ) 
     789               ENDIF 
     790  
     791               il_dim(1)=2*il_radius+1 
     792               IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1 
     793               il_dim(2)=2*il_radius+1 
     794               IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1 
     795               il_dim(3)=2*il_radius+1 
     796               IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1 
     797 
     798               ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) )  
     799 
     800               dl_coef(:,:,:)=extrap__3D_min_error_coef(dd_value( 1:il_dim(1), & 
     801               &                                                  1:il_dim(2), & 
     802               &                                                  1:il_dim(3), & 
     803               &                                                  jl )) 
     804 
     805               DO jk=1,il_shape(3) 
     806                  IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 
     807                  DO jj=1,il_shape(2) 
     808                     IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 
     809                     DO ji=1,il_shape(1) 
     810 
     811                        IF( il_detect(ji,jj,jk) == 1 )THEN 
     812                           
     813                           il_imin=MAX(ji-il_radius,1) 
     814                           il_imax=MIN(ji+il_radius,il_shape(1)) 
     815                           IF( il_dim(1) == 1 )THEN 
     816                              il_imin=ji 
     817                              il_imax=ji 
     818                           ENDIF 
     819 
     820                           il_jmin=MAX(jj-il_radius,1) 
     821                           il_jmax=MIN(jj+il_radius,il_shape(2)) 
     822                           IF( il_dim(2) == 1 )THEN 
     823                              il_jmin=jj 
     824                              il_jmax=jj 
     825                           ENDIF 
     826 
     827                           il_kmin=MAX(jk-il_radius,1) 
     828                           il_kmax=MIN(jk+il_radius,il_shape(3)) 
     829                           IF( il_dim(3) == 1 )THEN 
     830                              il_kmin=jk 
     831                              il_kmax=jk 
     832                           ENDIF 
     833 
     834                           dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( & 
     835                           &  dd_value( il_imin:il_imax, & 
     836                           &            il_jmin:il_jmax, & 
     837                           &            il_kmin:il_kmax,jl ), dd_fill, il_radius, & 
     838                           &  dl_dfdx(  il_imin:il_imax, & 
     839                           &            il_jmin:il_jmax, & 
     840                           &            il_kmin:il_kmax ), & 
     841                           &  dl_dfdy(  il_imin:il_imax, & 
     842                           &            il_jmin:il_jmax, & 
     843                           &            il_kmin:il_kmax ), & 
     844                           &  dl_dfdz(  il_imin:il_imax, & 
     845                           &            il_jmin:il_jmax, & 
     846                           &            il_kmin:il_kmax ), & 
     847                           &  dl_coef(:,:,:) ) 
     848 
     849                           IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 
     850                              il_detect(ji,jj,jk)= 0 
     851                           ENDIF 
     852 
    783853                        ENDIF 
    784854 
    785                         il_jmin=MAX(jj-id_ext,1) 
    786                         il_jmax=MIN(jj+id_ext,il_shape(2)) 
    787                         IF( il_dim(2) == 1 )THEN 
    788                            il_jmin=jj 
    789                            il_jmax=jj 
    790                         ENDIF 
    791  
    792                         il_kmin=MAX(jk-id_ext,1) 
    793                         il_kmax=MIN(jk+id_ext,il_shape(3)) 
    794                         IF( il_dim(3) == 1 )THEN 
    795                            il_kmin=jk 
    796                            il_kmax=jk 
    797                         ENDIF 
    798  
    799                         dd_value(ji,jj,jk)=extrap__3D_min_error_fill( & 
    800                         &  dd_value( il_imin:il_imax, & 
    801                         &            il_jmin:il_jmax, & 
    802                         &            il_kmin:il_kmax ), dd_fill, id_ext, & 
    803                         &  dl_dfdx(  il_imin:il_imax, & 
    804                         &            il_jmin:il_jmax, & 
    805                         &            il_kmin:il_kmax ), & 
    806                         &  dl_dfdy(  il_imin:il_imax, & 
    807                         &            il_jmin:il_jmax, & 
    808                         &            il_kmin:il_kmax ), & 
    809                         &  dl_dfdz(  il_imin:il_imax, & 
    810                         &            il_jmin:il_jmax, & 
    811                         &            il_kmin:il_kmax ), & 
    812                         &  dl_coef(:,:,:) ) 
    813  
    814                         IF( dd_value(ji,jj,jk) /= dd_fill )THEN 
    815                            id_detect(ji,jj,jk)= 0 
    816                         ENDIF 
    817  
    818                      ENDIF 
    819  
     855                     ENDDO 
    820856                  ENDDO 
    821857               ENDDO 
     858 
     859               DEALLOCATE( dl_dfdx ) 
     860               DEALLOCATE( dl_dfdy ) 
     861               DEALLOCATE( dl_dfdz ) 
     862               DEALLOCATE( dl_coef ) 
     863 
     864               il_iter=il_iter+1 
    822865            ENDDO 
    823  
    824             IF( ALLOCATED(dl_dfdx) ) DEALLOCATE( dl_dfdx ) 
    825             IF( ALLOCATED(dl_dfdy) ) DEALLOCATE( dl_dfdy ) 
    826             IF( ALLOCATED(dl_dfdz) ) DEALLOCATE( dl_dfdz ) 
    827             IF( ALLOCATED(dl_coef) ) DEALLOCATE( dl_coef ) 
    828  
    829          CASE DEFAULT ! 'dist_weight' 
    830  
    831             il_dim(1)=2*id_ext+1 
    832             IF( il_shape(1) < 2*id_ext+1 ) il_dim(1)=1 
    833             il_dim(2)=2*id_ext+1 
    834             IF( il_shape(2) < 2*id_ext+1 ) il_dim(2)=1 
    835             il_dim(3)=2*id_ext+1 
    836             IF( il_shape(3) < 2*id_ext+1 ) il_dim(3)=1 
    837              
    838             ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 
    839  
    840             dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), & 
    841             &                                                   1:il_dim(2), & 
    842             &                                                   1:il_dim(3)) ) 
    843  
    844             DO jk=1,il_shape(3) 
    845                IF( ALL(id_detect(:,:,jk) == 0) ) CYCLE 
    846                DO jj=1,il_shape(2) 
    847                   IF( ALL(id_detect(:,jj,jk) == 0) ) CYCLE 
    848                   DO ji=1,il_shape(1) 
    849  
    850                      IF( id_detect(ji,jj,jk) == 1 )THEN 
    851                          
    852                         il_imin=MAX(ji-id_ext,1) 
    853                         il_imax=MIN(ji+id_ext,il_shape(1)) 
    854                         IF( il_dim(1) == 1 )THEN 
    855                            il_imin=ji 
    856                            il_imax=ji 
     866         ENDDO 
     867 
     868      CASE DEFAULT ! 'dist_weight' 
     869         DO jl=1,il_shape(4) 
     870 
     871            ! intitialise table of poitn to be extrapolated 
     872            il_detect(:,:,:)=id_detect(:,:,:) 
     873 
     874            il_iter=1 
     875            DO WHILE( ANY(il_detect(:,:,:)==1) ) 
     876               ! change extend value to minimize number of iteration 
     877               il_radius=id_radius+(il_iter/id_maxiter) 
     878 
     879               il_dim(1)=2*il_radius+1 
     880               IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1 
     881               il_dim(2)=2*il_radius+1 
     882               IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1 
     883               il_dim(3)=2*il_radius+1 
     884               IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1 
     885                
     886               ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 
     887 
     888               dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), & 
     889               &                                                   1:il_dim(2), & 
     890               &                                                   1:il_dim(3), & 
     891               &                                                   jl ) ) 
     892 
     893               DO jk=1,il_shape(3) 
     894                  IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 
     895                  DO jj=1,il_shape(2) 
     896                     IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 
     897                     DO ji=1,il_shape(1) 
     898 
     899                        IF( il_detect(ji,jj,jk) == 1 )THEN 
     900                            
     901                           il_imin=MAX(ji-il_radius,1) 
     902                           il_imax=MIN(ji+il_radius,il_shape(1)) 
     903                           IF( il_dim(1) == 1 )THEN 
     904                              il_imin=ji 
     905                              il_imax=ji 
     906                           ENDIF 
     907 
     908                           il_jmin=MAX(jj-il_radius,1) 
     909                           il_jmax=MIN(jj+il_radius,il_shape(2)) 
     910                           IF( il_dim(2) == 1 )THEN 
     911                              il_jmin=jj 
     912                              il_jmax=jj 
     913                           ENDIF 
     914 
     915                           il_kmin=MAX(jk-il_radius,1) 
     916                           il_kmax=MIN(jk+il_radius,il_shape(3)) 
     917                           IF( il_dim(3) == 1 )THEN 
     918                              il_kmin=jk 
     919                              il_kmax=jk 
     920                           ENDIF 
     921 
     922                           dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( & 
     923                           &  dd_value( il_imin:il_imax, & 
     924                           &            il_jmin:il_jmax, & 
     925                           &            il_kmin:il_kmax, & 
     926                           &            jl), dd_fill, il_radius, & 
     927                           &  dl_coef(:,:,:) ) 
     928 
     929                           IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 
     930                              il_detect(ji,jj,jk)= 0 
     931                           ENDIF 
     932 
    857933                        ENDIF 
    858934 
    859                         il_jmin=MAX(jj-id_ext,1) 
    860                         il_jmax=MIN(jj+id_ext,il_shape(2)) 
    861                         IF( il_dim(2) == 1 )THEN 
    862                            il_jmin=jj 
    863                            il_jmax=jj 
    864                         ENDIF 
    865  
    866                         il_kmin=MAX(jk-id_ext,1) 
    867                         il_kmax=MIN(jk+id_ext,il_shape(3)) 
    868                         IF( il_dim(3) == 1 )THEN 
    869                            il_kmin=jk 
    870                            il_kmax=jk 
    871                         ENDIF 
    872  
    873                         dd_value(ji,jj,jk)=extrap__3D_dist_weight_fill( & 
    874                         &  dd_value( il_imin:il_imax, & 
    875                         &            il_jmin:il_jmax, & 
    876                         &            il_kmin:il_kmax ), dd_fill, id_ext, & 
    877                         &  dl_coef(:,:,:) ) 
    878  
    879                         IF( dd_value(ji,jj,jk) /= dd_fill )THEN 
    880                            id_detect(ji,jj,jk)= 0 
    881                         ENDIF 
    882  
    883                      ENDIF 
    884  
     935                     ENDDO 
    885936                  ENDDO 
    886937               ENDDO 
     938 
     939               DEALLOCATE( dl_coef ) 
     940               il_iter=il_iter+1 
    887941            ENDDO 
    888  
    889             IF( ALLOCATED(dl_coef) ) DEALLOCATE( dl_coef ) 
    890              
     942         ENDDO             
    891943      END SELECT 
    892944 
     945      DEALLOCATE( il_detect ) 
     946 
    893947   END SUBROUTINE extrap__3D 
    894    !> @endcode 
    895948   !------------------------------------------------------------------- 
    896949   !> @brief 
    897    !> This function compute derivative of a function in i- and 
    898    !> j-direction 
     950   !> This function compute derivative of 1D array. 
    899951   !>  
    900952   !> @details  
     953   !> optionaly you could specify to take into account east west discontinuity 
     954   !> (-180° 180° or 0° 360° for longitude variable) 
    901955   !> 
    902956   !> @author J.Paul 
    903    !> - Nov, 2013- Initial Version 
     957   !> - November, 2013- Initial Version 
    904958   ! 
    905    !> @param[in] dd_value : 1D table of variable to be extrapolated 
    906    !> @param[in] dd_fill : FillValue of variable 
    907    !> @param[in] cd_dim : compute derivative on first (I) or second (J) dimension  
    908    !------------------------------------------------------------------- 
    909    !> @code 
     959   !> @param[in] dd_value     1D array of variable to be extrapolated 
     960   !> @param[in] dd_fill      FillValue of variable 
     961   !> @param[in] ld_discont   logical to take into account east west discontinuity  
     962   !------------------------------------------------------------------- 
    910963   PURE FUNCTION extrap_deriv_1D( dd_value, dd_fill, ld_discont ) 
    911964 
     
    10031056 
    10041057   END FUNCTION extrap_deriv_1D 
    1005    !> @endcode 
    10061058   !------------------------------------------------------------------- 
    10071059   !> @brief 
    1008    !> This function compute derivative of a function in i- and 
    1009    !> j-direction 
    1010    !>  
     1060   !> This function compute derivative of 2D array. 
     1061   !> you have to specify in which direction derivative have to be computed: 
     1062   !> first (I) or second (J) dimension.  
     1063   !> 
    10111064   !> @details  
     1065   !> optionaly you could specify to take into account east west discontinuity 
     1066   !> (-180° 180° or 0° 360° for longitude variable) 
    10121067   !> 
    10131068   !> @author J.Paul 
    1014    !> - Nov, 2013- Initial Version 
     1069   !> - November, 2013- Initial Version 
    10151070   ! 
    1016    !> @param[in] dd_value : 2D table of variable to be extrapolated 
    1017    !> @param[in] dd_fill : FillValue of variable 
    1018    !> @param[in] cd_dim : compute derivative on first (I) or second (J) dimension  
    1019    !------------------------------------------------------------------- 
    1020    !> @code 
     1071   !> @param[in] dd_value     2D array of variable to be extrapolated 
     1072   !> @param[in] dd_fill     FillValue of variable 
     1073   !> @param[in] cd_dim       compute derivative on first (I) or second (J) dimension  
     1074   !> @param[in] ld_discont   logical to take into account east west discontinuity  
     1075   !------------------------------------------------------------------- 
    10211076   FUNCTION extrap_deriv_2D( dd_value, dd_fill, cd_dim, ld_discont ) 
    10221077 
     
    11231178            END WHERE 
    11241179 
    1125       ENDDO 
     1180         ENDDO 
    11261181 
    11271182      CASE('J') 
     
    11871242 
    11881243   END FUNCTION extrap_deriv_2D 
    1189    !> @endcode 
    11901244   !------------------------------------------------------------------- 
    11911245   !> @brief 
    1192    !> This function compute derivative of a function in i- and 
    1193    !> j-direction 
     1246   !> This function compute derivative of 3D array. 
     1247   !> you have to specify in which direction derivative have to be computed: 
     1248   !> first (I), second (J) or third (K) dimension. 
    11941249   !>  
    11951250   !> @details  
     1251   !> optionaly you could specify to take into account east west discontinuity 
     1252   !> (-180° 180° or 0° 360° for longitude variable) 
    11961253   !> 
    11971254   !> @author J.Paul 
    1198    !> - Nov, 2013- Initial Version 
     1255   !> - November, 2013- Initial Version 
    11991256   ! 
    1200    !> @param[inout] dd_value : 3D table of variable to be extrapolated 
    1201    !> @param[in] dd_fill : FillValue of variable 
    1202    !> @param[in] cd_dim : compute derivative on first (I) second (J) or third (K) dimension    
    1203    !------------------------------------------------------------------- 
    1204    !> @code 
     1257   !> @param[inout] dd_value  3D array of variable to be extrapolated 
     1258   !> @param[in] dd_fill     FillValue of variable 
     1259   !> @param[in] cd_dim       compute derivative on first (I) second (J) or third (K) dimension    
     1260   !> @param[in] ld_discont   logical to take into account east west discontinuity 
     1261   !------------------------------------------------------------------- 
    12051262   PURE FUNCTION extrap_deriv_3D( dd_value, dd_fill, cd_dim, ld_discont ) 
    12061263 
     
    14311488 
    14321489   END FUNCTION extrap_deriv_3D 
    1433    !> @endcode 
    14341490   !------------------------------------------------------------------- 
    14351491   !> @brief 
    1436    !> This function compute extrapolatd values by calculated minimum error using 
    1437    !> taylor expansion 
     1492   !> This function compute coefficient for min_error extrapolation. 
    14381493   !>  
    14391494   !> @details  
     1495   !> coefficients are  "grid distance" to the center of the box choosed to compute 
     1496   !> extrapolation. 
    14401497   !> 
    14411498   !> @author J.Paul 
    1442    !> - Nov, 2013- Initial Version 
     1499   !> - November, 2013- Initial Version 
    14431500   ! 
    1444    !> @param[in] dd_value : 3D table of variable to be extrapolated 
    1445    !> @param[in] dd_fill : FillValue of variable 
    1446    !> @param[in] dd_ideriv : derivative of function in i-direction 
    1447    !> @param[in] dd_jderiv : derivative of function in j-direction 
    1448    !> @param[in] dd_kderiv : derivative of function in k-direction 
    1449    !> @param[in] id_ji  : i-direction indices table 
    1450    !> @param[in] id_jj  : j-direction indices table  
    1451    !> @param[in] id_ii  : i-direction indices of the point to extrapolate  
    1452    !> @param[in] id_ij  : j-direction indices of the point to extrapolate 
    1453    !------------------------------------------------------------------- 
    1454    !> @code 
     1501   !> @param[in] dd_value  3D array of variable to be extrapolated 
     1502   !> @return 3D array of coefficient for minimum error extrapolation 
     1503   !------------------------------------------------------------------- 
    14551504   PURE FUNCTION extrap__3D_min_error_coef( dd_value ) 
    14561505 
     
    15141563 
    15151564   END FUNCTION extrap__3D_min_error_coef 
    1516    !> @endcode 
    15171565   !------------------------------------------------------------------- 
    15181566   !> @brief 
    1519    !> This function compute extrapolatd values by calculated minimum error using 
     1567   !> This function compute extrapolatd value by calculated minimum error using 
    15201568   !> taylor expansion 
    15211569   !>  
    1522    !> @details  
     1570   !> @author J.Paul 
     1571   !> - November, 2013- Initial Version 
    15231572   !> 
    1524    !> @author J.Paul 
    1525    !> - Nov, 2013- Initial Version 
    1526    ! 
    1527    !> @param[in] dd_value : 3D table of variable to be extrapolated 
    1528    !> @param[in] dd_fill : FillValue of variable 
    1529    !> @param[in] dd_dfdx : derivative of function in i-direction 
    1530    !> @param[in] dd_dfdy : derivative of function in j-direction 
    1531    !> @param[in] dd_dfdz : derivative of function in k-direction 
    1532    !> @param[in] dd_coef :  
    1533    !------------------------------------------------------------------- 
    1534    !> @code 
    1535    PURE FUNCTION extrap__3D_min_error_fill( dd_value, dd_fill, id_ext, & 
     1573   !> @param[in] dd_value  3D array of variable to be extrapolated 
     1574   !> @param[in] dd_fill   FillValue of variable 
     1575   !> @param[in] id_radius radius of the halo used to compute extrapolation 
     1576   !> @param[in] dd_dfdx   derivative of function in i-direction 
     1577   !> @param[in] dd_dfdy   derivative of function in j-direction 
     1578   !> @param[in] dd_dfdz   derivative of function in k-direction 
     1579   !> @param[in] dd_coef   array of coefficient for min_error extrapolation  
     1580   !> @return extrapolatd value 
     1581   !------------------------------------------------------------------- 
     1582   PURE FUNCTION extrap__3D_min_error_fill( dd_value, dd_fill, id_radius, & 
    15361583   &                                        dd_dfdx, dd_dfdy, dd_dfdz, & 
    15371584   &                                        dd_coef ) 
     
    15401587      REAL(dp)   , DIMENSION(:,:,:), INTENT(IN) :: dd_value 
    15411588      REAL(dp)   ,                   INTENT(IN) :: dd_fill 
    1542       INTEGER(i4),                   INTENT(IN) :: id_ext 
     1589      INTEGER(i4),                   INTENT(IN) :: id_radius 
    15431590      REAL(dp)   , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdx 
    15441591      REAL(dp)   , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdy 
     
    15641611      extrap__3D_min_error_fill=dd_fill 
    15651612 
    1566       il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_ext*2)) 
     1613      il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2)) 
    15671614 
    15681615      IF( COUNT(dd_value(:,:,:) /= dd_fill) >= il_min )THEN 
     
    16021649 
    16031650   END FUNCTION extrap__3D_min_error_fill 
    1604    !> @endcode 
    16051651   !------------------------------------------------------------------- 
    16061652   !> @brief 
    1607    !> This function compute extrapolatd values by calculated minimum error using 
    1608    !> taylor expansion 
     1653   !> This function compute coefficient for inverse distance weighted method  
    16091654   !>  
    16101655   !> @details  
     1656   !> coefficients are inverse "grid distance" to the center of the box choosed to compute 
     1657   !> extrapolation. 
    16111658   !> 
    16121659   !> @author J.Paul 
    1613    !> - Nov, 2013- Initial Version 
     1660   !> - November, 2013- Initial Version 
    16141661   ! 
    1615    !> @param[in] dd_value : 3D table of variable to be extrapolated 
    1616    !------------------------------------------------------------------- 
    1617    !> @code 
     1662   !> @param[in] dd_value  3D array of variable to be extrapolated 
     1663   !> @return 3D array of coefficient for inverse distance weighted extrapolation 
     1664   !------------------------------------------------------------------- 
    16181665   PURE FUNCTION extrap__3D_dist_weight_coef( dd_value ) 
    16191666 
     
    16771724 
    16781725   END FUNCTION extrap__3D_dist_weight_coef 
    1679    !> @endcode 
    16801726   !------------------------------------------------------------------- 
    16811727   !> @brief 
    1682    !> This function compute extrapolatd values by calculated minimum error using 
    1683    !> taylor expansion 
     1728   !> This function compute extrapolatd value using inverse distance weighted 
     1729   !> method 
    16841730   !>  
    16851731   !> @details  
    16861732   !> 
    16871733   !> @author J.Paul 
    1688    !> - Nov, 2013- Initial Version 
     1734   !> - November, 2013- Initial Version 
    16891735   ! 
    1690    !> @param[in] dd_value : 3D table of variable to be extrapolated 
    1691    !> @param[in] dd_fill : FillValue of variable 
    1692    !> @param[in] dd_coef :   
    1693    !------------------------------------------------------------------- 
    1694    !> @code 
    1695    FUNCTION extrap__3D_dist_weight_fill( dd_value, dd_fill, id_ext, & 
    1696    &                                          dd_coef ) 
     1736   !> @param[in] dd_value  3D array of variable to be extrapolated 
     1737   !> @param[in] dd_fill   FillValue of variable 
     1738   !> @param[in] id_radius radius of the halo used to compute extrapolation 
     1739   !> @param[in] dd_coef   3D array of coefficient for inverse distance weighted extrapolation  
     1740   !> @return extrapolatd value 
     1741   !------------------------------------------------------------------- 
     1742   FUNCTION extrap__3D_dist_weight_fill( dd_value, dd_fill, id_radius, & 
     1743   &                                     dd_coef ) 
    16971744      IMPLICIT NONE 
    16981745      ! Argument 
    16991746      REAL(dp)   , DIMENSION(:,:,:), INTENT(IN) :: dd_value 
    17001747      REAL(dp)   ,                   INTENT(IN) :: dd_fill 
    1701       INTEGER(i4),                   INTENT(IN) :: id_ext 
     1748      INTEGER(i4),                   INTENT(IN) :: id_radius 
    17021749      REAL(dp)   , DIMENSION(:,:,:), INTENT(IN) :: dd_coef 
    17031750 
     
    17221769      extrap__3D_dist_weight_fill=dd_fill 
    17231770 
    1724       il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_ext*2)) 
     1771      il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2)) 
    17251772 
    17261773      IF( COUNT(dd_value(:,:,:)/= dd_fill) >= il_min )THEN 
     
    17331780         dl_coef(:,:,:)=0 
    17341781 
    1735          FORALL( ji=1:il_shape(1), & 
    1736          &       jj=1:il_shape(2), & 
    1737          &       jk=1:il_shape(3), & 
    1738          &       dd_value(ji,jj,jk) /= dd_fill ) 
    1739  
    1740             ! compute factor 
    1741             dl_value(ji,jj,jk)=dd_coef(ji,jj,jk)*dd_value(ji,jj,jk) 
    1742             dl_coef(ji,jj,jk)=dd_coef(ji,jj,jk) 
    1743  
    1744          END FORALL 
     1782         DO jk=1,il_shape(3) 
     1783            DO jj=1,il_shape(2) 
     1784               DO ji=1,il_shape(1) 
     1785 
     1786                  IF( dd_value(ji,jj,jk) /= dd_fill )THEN 
     1787                     ! compute factor 
     1788                     dl_value(ji,jj,jk)=dd_coef(ji,jj,jk)*dd_value(ji,jj,jk) 
     1789                     dl_coef(ji,jj,jk)=dd_coef(ji,jj,jk) 
     1790                  ENDIF 
     1791 
     1792               ENDDO 
     1793            ENDDO 
     1794         ENDDO 
    17451795 
    17461796         ! return value 
    17471797         IF( SUM( dl_coef(:,:,:) ) /= 0 )THEN 
    1748             extrap__3D_dist_weight_fill=SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) ) 
     1798            extrap__3D_dist_weight_fill = & 
     1799            &  SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) ) 
    17491800         ENDIF 
    17501801 
     
    17551806 
    17561807   END FUNCTION extrap__3D_dist_weight_fill 
    1757    !> @endcode 
    17581808   !------------------------------------------------------------------- 
    17591809   !> @brief 
     
    17611811   !> extraband of N points at north,south,east and west boundaries. 
    17621812   !>  
     1813   !> @details 
     1814   !> optionaly you could specify size of extra bands in i- and j-direction 
     1815   !> 
    17631816   !> @author J.Paul 
    1764    !> - 2013-Initial version 
     1817   !> - November, 2013-Initial version 
    17651818   ! 
    1766    !> @param[inout] td_var : variable  
    1767    !> @param[in] id_isize : i-direction size of extra bands (default=im_minext) 
    1768    !> @param[in] id_jsize : j-direction size of extra bands (default=im_minext) 
     1819   !> @param[inout] td_var variable  
     1820   !> @param[in] id_isize i-direction size of extra bands (default=im_minext) 
     1821   !> @param[in] id_jsize j-direction size of extra bands (default=im_minext) 
    17691822   !> @todo 
    17701823   !> - invalid special case for grid with north fold 
    17711824   !------------------------------------------------------------------- 
    1772    !> @code 
    17731825   SUBROUTINE extrap_add_extrabands(td_var, id_isize, id_jsize ) 
    17741826      IMPLICIT NONE 
     
    18211873      dl_value(:,:,:,:)=td_var%d_value(:,:,:,:) 
    18221874 
     1875 
    18231876      td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len + 2*il_isize 
    18241877      td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len + 2*il_jsize 
     
    18551908 
    18561909   END SUBROUTINE extrap_add_extrabands 
    1857    !> @endcode 
    18581910   !------------------------------------------------------------------- 
    18591911   !> @brief 
     
    18611913   !> of N points at north,south,east and west boundaries. 
    18621914   !>  
     1915   !> @details 
     1916   !> optionaly you could specify size of extra bands in i- and j-direction 
     1917   !> 
    18631918   !> @author J.Paul 
    1864    !> - 2013-Initial version 
    1865    ! 
    1866    !> @param[inout] td_var : variable  
    1867    !> @param[in] id_isize : i-direction size of extra bands (default=im_minext) 
    1868    !> @param[in] id_jsize : j-direction size of extra bands (default=im_minext) 
    1869    !------------------------------------------------------------------- 
    1870    !> @code 
     1919   !> - November, 2013-Initial version 
     1920   !> 
     1921   !> @param[inout] td_var variable  
     1922   !> @param[in] id_isize  i-direction size of extra bands (default=im_minext) 
     1923   !> @param[in] id_jsize  j-direction size of extra bands (default=im_minext) 
     1924   !------------------------------------------------------------------- 
    18711925   SUBROUTINE extrap_del_extrabands(td_var, id_isize, id_jsize ) 
    18721926      IMPLICIT NONE 
    18731927      ! Argument 
    1874       TYPE(TVAR) , INTENT(INOUT)  :: td_var 
     1928      TYPE(TVAR) , INTENT(INOUT) :: td_var 
    18751929      INTEGER(i4), INTENT(IN   ), OPTIONAL :: id_isize 
    18761930      INTEGER(i4), INTENT(IN   ), OPTIONAL :: id_jsize 
     
    18811935      INTEGER(i4) :: il_isize 
    18821936      INTEGER(i4) :: il_jsize 
    1883        
     1937  
    18841938      INTEGER(i4) :: il_imin 
    18851939      INTEGER(i4) :: il_imax 
     
    19351989 
    19361990   END SUBROUTINE extrap_del_extrabands 
    1937    !> @endcode 
    1938 !   !------------------------------------------------------------------- 
    1939 !   !> @brief 
    1940 !   !> This function 
    1941 !   !>  
    1942 !   !> @details  
    1943 !   !> 
    1944 !   !> @author J.Paul 
    1945 !   !> - Nov, 2013- Initial Version 
    1946 !   ! 
    1947 !   !> @param[in] 
    1948 !   !> @param[out]  
    1949 !   !------------------------------------------------------------------- 
    1950 !   !> @code 
    1951 !   FUNCTION extrap_( ) 
    1952 !      IMPLICIT NONE 
    1953 !      ! Argument 
    1954 ! 
    1955 !      ! local variable 
    1956 ! 
    1957 !      ! loop indices 
    1958 !      !---------------------------------------------------------------- 
    1959 !   END FUNCTION extrap_ 
    1960 !   !> @endcode 
    1961 !   !------------------------------------------------------------------- 
    1962 !   !> @brief 
    1963 !   !> This subroutine 
    1964 !   !>  
    1965 !   !> @details  
    1966 !   !> 
    1967 !   !> @author J.Paul 
    1968 !   !> - Nov, 2013- Initial Version 
    1969 !   ! 
    1970 !   !> @param[in] 
    1971 !   !> @param[out]  
    1972 !   !------------------------------------------------------------------- 
    1973 !   !> @code 
    1974 !   SUBROUTINE extrap_( ) 
    1975 !      IMPLICIT NONE 
    1976 !      ! Argument 
    1977 ! 
    1978 !      ! local variable 
    1979 ! 
    1980 !      ! loop indices 
    1981 !      !---------------------------------------------------------------- 
    1982 !   END SUBROUTINE extrap_ 
    1983 !   !> @endcode 
    19841991END MODULE extrap 
Note: See TracChangeset for help on using the changeset viewer.