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 5989 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/TOOLS/SIREN/src/extrap.f90 – NEMO

Ignore:
Timestamp:
2015-12-03T09:10:32+01:00 (8 years ago)
Author:
deazer
Message:

Merging TMB and 25h diagnostics to head of trunk
added brief documentation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/TOOLS/SIREN/src/extrap.f90

    r5260 r5989  
    1919!>    defining string character _cn\_varinfo_. By default _dist_weight_.<br/> 
    2020!>    Example: 
    21 !>       - cn_varinfo='varname1:dist_weight', 'varname2:min_error' 
     21!>       - cn_varinfo='varname1:ext=dist_weight', 'varname2:ext=min_error' 
    2222!> 
    2323!>    to detect point to be extrapolated:<br/> 
    2424!> @code 
    25 !>    il_detect(:,:,:)=extrap_detect(td_var, [td_level], [id_offset,] [id_rho,] [id_ext])  
     25!>    il_detect(:,:,:)=extrap_detect(td_var) 
    2626!> @endcode 
    2727!>       - il_detect(:,:,:) is 3D array of point to be extrapolated 
    2828!>       - 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] 
    3329!> 
    3430!>    to extrapolate variable:<br/> 
    3531!> @code 
    36 !>    CALL extrap_fill_value( td_var, [td_level], [id_offset], [id_rho], [id_iext], [id_jext], [id_kext], [id_radius], [id_maxiter]) 
     32!>    CALL extrap_fill_value( td_var, [id_radius]) 
    3733!> @endcode 
    3834!>       - 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] 
    4535!>       - id_radius is radius of the halo used to compute extrapolation [optional] 
    46 !>       - id_maxiter is maximum number of iteration [optional] 
    4736!> 
    4837!>    to add extraband to the variable (to be extrapolated):<br/> 
     
    6251!>       - id_jsize : j-direction size of extra bands [optional] 
    6352!> 
    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 !> 
    9053!> @warning _FillValue must not be zero (use var_chg_FillValue()) 
    9154!> 
     
    9356!> J.Paul 
    9457! REVISION HISTORY: 
    95 !> @date Nov, 2013 - Initial Version 
     58!> @date November, 2013 - Initial Version 
    9659!> @date September, 2014 
    9760!> - add header 
     61!> @date June, 2015 
     62!> - extrapolate all land points (_FillValue) 
     63!> - move deriv function to math module 
     64!> @date July, 2015 
     65!> - compute extrapolation from north west to south east,  
     66!> and from south east to north west 
    9867!> 
    9968!> @todo 
    10069!> - create module for each extrapolation method 
     70!> - smooth extrapolated points 
    10171!> 
    10272!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    11080   USE date                            ! date manager 
    11181   USE logger                          ! log file manager 
     82   USE math                            ! mathematical function 
    11283   USE att                             ! attribute manager 
    11384   USE dim                             ! dimension manager 
     
    11889 
    11990   ! type and variable 
    120    PRIVATE :: im_maxiter   !< default maximum number of iteration  
    12191   PRIVATE :: im_minext    !< default minumum number of point to extrapolate 
    12292   PRIVATE :: im_mincubic  !< default minumum number of point to extrapolate for cubic interpolation 
     
    12797   PUBLIC :: extrap_add_extrabands !< add extraband to the variable (to be extrapolated)  
    12898   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 
    13299 
    133100   PRIVATE :: extrap__detect_wrapper      ! detected point to be extrapolated wrapper 
     
    141108   PRIVATE :: extrap__3D_dist_weight_fill !  
    142109 
    143    INTEGER(i4), PARAMETER :: im_maxiter = 10 !< default maximum number of iteration 
    144110   INTEGER(i4), PARAMETER :: im_minext  = 2  !< default minumum number of point to extrapolate 
    145111   INTEGER(i4), PARAMETER :: im_mincubic= 4  !< default minumum number of point to extrapolate for cubic interpolation 
     
    171137   !>  
    172138   !> @author J.Paul 
    173    !> - November, 2013- Initial Version 
     139   !> @date November, 2013 - Initial Version 
     140   !> @date June, 2015 
     141   !> - do not use level to select points to be extrapolated 
    174142   ! 
    175143   !> @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 
    180144   !> @return array of point to be extrapolated 
    181145   !------------------------------------------------------------------- 
    182    FUNCTION extrap__detect( td_var0, td_level1, & 
    183    &                        id_offset, id_rho, id_ext ) 
     146   FUNCTION extrap__detect( td_var0 )  
    184147      IMPLICIT NONE 
    185148      ! Argument 
    186149      TYPE(TVAR) ,                 INTENT(IN   ) :: td_var0 
    187       TYPE(TVAR) , DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: td_level1 
    188       INTEGER(i4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset 
    189       INTEGER(i4), DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_rho 
    190       INTEGER(i4), DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_ext 
    191150 
    192151      ! function 
     
    196155 
    197156      ! local variable 
    198       CHARACTER(LEN=lc)                                :: cl_level 
    199  
    200       INTEGER(i4)                                      :: il_ind 
    201       INTEGER(i4)      , DIMENSION(:,:,:), ALLOCATABLE :: il_detect 
    202       INTEGER(i4)      , DIMENSION(:,:,:), ALLOCATABLE :: il_tmp 
    203       INTEGER(i4)      , DIMENSION(:,:)  , ALLOCATABLE :: il_offset 
    204       INTEGER(i4)      , DIMENSION(:,:)  , ALLOCATABLE :: il_level1 
    205       INTEGER(i4)      , DIMENSION(:,:)  , ALLOCATABLE :: il_level1_G0 
    206       INTEGER(i4)      , DIMENSION(:,:)  , ALLOCATABLE :: il_extra 
    207       INTEGER(i4)      , DIMENSION(:)    , ALLOCATABLE :: il_ext 
    208       INTEGER(i4)      , DIMENSION(:)    , ALLOCATABLE :: il_rho 
    209       INTEGER(i4)      , DIMENSION(:)    , ALLOCATABLE :: il_dim0 
    210  
    211       TYPE(TVAR)                                       :: tl_var1 
    212  
    213157      ! loop indices 
    214158      INTEGER(i4) :: ji0 
    215159      INTEGER(i4) :: jj0 
    216160      INTEGER(i4) :: jk0 
    217       INTEGER(i4) :: ji1 
    218       INTEGER(i4) :: jj1 
    219       INTEGER(i4) :: ji1m 
    220       INTEGER(i4) :: jj1m 
    221       INTEGER(i4) :: ji1p 
    222       INTEGER(i4) :: jj1p 
    223161      !---------------------------------------------------------------- 
    224162 
    225       ! init 
    226       extrap__detect(:,:,:)=0 
    227  
    228       ALLOCATE( il_dim0(3) ) 
    229       il_dim0(:)=td_var0%t_dim(1:3)%i_len 
    230  
    231       ! optional argument 
    232       ALLOCATE( il_rho(ip_maxdim) ) 
    233       il_rho(:)=1 
    234       IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:) 
    235  
    236       ALLOCATE( il_offset(ip_maxdim,2) ) 
    237       il_offset(:,:)=0 
    238       IF( PRESENT(id_offset) )THEN 
    239          il_offset(1:SIZE(id_offset(:,:),DIM=1),& 
    240          &         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) 
    244       ENDIF 
    245  
    246       ALLOCATE( il_ext(ip_maxdim) ) 
    247       il_ext(:)=im_minext 
    248       IF( PRESENT(id_ext) ) il_ext(1:SIZE(id_ext(:)))=id_ext(:) 
    249  
    250       ALLOCATE( il_detect(il_dim0(1),& 
    251       &                   il_dim0(2),& 
    252       &                   il_dim0(3)) ) 
    253       il_detect(:,:,:)=0 
    254  
    255       ! select point already inform 
    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 
    263   
    264       IF( PRESENT(td_level1) )THEN 
    265          SELECT CASE(TRIM(td_var0%c_point)) 
    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' 
    274          END SELECT 
    275  
    276          il_ind=var_get_index(td_level1(:),TRIM(cl_level)) 
    277          IF( il_ind == 0 )THEN 
    278             CALL logger_error("EXTRAP DETECT: can not compute point to be "//& 
    279             &     "extrapolated for variable "//TRIM(td_var0%c_name)//& 
    280             &      ". can not find "//& 
    281             &     "level for variable point "//TRIM(TRIM(td_var0%c_point))) 
    282          ELSE 
    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  
    344          ENDIF 
    345       ENDIF 
    346  
    347       ! clean 
    348       DEALLOCATE( il_offset ) 
    349   
    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 
    371        
    372       ! clean 
    373       DEALLOCATE( il_tmp ) 
     163      ! force to extrapolated all points 
     164      extrap__detect(:,:,:)=1 
    374165 
    375166      ! do not compute grid point already inform 
     
    377168         DO jj0=1,td_var0%t_dim(2)%i_len 
    378169            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 
     170               IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill )THEN 
     171                  extrap__detect(ji0,jj0,jk0)=0 
     172               ENDIF 
    380173            ENDDO 
    381174         ENDDO 
    382175      ENDDO 
    383  
    384       ! save result 
    385       extrap__detect(:,:,:)=il_detect(:,:,:) 
    386  
    387       ! clean 
    388       DEALLOCATE( il_dim0 ) 
    389       DEALLOCATE( il_ext ) 
    390       DEALLOCATE( il_detect ) 
    391       DEALLOCATE( il_rho ) 
    392176 
    393177   END FUNCTION extrap__detect 
     
    398182   !>  
    399183   !> @author J.Paul 
    400    !> - November, 2013- Initial Version 
     184   !> @date November, 2013 - Initial Version 
     185   !> @date June, 2015 
     186   !> - select all land points for extrapolation 
    401187   !> 
    402188   !> @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 
    407189   !> @return 3D array of point to be extrapolated 
    408190   !------------------------------------------------------------------- 
    409    FUNCTION extrap__detect_wrapper( td_var, td_level, & 
    410    &                                id_offset, id_rho, id_ext ) 
     191   FUNCTION extrap__detect_wrapper( td_var ) 
    411192 
    412193      IMPLICIT NONE 
    413194      ! Argument 
    414195      TYPE(TVAR) ,                 INTENT(IN   ) :: td_var 
    415       TYPE(TVAR) , DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: td_level 
    416       INTEGER(i4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset 
    417       INTEGER(i4), DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_rho 
    418       INTEGER(i4), DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_ext 
    419196 
    420197      ! function 
     
    439216         &              " for variable "//TRIM(td_var%c_name) ) 
    440217          
    441          extrap__detect_wrapper(:,:,:)=extrap__detect( td_var, td_level, & 
    442          &                                             id_offset, & 
    443          &                                             id_rho,    & 
    444          &                                             id_ext     ) 
     218         extrap__detect_wrapper(:,:,:)=extrap__detect( td_var ) 
    445219 
    446220      ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN 
     
    450224         &              " for variable "//TRIM(td_var%c_name) ) 
    451225          
    452          extrap__detect_wrapper(:,:,1:1)=extrap__detect( td_var , td_level,& 
    453          &                                               id_offset, & 
    454          &                                               id_rho,    & 
    455          &                                               id_ext     ) 
     226         extrap__detect_wrapper(:,:,1:1)=extrap__detect( td_var ) 
    456227 
    457228      ELSE IF( td_var%t_dim(3)%l_use )THEN 
     
    461232         &              " for variable "//TRIM(td_var%c_name) ) 
    462233          
    463          extrap__detect_wrapper(1:1,1:1,:)=extrap__detect( td_var , td_level, & 
    464          &                                                 id_offset, & 
    465          &                                                 id_rho,    & 
    466          &                                                 id_ext     ) 
     234         extrap__detect_wrapper(1:1,1:1,:)=extrap__detect( td_var ) 
    467235 
    468236      ENDIF               
     
    489257   !> 
    490258   !> @author J.Paul 
    491    !> - Nov, 2013- Initial Version 
     259   !> @date November, 2013 - Initial Version 
     260   !> @date June, 2015 
     261   !> - select all land points for extrapolation 
    492262   ! 
    493263   !> @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 
    500264   !> @param[in] id_radius    radius of the halo used to compute extrapolation  
    501    !> @param[in] id_maxiter   maximum number of iteration 
    502    !------------------------------------------------------------------- 
    503    SUBROUTINE extrap__fill_value_wrapper( td_var, td_level, & 
    504    &                                      id_offset,        & 
    505    &                                      id_rho,           & 
    506    &                                      id_iext, id_jext, id_kext, & 
    507    &                                      id_radius, id_maxiter ) 
     265   !------------------------------------------------------------------- 
     266   SUBROUTINE extrap__fill_value_wrapper( td_var, &  
     267   &                                      id_radius ) 
    508268      IMPLICIT NONE 
    509269      ! Argument 
    510270      TYPE(TVAR) ,                  INTENT(INOUT) :: td_var 
    511       TYPE(TVAR) ,  DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: td_level 
    512       INTEGER(i4),  DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset 
    513       INTEGER(i4),  DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_rho 
    514       INTEGER(i4),                  INTENT(IN   ), OPTIONAL :: id_iext 
    515       INTEGER(i4),                  INTENT(IN   ), OPTIONAL :: id_jext 
    516       INTEGER(i4),                  INTENT(IN   ), OPTIONAL :: id_kext 
    517271      INTEGER(i4),                  INTENT(IN   ), OPTIONAL :: id_radius 
    518       INTEGER(i4),                  INTENT(IN   ), OPTIONAL :: id_maxiter 
    519272 
    520273      ! local variable 
    521       INTEGER(i4) :: il_iext 
    522       INTEGER(i4) :: il_jext 
    523       INTEGER(i4) :: il_kext 
    524274      INTEGER(i4) :: il_radius 
    525       INTEGER(i4) :: il_maxiter 
    526275 
    527276      CHARACTER(LEN=lc) :: cl_method 
     
    544293         END SELECT 
    545294 
    546          il_iext=im_minext 
    547          IF( PRESENT(id_iext) ) il_iext=id_iext 
    548          il_jext=im_minext 
    549          IF( PRESENT(id_jext) ) il_jext=id_jext 
    550          il_kext=0 
    551          IF( PRESENT(id_kext) ) il_kext=id_kext 
    552  
    553          IF( TRIM(td_var%c_interp(1)) == 'cubic')THEN 
    554             IF( il_iext > 0 .AND. il_iext < im_mincubic ) il_iext=im_mincubic 
    555             IF( il_jext > 0 .AND. il_jext < im_mincubic ) il_jext=im_mincubic 
     295         ! number of point use to compute box 
     296         il_radius=1 
     297         IF( PRESENT(id_radius) ) il_radius=id_radius 
     298         IF( il_radius < 0 )THEN 
     299            CALL logger_error("EXTRAP FILL VALUE: invalid "//& 
     300            &  " radius of the box used to compute extrapolation "//& 
     301            &  "("//TRIM(fct_str(il_radius))//")") 
    556302         ENDIF 
    557303 
    558          IF( il_iext < 0 )THEN 
    559             CALL logger_error("EXTRAP FILL VALUE: invalid "//& 
    560             &  " number of points to be extrapolated in i-direction "//& 
    561             &  "("//TRIM(fct_str(il_iext))//")") 
    562          ENDIF 
    563  
    564          IF( il_jext < 0 )THEN 
    565             CALL logger_error("EXTRAP FILL VALUE: invalid "//& 
    566             &  " number of points to be extrapolated in j-direction "//& 
    567             &  "("//TRIM(fct_str(il_jext))//")") 
    568          ENDIF 
    569  
    570          IF( il_kext < 0 )THEN 
    571             CALL logger_error("EXTRAP FILL VALUE: invalid "//& 
    572             &  " number of points to be extrapolated in k-direction "//& 
    573             &  "("//TRIM(fct_str(il_kext))//")") 
    574          ENDIF 
    575  
    576          IF( (il_iext /= 0 .AND. td_var%t_dim(1)%l_use) .OR. & 
    577          &   (il_jext /= 0 .AND. td_var%t_dim(2)%l_use) .OR. & 
    578          &   (il_kext /= 0 .AND. td_var%t_dim(3)%l_use) )THEN 
    579  
    580             ! number of point use to compute box 
    581             il_radius=1 
    582             IF( PRESENT(id_radius) ) il_radius=id_radius 
    583             IF( il_radius < 0 )THEN 
    584                CALL logger_error("EXTRAP FILL VALUE: invalid "//& 
    585                &  " radius of the box used to compute extrapolation "//& 
    586                &  "("//TRIM(fct_str(il_radius))//")") 
    587             ENDIF 
    588  
    589             ! maximum number of iteration 
    590             il_maxiter=im_maxiter 
    591             IF( PRESENT(id_maxiter) ) il_maxiter=id_maxiter 
    592             IF( il_maxiter < 0 )THEN 
    593                CALL logger_error("EXTRAP FILL VALUE: invalid "//& 
    594                &  " maximum nuber of iteration "//& 
    595                &  "("//TRIM(fct_str(il_maxiter))//")") 
    596             ENDIF 
    597  
    598             CALL logger_info("EXTRAP FILL: extrapolate "//TRIM(td_var%c_name)//& 
    599             &  " using "//TRIM(cl_method)//" method." ) 
    600  
    601             CALL extrap__fill_value( td_var, cl_method, & 
    602             &                        il_iext, il_jext, il_kext,   & 
    603             &                        il_radius, il_maxiter,       & 
    604             &                        td_level,                    & 
    605             &                        id_offset, id_rho ) 
    606   
    607          ENDIF 
     304         CALL logger_info("EXTRAP FILL: extrapolate "//TRIM(td_var%c_name)//& 
     305         &  " using "//TRIM(cl_method)//" method." ) 
     306 
     307         CALL extrap__fill_value( td_var, cl_method, & 
     308         &                        il_radius ) 
    608309  
    609310      ENDIF 
     
    621322   !> 
    622323   !> @author J.Paul 
    623    !> - November, 2013- Initial Version 
     324   !> @date November, 2013 - Initial Version 
     325   !> @date June, 2015 
     326   !> - select all land points for extrapolation 
    624327   ! 
    625328   !> @param[inout] td_var    variable structure 
    626329   !> @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 
    630330   !> @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 
    635331   !------------------------------------------------------------------- 
    636332   SUBROUTINE extrap__fill_value( td_var, cd_method, & 
    637    &                              id_iext, id_jext, id_kext, & 
    638    &                              id_radius, id_maxiter, & 
    639    &                              td_level,          & 
    640    &                              id_offset,         & 
    641    &                              id_rho ) 
     333   &                              id_radius ) 
    642334      IMPLICIT NONE 
    643335      ! Argument 
    644336      TYPE(TVAR)      ,                 INTENT(INOUT) :: td_var 
    645337      CHARACTER(LEN=*),                 INTENT(IN   ) :: cd_method 
    646       INTEGER(i4)     ,                 INTENT(IN   ) :: id_iext 
    647       INTEGER(i4)     ,                 INTENT(IN   ) :: id_jext 
    648       INTEGER(i4)     ,                 INTENT(IN   ) :: id_kext 
    649338      INTEGER(i4)     ,                 INTENT(IN   ) :: id_radius 
    650       INTEGER(i4)     ,                 INTENT(IN   ) :: id_maxiter 
    651       TYPE(TVAR)      , DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: td_level 
    652       INTEGER(i4)     , DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset 
    653       INTEGER(i4)     , DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_rho 
    654339 
    655340      ! local variable 
     
    668353      &                    td_var%t_dim(3)%i_len) ) 
    669354 
    670       il_detect(:,:,:) = extrap_detect( td_var, td_level, & 
    671       &                                 id_offset,        & 
    672       &                                 id_rho,           & 
    673       &                                 id_ext=(/id_iext, id_jext, id_kext/) ) 
     355      il_detect(:,:,:) = extrap_detect( td_var ) 
     356 
    674357      !2- add attribute to variable 
    675358      cl_extrap=fct_concat(td_var%c_extrap(:)) 
     
    679362      CALL att_clean(tl_att) 
    680363 
    681       CALL logger_info(" EXTRAP FILL: "//& 
    682          &              TRIM(fct_str(SUM(il_detect(:,:,:))))//& 
    683          &              " point(s) to extrapolate " ) 
    684  
    685       !3- extrapolate 
    686       CALL extrap__3D(td_var%d_value(:,:,:,:), td_var%d_fill,    & 
    687       &               il_detect(:,:,:),                           & 
    688       &               cd_method, id_radius, id_maxiter  ) 
     364      IF( ALL(il_detect(:,:,:)==1) )THEN 
     365         CALL logger_warn(" EXTRAP FILL: "//& 
     366            &  " can not extrapolate "//TRIM(td_var%c_name)//& 
     367            &  ". no value inform." ) 
     368      ELSE 
     369         CALL logger_info(" EXTRAP FILL: "//& 
     370            &              TRIM(fct_str(SUM(il_detect(:,:,:))))//& 
     371            &              " point(s) to extrapolate " ) 
     372 
     373         CALL logger_info(" EXTRAP FILL: method "//& 
     374            &  TRIM(cd_method) ) 
     375 
     376         !3- extrapolate 
     377         CALL extrap__3D(td_var%d_value(:,:,:,:), td_var%d_fill, & 
     378         &               il_detect(:,:,:),                       & 
     379         &               cd_method, id_radius ) 
     380      ENDIF 
    689381 
    690382      DEALLOCATE(il_detect) 
     
    705397   !> 
    706398   !> @author J.Paul 
    707    !> - Nov, 2013- Initial Version 
     399   !> @date November, 2013 - Initial Version 
     400   !> @date July, 2015  
     401   !> - compute coef indices to be used 
     402   !> - bug fix: force coef indice to 1, for dimension lenth equal to 1 
    708403   ! 
    709404   !> @param[inout] dd_value  3D array of variable to be extrapolated 
     
    714409   !------------------------------------------------------------------- 
    715410   SUBROUTINE extrap__3D( dd_value, dd_fill, id_detect,& 
    716    &                      cd_method, id_radius, id_maxiter ) 
     411   &                      cd_method, id_radius ) 
    717412      IMPLICIT NONE 
    718413      ! Argument 
    719414      REAL(dp)   , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_value 
    720       REAL(dp)   ,                   INTENT(IN   ) :: dd_fill 
    721       INTEGER(i4), DIMENSION(:,:,:), INTENT(INOUT) :: id_detect 
    722       CHARACTER(LEN=*),              INTENT(IN   ) :: cd_method 
    723       INTEGER(i4),                   INTENT(IN   ) :: id_radius 
    724       INTEGER(i4),                   INTENT(IN   ) :: id_maxiter 
     415      REAL(dp)   ,                     INTENT(IN   ) :: dd_fill 
     416      INTEGER(i4), DIMENSION(:,:,:)  , INTENT(INOUT) :: id_detect 
     417      CHARACTER(LEN=*),                INTENT(IN   ) :: cd_method 
     418      INTEGER(i4),                     INTENT(IN   ) :: id_radius 
    725419 
    726420      ! local variable 
    727       INTEGER(i4) :: il_imin 
    728       INTEGER(i4) :: il_imax 
    729       INTEGER(i4) :: il_jmin 
    730       INTEGER(i4) :: il_jmax 
    731       INTEGER(i4) :: il_kmin 
    732       INTEGER(i4) :: il_kmax 
    733       INTEGER(i4) :: il_iter 
    734       INTEGER(i4) :: il_radius 
    735  
    736       INTEGER(i4), DIMENSION(4) :: il_shape 
    737       INTEGER(i4), DIMENSION(3) :: il_dim 
     421      INTEGER(i4)                                :: il_imin 
     422      INTEGER(i4)                                :: il_imax 
     423      INTEGER(i4)                                :: il_jmin 
     424      INTEGER(i4)                                :: il_jmax 
     425      INTEGER(i4)                                :: il_kmin 
     426      INTEGER(i4)                                :: il_kmax 
     427      INTEGER(i4)                                :: il_iter 
     428      INTEGER(i4)                                :: il_radius 
     429      INTEGER(i4)                                :: il_i1 
     430      INTEGER(i4)                                :: il_i2 
     431      INTEGER(i4)                                :: il_j1 
     432      INTEGER(i4)                                :: il_j2 
     433      INTEGER(i4)                                :: il_k1 
     434      INTEGER(i4)                                :: il_k2 
     435 
     436      INTEGER(i4), DIMENSION(4)                  :: il_shape 
     437      INTEGER(i4), DIMENSION(3)                  :: il_dim 
    738438 
    739439      INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect 
     
    743443      REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdz 
    744444      REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef 
     445 
     446      LOGICAL                                    :: ll_iter 
    745447 
    746448      ! loop indices 
     
    765467            DO WHILE( ANY(il_detect(:,:,:)==1) ) 
    766468               ! change extend value to minimize number of iteration 
    767                il_radius=id_radius+(il_iter/id_maxiter) 
     469               il_radius=id_radius+(il_iter-1) 
     470               ll_iter=.TRUE. 
    768471 
    769472               ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) )  
     
    774477               dl_dfdx(:,:,:)=dd_fill 
    775478               IF( il_shape(1) > 1 )THEN 
    776                   dl_dfdx(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'I' ) 
     479                  dl_dfdx(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 
     480                     &                          dd_fill, 'I' ) 
    777481               ENDIF 
    778482 
     
    780484               dl_dfdy(:,:,:)=dd_fill 
    781485               IF( il_shape(2) > 1 )THEN 
    782                   dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'J' ) 
     486                  dl_dfdy(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 
     487                     &                          dd_fill, 'J' ) 
    783488               ENDIF 
    784489  
     
    786491               dl_dfdz(:,:,:)=dd_fill 
    787492               IF( il_shape(3) > 1 )THEN 
    788                   dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'K' ) 
     493                  dl_dfdz(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), & 
     494                     &                          dd_fill, 'K' ) 
    789495               ENDIF 
    790496  
     
    804510 
    805511               DO jk=1,il_shape(3) 
     512                  ! from North West(1,1) to South East(il_shape(1),il_shape(2)) 
    806513                  IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 
    807514                  DO jj=1,il_shape(2) 
     
    813520                           il_imin=MAX(ji-il_radius,1) 
    814521                           il_imax=MIN(ji+il_radius,il_shape(1)) 
     522                           ! coef indices to be used 
     523                           il_i1 = il_radius-(ji-il_imin)+1 
     524                           il_i2 = il_radius+(il_imax-ji)+1 
    815525                           IF( il_dim(1) == 1 )THEN 
    816526                              il_imin=ji 
    817527                              il_imax=ji 
    818                            ENDIF 
     528                              ! coef indices to be used 
     529                              il_i1 = 1 
     530                              il_i2 = 1 
     531                           ENDIF 
     532 
    819533 
    820534                           il_jmin=MAX(jj-il_radius,1) 
    821535                           il_jmax=MIN(jj+il_radius,il_shape(2)) 
     536                           ! coef indices to be used 
     537                           il_j1 = il_radius-(jj-il_jmin)+1 
     538                           il_j2 = il_radius+(il_jmax-jj)+1 
    822539                           IF( il_dim(2) == 1 )THEN 
    823540                              il_jmin=jj 
    824541                              il_jmax=jj 
     542                              ! coef indices to be used 
     543                              il_j1 = 1 
     544                              il_j2 = 1 
    825545                           ENDIF 
    826546 
    827547                           il_kmin=MAX(jk-il_radius,1) 
    828548                           il_kmax=MIN(jk+il_radius,il_shape(3)) 
     549                           ! coef indices to be used 
     550                           il_k1 = il_radius-(jk-il_kmin)+1 
     551                           il_k2 = il_radius+(il_kmax-jk)+1 
    829552                           IF( il_dim(3) == 1 )THEN 
    830553                              il_kmin=jk 
    831554                              il_kmax=jk 
     555                              ! coef indices to be used 
     556                              il_k1 = 1 
     557                              il_k2 = 1 
    832558                           ENDIF 
    833559 
     
    845571                           &            il_jmin:il_jmax, & 
    846572                           &            il_kmin:il_kmax ), & 
    847                            &  dl_coef(:,:,:) ) 
     573                           &  dl_coef(il_i1:il_i2, & 
     574                           &          il_j1:il_j2, & 
     575                           &          il_k1:il_k2) ) 
    848576 
    849577                           IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 
    850578                              il_detect(ji,jj,jk)= 0 
     579                              ll_iter=.FALSE. 
     580                           ENDIF 
     581 
     582                        ENDIF 
     583 
     584                     ENDDO 
     585                  ENDDO 
     586                  ! from South East(il_shape(1),il_shape(2)) to North West(1,1) 
     587                  IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 
     588                  DO jj=il_shape(2),1,-1 
     589                     IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 
     590                     DO ji=il_shape(1),1,-1 
     591 
     592                        IF( il_detect(ji,jj,jk) == 1 )THEN 
     593                           
     594                           il_imin=MAX(ji-il_radius,1) 
     595                           il_imax=MIN(ji+il_radius,il_shape(1)) 
     596                           ! coef indices to be used 
     597                           il_i1 = il_radius-(ji-il_imin)+1 
     598                           il_i2 = il_radius+(il_imax-ji)+1 
     599                           IF( il_dim(1) == 1 )THEN 
     600                              il_imin=ji 
     601                              il_imax=ji 
     602                              ! coef indices to be used 
     603                              il_i1 = 1 
     604                              il_i2 = 1 
     605                           ENDIF 
     606 
     607 
     608                           il_jmin=MAX(jj-il_radius,1) 
     609                           il_jmax=MIN(jj+il_radius,il_shape(2)) 
     610                           ! coef indices to be used 
     611                           il_j1 = il_radius-(jj-il_jmin)+1 
     612                           il_j2 = il_radius+(il_jmax-jj)+1 
     613                           IF( il_dim(2) == 1 )THEN 
     614                              il_jmin=jj 
     615                              il_jmax=jj 
     616                              ! coef indices to be used 
     617                              il_j1 = 1 
     618                              il_j2 = 1 
     619                           ENDIF 
     620 
     621                           il_kmin=MAX(jk-il_radius,1) 
     622                           il_kmax=MIN(jk+il_radius,il_shape(3)) 
     623                           ! coef indices to be used 
     624                           il_k1 = il_radius-(jk-il_kmin)+1 
     625                           il_k2 = il_radius+(il_kmax-jk)+1 
     626                           IF( il_dim(3) == 1 )THEN 
     627                              il_kmin=jk 
     628                              il_kmax=jk 
     629                              ! coef indices to be used 
     630                              il_k1 = 1 
     631                              il_k2 = 1 
     632                           ENDIF 
     633 
     634                           dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( & 
     635                           &  dd_value( il_imin:il_imax, & 
     636                           &            il_jmin:il_jmax, & 
     637                           &            il_kmin:il_kmax,jl ), dd_fill, il_radius, & 
     638                           &  dl_dfdx(  il_imin:il_imax, & 
     639                           &            il_jmin:il_jmax, & 
     640                           &            il_kmin:il_kmax ), & 
     641                           &  dl_dfdy(  il_imin:il_imax, & 
     642                           &            il_jmin:il_jmax, & 
     643                           &            il_kmin:il_kmax ), & 
     644                           &  dl_dfdz(  il_imin:il_imax, & 
     645                           &            il_jmin:il_jmax, & 
     646                           &            il_kmin:il_kmax ), & 
     647                           &  dl_coef(il_i1:il_i2, & 
     648                           &          il_j1:il_j2, & 
     649                           &          il_k1:il_k2) ) 
     650 
     651                           IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 
     652                              il_detect(ji,jj,jk)= 0 
     653                              ll_iter=.FALSE. 
    851654                           ENDIF 
    852655 
     
    862665               DEALLOCATE( dl_coef ) 
    863666 
    864                il_iter=il_iter+1 
     667               IF( ll_iter ) il_iter=il_iter+1 
    865668            ENDDO 
    866669         ENDDO 
     
    875678            DO WHILE( ANY(il_detect(:,:,:)==1) ) 
    876679               ! change extend value to minimize number of iteration 
    877                il_radius=id_radius+(il_iter/id_maxiter) 
     680               il_radius=id_radius+(il_iter-1) 
     681               ll_iter=.TRUE. 
    878682 
    879683               il_dim(1)=2*il_radius+1 
     
    886690               ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) 
    887691 
    888                dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), & 
    889                &                                                   1:il_dim(2), & 
    890                &                                                   1:il_dim(3), & 
     692               dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1),& 
     693               &                                                   1:il_dim(2),& 
     694               &                                                   1:il_dim(3),& 
    891695               &                                                   jl ) ) 
    892  
     696                
    893697               DO jk=1,il_shape(3) 
     698                  ! from North West(1,1) to South East(il_shape(1),il_shape(2)) 
    894699                  IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 
    895700                  DO jj=1,il_shape(2) 
     
    901706                           il_imin=MAX(ji-il_radius,1) 
    902707                           il_imax=MIN(ji+il_radius,il_shape(1)) 
     708                           ! coef indices to be used 
     709                           il_i1 = il_radius-(ji-il_imin)+1 
     710                           il_i2 = il_radius+(il_imax-ji)+1 
    903711                           IF( il_dim(1) == 1 )THEN 
    904712                              il_imin=ji 
    905713                              il_imax=ji 
     714                              ! coef indices to be used 
     715                              il_i1 = 1 
     716                              il_i2 = 1 
    906717                           ENDIF 
    907718 
    908719                           il_jmin=MAX(jj-il_radius,1) 
    909720                           il_jmax=MIN(jj+il_radius,il_shape(2)) 
     721                           ! coef indices to be used 
     722                           il_j1 = il_radius-(jj-il_jmin)+1 
     723                           il_j2 = il_radius+(il_jmax-jj)+1 
    910724                           IF( il_dim(2) == 1 )THEN 
    911725                              il_jmin=jj 
    912726                              il_jmax=jj 
     727                              ! coef indices to be used 
     728                              il_j1 = 1 
     729                              il_j2 = 1 
    913730                           ENDIF 
    914731 
    915732                           il_kmin=MAX(jk-il_radius,1) 
    916733                           il_kmax=MIN(jk+il_radius,il_shape(3)) 
     734                           ! coef indices to be used 
     735                           il_k1 = il_radius-(jk-il_kmin)+1 
     736                           il_k2 = il_radius+(il_kmax-jk)+1 
    917737                           IF( il_dim(3) == 1 )THEN 
    918738                              il_kmin=jk 
    919739                              il_kmax=jk 
     740                              ! coef indices to be used 
     741                              il_k1 = 1 
     742                              il_k2 = 1 
    920743                           ENDIF 
    921744 
     
    925748                           &            il_kmin:il_kmax, & 
    926749                           &            jl), dd_fill, il_radius, & 
    927                            &  dl_coef(:,:,:) ) 
     750                           &  dl_coef(il_i1:il_i2, & 
     751                           &          il_j1:il_j2, & 
     752                           &          il_k1:il_k2) ) 
    928753 
    929754                           IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 
    930755                              il_detect(ji,jj,jk)= 0 
     756                              ll_iter=.FALSE. 
     757                           ENDIF 
     758 
     759                        ENDIF 
     760 
     761                     ENDDO 
     762                  ENDDO 
     763                  ! from South East(il_shape(1),il_shape(2)) to North West(1,1) 
     764                  IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE 
     765                  DO jj=il_shape(2),1,-1 
     766                     IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE 
     767                     DO ji=il_shape(1),1,-1 
     768 
     769                        IF( il_detect(ji,jj,jk) == 1 )THEN 
     770                            
     771                           il_imin=MAX(ji-il_radius,1) 
     772                           il_imax=MIN(ji+il_radius,il_shape(1)) 
     773                           ! coef indices to be used 
     774                           il_i1 = il_radius-(ji-il_imin)+1 
     775                           il_i2 = il_radius+(il_imax-ji)+1 
     776                           IF( il_dim(1) == 1 )THEN 
     777                              il_imin=ji 
     778                              il_imax=ji 
     779                              ! coef indices to be used 
     780                              il_i1 = 1 
     781                              il_i2 = 1 
     782                           ENDIF 
     783 
     784                           il_jmin=MAX(jj-il_radius,1) 
     785                           il_jmax=MIN(jj+il_radius,il_shape(2)) 
     786                           ! coef indices to be used 
     787                           il_j1 = il_radius-(jj-il_jmin)+1 
     788                           il_j2 = il_radius+(il_jmax-jj)+1 
     789                           IF( il_dim(2) == 1 )THEN 
     790                              il_jmin=jj 
     791                              il_jmax=jj 
     792                              ! coef indices to be used 
     793                              il_j1 = 1 
     794                              il_j2 = 1 
     795                           ENDIF 
     796 
     797                           il_kmin=MAX(jk-il_radius,1) 
     798                           il_kmax=MIN(jk+il_radius,il_shape(3)) 
     799                           ! coef indices to be used 
     800                           il_k1 = il_radius-(jk-il_kmin)+1 
     801                           il_k2 = il_radius+(il_kmax-jk)+1 
     802                           IF( il_dim(3) == 1 )THEN 
     803                              il_kmin=jk 
     804                              il_kmax=jk 
     805                              ! coef indices to be used 
     806                              il_k1 = 1 
     807                              il_k2 = 1 
     808                           ENDIF 
     809 
     810                           dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( & 
     811                           &  dd_value( il_imin:il_imax, & 
     812                           &            il_jmin:il_jmax, & 
     813                           &            il_kmin:il_kmax, & 
     814                           &            jl), dd_fill, il_radius, & 
     815                           &  dl_coef(il_i1:il_i2, & 
     816                           &          il_j1:il_j2, & 
     817                           &          il_k1:il_k2) ) 
     818 
     819                           IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN 
     820                              il_detect(ji,jj,jk)= 0 
     821                              ll_iter=.FALSE. 
    931822                           ENDIF 
    932823 
     
    936827                  ENDDO 
    937828               ENDDO 
    938  
     829               CALL logger_info(" EXTRAP 3D: "//& 
     830               &              TRIM(fct_str(SUM(il_detect(:,:,:))))//& 
     831               &              " point(s) to extrapolate " ) 
     832             
    939833               DEALLOCATE( dl_coef ) 
    940                il_iter=il_iter+1 
     834               IF( ll_iter ) il_iter=il_iter+1 
    941835            ENDDO 
    942836         ENDDO             
     
    946840 
    947841   END SUBROUTINE extrap__3D 
    948    !------------------------------------------------------------------- 
    949    !> @brief 
    950    !> This function compute derivative of 1D array. 
    951    !>  
    952    !> @details  
    953    !> optionaly you could specify to take into account east west discontinuity 
    954    !> (-180° 180° or 0° 360° for longitude variable) 
    955    !> 
    956    !> @author J.Paul 
    957    !> - November, 2013- Initial Version 
    958    ! 
    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    !------------------------------------------------------------------- 
    963    PURE FUNCTION extrap_deriv_1D( dd_value, dd_fill, ld_discont ) 
    964  
    965       IMPLICIT NONE 
    966       ! Argument 
    967       REAL(dp)   , DIMENSION(:), INTENT(IN) :: dd_value 
    968       REAL(dp)                 , INTENT(IN) :: dd_fill 
    969       LOGICAL                  , INTENT(IN), OPTIONAL :: ld_discont 
    970  
    971       ! function 
    972       REAL(dp), DIMENSION(SIZE(dd_value,DIM=1) ) :: extrap_deriv_1D 
    973  
    974       ! local variable 
    975       INTEGER(i4)                            :: il_imin 
    976       INTEGER(i4)                            :: il_imax 
    977       INTEGER(i4), DIMENSION(1)              :: il_shape 
    978  
    979       REAL(dp)                               :: dl_min 
    980       REAL(dp)                               :: dl_max 
    981       REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_value 
    982  
    983       LOGICAL                                :: ll_discont 
    984  
    985       ! loop indices 
    986       INTEGER(i4) :: ji 
    987  
    988       INTEGER(i4) :: i1 
    989       INTEGER(i4) :: i2 
    990       !---------------------------------------------------------------- 
    991       ! init 
    992       extrap_deriv_1D(:)=dd_fill 
    993  
    994       ll_discont=.FALSE. 
    995       IF( PRESENT(ld_discont) ) ll_discont=ld_discont 
    996  
    997       il_shape(:)=SHAPE(dd_value(:)) 
    998  
    999       ALLOCATE( dl_value(3)) 
    1000  
    1001       ! compute derivative in i-direction 
    1002       DO ji=1,il_shape(1) 
    1003           
    1004             il_imin=MAX(ji-1,1) 
    1005             il_imax=MIN(ji+1,il_shape(1)) 
    1006  
    1007             IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN 
    1008                i1=1  ; i2=3 
    1009             ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN 
    1010                i1=1  ; i2=2 
    1011             ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN 
    1012                i1=2  ; i2=3 
    1013             ENDIF 
    1014  
    1015             dl_value(i1:i2)=dd_value(il_imin:il_imax) 
    1016             IF( il_imin == 1 )THEN 
    1017                dl_value(:)=EOSHIFT( dl_value(:), & 
    1018                &                    DIM=1,         & 
    1019                &                    SHIFT=-1,      & 
    1020                &                    BOUNDARY=dl_value(1) ) 
    1021             ENDIF 
    1022             IF( il_imax == il_shape(1) )THEN 
    1023                dl_value(:)=EOSHIFT( dl_value(:), & 
    1024                &                    DIM=1,         & 
    1025                &                    SHIFT=1,       & 
    1026                &                    BOUNDARY=dl_value(3)) 
    1027             ENDIF 
    1028  
    1029             IF( ll_discont )THEN 
    1030                dl_min=MINVAL( dl_value(:), dl_value(:)/=dd_fill ) 
    1031                dl_max=MAXVAL( dl_value(:), dl_value(:)/=dd_fill ) 
    1032                IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 
    1033                   WHERE( dl_value(:) < 0._dp )  
    1034                      dl_value(:) = dl_value(:)+360._dp 
    1035                   END WHERE 
    1036                ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 
    1037                   WHERE( dl_value(:) > 180._dp )  
    1038                      dl_value(:) = dl_value(:)-180._dp 
    1039                   END WHERE 
    1040                ENDIF 
    1041             ENDIF 
    1042  
    1043          IF( dl_value( 2) /= dd_fill .AND. & ! ji 
    1044          &   dl_value( 3) /= dd_fill .AND. & ! ji+1 
    1045          &   dl_value( 1) /= dd_fill )THEN   ! ji-1 
    1046  
    1047             extrap_deriv_1D(ji)=& 
    1048             &  ( dl_value(3) - dl_value(1) ) / & 
    1049             &  REAL( il_imax-il_imin ,dp) 
    1050  
    1051          ENDIF 
    1052  
    1053       ENDDO 
    1054  
    1055       DEALLOCATE( dl_value ) 
    1056  
    1057    END FUNCTION extrap_deriv_1D 
    1058    !------------------------------------------------------------------- 
    1059    !> @brief 
    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    !> 
    1064    !> @details  
    1065    !> optionaly you could specify to take into account east west discontinuity 
    1066    !> (-180° 180° or 0° 360° for longitude variable) 
    1067    !> 
    1068    !> @author J.Paul 
    1069    !> - November, 2013- Initial Version 
    1070    ! 
    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    !------------------------------------------------------------------- 
    1076    FUNCTION extrap_deriv_2D( dd_value, dd_fill, cd_dim, ld_discont ) 
    1077  
    1078       IMPLICIT NONE 
    1079       ! Argument 
    1080       REAL(dp)   , DIMENSION(:,:), INTENT(IN) :: dd_value 
    1081       REAL(dp)                   , INTENT(IN) :: dd_fill 
    1082       CHARACTER(LEN=*)           , INTENT(IN) :: cd_dim 
    1083       LOGICAL                    , INTENT(IN), OPTIONAL :: ld_discont 
    1084  
    1085       ! function 
    1086       REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), & 
    1087       &                   SIZE(dd_value,DIM=2) ) :: extrap_deriv_2D 
    1088  
    1089       ! local variable 
    1090       INTEGER(i4)                              :: il_imin 
    1091       INTEGER(i4)                              :: il_imax 
    1092       INTEGER(i4)                              :: il_jmin 
    1093       INTEGER(i4)                              :: il_jmax 
    1094       INTEGER(i4), DIMENSION(2)                :: il_shape 
    1095  
    1096       REAL(dp)                                 :: dl_min 
    1097       REAL(dp)                                 :: dl_max 
    1098       REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_value 
    1099  
    1100       LOGICAL                                  :: ll_discont 
    1101  
    1102       ! loop indices 
    1103       INTEGER(i4) :: ji 
    1104       INTEGER(i4) :: jj 
    1105  
    1106       INTEGER(i4) :: i1 
    1107       INTEGER(i4) :: i2 
    1108  
    1109       INTEGER(i4) :: j1 
    1110       INTEGER(i4) :: j2 
    1111       !---------------------------------------------------------------- 
    1112       ! init 
    1113       extrap_deriv_2D(:,:)=dd_fill 
    1114  
    1115       ll_discont=.FALSE. 
    1116       IF( PRESENT(ld_discont) ) ll_discont=ld_discont 
    1117  
    1118       il_shape(:)=SHAPE(dd_value(:,:)) 
    1119  
    1120       SELECT CASE(TRIM(fct_upper(cd_dim))) 
    1121  
    1122       CASE('I') 
    1123  
    1124          ALLOCATE( dl_value(3,il_shape(2)) ) 
    1125          ! compute derivative in i-direction 
    1126          DO ji=1,il_shape(1) 
    1127  
    1128             ! init 
    1129             dl_value(:,:)=dd_fill 
    1130              
    1131             il_imin=MAX(ji-1,1) 
    1132             il_imax=MIN(ji+1,il_shape(1)) 
    1133  
    1134             IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN 
    1135                i1=1  ; i2=3 
    1136             ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN 
    1137                i1=1  ; i2=2 
    1138             ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN 
    1139                i1=2  ; i2=3 
    1140             ENDIF 
    1141  
    1142             dl_value(i1:i2,:)=dd_value(il_imin:il_imax,:) 
    1143             IF( il_imin == 1 )THEN 
    1144                dl_value(:,:)=EOSHIFT( dl_value(:,:), & 
    1145                &                      DIM=1,         & 
    1146                &                      SHIFT=-1,      & 
    1147                &                      BOUNDARY=dl_value(1,:) ) 
    1148             ENDIF 
    1149             IF( il_imax == il_shape(1) )THEN 
    1150                dl_value(:,:)=EOSHIFT( dl_value(:,:), & 
    1151                &                      DIM=1,         & 
    1152                &                      SHIFT=1,       & 
    1153                &                      BOUNDARY=dl_value(3,:)) 
    1154             ENDIF 
    1155  
    1156             IF( ll_discont )THEN 
    1157                dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 
    1158                dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 
    1159                IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 
    1160                   WHERE( dl_value(:,:) < 0_dp )  
    1161                      dl_value(:,:) = dl_value(:,:)+360._dp 
    1162                   END WHERE 
    1163                ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 
    1164                   WHERE( dl_value(:,:) > 180 )  
    1165                      dl_value(:,:) = dl_value(:,:)-180._dp 
    1166                   END WHERE 
    1167                ENDIF 
    1168             ENDIF 
    1169              
    1170             WHERE( dl_value(2,:) /= dd_fill .AND. &  ! ji 
    1171             &      dl_value(3,:) /= dd_fill .AND. &  ! ji+1 
    1172             &      dl_value(1,:) /= dd_fill )        ! ji-1 
    1173  
    1174                extrap_deriv_2D(ji,:)=& 
    1175                &  ( dl_value(3,:) - dl_value(1,:) ) / & 
    1176                &    REAL( il_imax-il_imin,dp) 
    1177  
    1178             END WHERE 
    1179  
    1180          ENDDO 
    1181  
    1182       CASE('J') 
    1183  
    1184          ALLOCATE( dl_value(il_shape(1),3) ) 
    1185          ! compute derivative in j-direction 
    1186          DO jj=1,il_shape(2) 
    1187           
    1188             il_jmin=MAX(jj-1,1) 
    1189             il_jmax=MIN(jj+1,il_shape(2)) 
    1190  
    1191             IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN 
    1192                j1=1  ; j2=3 
    1193             ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN 
    1194                j1=1  ; j2=2 
    1195             ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN 
    1196                j1=2  ; j2=3 
    1197             ENDIF 
    1198  
    1199             dl_value(:,j1:j2)=dd_value(:,il_jmin:il_jmax) 
    1200             IF( il_jmin == 1 )THEN 
    1201                dl_value(:,:)=EOSHIFT( dl_value(:,:), & 
    1202                &                      DIM=2,         & 
    1203                &                      SHIFT=-1,      & 
    1204                &                      BOUNDARY=dl_value(:,1)) 
    1205             ENDIF 
    1206             IF( il_jmax == il_shape(2) )THEN 
    1207                dl_value(:,:)=EOSHIFT( dl_value(:,:), & 
    1208                &                      DIM=2,         & 
    1209                &                      SHIFT=1,       & 
    1210                &                      BOUNDARY=dl_value(:,3)) 
    1211             ENDIF 
    1212  
    1213             IF( ll_discont )THEN 
    1214                dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 
    1215                dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) 
    1216                IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 
    1217                   WHERE( dl_value(:,:) < 0_dp )  
    1218                      dl_value(:,:) = dl_value(:,:)+360._dp 
    1219                   END WHERE 
    1220                ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 
    1221                   WHERE( dl_value(:,:) > 180 )  
    1222                      dl_value(:,:) = dl_value(:,:)-180._dp 
    1223                   END WHERE 
    1224                ENDIF 
    1225             ENDIF 
    1226  
    1227             WHERE( dl_value(:, 2) /= dd_fill .AND. & ! jj 
    1228             &      dl_value(:, 3) /= dd_fill .AND. & ! jj+1 
    1229             &      dl_value(:, 1) /= dd_fill )       ! jj-1 
    1230  
    1231                extrap_deriv_2D(:,jj)=& 
    1232                &  ( dl_value(:,3) - dl_value(:,1) ) / & 
    1233                &   REAL(il_jmax-il_jmin,dp)          
    1234  
    1235             END WHERE 
    1236  
    1237          ENDDO 
    1238           
    1239       END SELECT 
    1240  
    1241       DEALLOCATE( dl_value ) 
    1242  
    1243    END FUNCTION extrap_deriv_2D 
    1244    !------------------------------------------------------------------- 
    1245    !> @brief 
    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. 
    1249    !>  
    1250    !> @details  
    1251    !> optionaly you could specify to take into account east west discontinuity 
    1252    !> (-180° 180° or 0° 360° for longitude variable) 
    1253    !> 
    1254    !> @author J.Paul 
    1255    !> - November, 2013- Initial Version 
    1256    ! 
    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    !------------------------------------------------------------------- 
    1262    PURE FUNCTION extrap_deriv_3D( dd_value, dd_fill, cd_dim, ld_discont ) 
    1263  
    1264       IMPLICIT NONE 
    1265       ! Argument 
    1266       REAL(dp)        , DIMENSION(:,:,:), INTENT(IN) :: dd_value 
    1267       REAL(dp)                          , INTENT(IN) :: dd_fill 
    1268       CHARACTER(LEN=*)                  , INTENT(IN) :: cd_dim 
    1269       LOGICAL                           , INTENT(IN), OPTIONAL :: ld_discont 
    1270  
    1271       ! function 
    1272       REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), & 
    1273       &                   SIZE(dd_value,DIM=2), & 
    1274       &                   SIZE(dd_value,DIM=3)) :: extrap_deriv_3D 
    1275  
    1276       ! local variable 
    1277       INTEGER(i4)                                :: il_imin 
    1278       INTEGER(i4)                                :: il_imax 
    1279       INTEGER(i4)                                :: il_jmin 
    1280       INTEGER(i4)                                :: il_jmax 
    1281       INTEGER(i4)                                :: il_kmin 
    1282       INTEGER(i4)                                :: il_kmax 
    1283       INTEGER(i4), DIMENSION(3)                  :: il_shape 
    1284  
    1285       REAL(dp)                                   :: dl_min 
    1286       REAL(dp)                                   :: dl_max 
    1287       REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_value 
    1288  
    1289       LOGICAL                                    :: ll_discont 
    1290  
    1291       ! loop indices 
    1292       INTEGER(i4) :: ji 
    1293       INTEGER(i4) :: jj 
    1294       INTEGER(i4) :: jk 
    1295  
    1296       INTEGER(i4) :: i1 
    1297       INTEGER(i4) :: i2 
    1298  
    1299       INTEGER(i4) :: j1 
    1300       INTEGER(i4) :: j2 
    1301        
    1302       INTEGER(i4) :: k1 
    1303       INTEGER(i4) :: k2       
    1304       !---------------------------------------------------------------- 
    1305       ! init 
    1306       extrap_deriv_3D(:,:,:)=dd_fill 
    1307  
    1308       ll_discont=.FALSE. 
    1309       IF( PRESENT(ld_discont) ) ll_discont=ld_discont 
    1310  
    1311       il_shape(:)=SHAPE(dd_value(:,:,:)) 
    1312  
    1313  
    1314       SELECT CASE(TRIM(fct_upper(cd_dim))) 
    1315  
    1316       CASE('I') 
    1317  
    1318          ALLOCATE( dl_value(3,il_shape(2),il_shape(3)) ) 
    1319          ! compute derivative in i-direction 
    1320          DO ji=1,il_shape(1) 
    1321              
    1322             il_imin=MAX(ji-1,1) 
    1323             il_imax=MIN(ji+1,il_shape(1)) 
    1324  
    1325             IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN 
    1326                i1=1  ; i2=3 
    1327             ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN 
    1328                i1=1  ; i2=2 
    1329             ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN 
    1330                i1=2  ; i2=3 
    1331             ENDIF 
    1332  
    1333             dl_value(i1:i2,:,:)=dd_value(il_imin:il_imax,:,:) 
    1334             IF( il_imin == 1 )THEN 
    1335                dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 
    1336                &                      DIM=1,         & 
    1337                &                      SHIFT=-1,      & 
    1338                &                      BOUNDARY=dl_value(1,:,:) ) 
    1339             ENDIF 
    1340             IF( il_imax == il_shape(1) )THEN 
    1341                dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 
    1342                &                      DIM=1,         & 
    1343                &                      SHIFT=1,       & 
    1344                &                      BOUNDARY=dl_value(3,:,:)) 
    1345             ENDIF 
    1346  
    1347             IF( ll_discont )THEN 
    1348                dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 
    1349                dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 
    1350                IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 
    1351                   WHERE( dl_value(:,:,:) < 0_dp )  
    1352                      dl_value(:,:,:) = dl_value(:,:,:)+360._dp 
    1353                   END WHERE 
    1354                ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 
    1355                   WHERE( dl_value(:,:,:) > 180 )  
    1356                      dl_value(:,:,:) = dl_value(:,:,:)-180._dp 
    1357                   END WHERE 
    1358                ENDIF 
    1359             ENDIF 
    1360  
    1361             WHERE( dl_value(2,:,:) /= dd_fill .AND. & ! ji 
    1362             &      dl_value(3,:,:) /= dd_fill .AND. & !ji+1  
    1363             &      dl_value(1,:,:) /= dd_fill )       !ji-1 
    1364  
    1365                extrap_deriv_3D(ji,:,:)= & 
    1366                &  ( dl_value(3,:,:) - dl_value(1,:,:) ) / & 
    1367                &  REAL( il_imax-il_imin ,dp) 
    1368  
    1369             END WHERE 
    1370  
    1371          ENDDO 
    1372  
    1373       CASE('J') 
    1374  
    1375          ALLOCATE( dl_value(il_shape(1),3,il_shape(3)) ) 
    1376          ! compute derivative in j-direction 
    1377          DO jj=1,il_shape(2) 
    1378           
    1379             il_jmin=MAX(jj-1,1) 
    1380             il_jmax=MIN(jj+1,il_shape(2)) 
    1381  
    1382             IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN 
    1383                j1=1  ; j2=3 
    1384             ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN 
    1385                j1=1  ; j2=2 
    1386             ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN 
    1387                j1=2  ; j2=3 
    1388             ENDIF 
    1389  
    1390             dl_value(:,j1:j2,:)=dd_value(:,il_jmin:il_jmax,:) 
    1391             IF( il_jmin == 1 )THEN 
    1392                dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 
    1393                &                      DIM=2,         & 
    1394                &                      SHIFT=-1,      & 
    1395                &                      BOUNDARY=dl_value(:,1,:) ) 
    1396             ENDIF 
    1397             IF( il_jmax == il_shape(2) )THEN 
    1398                dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 
    1399                &                      DIM=2,         & 
    1400                &                      SHIFT=1,       & 
    1401                &                      BOUNDARY=dl_value(:,3,:)) 
    1402             ENDIF 
    1403  
    1404             IF( ll_discont )THEN 
    1405                dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 
    1406                dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 
    1407                IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 
    1408                   WHERE( dl_value(:,:,:) < 0_dp )  
    1409                      dl_value(:,:,:) = dl_value(:,:,:)+360._dp 
    1410                   END WHERE 
    1411                ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 
    1412                   WHERE( dl_value(:,:,:) > 180 )  
    1413                      dl_value(:,:,:) = dl_value(:,:,:)-180._dp 
    1414                   END WHERE 
    1415                ENDIF 
    1416             ENDIF 
    1417  
    1418             WHERE( dl_value(:, 2,:) /= dd_fill .AND. & ! jj 
    1419                &   dl_value(:, 3,:) /= dd_fill .AND. & ! jj+1 
    1420             &      dl_value(:, 1,:) /= dd_fill )       ! jj-1 
    1421  
    1422                extrap_deriv_3D(:,jj,:)=& 
    1423                &  ( dl_value(:,3,:) - dl_value(:,1,:) ) / & 
    1424                &  REAL( il_jmax - il_jmin ,dp)          
    1425  
    1426             END WHERE 
    1427  
    1428          ENDDO 
    1429           
    1430       CASE('K') 
    1431          ! compute derivative in k-direction 
    1432          DO jk=1,il_shape(3) 
    1433  
    1434             il_kmin=MAX(jk-1,1) 
    1435             il_kmax=MIN(jk+1,il_shape(3)) 
    1436  
    1437             IF( il_kmin==jk-1 .AND. il_kmax==jk+1 )THEN 
    1438                k1=1  ; k2=3 
    1439             ELSEIF( il_kmin==jk .AND. il_kmax==jk+1 )THEN 
    1440                k1=1  ; k2=2 
    1441             ELSEIF( il_kmin==jk-1 .AND. il_kmax==jk )THEN 
    1442                k1=2  ; k2=3 
    1443             ENDIF 
    1444  
    1445             dl_value(:,:,k1:k2)=dd_value(:,:,il_kmin:il_kmax) 
    1446             IF( il_kmin == 1 )THEN 
    1447                dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 
    1448                &                      DIM=3,         & 
    1449                &                      SHIFT=-1,      & 
    1450                &                      BOUNDARY=dl_value(:,:,1) ) 
    1451             ENDIF 
    1452             IF( il_kmax == il_shape(3) )THEN 
    1453                dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & 
    1454                &                        DIM=3,         & 
    1455                &                        SHIFT=1,       & 
    1456                &                        BOUNDARY=dl_value(:,:,3)) 
    1457             ENDIF 
    1458  
    1459             IF( ll_discont )THEN 
    1460                dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 
    1461                dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) 
    1462                IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 
    1463                   WHERE( dl_value(:,:,:) < 0_dp )  
    1464                      dl_value(:,:,:) = dl_value(:,:,:)+360._dp 
    1465                   END WHERE 
    1466                ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 
    1467                   WHERE( dl_value(:,:,:) > 180 )  
    1468                      dl_value(:,:,:) = dl_value(:,:,:)-180._dp 
    1469                   END WHERE 
    1470                ENDIF 
    1471             ENDIF          
    1472  
    1473             WHERE( dl_value(:,:, 2) /= dd_fill .AND. & ! jk 
    1474                &   dl_value(:,:, 3) /= dd_fill .AND. & ! jk+1 
    1475                &   dl_value(:,:, 1) /= dd_fill )       ! jk-1 
    1476  
    1477                extrap_deriv_3D(:,:,jk)=& 
    1478                &  ( dl_value(:,:,3) - dl_value(:,:,1) ) / & 
    1479                &  REAL( il_kmax-il_kmin,dp)          
    1480  
    1481             END WHERE 
    1482  
    1483          ENDDO 
    1484  
    1485       END SELECT 
    1486  
    1487       DEALLOCATE( dl_value ) 
    1488  
    1489    END FUNCTION extrap_deriv_3D 
    1490842   !------------------------------------------------------------------- 
    1491843   !> @brief 
     
    1493845   !>  
    1494846   !> @details  
    1495    !> coefficients are  "grid distance" to the center of the box choosed to compute 
    1496    !> extrapolation. 
     847   !> coefficients are  "grid distance" to the center of the box  
     848   !> choosed to compute extrapolation. 
    1497849   !> 
    1498850   !> @author J.Paul 
    1499    !> - November, 2013- Initial Version 
     851   !> @date November, 2013 - Initial Version 
     852   !> @date July, 2015  
     853   !> - decrease weight of third dimension 
    1500854   ! 
    1501855   !> @param[in] dd_value  3D array of variable to be extrapolated 
     
    1544898 
    1545899               ! compute distance 
     900               ! "vertical weight" is lower than horizontal  
    1546901               dl_dist(ji,jj,jk) = (ji-il_imid)**2 + & 
    1547902               &                   (jj-il_jmid)**2 + & 
    1548                &                   (jk-il_kmid)**2 
     903               &                 3*(jk-il_kmid)**2 
    1549904 
    1550905               IF( dl_dist(ji,jj,jk) /= 0 )THEN 
     
    1569924   !>  
    1570925   !> @author J.Paul 
    1571    !> - November, 2013- Initial Version 
     926   !> @date November, 2013 - Initial Version 
    1572927   !> 
    1573928   !> @param[in] dd_value  3D array of variable to be extrapolated 
     
    16581013   !> 
    16591014   !> @author J.Paul 
    1660    !> - November, 2013- Initial Version 
     1015   !> @date November, 2013 - Initial Version 
     1016   !> @date July, 2015  
     1017   !> - decrease weight of third dimension 
    16611018   ! 
    16621019   !> @param[in] dd_value  3D array of variable to be extrapolated 
     
    17051062 
    17061063               ! compute distance 
     1064               ! "vertical weight" is lower than horizontal  
    17071065               dl_dist(ji,jj,jk) = (ji-il_imid)**2 + & 
    17081066               &                   (jj-il_jmid)**2 + & 
    1709                &                   (jk-il_kmid)**2 
     1067               &                 3*(jk-il_kmid)**2 
    17101068 
    17111069               IF( dl_dist(ji,jj,jk) /= 0 )THEN 
     
    17321090   !> 
    17331091   !> @author J.Paul 
    1734    !> - November, 2013- Initial Version 
     1092   !> @date November, 2013 - Initial Version 
    17351093   ! 
    17361094   !> @param[in] dd_value  3D array of variable to be extrapolated 
     
    17631121      INTEGER(i4) :: jj 
    17641122      INTEGER(i4) :: jk 
    1765  
    17661123      !---------------------------------------------------------------- 
    17671124 
     
    17931150            ENDDO 
    17941151         ENDDO 
     1152 
    17951153 
    17961154         ! return value 
     
    18151173   !> 
    18161174   !> @author J.Paul 
    1817    !> - November, 2013-Initial version 
     1175   !> @date November, 2013 - Initial version 
    18181176   ! 
    18191177   !> @param[inout] td_var variable  
     
    19171275   !> 
    19181276   !> @author J.Paul 
    1919    !> - November, 2013-Initial version 
     1277   !> @date November, 2013 - Initial version 
    19201278   !> 
    19211279   !> @param[inout] td_var variable  
Note: See TracChangeset for help on using the changeset viewer.