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 5602 for branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/TOOLS/SIREN/src/interp.f90 – NEMO

Ignore:
Timestamp:
2015-07-16T13:55:15+02:00 (9 years ago)
Author:
cbricaud
Message:

merge change from trunk rev 5003 to 5519 ( rev where branche 3.6_stable were created )

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/TOOLS/SIREN/src/interp.f90

    r4213 r5602  
    88!> @brief  
    99!> This module manage interpolation on regular grid. 
    10 !> @note It is used to work on ORCA grid, as we work only with grid indices. 
     10!> 
     11!> @details Interpolation method to be used is specify inside variable  
     12!>    strcuture, as array of string character.<br/>  
     13!>    - td_var\%c_interp(1) string character is the interpolation name choose between: 
     14!>       - 'nearest' 
     15!>       - 'cubic  ' 
     16!>       - 'linear ' 
     17!>    - td_var\%c_interp(2) string character is an operation to be used  
     18!>    on interpolated value.<br/> 
     19!>          operation have to be mulitplication '*' or division '/'.<br/> 
     20!>          coefficient have to be refinement factor following i-direction 'rhoi',  
     21!>          j-direction 'rhoj', or k-direction 'rhok'.<br/> 
     22!> 
     23!>          Examples: '*rhoi', '/rhoj'. 
     24!>  
     25!>    @note Those informations are read from namelist or variable configuration file (default).<br/> 
     26!>    Interplation method could be specify for each variable in namelist _namvar_, 
     27!>    defining string character _cn\_varinfo_.<br/> 
     28!>    Example: 
     29!>       - cn_varinfo='varname1:cubic/rhoi', 'varname2:linear'  
     30!> 
     31!>    to create mixed grid (with coarse grid point needed to compute 
     32!> interpolation):<br/> 
     33!> @code 
     34!>    CALL interp_create_mixed_grid( td_var, td_mix [,id_rho] ) 
     35!> @endcode 
     36!>       - td_var is coarse grid variable (should be extrapolated) 
     37!>       - td_mix is mixed grid variable structure [output] 
     38!>       - id_rho is array of refinment factor [optional] 
     39!> 
     40!>    to detected point to be interpolated:<br/> 
     41!> @code 
     42!>    il_detect(:,:,:)=interp_detect( td_mix [,id_rho] ) 
     43!> @endcode 
     44!>       - il_detect(:,:,:) is 3D array of detected point to be interpolated 
     45!>       - td_mix is mixed grid variable 
     46!>       - id_rho is array of refinement factor [optional] 
     47!> 
     48!>    to interpolate variable value:<br/> 
     49!> @code 
     50!>    CALL  interp_fill_value( td_var [,id_rho] [,id_offset] ) 
     51!> @endcode 
     52!>       - td_var is variable structure 
     53!>       - id_rho is array of refinement factor [optional] 
     54!>       - id_offset is array of offset between fine and coarse grid [optional] 
     55!> 
     56!>    to clean mixed grid (remove points added on mixed grid to compute interpolation):<br/> 
     57!> @code 
     58!>    CALL interp_clean_mixed_grid( td_mix, td_var, id_rho ) 
     59!> @endcode 
     60!>       - td_mix is mixed grid variable structure 
     61!>       - td_var is variable structure [output] 
     62!>       - id_rho is array of refinement factor [optional] 
     63!>       - id_offset is array of offset between fine and coarse grid [optional] 
     64!> 
     65!> @note It use to work on ORCA grid, as we work only with grid indices. 
    1166!> 
    1267!> @warning due to the use of second derivative when using cubic interpolation 
    13 !> you should add 2 extrabands 
     68!> you should add at least 2 extrabands. 
    1469!> 
    15 !> @details 
    1670!> @author 
    1771!> J.Paul 
    1872! REVISION HISTORY: 
    19 !> @date Nov, 2013 - Initial Version 
    20 ! 
     73!> @date November, 2013 - Initial Version 
     74!> @date September, 2014 
     75!> - add header 
     76!> - use interpolation method modules 
     77!> 
    2178!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    22 !> @todo 
    23 !> - interp 3D 
    24 !> - see issue when fill value is zero for check cubic_fill.. 
    2579!---------------------------------------------------------------------- 
    2680MODULE interp 
     
    2983   USE global                          ! global variable 
    3084   USE kind                            ! F90 kind parameter 
    31    USE logger                             ! log file manager 
     85   USE logger                          ! log file manager 
    3286   USE fct                             ! basic useful function 
    33 !   USE date                            ! date manager 
     87   USE date                            ! date manager 
    3488   USE att                             ! attribute manager 
    3589   USE dim                             ! dimension manager 
    3690   USE var                             ! variable manager 
    37 !   USE file                            ! file manager 
    38 !   USE iom                             ! I/O manager 
    39 !   USE dom                             ! domain manager 
    4091   USE grid                            ! grid manager 
    4192   USE extrap                          ! extrapolation manager 
    42 !   USE interp                          ! interpolation manager 
    43 !   USE filter                          ! filter manager 
    44 !   USE mpp                             ! MPP manager 
    45 !   USE iom_mpp                         ! MPP I/O manager 
     93   USE interp_cubic                    !   cubic interpolation manager 
     94   USE interp_linear                   !  linear interpolation manager 
     95   USE interp_nearest                  ! nearest interpolation manager 
    4696 
    4797   IMPLICIT NONE 
    48    PRIVATE 
    4998   ! NOTE_avoid_public_variables_if_possible 
    5099 
     
    53102   ! function and subroutine 
    54103   PUBLIC :: interp_detect             !< detected point to be interpolated 
    55    PUBLIC :: interp_fill_value         !< interpolate value over detectected point 
     104   PUBLIC :: interp_fill_value         !< interpolate value 
    56105   PUBLIC :: interp_create_mixed_grid  !< create mixed grid 
    57106   PUBLIC :: interp_clean_mixed_grid   !< clean mixed grid 
    58107 
    59    PRIVATE :: interp__detect           !< detected point to be interpolated 
    60    PRIVATE :: interp__detect_wrapper   !< detected point to be interpolated 
    61    PRIVATE :: interp__fill_value_wrapper !< interpolate value over detectected point 
    62    PRIVATE :: interp__fill_value       !< interpolate value over detectected point 
    63    PRIVATE :: interp__clean_even_grid  !< clean even mixed grid 
    64    PRIVATE :: interp__del_offset       !< remove offset from interpolated grid 
    65    PRIVATE :: interp__check_method     !< check if interpolation method available 
    66 !   PRIVATE :: interp__3D               !< interpolate 3D grid 
    67    PRIVATE :: interp__2D               !< interpolate 2D grid 
    68    PRIVATE :: interp__1D               !< interpolate 1D grid  
    69    PRIVATE :: interp__2D_cubic_coef    !< compute coefficient for bicubic interpolation  
    70    PRIVATE :: interp__2D_cubic_fill    !< compute bicubic interpolation  
    71    PRIVATE :: interp__2D_linear_coef   !< compute coefficient for bilinear interpolation 
    72    PRIVATE :: interp__2D_linear_fill   !< compute bilinear interpolation  
    73    PRIVATE :: interp__2D_nearest_fill  !< compute nearest interpolation  
    74    PRIVATE :: interp__1D_cubic_coef    !< compute coefficient for cubic interpolation  
    75    PRIVATE :: interp__1D_cubic_fill    !< compute cubic interpolation  
    76    PRIVATE :: interp__1D_linear_coef   !< compute coefficient for linear interpolation 
    77    PRIVATE :: interp__1D_linear_fill   !< compute linear interpolation  
    78    PRIVATE :: interp__1D_nearest_fill  !< compute nearest interpolation  
    79 !   PRIVATE :: interp__longitude 
     108   PRIVATE :: interp__detect              ! detected point to be interpolated 
     109   PRIVATE :: interp__detect_wrapper      ! detected point to be interpolated 
     110   PRIVATE :: interp__fill_value_wrapper  ! interpolate value over detectected point 
     111   PRIVATE :: interp__fill_value          ! interpolate value over detectected point 
     112   PRIVATE :: interp__clean_even_grid     ! clean even mixed grid 
     113   PRIVATE :: interp__check_method        ! check if interpolation method available 
    80114    
    81115   TYPE TINTERP 
    82       !CHARACTER(LEN=lc) :: c_name = 'unknown'   !< interpolation method name 
    83       !CHARACTER(LEN=lc) :: c_factor = 'unknown' !< interpolation factor 
    84       !CHARACTER(LEN=lc) :: c_divisor = 'unknown'   !< interpolation divisor 
    85116      CHARACTER(LEN=lc) :: c_name   = '' !< interpolation method name 
    86117      CHARACTER(LEN=lc) :: c_factor = '' !< interpolation factor 
     
    89120 
    90121   INTERFACE interp_detect 
    91       MODULE PROCEDURE interp__detect_wrapper !< detected point to be interpolated 
     122      MODULE PROCEDURE interp__detect_wrapper  
    92123   END INTERFACE interp_detect 
    93124 
    94125   INTERFACE interp_fill_value 
    95       MODULE PROCEDURE interp__fill_value_wrapper !< detected point to be interpolated 
     126      MODULE PROCEDURE interp__fill_value_wrapper  
    96127   END INTERFACE interp_fill_value 
    97128 
     
    99130   !------------------------------------------------------------------- 
    100131   !> @brief 
    101    !> This function check if interpolation method available. 
     132   !> This function check if interpolation method is available. 
    102133   !>  
    103134   !> @details  
     135   !> check if name of interpolation method is present in global list of string 
     136   !> character cp_interp_list (see global.f90). 
    104137   !> 
    105138   !> @author J.Paul 
    106    !> - Nov, 2013- Initial Version 
     139   !> - November, 2013- Initial Version 
    107140   ! 
    108    !> @param[in] cd_method : interpolation method 
     141   !> @param[in] cd_method interpolation method 
    109142   !> @return  
    110    !> @todo see extrap_detect 
    111    !------------------------------------------------------------------- 
    112    !> @code 
     143   !------------------------------------------------------------------- 
    113144   FUNCTION interp__check_method( cd_method ) 
    114145      IMPLICIT NONE 
     
    130161 
    131162      interp__check_method=.FALSE. 
    132       DO ji=1,ig_ninterp 
    133          cl_interp=fct_lower(cg_interp_list(ji)) 
     163      DO ji=1,ip_ninterp 
     164         cl_interp=fct_lower(cp_interp_list(ji)) 
    134165         IF( TRIM(cl_interp) == TRIM(cl_method) )THEN 
    135166            interp__check_method=.TRUE. 
     
    139170 
    140171   END FUNCTION interp__check_method 
    141    !> @endcode 
    142172   !------------------------------------------------------------------- 
    143173   !> @brief 
     
    149179   !> 
    150180   !> @author J.Paul 
    151    !> - Nov, 2013- Initial Version 
    152    ! 
    153    !> @param[in] td_mix : mixed grid variable (to interpolate) 
    154    !> @param[in] id_rho : table of refinement factor 
    155    !> @return table of detected point to be interpolated  
    156    !------------------------------------------------------------------- 
    157    !> @code 
     181   !> - November, 2013- Initial Version 
     182   !> 
     183   !> @param[in] td_mix mixed grid variable (to interpolate) 
     184   !> @param[in] id_rho array of refinement factor 
     185   !> @return 3D array of detected point to be interpolated  
     186   !------------------------------------------------------------------- 
    158187   FUNCTION interp__detect_wrapper( td_mix, id_rho ) 
    159188      IMPLICIT NONE 
     
    168197 
    169198      ! local variable 
    170  
    171199      ! loop indices 
    172200      !---------------------------------------------------------------- 
     
    208236 
    209237   END FUNCTION interp__detect_wrapper 
    210    !> @endcode 
    211238   !------------------------------------------------------------------- 
    212239   !> @brief 
     
    217244   !> 
    218245   !> @author J.Paul 
    219    !> - Nov, 2013- Initial Version 
     246   !> - November, 2013- Initial Version 
    220247   ! 
    221    !> @param[in] td_mix : mixed grid variable (to interpolate) 
    222    !> @param[in] id_rho : table of refinement factor 
    223    !> @return table of detected point to be interpolated  
    224    !------------------------------------------------------------------- 
    225    !> @code 
     248   !> @param[in] td_mix mixed grid variable (to interpolate) 
     249   !> @param[in] id_rho array of refinement factor 
     250   !> @return 3D array of detected point to be interpolated  
     251   !------------------------------------------------------------------- 
    226252   FUNCTION interp__detect( td_mix, id_rho ) 
    227253      IMPLICIT NONE 
     
    253279      ALLOCATE( il_rho(ip_maxdim) ) 
    254280      il_rho(:)=1 
    255       IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:) 
     281      IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:) 
    256282 
    257283      ! special case for even refinement on ARAKAWA-C grid 
     
    289315      il_dim(:)=td_mix%t_dim(1:3)%i_len 
    290316 
     317      ! init 
    291318      interp__detect(:,:,:)=1 
    292319 
     
    297324 
    298325      ! do not compute point near fill value 
    299       FORALL( ji=1:il_dim(1):il_rho(jp_I), & 
    300       &       jj=1:il_dim(2):il_rho(jp_J), & 
    301       &       jk=1:il_dim(3):il_rho(jp_K), & 
    302       &       td_mix%d_value(ji,jj,jk,1) == td_mix%d_fill ) 
    303  
    304          ! i-direction 
    305          interp__detect(MAX(1,ji-il_xextra):MIN(ji+il_xextra,il_dim(1)),& 
    306          &              MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),& 
    307          &              MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0  
    308          ! j-direction 
    309          interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),& 
    310          &              MAX(1,jj-il_yextra):MIN(jj+il_yextra,il_dim(2)),& 
    311          &              MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0 
    312          ! k-direction 
    313          interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),& 
    314          &              MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),& 
    315          &              MAX(1,jk-il_zextra):MIN(jk+il_zextra,il_dim(3)) )=0          
    316  
    317       END FORALL 
     326      DO jk=1,il_dim(3),il_rho(jp_K) 
     327         DO jj=1,il_dim(2),il_rho(jp_J) 
     328            DO ji=1,il_dim(1),il_rho(jp_I) 
     329 
     330               IF( td_mix%d_value(ji,jj,jk,1) == td_mix%d_fill )THEN 
     331 
     332                  ! i-direction 
     333                  interp__detect(MAX(1,ji-il_xextra):MIN(ji+il_xextra,il_dim(1)),& 
     334                  &              MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),& 
     335                  &              MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0  
     336                  ! j-direction 
     337                  interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),& 
     338                  &              MAX(1,jj-il_yextra):MIN(jj+il_yextra,il_dim(2)),& 
     339                  &              MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0 
     340                  ! k-direction 
     341                  interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),& 
     342                  &              MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),& 
     343                  &              MAX(1,jk-il_zextra):MIN(jk+il_zextra,il_dim(3)) )=0          
     344 
     345               ENDIF 
     346 
     347            ENDDO 
     348         ENDDO 
     349      ENDDO 
     350 
     351      DEALLOCATE( il_rho ) 
    318352 
    319353   END FUNCTION interp__detect 
    320    !> @endcode 
    321354   !------------------------------------------------------------------- 
    322355   !> @brief 
     
    330363   !> 
    331364   !> @author J.Paul 
    332    !> - Nov, 2013- Initial Version 
    333    ! 
    334    !> @param[in] td_var : coarse grid variable (should be extrapolated) 
    335    !> @param[out] td_mix : mixed grid variable 
    336    !> @param[in] id_rho  : table of refinment factor  
    337    !> 
    338    !> @todo 
    339    !------------------------------------------------------------------- 
    340    !> @code 
     365   !> - November, 2013- Initial Version 
     366   !> 
     367   !> @param[in] td_var    coarse grid variable (should be extrapolated) 
     368   !> @param[out] td_mix   mixed grid variable 
     369   !> @param[in] id_rho    array of refinment factor (default 1)  
     370   !------------------------------------------------------------------- 
    341371   SUBROUTINE interp_create_mixed_grid( td_var, td_mix, id_rho ) 
    342372      IMPLICIT NONE 
     
    359389      ALLOCATE(il_rho(ip_maxdim)) 
    360390      il_rho(:)=1 
    361       IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:) 
     391      IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:) 
    362392 
    363393      ! special case for even refinement on ARAKAWA-C grid 
     
    381411 
    382412      ! copy variable 
    383       td_mix=td_var 
     413      td_mix=var_copy(td_var) 
    384414 
    385415      ! compute new dimension length 
     
    408438      &     td_var%d_value(:,:,:,:) 
    409439 
     440      DEALLOCATE(il_rho) 
     441 
    410442   END SUBROUTINE interp_create_mixed_grid 
    411    !> @endcode 
    412443   !------------------------------------------------------------------- 
    413444   !> @brief 
     
    418449   !> 
    419450   !> @author J.Paul 
    420    !> - Nov, 2013- Initial Version 
    421    ! 
    422    !> @param[in]  td_mix : mixed grid variable 
    423    !> @param[in] id_rho  : table of refinment factor  
    424    !> 
    425    !> @todo 
    426    !------------------------------------------------------------------- 
    427    !> @code 
     451   !> - November, 2013- Initial Version 
     452   !> 
     453   !> @param[inout] td_mix mixed grid variable 
     454   !> @param[in] id_rho    array of refinment factor  
     455   !------------------------------------------------------------------- 
    428456   SUBROUTINE interp__clean_even_grid( td_mix, id_rho ) 
    429457      IMPLICIT NONE 
    430458      ! Argument 
    431       TYPE(TVAR) , INTENT(INOUT) :: td_mix  
     459      TYPE(TVAR) ,               INTENT(INOUT) :: td_mix  
    432460      INTEGER(I4), DIMENSION(:), INTENT(IN   ), OPTIONAL :: id_rho 
    433461 
     
    468496      END SELECT 
    469497 
    470       ! copy variable 
    471       tl_mix=td_mix 
    472  
    473498      ! remove some point only if refinement in some direction is even 
    474499      IF( ANY(ll_even(:)) )THEN 
     500 
     501         ! copy variable 
     502         tl_mix=var_copy(td_mix) 
    475503 
    476504         ALLOCATE( ll_mask( tl_mix%t_dim(1)%i_len, & 
     
    554582 
    555583            td_mix%d_value(:,:,:,:)=RESHAPE( dl_vect(:), & 
    556             &                                SHAPE=td_mix%t_dim(:)%i_len ) 
     584            &                                SHAPE=(/td_mix%t_dim(1)%i_len, & 
     585            &                                        td_mix%t_dim(2)%i_len, & 
     586            &                                        td_mix%t_dim(3)%i_len, & 
     587            &                                        td_mix%t_dim(4)%i_len/) ) 
    557588 
    558589            DEALLOCATE( dl_vect ) 
     
    562593         DEALLOCATE( ll_mask ) 
    563594 
     595         ! clean 
     596         CALL var_clean(tl_mix) 
     597 
    564598      ENDIF 
    565599 
    566       CALL var_clean(tl_mix) 
     600      ! clean 
     601      DEALLOCATE(il_rho) 
    567602 
    568603   END SUBROUTINE interp__clean_even_grid 
    569    !> @endcode 
    570604   !------------------------------------------------------------------- 
    571605   !> @brief 
    572    !> This subroutine save interpolated value over domain.  
    573    !> And so remove points added on mixed grid  
    574    !> to compute interpolation 
     606   !> This subroutine remove points added on mixed grid  
     607   !> to compute interpolation. And save interpolated value over domain.  
    575608   !>  
    576609   !> @details  
    577610   !> 
    578611   !> @author J.Paul 
    579    !> - Nov, 2013- Initial Version 
    580    ! 
    581    !> @param[in] td_var : table of mixed grid variable (to interpolate) 
    582    !> @param[in] id_rho : table of refinement factor 
    583    !> @return  
    584    !------------------------------------------------------------------- 
    585    !> @code 
     612   !> - November, 2013- Initial Version 
     613   !> @date September, 2014 
     614   !> - use offset to save useful domain 
     615   !> 
     616   !> @param[in] td_mix    mixed grid variable structure 
     617   !> @param[out] td_var   variable structure 
     618   !> @param[in] id_rho    array of refinement factor (default 1) 
     619   !> @param[in] id_offset 2D array of offset between fine and coarse grid 
     620   !------------------------------------------------------------------- 
    586621   SUBROUTINE interp_clean_mixed_grid( td_mix, td_var, & 
    587    &                                   id_rho, & 
    588    &                                   id_offset ) 
     622   &                                   id_rho, id_offset ) 
    589623      IMPLICIT NONE 
    590624      ! Argument 
    591625      TYPE(TVAR)                 , INTENT(IN   ) :: td_mix 
    592626      TYPE(TVAR)                 , INTENT(  OUT) :: td_var 
    593       INTEGER(I4), DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_rho 
    594       INTEGER(I4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset 
     627      INTEGER(I4), DIMENSION(:)  , INTENT(IN   ) :: id_rho 
     628      INTEGER(I4), DIMENSION(2,2), INTENT(IN   ) :: id_offset 
    595629 
    596630      ! local variable 
     
    605639      INTEGER(i4) :: il_jmax1 
    606640 
    607       INTEGER(i4), DIMENSION(2,2) :: il_offset 
    608  
    609641      REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value 
    610642 
     
    613645      ! loop indices 
    614646      !---------------------------------------------------------------- 
    615       il_offset(:,:)=0  
    616       IF( PRESENT(id_offset) )THEN 
    617          IF( ANY( SHAPE(id_offset(:,:)) /= SHAPE(il_offset(:,:)) ) )THEN 
    618             CALL logger_error("INTERP CLEAN MIXED GRID: invalid dimension of"//& 
    619             &                 " offset table") 
    620          ELSE 
    621             il_offset(:,:)=id_offset(:,:) 
    622          ENDIF 
    623       ENDIF 
    624647      ! copy mixed variable in temporary structure 
    625       tl_mix=td_mix 
     648      tl_mix=var_copy(td_mix) 
    626649 
    627650      ! remove unusefull points over mixed grid for even refinement 
     
    629652 
    630653      ! copy cleaned mixed variable 
    631       td_var=tl_mix 
    632       IF( ALL(td_var%t_dim(1:2)%l_use) )THEN 
    633  
    634          ! delete table of value 
    635          CALL var_del_value(td_var) 
    636  
    637          ! compute domain indices 
    638          il_imin0=1  ; il_imax0=td_var%t_dim(1)%i_len 
    639          il_jmin0=1  ; il_jmax0=td_var%t_dim(2)%i_len 
    640     
    641          il_imin1=il_imin0+il_offset(1,1) 
    642          il_jmin1=il_jmin0+il_offset(2,1) 
    643     
    644          il_imax1=il_imax0-il_offset(1,2) 
    645          il_jmax1=il_jmax0-il_offset(2,2) 
    646  
    647          SELECT CASE(TRIM(td_var%c_point)) 
    648          CASE('U') 
    649             il_imin1=il_imin0+(il_offset(1,1)-1) 
    650             il_imax1=il_imax0-(il_offset(1,2)+1) 
    651          CASE('V') 
    652             il_jmin1=il_jmin0+(il_offset(2,1)-1) 
    653             il_jmax1=il_jmax0-(il_offset(2,2)+1) 
    654          CASE('F') 
    655             il_imin1=il_imin0+(il_offset(1,1)-1) 
    656             il_imax1=il_imax0-(il_offset(1,2)+1) 
    657  
    658             il_jmin1=il_jmin0+(il_offset(2,1)-1) 
    659             il_jmax1=il_jmax0-(il_offset(2,2)+1) 
    660          END SELECT 
    661  
    662          ! compute new dimension 
    663          td_var%t_dim(1)%i_len=il_imax1-il_imin1+1 
    664          td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1 
     654      td_var=var_copy(tl_mix) 
     655 
     656      ! delete array of value 
     657      CALL var_del_value(td_var) 
     658 
     659      ! compute domain indices in i-direction 
     660      il_imin0=1  ; il_imax0=td_var%t_dim(1)%i_len 
    665661  
    666          ALLOCATE(dl_value(td_var%t_dim(1)%i_len, & 
    667          &                 td_var%t_dim(2)%i_len, & 
    668          &                 td_var%t_dim(3)%i_len, & 
    669          &                 td_var%t_dim(4)%i_len) ) 
    670  
    671          dl_value( 1:td_var%t_dim(1)%i_len, & 
    672          &         :,:,:) = tl_mix%d_value( il_imin1:il_imax1, & 
    673          &                                  il_jmin1:il_jmax1, & 
    674          &                                      :, : ) 
    675  
    676          ! add variable value 
    677          CALL var_add_value(td_var,dl_value(:,:,:,:)) 
    678  
    679          ! save variable type 
    680          td_var%i_type=tl_mix%i_type 
    681  
    682          DEALLOCATE(dl_value) 
     662      IF( td_var%t_dim(1)%l_use )THEN 
     663         il_imin1=il_imin0+id_offset(Jp_I,1) 
     664         il_imax1=il_imax0-id_offset(Jp_I,2) 
     665      ELSE 
     666 
     667         il_imin1=il_imin0 
     668         il_imax1=il_imax0 
    683669 
    684670      ENDIF 
    685        
     671 
     672      ! compute domain indices in j-direction 
     673      il_jmin0=1  ; il_jmax0=td_var%t_dim(2)%i_len 
     674  
     675      IF( td_var%t_dim(2)%l_use )THEN 
     676         il_jmin1=il_jmin0+id_offset(Jp_J,1) 
     677         il_jmax1=il_jmax0-id_offset(Jp_J,2) 
     678      ELSE 
     679 
     680         il_jmin1=il_jmin0 
     681         il_jmax1=il_jmax0 
     682  
     683      ENDIF 
     684 
     685      ! compute new dimension 
     686      td_var%t_dim(1)%i_len=il_imax1-il_imin1+1 
     687      td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1 
     688  
     689      ALLOCATE(dl_value(td_var%t_dim(1)%i_len, & 
     690      &                 td_var%t_dim(2)%i_len, & 
     691      &                 td_var%t_dim(3)%i_len, & 
     692      &                 td_var%t_dim(4)%i_len) ) 
     693 
     694      dl_value( 1:td_var%t_dim(1)%i_len, & 
     695      &         1:td_var%t_dim(2)%i_len, & 
     696      &         :,:) = tl_mix%d_value( il_imin1:il_imax1, & 
     697      &                                il_jmin1:il_jmax1, & 
     698      &                                     :, : ) 
     699 
     700      ! add variable value 
     701      CALL var_add_value(td_var,dl_value(:,:,:,:)) 
     702 
     703      DEALLOCATE(dl_value) 
     704  
     705      ! clean 
    686706      CALL var_clean(tl_mix) 
    687707 
    688708   END SUBROUTINE interp_clean_mixed_grid 
    689    !> @endcode 
    690709   !------------------------------------------------------------------- 
    691710   !> @brief 
    692    !> This subroutine save interpolated value over domain.  
    693    !> And so remove points added on mixed grid  
    694    !> to compute interpolation 
    695    !>  
    696    !> @details  
    697    !> 
    698    !> @author J.Paul 
    699    !> - Nov, 2013- Initial Version 
    700    ! 
    701    !> @param[in] td_var : table of mixed grid variable (to interpolate) 
    702    !> @param[in] id_offset : table of offset 
    703    !> @return  
    704    !------------------------------------------------------------------- 
    705    !> @code 
    706    SUBROUTINE interp__del_offset( td_var, id_offset ) 
    707       IMPLICIT NONE 
    708       ! Argument 
    709       TYPE(TVAR)                 , INTENT(INOUT) :: td_var 
    710       INTEGER(i4), DIMENSION(:,:), INTENT(IN   ) :: id_offset 
    711  
    712       ! local variable 
    713       INTEGER(i4) :: il_imin1 
    714       INTEGER(i4) :: il_jmin1 
    715       INTEGER(i4) :: il_imax1 
    716       INTEGER(i4) :: il_jmax1 
    717  
    718       REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value 
    719  
    720       ! loop indices 
    721       !---------------------------------------------------------------- 
    722     
    723       IF( ALL(td_var%t_dim(1:2)%l_use) )THEN 
    724  
    725          il_imin1=1+id_offset(1,1) 
    726          il_jmin1=1+id_offset(2,1) 
    727  
    728          il_imax1=td_var%t_dim(1)%i_len-id_offset(2,1) 
    729          il_jmax1=td_var%t_dim(2)%i_len-id_offset(2,2) 
    730  
    731          ! compute new dimension 
    732          td_var%t_dim(1)%i_len=il_imax1-il_imin1+1 
    733          td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1 
    734          ALLOCATE( dl_value( td_var%t_dim(1)%i_len, & 
    735          &                   td_var%t_dim(2)%i_len, & 
    736          &                   td_var%t_dim(3)%i_len, & 
    737          &                   td_var%t_dim(4)%i_len) ) 
    738  
    739          dl_value(:,:,:,:)=td_var%d_value( il_imin1:il_imax1, & 
    740          &                                 il_jmin1:il_jmax1, & 
    741          &                                 :,: ) 
    742  
    743          ! delete table of value 
    744          CALL var_del_value(td_var) 
    745  
    746          ! add variable value 
    747          CALL var_add_value(td_var,dl_value(:,:,:,:)) 
    748  
    749          DEALLOCATE(dl_value) 
    750  
    751       ENDIF 
    752  
    753    END SUBROUTINE interp__del_offset 
    754    !> @endcode 
    755    !------------------------------------------------------------------- 
    756    !> @brief 
    757    !> This subroutine interpolate detected point. 
     711   !> This subroutine interpolate variable value. 
    758712   !>  
    759713   !> @details  
     
    762716   !> 
    763717   !> @author J.Paul 
    764    !> - Nov, 2013- Initial Version 
    765    ! 
    766    !> @param[in] td_var : variable to be interpolated 
    767    !> @param[in] id_rho : table of refinement factor 
    768    !------------------------------------------------------------------- 
    769    !> @code 
     718   !> - November, 2013- Initial Version 
     719   !> 
     720   !> @param[inout] td_var variable structure 
     721   !> @param[in] id_rho    array of refinement factor 
     722   !> @param[in] id_offset 2D array of offset between fine and coarse grid 
     723   !------------------------------------------------------------------- 
    770724   SUBROUTINE interp__fill_value_wrapper( td_var, &  
    771725   &                                      id_rho, & 
     
    774728      ! Argument 
    775729      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var 
    776       INTEGER(I4), DIMENSION(:) INTENT(IN   ), OPTIONAL :: id_rho 
     730      INTEGER(I4), DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_rho 
    777731      INTEGER(I4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset 
    778732 
    779733      ! local variable 
    780       CHARACTER(LEN=lc)                 :: cl_method 
    781       INTEGER(i4)      , DIMENSION(:), ALLOCATABLE  :: il_rho 
    782       INTEGER(i4)      , DIMENSION(2,2) :: il_offset 
     734      CHARACTER(LEN=lc)                            :: cl_method 
     735      INTEGER(i4)      , DIMENSION(:), ALLOCATABLE :: il_rho 
     736      INTEGER(i4)      , DIMENSION(2,2)            :: il_offset 
    783737 
    784738      ! loop indices 
    785739      !---------------------------------------------------------------- 
    786740 
    787       SELECT CASE(TRIM(td_var%c_interp(1))) 
    788          CASE('cubic','linear','nearest') 
    789             cl_method=TRIM(td_var%c_interp(1)) 
    790          CASE DEFAULT 
    791             CALL logger_warn("INTERP FILL: interpolation method unknown."//& 
    792             &  " use linear interpolation") 
    793             cl_method='linear' 
    794  
    795             ! update variable structure value 
    796             td_var%c_interp(1)='linear' 
    797       END SELECT 
    798  
    799741      ALLOCATE( il_rho(ip_maxdim) ) 
    800742      il_rho(:)=1 
    801       IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:) 
     743      IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:) 
    802744      IF( ANY(il_rho(:) < 0) )THEN 
    803745         CALL logger_error("INTERP FILL VALUE: invalid "//& 
     
    808750      IF( PRESENT(id_offset) )THEN 
    809751         IF( ANY(SHAPE(id_offset(:,:)) /= (/2,2/)) )THEN 
    810             CALL logger_error("INTERP FILL VALUE: invalid table of offset") 
     752            CALL logger_error("INTERP FILL VALUE: invalid array of offset") 
    811753         ELSE 
    812754            il_offset(:,:)=id_offset(:,:) 
     
    817759      &   (il_rho(jp_J) /= 1 .AND. td_var%t_dim(2)%l_use) .OR. & 
    818760      &   (il_rho(jp_K) /= 1 .AND. td_var%t_dim(3)%l_use) )THEN 
     761 
     762         SELECT CASE(TRIM(td_var%c_interp(1))) 
     763            CASE('cubic','linear','nearest') 
     764               cl_method=TRIM(td_var%c_interp(1)) 
     765            CASE DEFAULT 
     766               CALL logger_warn("INTERP FILL VALUE: interpolation method unknown."//& 
     767               &  " use linear interpolation") 
     768               cl_method='linear' 
     769               ! update variable structure value 
     770               td_var%c_interp(1)='linear' 
     771         END SELECT 
    819772 
    820773         CALL logger_info("INTERP FILL: interpolate "//TRIM(td_var%c_name)//& 
     
    826779  
    827780         CALL interp__fill_value( td_var, cl_method, &  
    828          &                        il_rho(:), & 
    829          &                        il_offset(:,:) ) 
     781         &                        il_rho(:), il_offset(:,:) ) 
    830782 
    831783         SELECT CASE(TRIM(td_var%c_interp(2))) 
     
    864816         END SELECT 
    865817 
     818      ELSE 
     819         td_var%c_interp(:)='' 
    866820      ENDIF 
    867821 
     822      DEALLOCATE(il_rho) 
     823 
    868824   END SUBROUTINE interp__fill_value_wrapper 
    869    !> @endcode 
    870825   !------------------------------------------------------------------- 
    871826   !> @brief 
    872827   !> This subroutine interpolate value over mixed grid. 
    873828   !>  
    874    !> @details  
    875    !>  
    876    !> 
    877829   !> @author J.Paul 
    878    !> - Nov, 2013- Initial Version 
    879    ! 
    880    !> @param[inout] td_var : variable 
    881    !> @param[in] cd_method : interpolation method  
    882    !> @param[in] id_rho : table of refinment factor  
    883    !------------------------------------------------------------------- 
    884    !> @code 
     830   !> - November, 2013- Initial Version 
     831   !> @date September, 2014 
     832   !> - use interpolation method modules 
     833   !> 
     834   !> @param[inout] td_var variable structure 
     835   !> @param[in] cd_method interpolation method  
     836   !> @param[in] id_rho    array of refinment factor  
     837   !> @param[in] id_offset 2D array of offset between fine and coarse grid 
     838   !------------------------------------------------------------------- 
    885839   SUBROUTINE interp__fill_value( td_var, cd_method, & 
    886    &                              id_rho, & 
    887    &                              id_offset ) 
     840   &                              id_rho, id_offset ) 
    888841      IMPLICIT NONE 
    889842      ! Argument 
     
    891844      CHARACTER(LEN=*)                , INTENT(IN   ) :: cd_method 
    892845      INTEGER(I4)     , DIMENSION(:)  , INTENT(IN   ) :: id_rho 
    893       INTEGER(I4)     , DIMENSION(:,:), INTENT(IN   ) :: id_offset 
     846      INTEGER(I4)     , DIMENSION(2,2), INTENT(IN   ) :: id_offset 
    894847 
    895848      ! local variable 
     
    908861 
    909862      TYPE(TATT)                                       :: tl_att 
    910  
     863  
    911864      ! loop indices 
    912       INTEGER(i4) :: ji 
    913       INTEGER(i4) :: jj 
    914       INTEGER(i4) :: jk 
    915       INTEGER(i4) :: jl 
    916865      !---------------------------------------------------------------- 
    917866 
     
    935884      tl_att=att_init('interpolation',cl_interp) 
    936885      CALL var_move_att(tl_mix, tl_att) 
     886 
     887      ! clean 
     888      CALL att_clean(tl_att) 
    937889 
    938890      ! special case for even refinement on ARAKAWA-C grid 
     
    972924 
    973925      !3- interpolate 
    974       DO jl=1,tl_mix%t_dim(4)%i_len 
    975          IF( il_rho(jp_K) /= 1 )THEN 
    976 !            CALL interp__3D(tl_mix%d_value(:,:,:,jl), tl_mix%d_fill, & 
    977 !              &             il_detect(:,:,:), cd_method,             & 
    978 !              &             il_rhoi, il_rhoj, il_rhok,               & 
    979 !              &             ll_even(:), ll_discont ) 
    980              CALL logger_error("INTERP FILL: can not interpolate "//& 
    981              &                "vertically for now ") 
    982          ENDIF 
    983  
    984          IF( ANY(il_detect(:,:,:)==1) )THEN 
    985             ! I-J plan 
    986             DO jk=1,tl_mix%t_dim(3)%i_len 
    987               CALL interp__2D(tl_mix%d_value(:,:,jk,jl), tl_mix%d_fill, & 
    988               &               il_detect(:,:,jk), cd_method,             & 
    989               &               il_rho(jp_I), il_rho(jp_J),               & 
    990               &               ll_even(1:2), ll_discont) 
    991             ENDDO 
    992             IF( ALL(il_detect(:,:,:)==0) ) CYCLE 
    993             IF( il_rho(jp_K) /= 1 )THEN 
    994                ! I-K plan 
    995                DO jj=1,tl_mix%t_dim(2)%i_len 
    996                   CALL interp__2D(tl_mix%d_value(:,jj,:,jl), tl_mix%d_fill, & 
    997                   &               il_detect(:,jj,:), cd_method,             & 
    998                   &               il_rho(jp_J), il_rho(jp_K),               & 
    999                   &               ll_even(1:3:2), ll_discont  ) 
    1000                ENDDO 
    1001                IF( ALL(il_detect(:,:,:)==0) ) CYCLE 
    1002                ! J-K plan 
    1003                DO ji=1,tl_mix%t_dim(1)%i_len 
    1004                   CALL interp__2D(tl_mix%d_value(ji,:,:,jl), tl_mix%d_fill, & 
    1005                   &               il_detect(ji,:,:), cd_method,             & 
    1006                   &               il_rho(jp_J), il_rho(jp_K),               & 
    1007                   &               ll_even(2:3), ll_discont ) 
    1008                ENDDO 
    1009  
    1010             ENDIF 
    1011             IF( ANY(il_detect(:,:,:)==1) )THEN 
    1012                ! I direction 
    1013                DO jk=1,tl_mix%t_dim(3)%i_len 
    1014                   DO jj=1,tl_mix%t_dim(2)%i_len 
    1015                      CALL interp__1D( tl_mix%d_value(:,jj,jk,jl),    & 
    1016                      &                tl_mix%d_fill,                 & 
    1017                      &                il_detect(:,jj,jk), cd_method, & 
    1018                      &                il_rho(jp_I),                  & 
    1019                      &                ll_even(1), ll_discont ) 
    1020                   ENDDO 
    1021                ENDDO 
    1022                IF( ALL(il_detect(:,:,:)==0) ) CYCLE 
    1023                ! J direction 
    1024                DO jk=1,tl_mix%t_dim(3)%i_len 
    1025                   DO ji=1,tl_mix%t_dim(1)%i_len 
    1026                      CALL interp__1D( tl_mix%d_value(ji,:,jk,jl),    & 
    1027                      &                tl_mix%d_fill,                 & 
    1028                      &                il_detect(ji,:,jk), cd_method, & 
    1029                      &                il_rho(jp_J),                  & 
    1030                      &                ll_even(2), ll_discont ) 
    1031                   ENDDO 
    1032                ENDDO 
    1033                IF( il_rho(jp_K) /= 1 )THEN 
    1034                   IF( ALL(il_detect(:,:,:)==0) ) CYCLE 
    1035                   ! K direction 
    1036                   DO jj=1,tl_mix%t_dim(2)%i_len 
    1037                      DO ji=1,tl_mix%t_dim(1)%i_len 
    1038                         CALL interp__1D( tl_mix%d_value(ji,jj,:,jl),    & 
    1039                         &                tl_mix%d_fill,                 &  
    1040                         &                il_detect(ji,jj,:), cd_method, & 
    1041                         &                il_rho(jp_K),                  & 
    1042                         &                ll_even(3), ll_discont  ) 
    1043                      ENDDO 
    1044                   ENDDO 
    1045                ENDIF 
    1046             ENDIF 
    1047          ENDIF 
    1048       ENDDO 
     926      CALL logger_debug("INTERP 2D: interpolation method "//TRIM(cd_method)//& 
     927      &  " discont "//TRIM(fct_str(ll_discont)) ) 
     928      SELECT CASE(TRIM(cd_method)) 
     929      CASE('cubic') 
     930         CALL interp_cubic_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, & 
     931           &                    il_detect(:,:,:),                        & 
     932           &                    il_rho(:), ll_even(:), ll_discont ) 
     933      CASE('nearest') 
     934         CALL interp_nearest_fill(tl_mix%d_value(:,:,:,:), & 
     935              &                   il_detect(:,:,:),         & 
     936              &                   il_rho(:) ) 
     937      CASE DEFAULT ! linear 
     938         CALL interp_linear_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, & 
     939              &                  il_detect(:,:,:),                        & 
     940              &                  il_rho(:), ll_even(:), ll_discont ) 
     941      END SELECT          
    1049942 
    1050943      IF( ANY(il_detect(:,:,:)==1) )THEN 
     
    1056949      !4- save useful domain (remove offset) 
    1057950      CALL interp_clean_mixed_grid( tl_mix, td_var, & 
    1058       &                             id_rho(:),      & 
    1059       &                             id_offset(:,:) )       
     951      &                             id_rho(:), id_offset(:,:)  )   
    1060952 
    1061953      ! clean variable structure 
     954      DEALLOCATE(il_rho) 
    1062955      CALL var_clean(tl_mix) 
    1063956 
    1064957   END SUBROUTINE interp__fill_value 
    1065    !> @endcode 
    1066 !   !------------------------------------------------------------------- 
    1067 !   !> @brief 
    1068 !   !> This subroutine  
    1069 !   !>  
    1070 !   !> @details  
    1071 !   !> 
    1072 !   !> @author J.Paul 
    1073 !   !> - Nov, 2013- Initial Version 
    1074 !   ! 
    1075 !   !> @param[in] 
    1076 !   !> @param[out]  
    1077 !   !------------------------------------------------------------------- 
    1078 !   !> @code 
    1079 !   SUBROUTINE interp__3D( dd_value,  dd_fill,   & 
    1080 !      &                   id_detect, cd_method, & 
    1081 !      &                   id_rhoi, id_rhoj, id_rhok, & 
    1082 !      &                   ld_even,  ld_discont ) 
    1083 !      IMPLICIT NONE 
    1084 !      ! Argument 
    1085 !      REAL(dp)        , DIMENSION(:,:,:), INTENT(INOUT) :: dd_value  
    1086 !      REAL(dp)                          , INTENT(IN   ) :: dd_fill  
    1087 !      INTEGER(I4)     , DIMENSION(:,:,:), INTENT(INOUT) :: id_detect 
    1088 !      CHARACTER(LEN=*)                  , INTENT(IN   ) :: cd_method 
    1089 !      INTEGER(I4)                       , INTENT(IN   ) :: id_rhoi 
    1090 !      INTEGER(I4)                       , INTENT(IN   ) :: id_rhoj 
    1091 !      INTEGER(I4)                       , INTENT(IN   ) :: id_rhok 
    1092 !      LOGICAL         , DIMENSION(:)    , INTENT(IN   ) :: ld_even 
    1093 !      LOGICAL                           , INTENT(IN   ), OPTIONAL :: ld_discont 
    1094 ! 
    1095 !      ! local variable 
    1096 !      INTEGER(i4), DIMENSION(3)                  :: il_shape 
    1097 !      INTEGER(i4), DIMENSION(3)                  :: il_dim 
    1098 ! 
    1099 !      REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_coarse 
    1100 !      REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx 
    1101 !      REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdy 
    1102 !      REAL(dp)   , DIMENSION(:,:,:), ALLOCATABLE :: dl_d2fdxy 
    1103 ! 
    1104 !      LOGICAL                                    :: ll_discont 
    1105 ! 
    1106 !      ! loop indices 
    1107 !!      INTEGER(i4) :: ji 
    1108 !!      INTEGER(i4) :: jj 
    1109 !!      INTEGER(i4) :: ii 
    1110 !!      INTEGER(i4) :: ij 
    1111 !      !---------------------------------------------------------------- 
    1112 !      ll_discont=.FALSE. 
    1113 !      IF( PRESENT(ld_discont) ) ll_discont=ld_discont 
    1114 ! 
    1115 !      il_shape(:)=SHAPE(dd_value) 
    1116 !      il_dim(:)=(il_shape(:)-1)/2 
    1117 ! 
    1118 !      ALLOCATE( dl_coarse(il_dim(1),il_dim(2),il_dim(3)) ) 
    1119 !      ! value on coarse grid 
    1120 !      dl_coarse(:,:,:)=dd_value( 1:il_shape(1):id_rhoi, & 
    1121 !      &                        1:il_shape(2):id_rhoj, & 
    1122 !      &                        1:il_shape(3):id_rhok ) 
    1123 !      SELECT CASE(TRIM(cd_method)) 
    1124 ! 
    1125 !         CASE('cubic') 
    1126 !             
    1127 !            ALLOCATE( dl_dfdx(  il_dim(1),il_dim(2),il_dim(3)), & 
    1128 !            &         dl_dfdy(  il_dim(1),il_dim(2),il_dim(3)), & 
    1129 !            &         dl_d2fdxy(il_dim(1),il_dim(2),il_dim(3)) ) 
    1130 ! 
    1131 !!            ! compute derivative on coarse grid 
    1132 !!            dl_dfdx(:,:)=extrap_deriv_2D(dl_coarse(:,:,:), dd_fill, 'I') 
    1133 !!            dl_dfdy(:,:)=extrap_deriv_2D(dl_coarse(:,:,:), dd_fill, 'J') 
    1134 !! 
    1135 !!            ! compute cross derivative on coarse grid 
    1136 !!            dl_d2fdxy(:,:)=extrap_deriv_2D(dl_dfdx(:,:,:), dd_fill, 'J') 
    1137 !! 
    1138 !!            DO jj=1,il_shape(2)-1,id_rhoj 
    1139 !!               ij=((jj-1)/id_rhoj)+1 
    1140 !!               DO ji=1,il_shape(1)-1,id_rhoi 
    1141 !!                  ii=((ji-1)/id_rhoi)+1 
    1142 !!             
    1143 !!                  IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE 
    1144 !!                  ! compute bicubic coefficient 
    1145 !!                  dl_coef(:)=interp__2D_cubic_coef(dl_coarse(ii:ii+1,ij:ij+1),& 
    1146 !!                  &                                dl_dfdx(  ii:ii+1,ij:ij+1),& 
    1147 !!                  &                                dl_dfdy(  ii:ii+1,ij:ij+1),& 
    1148 !!                  &                                dl_d2fdxy(ii:ii+1,ij:ij+1) ) 
    1149 !! 
    1150 !!                  ! compute value on detetected point 
    1151 !!                  CALL interp__2D_cubic_fill(dl_coef(:), & 
    1152 !!                  &                          dd_value( ji:ji+id_rhoi,   & 
    1153 !!                  &                                    jj:jj+id_rhoj ), & 
    1154 !!                  &                          id_detect(ji:ji+id_rhoi,   & 
    1155 !!                  &                                    jj:jj+id_rhoj )  ) 
    1156 !! 
    1157 !!               ENDDO 
    1158 !!            ENDDO 
    1159 ! 
    1160 !            DEALLOCATE( dl_dfdx, & 
    1161 !            &           dl_dfdy, & 
    1162 !            &           dl_d2fdxy ) 
    1163 ! 
    1164 !         CASE('nearest') 
    1165 ! 
    1166 !         CASE DEFAULT ! linear  
    1167 ! 
    1168 !!            DO jj=1,il_shape(2)-1,id_rhoj 
    1169 !!               ij=((jj-1)/id_rhoj)+1 
    1170 !!               DO ji=1,il_shape(1)-1,id_rhoi 
    1171 !!                  ii=((ji-1)/id_rhoi)+1             
    1172 !! 
    1173 !!                  IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE 
    1174 !!                  ! compute bilinear coefficient 
    1175 !!                  dl_coef(:)=interp__2D_linear_coef(dl_coarse(ii:ii+1,ij:ij+1)) 
    1176 !! 
    1177 !!                  ! compute value on detetected point 
    1178 !!                  CALL interp__2D_linear_fill(dl_coef(:), & 
    1179 !!                  &                          dd_value( ji:ji+id_rhoi,   & 
    1180 !!                  &                                    jj:jj+id_rhoj ), & 
    1181 !!                  &                          id_detect(ji:ji+id_rhoi,   & 
    1182 !!                  &                                    jj:jj+id_rhoj )  ) 
    1183 !! 
    1184 !!               ENDDO 
    1185 !!            ENDDO 
    1186 ! 
    1187 !      END SELECT 
    1188 ! 
    1189 !      DEALLOCATE( dl_coarse ) 
    1190 ! 
    1191 !   END SUBROUTINE interp__3D 
    1192 !   !> @endcode 
    1193    !------------------------------------------------------------------- 
    1194    !> @brief 
    1195    !> This subroutine  
    1196    !>  
    1197    !> @details  
    1198    !> 
    1199    !> @author J.Paul 
    1200    !> - Nov, 2013- Initial Version 
    1201    ! 
    1202    !> @param[in] 
    1203    !> @param[out]  
    1204    !------------------------------------------------------------------- 
    1205    !> @code 
    1206    SUBROUTINE interp__2D( dd_value,  dd_fill,   & 
    1207       &                   id_detect, cd_method, & 
    1208       &                   id_rhoi, id_rhoj,     & 
    1209       &                   ld_even,  ld_discont ) 
    1210  
    1211       IMPLICIT NONE 
    1212       ! Argument 
    1213       REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value  
    1214       REAL(dp)                        , INTENT(IN   ) :: dd_fill  
    1215       INTEGER(I4)     , DIMENSION(:,:), INTENT(INOUT) :: id_detect 
    1216       CHARACTER(LEN=*)                , INTENT(IN   ) :: cd_method 
    1217       INTEGER(I4)                     , INTENT(IN   ) :: id_rhoi 
    1218       INTEGER(I4)                     , INTENT(IN   ) :: id_rhoj 
    1219       LOGICAL         , DIMENSION(:)  , INTENT(IN   ) :: ld_even 
    1220       LOGICAL                         , INTENT(IN   ), OPTIONAL :: ld_discont 
    1221  
    1222       ! local variable 
    1223       INTEGER(I4)                              :: il_xextra 
    1224       INTEGER(I4)                              :: il_yextra 
    1225       INTEGER(i4), DIMENSION(2)                :: il_shape 
    1226       INTEGER(i4), DIMENSION(2)                :: il_dim 
    1227  
    1228       REAL(dp)                                 :: dl_min  
    1229       REAL(dp)                                 :: dl_max  
    1230       REAL(dp)   , DIMENSION(:)  , ALLOCATABLE :: dl_coef 
    1231       REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_coarse 
    1232       REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_tmp 
    1233       REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_dfdx 
    1234       REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_dfdy 
    1235       REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_d2fdxy 
    1236  
    1237       LOGICAL                                  :: ll_discont 
    1238       ! loop indices 
    1239       INTEGER(i4) :: ji 
    1240       INTEGER(i4) :: jj 
    1241       INTEGER(i4) :: ii 
    1242       INTEGER(i4) :: ij 
    1243  
    1244       !---------------------------------------------------------------- 
    1245       ll_discont=.FALSE. 
    1246       IF( PRESENT(ld_discont) ) ll_discont=ld_discont 
    1247  
    1248       CALL logger_debug("INTERP 2D: interpolation method "//TRIM(cd_method)//& 
    1249       &  " discont "//TRIM(fct_str(ll_discont)) ) 
    1250  
    1251       il_shape(:)=SHAPE(dd_value) 
    1252  
    1253       ! compute coarse grid dimension 
    1254       il_xextra=id_rhoi-1 
    1255       il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi 
    1256  
    1257       il_yextra=id_rhoj-1 
    1258       il_dim(2)=(il_shape(2)+il_yextra)/id_rhoj 
    1259  
    1260       ALLOCATE( dl_coarse(il_dim(1),il_dim(2)) ) 
    1261  
    1262       ! value on coarse grid 
    1263       dl_coarse(:,:)=dd_value( 1:il_shape(1):id_rhoi, & 
    1264       &                        1:il_shape(2):id_rhoj ) 
    1265  
    1266       SELECT CASE(TRIM(cd_method)) 
    1267  
    1268          CASE('cubic') 
    1269              
    1270             ALLOCATE( dl_dfdx(  il_dim(1),il_dim(2)), & 
    1271             &         dl_dfdy(  il_dim(1),il_dim(2)), & 
    1272             &         dl_d2fdxy(il_dim(1),il_dim(2)) ) 
    1273  
    1274             ! compute derivative on coarse grid 
    1275             dl_dfdx(:,:)=extrap_deriv_2D(dl_coarse(:,:), dd_fill, 'I', ll_discont) 
    1276             dl_dfdy(:,:)=extrap_deriv_2D(dl_coarse(:,:), dd_fill, 'J', ll_discont) 
    1277  
    1278             ! compute cross derivative on coarse grid 
    1279             dl_d2fdxy(:,:)=extrap_deriv_2D(dl_dfdx(:,:), dd_fill, 'J', ll_discont) 
    1280  
    1281             ALLOCATE( dl_tmp(2,2) ) 
    1282             ALLOCATE( dl_coef(16) ) 
    1283  
    1284             DO jj=1,il_shape(2)-1,id_rhoj 
    1285                ij=((jj-1)/id_rhoj)+1 
    1286                DO ji=1,il_shape(1)-1,id_rhoi 
    1287                   ii=((ji-1)/id_rhoi)+1 
    1288              
    1289                   IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill) .OR. & 
    1290                   &   ANY(  dl_dfdx(ii:ii+1,ij:ij+1)==dd_fill) .OR. & 
    1291                   &   ANY(  dl_dfdy(ii:ii+1,ij:ij+1)==dd_fill) .OR. & 
    1292                   &   ANY(dl_d2fdxy(ii:ii+1,ij:ij+1)==dd_fill) ) CYCLE 
    1293  
    1294                   dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1) 
    1295                   IF( ll_discont )THEN 
    1296  
    1297                      dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill ) 
    1298                      dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill ) 
    1299                      IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 
    1300                         WHERE( dl_tmp(:,:) < 0_dp )  
    1301                            dl_tmp(:,:) = dl_tmp(:,:)+360._dp 
    1302                         END WHERE 
    1303                      ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 
    1304                         WHERE( dl_tmp(:,:) > 180_dp )  
    1305                            dl_tmp(:,:) = dl_tmp(:,:)-180._dp 
    1306                         END WHERE 
    1307                      ENDIF 
    1308                   ENDIF 
    1309  
    1310                   ! compute bicubic coefficient 
    1311                   dl_coef(:)=interp__2D_cubic_coef(dl_tmp(:,:),& 
    1312                   &                                dl_dfdx(   ii:ii+1,ij:ij+1),& 
    1313                   &                                dl_dfdy(   ii:ii+1,ij:ij+1),& 
    1314                   &                                dl_d2fdxy( ii:ii+1,ij:ij+1),& 
    1315                   &                                dd_fill ) 
    1316  
    1317                   ! compute value on detetected point 
    1318                   CALL interp__2D_cubic_fill(dd_value( ji:ji+id_rhoi,   & 
    1319                   &                                    jj:jj+id_rhoj ), & 
    1320                   &                          id_detect(ji:ji+id_rhoi,   & 
    1321                   &                                    jj:jj+id_rhoj ), & 
    1322                   &                          dl_coef(:), dd_fill,       & 
    1323                   &                          ld_even(:), id_rhoi, id_rhoj ) 
    1324  
    1325                   IF( ll_discont )THEN 
    1326                      WHERE( dd_value( ji:ji+id_rhoi, & 
    1327                         &             jj:jj+id_rhoj ) >= 180._dp .AND. & 
    1328                         &   dd_value( ji:ji+id_rhoi, & 
    1329                         &             jj:jj+id_rhoj ) /= dd_fill ) 
    1330                         dd_value( ji:ji+id_rhoi, & 
    1331                         &         jj:jj+id_rhoj ) = & 
    1332                         &           dd_value( ji:ji+id_rhoi, & 
    1333                         &                     jj:jj+id_rhoj ) - 360._dp 
    1334                      END WHERE 
    1335                   ENDIF 
    1336  
    1337                ENDDO 
    1338             ENDDO 
    1339  
    1340             DEALLOCATE(dl_coef) 
    1341             DEALLOCATE(dl_tmp ) 
    1342  
    1343             DEALLOCATE(dl_dfdx, & 
    1344             &          dl_dfdy, & 
    1345             &          dl_d2fdxy ) 
    1346  
    1347          CASE('nearest')  
    1348  
    1349             DO jj=1,il_shape(2)-1,id_rhoj 
    1350                ij=((jj-1)/id_rhoj)+1 
    1351                DO ji=1,il_shape(1)-1,id_rhoi 
    1352                   ii=((ji-1)/id_rhoi)+1             
    1353  
    1354                   ! compute value on detetected point 
    1355                   CALL interp__2D_nearest_fill(dd_value( ji:ji+id_rhoi,   & 
    1356                   &                                      jj:jj+id_rhoj ), & 
    1357                   &                            id_detect(ji:ji+id_rhoi,   & 
    1358                   &                                      jj:jj+id_rhoj )  ) 
    1359  
    1360                ENDDO 
    1361             ENDDO 
    1362  
    1363          CASE DEFAULT ! linear  
    1364  
    1365             ALLOCATE( dl_coef(4)  ) 
    1366             ALLOCATE( dl_tmp(2,2) ) 
    1367  
    1368             DO jj=1,il_shape(2)-1,id_rhoj 
    1369                ij=((jj-1)/id_rhoj)+1 
    1370                DO ji=1,il_shape(1)-1,id_rhoi 
    1371                   ii=((ji-1)/id_rhoi)+1             
    1372  
    1373                   IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE 
    1374  
    1375                   dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1) 
    1376                   IF( ll_discont )THEN 
    1377  
    1378                      dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill ) 
    1379                      dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill ) 
    1380                      IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 
    1381                         WHERE( dl_tmp(:,:) < 0_dp )  
    1382                            dl_tmp(:,:) = dl_tmp(:,:)+360._dp 
    1383                         END WHERE 
    1384                      ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 
    1385                         WHERE( dl_tmp(:,:) > 180_dp )  
    1386                            dl_tmp(:,:) = dl_tmp(:,:)-180._dp 
    1387                         END WHERE 
    1388                      ENDIF 
    1389                   ENDIF 
    1390  
    1391                   ! compute bilinear coefficient 
    1392                   dl_coef(:)=interp__2D_linear_coef(dl_tmp(:,:), dd_fill) 
    1393  
    1394                   ! compute value on detetected point 
    1395                   CALL interp__2D_linear_fill(dd_value( ji:ji+id_rhoi,   & 
    1396                   &                                     jj:jj+id_rhoj ), & 
    1397                   &                           id_detect(ji:ji+id_rhoi,   & 
    1398                   &                                     jj:jj+id_rhoj ), & 
    1399                   &                           dl_coef(:), dd_fill,       & 
    1400                   &                           ld_even(:), id_rhoi, id_rhoj ) 
    1401  
    1402                   IF( ll_discont )THEN 
    1403                      WHERE( dd_value( ji:ji+id_rhoi, & 
    1404                         &             jj:jj+id_rhoj ) >= 180._dp .AND. & 
    1405                         &   dd_value( ji:ji+id_rhoi, & 
    1406                         &             jj:jj+id_rhoj ) /= dd_fill ) 
    1407                         dd_value( ji:ji+id_rhoi, & 
    1408                         &         jj:jj+id_rhoj ) = & 
    1409                         &           dd_value( ji:ji+id_rhoi, & 
    1410                         &                     jj:jj+id_rhoj ) - 360._dp 
    1411                      END WHERE 
    1412                   ENDIF 
    1413  
    1414                ENDDO 
    1415             ENDDO 
    1416  
    1417             DEALLOCATE(dl_coef) 
    1418             DEALLOCATE(dl_tmp ) 
    1419  
    1420       END SELECT 
    1421  
    1422       DEALLOCATE( dl_coarse ) 
    1423  
    1424    END SUBROUTINE interp__2D 
    1425    !> @endcode 
    1426    !------------------------------------------------------------------- 
    1427    !> @brief 
    1428    !> This subroutine  
    1429    !>  
    1430    !> @details  
    1431    !> 
    1432    !> @author J.Paul 
    1433    !> - Nov, 2013- Initial Version 
    1434    ! 
    1435    !> @param[in] 
    1436    !> @param[out]  
    1437    !------------------------------------------------------------------- 
    1438    !> @code 
    1439    SUBROUTINE interp__1D( dd_value,  dd_fill,   & 
    1440       &                   id_detect, cd_method, & 
    1441       &                   id_rhoi,              & 
    1442       &                   ld_even, ld_discont ) 
    1443       IMPLICIT NONE 
    1444       ! Argument 
    1445       REAL(dp)        , DIMENSION(:), INTENT(INOUT) :: dd_value  
    1446       REAL(dp)                      , INTENT(IN   ) :: dd_fill  
    1447       INTEGER(I4)     , DIMENSION(:), INTENT(INOUT) :: id_detect 
    1448       CHARACTER(LEN=*)              , INTENT(IN   ) :: cd_method 
    1449       INTEGER(I4)                   , INTENT(IN   ) :: id_rhoi 
    1450       LOGICAL                       , INTENT(IN   ) :: ld_even 
    1451       LOGICAL                       , INTENT(IN   ), OPTIONAL :: ld_discont 
    1452  
    1453       ! local variable 
    1454       INTEGER(i4), DIMENSION(1)              :: il_shape 
    1455       INTEGER(i4), DIMENSION(1)              :: il_dim 
    1456       INTEGER(I4)                            :: il_xextra 
    1457  
    1458       REAL(dp)                               :: dl_min  
    1459       REAL(dp)                               :: dl_max  
    1460       REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_coarse 
    1461       REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_tmp 
    1462       REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_dfdx 
    1463  
    1464       REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_coef 
    1465  
    1466       LOGICAL                                :: ll_discont 
    1467  
    1468       ! loop indices 
    1469       INTEGER(i4) :: ji 
    1470       INTEGER(i4) :: ii 
    1471       !---------------------------------------------------------------- 
    1472       ll_discont=.FALSE. 
    1473       IF( PRESENT(ld_discont) ) ll_discont=ld_discont 
    1474  
    1475       CALL logger_debug("INTERP 1D: interpolation method "//TRIM(cd_method)//& 
    1476       &  " discont "//TRIM(fct_str(ll_discont)) ) 
    1477  
    1478       il_shape(:)=SHAPE(dd_value) 
    1479  
    1480       il_xextra=id_rhoi-1 
    1481       il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi 
    1482  
    1483       ALLOCATE( dl_coarse(il_dim(1)) ) 
    1484       ! value on coarse grid 
    1485       dl_coarse(:)=dd_value( 1:il_shape(1):id_rhoi ) 
    1486  
    1487       SELECT CASE(TRIM(cd_method)) 
    1488  
    1489          CASE('cubic') 
    1490              
    1491             ALLOCATE( dl_dfdx( il_dim(1)) ) 
    1492  
    1493             ! compute derivative on coarse grid 
    1494             dl_dfdx(:)=extrap_deriv_1D(dl_coarse(:), dd_fill, ll_discont) 
    1495  
    1496             ALLOCATE( dl_coef(4)) 
    1497             ALLOCATE( dl_tmp(2) ) 
    1498  
    1499             DO ji=1,il_shape(1)-1,id_rhoi 
    1500                ii=((ji-1)/id_rhoi)+1 
    1501              
    1502                IF( ANY( dl_tmp(:)==dd_fill ) .OR. & 
    1503                &   ANY(dl_dfdx(:)==dd_fill ) ) CYCLE 
    1504  
    1505                dl_tmp(:)=dl_coarse(ii:ii+1) 
    1506                IF( ll_discont )THEN 
    1507  
    1508                   dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill ) 
    1509                   dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill ) 
    1510                   IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 
    1511                      WHERE( dl_tmp(:) < 0_dp )  
    1512                         dl_tmp(:) = dl_tmp(:)+360._dp 
    1513                      END WHERE 
    1514                   ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 
    1515                      WHERE( dl_tmp(:) > 180_dp )  
    1516                         dl_tmp(:) = dl_tmp(:)-180._dp 
    1517                      END WHERE 
    1518                   ENDIF 
    1519                ENDIF 
    1520  
    1521                ! compute bicubic coefficient 
    1522                dl_coef(:)=interp__1D_cubic_coef(dl_tmp(:),& 
    1523                &                                dl_dfdx(   ii:ii+1), & 
    1524                &                                dd_fill ) 
    1525  
    1526                ! compute value on detetected point 
    1527                CALL interp__1D_cubic_fill(dd_value( ji:ji+id_rhoi),& 
    1528                &                          id_detect(ji:ji+id_rhoi),& 
    1529                &                          dl_coef(:), dd_fill,     & 
    1530                &                          ld_even, id_rhoi ) 
    1531  
    1532                IF( ll_discont )THEN 
    1533                   WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. & 
    1534                      &   dd_value( ji:ji+id_rhoi ) /= dd_fill ) 
    1535                      dd_value( ji:ji+id_rhoi ) = & 
    1536                      &           dd_value( ji:ji+id_rhoi ) - 360._dp 
    1537                   END WHERE 
    1538                ENDIF 
    1539  
    1540             ENDDO 
    1541  
    1542             DEALLOCATE(dl_coef) 
    1543             DEALLOCATE(dl_tmp ) 
    1544  
    1545          CASE('nearest')  
    1546  
    1547             DO ji=1,il_shape(1)-1,id_rhoi 
    1548                ii=((ji-1)/id_rhoi)+1             
    1549  
    1550                ! compute value on detetected point 
    1551                CALL interp__1D_nearest_fill(dd_value( ji:ji+id_rhoi), & 
    1552                &                            id_detect(ji:ji+id_rhoi)  ) 
    1553  
    1554             ENDDO 
    1555  
    1556          CASE DEFAULT ! linear  
    1557  
    1558             ALLOCATE(dl_coef(2)) 
    1559             ALLOCATE( dl_tmp(2) ) 
    1560  
    1561             DO ji=1,il_shape(1)-1,id_rhoi 
    1562                ii=((ji-1)/id_rhoi)+1             
    1563  
    1564                IF( ANY(dl_coarse(ii:ii+1)==dd_fill ) ) CYCLE 
    1565  
    1566                dl_tmp(:)=dl_coarse(ii:ii+1) 
    1567                IF( ll_discont )THEN 
    1568  
    1569                   dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill ) 
    1570                   dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill ) 
    1571                   IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN 
    1572                      WHERE( dl_tmp(:) < 0_dp )  
    1573                         dl_tmp(:) = dl_tmp(:)+360._dp 
    1574                      END WHERE 
    1575                   ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN 
    1576                      WHERE( dl_tmp(:) > 180_dp )  
    1577                         dl_tmp(:) = dl_tmp(:)-180._dp 
    1578                      END WHERE 
    1579                   ENDIF 
    1580                ENDIF 
    1581  
    1582                ! compute bilinear coefficient 
    1583                dl_coef(:)=interp__1D_linear_coef(dl_tmp(:), dd_fill) 
    1584  
    1585                ! compute value on detetected point 
    1586                CALL interp__1D_linear_fill( dd_value( ji:ji+id_rhoi),& 
    1587                &                            id_detect(ji:ji+id_rhoi),& 
    1588                &                            dl_coef(:), dd_fill,     & 
    1589                &                            ld_even, id_rhoi  ) 
    1590  
    1591                IF( ll_discont )THEN 
    1592                   WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. & 
    1593                      &   dd_value( ji:ji+id_rhoi ) /= dd_fill ) 
    1594                      dd_value( ji:ji+id_rhoi ) = & 
    1595                      &           dd_value( ji:ji+id_rhoi ) - 360._dp 
    1596                   END WHERE 
    1597                ENDIF 
    1598  
    1599             ENDDO 
    1600  
    1601             DEALLOCATE(dl_coef) 
    1602  
    1603       END SELECT 
    1604  
    1605       DEALLOCATE( dl_coarse ) 
    1606  
    1607    END SUBROUTINE interp__1D 
    1608    !> @endcode 
    1609    !------------------------------------------------------------------- 
    1610    !> @brief 
    1611    !> This subroutine  
    1612    !>  
    1613    !> @details  
    1614    !> 
    1615    !> @author J.Paul 
    1616    !> - Nov, 2013- Initial Version 
    1617    ! 
    1618    !> @param[in] 
    1619    !> @param[out]  
    1620    !------------------------------------------------------------------- 
    1621    !> @code 
    1622    FUNCTION interp__2D_cubic_coef( dd_value, & 
    1623       &                            dd_dfdx,  & 
    1624       &                            dd_dfdy,  & 
    1625       &                            dd_d2fdxy,& 
    1626       &                            dd_fill ) 
    1627       IMPLICIT NONE 
    1628       ! Argument 
    1629       REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_value  
    1630       REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdx   
    1631       REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdy   
    1632       REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_d2fdxy 
    1633       REAL(dp)                , INTENT(IN) :: dd_fill 
    1634  
    1635       ! function 
    1636       REAL(dp), DIMENSION(16) :: interp__2D_cubic_coef 
    1637  
    1638       ! local variable 
    1639       REAL(dp), DIMENSION(16,16), PARAMETER :: dp_matrix = RESHAPE( & 
    1640       & (/ 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&  
    1641            0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&  
    1642           -3 , 3 , 0 , 0 ,-2 ,-1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&  
    1643            2 ,-2 , 0 , 0 , 1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&  
    1644            0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&  
    1645            0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 ,&  
    1646            0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,-3 , 3 , 0 , 0 ,-2 ,-1 , 0 , 0 ,&  
    1647            0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 ,-2 , 0 , 0 , 1 , 1 , 0 , 0 ,&  
    1648           -3 , 0 , 3 , 0 , 0 , 0 , 0 , 0 ,-2 , 0 ,-1 , 0 , 0 , 0 , 0 , 0 ,&  
    1649            0 , 0 , 0 , 0 ,-3 , 0 , 3 , 0 , 0 , 0 , 0 , 0 ,-2 , 0 ,-1 , 0 ,&  
    1650            9 ,-9 ,-9 , 9 , 6 , 3 ,-6 ,-3 , 6 ,-6 , 3 ,-3 , 4 , 2 , 2 , 1 ,&  
    1651           -6 , 6 , 6 ,-6 ,-3 ,-3 , 3 , 3 ,-4 , 4 ,-2 , 2 ,-2 ,-2 ,-1 ,-1 ,&  
    1652            2 , 0 ,-2 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 ,&  
    1653            0 , 0 , 0 , 0 , 2 , 0 ,-2 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 ,&  
    1654           -6 , 6 , 6 ,-6 ,-4 ,-2 , 4 , 2 ,-3 , 3 ,-3 , 3 ,-2 ,-1 ,-2 ,-1 ,&  
    1655            4 ,-4 ,-4 , 4 , 2 , 2 ,-2 ,-2 , 2 ,-2 , 2 ,-2 , 1 , 1 , 1 , 1 /), & 
    1656       & (/ 16, 16 /) ) 
    1657  
    1658       REAL(dp), DIMENSION(16) :: dl_vect 
    1659  
    1660       !---------------------------------------------------------------- 
    1661       ! init 
    1662       interp__2D_cubic_coef(:)=dd_fill 
    1663  
    1664       IF( ANY(SHAPE( dd_value(:,:))/= 2) .OR. & 
    1665       &   ANY(SHAPE(  dd_dfdx(:,:))/= 2) .OR. & 
    1666       &   ANY(SHAPE(  dd_dfdy(:,:))/= 2) .OR. & 
    1667       &   ANY(SHAPE(dd_d2fdxy(:,:))/= 2) )THEN 
    1668     
    1669          CALL logger_error("INTERP CUBIC COEF: invalid dimension of "//& 
    1670          &              "input tables. shape should be (/2,2/)") 
    1671  
    1672       ELSEIF( ANY( dd_value(:,:) == dd_fill) .OR. & 
    1673       &       ANY(  dd_dfdx(:,:) == dd_fill) .OR. & 
    1674       &       ANY(  dd_dfdy(:,:) == dd_fill) .OR. & 
    1675       &       ANY(dd_d2fdxy(:,:) == dd_fill) )THEN 
    1676        
    1677          CALL logger_warn("INTERP CUBIC COEF: fill value detected. "//& 
    1678          &  "can not compute coefficient ") 
    1679  
    1680       ELSE 
    1681  
    1682          dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. ) 
    1683          dl_vect( 5: 8)=PACK(dd_dfdx(:,:),.TRUE. ) 
    1684          dl_vect( 9:12)=PACK(dd_dfdy(:,:),.TRUE. ) 
    1685          dl_vect(13:16)=PACK(dd_d2fdxy(:,:),.TRUE. ) 
    1686  
    1687          interp__2D_cubic_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:)) 
    1688  
    1689       ENDIF 
    1690    END FUNCTION interp__2D_cubic_coef 
    1691    !> @endcode 
    1692    !------------------------------------------------------------------- 
    1693    !> @brief 
    1694    !> This subroutine  
    1695    !>  
    1696    !> @details  
    1697    !> 
    1698    !> @author J.Paul 
    1699    !> - Nov, 2013- Initial Version 
    1700    ! 
    1701    !> @param[in] 
    1702    !> @param[out]  
    1703    !------------------------------------------------------------------- 
    1704    !> @code 
    1705    SUBROUTINE interp__2D_cubic_fill( dd_value, id_detect, dd_coef, dd_fill, & 
    1706    &                                 ld_even, id_rhoi, id_rhoj ) 
    1707       IMPLICIT NONE 
    1708       ! Argument 
    1709       REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value  
    1710       INTEGER(i4)     , DIMENSION(:,:), INTENT(INOUT) :: id_detect 
    1711       REAL(dp)        , DIMENSION(:)  , INTENT(IN   ) :: dd_coef  
    1712       REAL(dp)                        , INTENT(IN   ) :: dd_fill  
    1713       LOGICAL         , DIMENSION(:),   INTENT(IN   ) :: ld_even 
    1714       INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoi 
    1715       INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoj 
    1716  
    1717       ! local variable 
    1718       INTEGER(i4), DIMENSION(2)  :: il_shape 
    1719  
    1720       REAL(dp)   , DIMENSION(16) :: dl_vect 
    1721       REAL(dp)                   :: dl_dx 
    1722       REAL(dp)                   :: dl_dy 
    1723       REAL(dp)                   :: dl_x 
    1724       REAL(dp)                   :: dl_x2 
    1725       REAL(dp)                   :: dl_x3 
    1726       REAL(dp)                   :: dl_y 
    1727       REAL(dp)                   :: dl_y2 
    1728       REAL(dp)                   :: dl_y3 
    1729  
    1730       ! loop indices 
    1731       INTEGER(i4) :: ji 
    1732       INTEGER(i4) :: jj 
    1733       !---------------------------------------------------------------- 
    1734  
    1735       IF( SIZE(dd_coef(:)) /= 16 )THEN 
    1736          CALL logger_error("INTERP CUBIC FILL: invalid dimension of "//& 
    1737          &              "coef table. shape should be (/16/)") 
    1738       ELSEIF( ANY(   dd_coef(:)==dd_fill ) )THEN 
    1739          CALL logger_error("INTERP CUBIC FILL: fill value detected in coef . "//& 
    1740          &              "can not compute interpolation.") 
    1741       ELSE 
    1742  
    1743          il_shape(:)=SHAPE(dd_value(:,:)) 
    1744  
    1745          dl_dx=1./REAL(id_rhoi) 
    1746          dl_dy=1./REAL(id_rhoj) 
    1747  
    1748          DO jj=1,il_shape(2) 
    1749  
    1750             IF( ld_even(2) )THEN 
    1751                dl_y=(jj-1)*dl_dy - dl_dy*0.5 
    1752             ELSE ! odd refinement 
    1753                dl_y=(jj-1)*dl_dy 
    1754             ENDIF 
    1755             dl_y2=dl_y*dl_y 
    1756             dl_y3=dl_y2*dl_y 
    1757  
    1758             DO ji=1,il_shape(1) 
    1759                 
    1760                IF(id_detect(ji,jj)==1)THEN 
    1761                    
    1762                   IF( ld_even(1) )THEN 
    1763                      dl_x=(ji-1)*dl_dx - dl_dx*0.5  
    1764                   ELSE ! odd refinement 
    1765                      dl_x=(ji-1)*dl_dx  
    1766                   ENDIF 
    1767                   dl_x2=dl_x*dl_x 
    1768                   dl_x3=dl_x2*dl_x 
    1769  
    1770                   dl_vect(:)=(/1._dp, dl_x      , dl_x2      , dl_x3      , & 
    1771                   &            dl_y , dl_x*dl_y , dl_x2*dl_y , dl_x3*dl_y , & 
    1772                   &            dl_y2, dl_x*dl_y2, dl_x2*dl_y2, dl_x3*dl_y2, & 
    1773                   &            dl_y3, dl_x*dl_y3, dl_x2*dl_y3, dl_x3*dl_y3 /) 
    1774  
    1775                   dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dl_vect(:)) 
    1776                   id_detect(ji,jj)=0 
    1777  
    1778                ENDIF 
    1779  
    1780             ENDDO 
    1781          ENDDO 
    1782  
    1783       ENDIF 
    1784  
    1785    END SUBROUTINE interp__2D_cubic_fill 
    1786    !> @endcode 
    1787    !------------------------------------------------------------------- 
    1788    !> @brief 
    1789    !> This subroutine  
    1790    !>  
    1791    !> @details  
    1792    !> 
    1793    !> @author J.Paul 
    1794    !> - Nov, 2013- Initial Version 
    1795    ! 
    1796    !> @param[in] 
    1797    !> @param[out]  
    1798    !------------------------------------------------------------------- 
    1799    !> @code 
    1800    FUNCTION interp__2D_linear_coef( dd_value, dd_fill ) 
    1801       IMPLICIT NONE 
    1802       ! Argument 
    1803       REAL(dp), DIMENSION(:,:)  , INTENT(IN) :: dd_value  
    1804       REAL(dp)                  , INTENT(IN) :: dd_fill  
    1805  
    1806       ! function 
    1807       REAL(dp), DIMENSION(4) :: interp__2D_linear_coef 
    1808  
    1809       ! local variable 
    1810  
    1811       REAL(dp), DIMENSION(4,4), PARAMETER :: dp_matrix = RESHAPE( & 
    1812       & (/  1 , 0 , 0 , 0 ,&  
    1813            -1 , 1 , 0 , 0 ,&  
    1814            -1 , 0 , 1 , 0 ,&  
    1815             1 ,-1 ,-1 , 1 /), & 
    1816       & (/ 4, 4 /) ) 
    1817  
    1818       REAL(dp), DIMENSION(4) :: dl_vect 
    1819  
    1820       !---------------------------------------------------------------- 
    1821  
    1822       IF( ANY(SHAPE(dd_value(:,:))/= 2) )THEN 
    1823          CALL logger_error("INTERP LINEAR COEF: invalid dimension of "//& 
    1824          &              "input tables. shape should be (/2,2/)") 
    1825       ELSEIF( ANY(dd_value(:,:)==dd_fill) )THEN 
    1826          CALL logger_error("INTERP LINEAR COEF: fill value detected. "//& 
    1827          &              "can not compute coefficient.") 
    1828       ELSE 
    1829  
    1830          dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. ) 
    1831  
    1832          interp__2D_linear_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:)) 
    1833  
    1834       ENDIF 
    1835    END FUNCTION interp__2D_linear_coef 
    1836    !> @endcode 
    1837    !------------------------------------------------------------------- 
    1838    !> @brief 
    1839    !> This subroutine  
    1840    !>  
    1841    !> @details  
    1842    !> 
    1843    !> @author J.Paul 
    1844    !> - Nov, 2013- Initial Version 
    1845    ! 
    1846    !> @param[in] 
    1847    !> @param[out]  
    1848    !------------------------------------------------------------------- 
    1849    !> @code 
    1850    SUBROUTINE interp__2D_linear_fill( dd_value, id_detect, dd_coef, dd_fill, & 
    1851    &                                  ld_even, id_rhoi, id_rhoj ) 
    1852       IMPLICIT NONE 
    1853       ! Argument 
    1854       REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value  
    1855       INTEGER(i4)     , DIMENSION(:,:), INTENT(INOUT) :: id_detect 
    1856       REAL(dp)        , DIMENSION(:)  , INTENT(IN   ) :: dd_coef  
    1857       REAL(dp)                        , INTENT(IN   ) :: dd_fill  
    1858       LOGICAL         , DIMENSION(:),   INTENT(IN   ) :: ld_even 
    1859       INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoi 
    1860       INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoj 
    1861  
    1862       ! local variable 
    1863       INTEGER(i4), DIMENSION(2) :: il_shape 
    1864  
    1865       REAL(dp)   , DIMENSION(4) :: dl_vect 
    1866       REAL(dp)                  :: dl_dx 
    1867       REAL(dp)                  :: dl_dy 
    1868       REAL(dp)                  :: dl_x 
    1869       REAL(dp)                  :: dl_y 
    1870  
    1871       ! loop indices 
    1872       INTEGER(i4) :: ji 
    1873       INTEGER(i4) :: jj 
    1874       !---------------------------------------------------------------- 
    1875  
    1876       IF( SIZE(dd_coef(:)) /= 4 )THEN 
    1877          CALL logger_error("INTERP LINEAR FILL: invalid dimension of "//& 
    1878          &              "coef table. shape should be (/4/)") 
    1879       ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN 
    1880          CALL logger_error("INTERP LINEAR FILL: fill value detected in coef. "//& 
    1881          &              "can not compute interpolation.") 
    1882       ELSE 
    1883  
    1884          il_shape(:)=SHAPE(dd_value(:,:)) 
    1885  
    1886          dl_dx=1./REAL(id_rhoi) 
    1887          dl_dy=1./REAL(id_rhoj) 
    1888  
    1889          DO jj=1,il_shape(2) 
    1890  
    1891             IF( ld_even(2) )THEN 
    1892                dl_y=(jj-1)*dl_dy - dl_dy*0.5 
    1893             ELSE ! odd refinement 
    1894                dl_y=(jj-1)*dl_dy 
    1895             ENDIF 
    1896  
    1897             DO ji=1,il_shape(1) 
    1898                 
    1899                IF(id_detect(ji,jj)==1)THEN 
    1900  
    1901                   IF( ld_even(1) )THEN 
    1902                      dl_x=(ji-1)*dl_dx - dl_dx*0.5  
    1903                   ELSE ! odd refinement 
    1904                      dl_x=(ji-1)*dl_dx  
    1905                   ENDIF 
    1906  
    1907                   dl_vect(:)=(/1._dp, dl_x, dl_y, dl_x*dl_y /) 
    1908  
    1909                   dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dl_vect(:)) 
    1910                   id_detect(ji,jj)=0 
    1911  
    1912                ENDIF 
    1913  
    1914             ENDDO 
    1915          ENDDO 
    1916  
    1917       ENDIF 
    1918  
    1919    END SUBROUTINE interp__2D_linear_fill 
    1920    !> @endcode 
    1921    !------------------------------------------------------------------- 
    1922    !> @brief 
    1923    !> This subroutine  
    1924    !>  
    1925    !> @details  
    1926    !> 
    1927    !> @author J.Paul 
    1928    !> - Nov, 2013- Initial Version 
    1929    ! 
    1930    !> @param[in] 
    1931    !> @param[out]  
    1932    !------------------------------------------------------------------- 
    1933    !> @code 
    1934    SUBROUTINE interp__2D_nearest_fill( dd_value, id_detect ) 
    1935       IMPLICIT NONE 
    1936       ! Argument 
    1937       REAL(dp)   , DIMENSION(:,:)  , INTENT(INOUT) :: dd_value  
    1938       INTEGER(i4), DIMENSION(:,:)  , INTENT(INOUT) :: id_detect 
    1939  
    1940       ! local variable 
    1941       INTEGER(i4), DIMENSION(2) :: il_shape 
    1942  
    1943       INTEGER(i4) :: il_i1 
    1944       INTEGER(i4) :: il_i2 
    1945       INTEGER(i4) :: il_j1 
    1946       INTEGER(i4) :: il_j2 
    1947  
    1948       INTEGER(i4) :: il_half1 
    1949       INTEGER(i4) :: il_half2 
    1950  
    1951       ! loop indices 
    1952       INTEGER(i4) :: ji 
    1953       INTEGER(i4) :: jj 
    1954       !---------------------------------------------------------------- 
    1955  
    1956       il_shape(:)=SHAPE(dd_value(:,:)) 
    1957  
    1958       il_i1=1 
    1959       il_i2=il_shape(1) 
    1960  
    1961       il_j1=1 
    1962       il_j2=il_shape(2) 
    1963  
    1964       il_half1=CEILING(il_shape(1)*0.5) 
    1965       il_half2=CEILING(il_shape(2)*0.5) 
    1966  
    1967       DO jj=1,il_half2 
    1968  
    1969          DO ji=1,il_half1 
    1970              
    1971             ! lower left point 
    1972             IF(id_detect(ji,jj)==1)THEN 
    1973  
    1974                dd_value( ji,jj)=dd_value(il_i1,il_j1) 
    1975                id_detect(ji,jj)=0 
    1976  
    1977             ENDIF 
    1978  
    1979             ! lower right point 
    1980             IF(id_detect(il_shape(1)-ji+1,jj)==1)THEN 
    1981  
    1982                dd_value( il_shape(1)-ji+1,jj)=dd_value(il_i2,il_j1) 
    1983                id_detect(il_shape(1)-ji+1,jj)=0 
    1984  
    1985             ENDIF 
    1986  
    1987             ! upper left point 
    1988             IF(id_detect(ji,il_shape(2)-jj+1)==1)THEN 
    1989  
    1990                dd_value( ji,il_shape(2)-jj+1)=dd_value(il_i1,il_j2) 
    1991                id_detect(ji,il_shape(2)-jj+1)=0 
    1992  
    1993             ENDIF             
    1994  
    1995             ! upper right point 
    1996             IF(id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)==1)THEN 
    1997  
    1998                dd_value( il_shape(1)-ji+1,il_shape(2)-jj+1)=dd_value(il_i2,il_j2) 
    1999                id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)=0 
    2000  
    2001             ENDIF             
    2002  
    2003          ENDDO 
    2004  
    2005       ENDDO 
    2006  
    2007    END SUBROUTINE interp__2D_nearest_fill 
    2008    !> @endcode 
    2009    !------------------------------------------------------------------- 
    2010    !> @brief 
    2011    !> This subroutine  
    2012    !>  
    2013    !> @details  
    2014    !> 
    2015    !> @author J.Paul 
    2016    !> - Nov, 2013- Initial Version 
    2017    ! 
    2018    !> @param[in] 
    2019    !> @param[out]  
    2020    !------------------------------------------------------------------- 
    2021    !> @code 
    2022    FUNCTION interp__1D_cubic_coef( dd_value, & 
    2023       &                            dd_dfdx,  & 
    2024       &                            dd_fill ) 
    2025       IMPLICIT NONE 
    2026       ! Argument 
    2027       REAL(dp), DIMENSION(:)  , INTENT(IN) :: dd_value  
    2028       REAL(dp), DIMENSION(:)  , INTENT(IN) :: dd_dfdx   
    2029       REAL(dp)                , INTENT(IN) :: dd_fill   
    2030  
    2031       ! function 
    2032       REAL(dp), DIMENSION(4) :: interp__1D_cubic_coef 
    2033  
    2034       ! local variable 
    2035  
    2036       REAL(dp), DIMENSION(4,4), PARAMETER :: dp_matrix = RESHAPE( & 
    2037       & (/  1 , 0 , 0 , 0 ,&  
    2038            -1 , 1 , 0 , 0 ,&  
    2039            -3 , 3 ,-2 ,-1 ,&  
    2040             2 ,-2 , 1 , 1  /), & 
    2041       & (/ 4, 4 /) ) 
    2042  
    2043       REAL(dp), DIMENSION(4) :: dl_vect 
    2044  
    2045       !---------------------------------------------------------------- 
    2046       IF( SIZE(dd_value(:))/= 2 .OR. & 
    2047       &   SIZE(dd_dfdx(:) )/= 2  )THEN 
    2048     
    2049          CALL logger_error("INTERP CUBIC COEF: invalid dimension of "//& 
    2050          &              "input tables. shape should be (/2,2/)") 
    2051  
    2052       ELSEIF( ANY(dd_value(:)==dd_fill) .OR. & 
    2053       &       ANY( dd_dfdx(:)==dd_fill) )THEN 
    2054          CALL logger_error("INTERP CUBIC COEF: fill value detected. "//& 
    2055          &              "can not compute coefficient.") 
    2056       ELSE 
    2057  
    2058          dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. ) 
    2059          dl_vect( 3: 4)=PACK(dd_dfdx(:),.TRUE. ) 
    2060  
    2061          interp__1D_cubic_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:)) 
    2062  
    2063       ENDIF 
    2064  
    2065    END FUNCTION interp__1D_cubic_coef 
    2066    !> @endcode 
    2067    !------------------------------------------------------------------- 
    2068    !> @brief 
    2069    !> This subroutine  
    2070    !>  
    2071    !> @details  
    2072    !> 
    2073    !> @author J.Paul 
    2074    !> - Nov, 2013- Initial Version 
    2075    ! 
    2076    !> @param[in] 
    2077    !> @param[out]  
    2078    !------------------------------------------------------------------- 
    2079    !> @code 
    2080    SUBROUTINE interp__1D_cubic_fill( dd_value, id_detect, dd_coef, dd_fill, & 
    2081    &                                 ld_even, id_rhoi ) 
    2082       IMPLICIT NONE 
    2083       ! Argument 
    2084       REAL(dp)        , DIMENSION(:), INTENT(INOUT) :: dd_value  
    2085       INTEGER(i4)     , DIMENSION(:), INTENT(INOUT) :: id_detect 
    2086       REAL(dp)        , DIMENSION(:), INTENT(IN   ) :: dd_coef  
    2087       REAL(dp)                      , INTENT(IN   ) :: dd_fill  
    2088       LOGICAL         ,               INTENT(IN   ) :: ld_even 
    2089       INTEGER(I4)     ,               INTENT(IN   ) :: id_rhoi 
    2090  
    2091       ! local variable 
    2092       INTEGER(i4), DIMENSION(1)  :: il_shape 
    2093  
    2094       REAL(dp)   , DIMENSION(4) :: dl_vect 
    2095       REAL(dp)                   :: dl_dx 
    2096       REAL(dp)                   :: dl_x 
    2097       REAL(dp)                   :: dl_x2 
    2098       REAL(dp)                   :: dl_x3 
    2099  
    2100       ! loop indices 
    2101       INTEGER(i4) :: ji 
    2102       !---------------------------------------------------------------- 
    2103  
    2104       IF( SIZE(dd_coef(:)) /= 4 )THEN 
    2105          CALL logger_error("INTERP CUBIC FILL: invalid dimension of "//& 
    2106          &              "coef table. shape should be (/4/)") 
    2107       !ELSEIF( ANY(dd_value(:)==dd_fill .AND. id_detect(:)==0 ) .OR. & 
    2108       ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN 
    2109          CALL logger_error("INTERP CUBIC FILL: fill value detected. "//& 
    2110          &              "can not compute interpolation") 
    2111       ELSE 
    2112  
    2113          il_shape(:)=SHAPE(dd_value(:)) 
    2114  
    2115          dl_dx=1./REAL(id_rhoi) 
    2116           
    2117          DO ji=1,il_shape(1) 
    2118              
    2119             IF(id_detect(ji)==1)THEN 
    2120  
    2121                IF( ld_even )THEN 
    2122                   dl_x=(ji-1)*dl_dx - dl_dx*0.5  
    2123                ELSE ! odd refinement 
    2124                   dl_x=(ji-1)*dl_dx  
    2125                ENDIF 
    2126                dl_x2=dl_x*dl_x 
    2127                dl_x3=dl_x2*dl_x 
    2128  
    2129                dl_vect(:)=(/1._dp, dl_x, dl_x2, dl_x3 /) 
    2130  
    2131                dd_value(ji)=DOT_PRODUCT(dd_coef(:),dl_vect(:)) 
    2132                id_detect(ji)=0 
    2133  
    2134             ENDIF 
    2135  
    2136          ENDDO 
    2137  
    2138       ENDIF 
    2139  
    2140    END SUBROUTINE interp__1D_cubic_fill 
    2141    !> @endcode 
    2142    !------------------------------------------------------------------- 
    2143    !> @brief 
    2144    !> This subroutine  
    2145    !>  
    2146    !> @details  
    2147    !> 
    2148    !> @author J.Paul 
    2149    !> - Nov, 2013- Initial Version 
    2150    ! 
    2151    !> @param[in] 
    2152    !> @param[out]  
    2153    !------------------------------------------------------------------- 
    2154    !> @code 
    2155    FUNCTION interp__1D_linear_coef( dd_value, dd_fill ) 
    2156       IMPLICIT NONE 
    2157       ! Argument 
    2158       REAL(dp), DIMENSION(:)  , INTENT(IN) :: dd_value  
    2159       REAL(dp)                , INTENT(IN) :: dd_fill   
    2160  
    2161       ! function 
    2162       REAL(dp), DIMENSION(2) :: interp__1D_linear_coef 
    2163  
    2164       ! local variable 
    2165  
    2166       REAL(dp), DIMENSION(2,2), PARAMETER :: dp_matrix = RESHAPE( & 
    2167       & (/  1 , 0 ,&  
    2168            -1 , 1  /), & 
    2169       &  (/ 2, 2 /) ) 
    2170  
    2171       REAL(dp), DIMENSION(2) :: dl_vect 
    2172  
    2173       !---------------------------------------------------------------- 
    2174  
    2175       IF( ANY(SHAPE(dd_value(:))/= 2) )THEN 
    2176          CALL logger_error("INTERP LINEAR COEF: invalid dimension of "//& 
    2177          &              "input tables. shape should be (/2/)") 
    2178       ELSEIF( ANY(dd_value(:)==dd_fill) )THEN 
    2179          CALL logger_error("INTERP LINEAR COEF: fill value detected. "//& 
    2180          &              "can not compute coefficient.") 
    2181       ELSE 
    2182           
    2183          dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. ) 
    2184  
    2185          interp__1D_linear_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:)) 
    2186  
    2187       ENDIF 
    2188  
    2189    END FUNCTION interp__1D_linear_coef 
    2190    !> @endcode 
    2191    !------------------------------------------------------------------- 
    2192    !> @brief 
    2193    !> This subroutine  
    2194    !>  
    2195    !> @details  
    2196    !> 
    2197    !> @author J.Paul 
    2198    !> - Nov, 2013- Initial Version 
    2199    ! 
    2200    !> @param[in] 
    2201    !> @param[out]  
    2202    !------------------------------------------------------------------- 
    2203    !> @code 
    2204    SUBROUTINE interp__1D_linear_fill( dd_value, id_detect, dd_coef, dd_fill, & 
    2205    &                                  ld_even, id_rhoi ) 
    2206       IMPLICIT NONE 
    2207       ! Argument 
    2208       REAL(dp)        , DIMENSION(:), INTENT(INOUT) :: dd_value  
    2209       INTEGER(i4)     , DIMENSION(:), INTENT(INOUT) :: id_detect 
    2210       REAL(dp)        , DIMENSION(:), INTENT(IN   ) :: dd_coef  
    2211       REAL(dp)                      , INTENT(IN   ) :: dd_fill  
    2212       LOGICAL         ,               INTENT(IN   ) :: ld_even 
    2213       INTEGER(I4)     ,               INTENT(IN   ) :: id_rhoi       
    2214  
    2215       ! local variable 
    2216       INTEGER(i4), DIMENSION(1) :: il_shape 
    2217  
    2218       REAL(dp)   , DIMENSION(2) :: dl_vect 
    2219       REAL(dp)                  :: dl_dx 
    2220       REAL(dp)                  :: dl_x 
    2221  
    2222       ! loop indices 
    2223       INTEGER(i4) :: ji 
    2224       !---------------------------------------------------------------- 
    2225  
    2226       IF( SIZE(dd_coef(:)) /= 2 )THEN 
    2227          CALL logger_error("INTERP LINEAR FILL: invalid dimension of "//& 
    2228          &              "coef table. shape should be (/2/)") 
    2229       !ELSEIF( ANY(dd_value(:)==dd_fill .AND. id_detect(:)==0 ) .OR. & 
    2230       ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN 
    2231          CALL logger_error("INTERP LINEAR FILL: fill value detected. "//& 
    2232          &              "can not compute interpolation") 
    2233       ELSE 
    2234  
    2235          il_shape(:)=SHAPE(dd_value) 
    2236  
    2237          dl_dx=1./REAL(id_rhoi) 
    2238  
    2239          DO ji=1,il_shape(1) 
    2240              
    2241             IF(id_detect(ji)==1)THEN 
    2242  
    2243                IF( ld_even )THEN 
    2244                   dl_x=(ji-1)*dl_dx - dl_dx*0.5  
    2245                ELSE ! odd refinement 
    2246                   dl_x=(ji-1)*dl_dx  
    2247                ENDIF 
    2248  
    2249                dl_vect(:)=(/1._dp, dl_x /) 
    2250  
    2251                dd_value(ji)=DOT_PRODUCT(dd_coef(:),dl_vect(:)) 
    2252                id_detect(ji)=0 
    2253  
    2254             ENDIF 
    2255  
    2256          ENDDO 
    2257  
    2258       ENDIF 
    2259  
    2260    END SUBROUTINE interp__1D_linear_fill 
    2261    !> @endcode 
    2262    !------------------------------------------------------------------- 
    2263    !> @brief 
    2264    !> This subroutine  
    2265    !>  
    2266    !> @details  
    2267    !> 
    2268    !> @author J.Paul 
    2269    !> - Nov, 2013- Initial Version 
    2270    ! 
    2271    !> @param[in] 
    2272    !> @param[out]  
    2273    !------------------------------------------------------------------- 
    2274    !> @code 
    2275    SUBROUTINE interp__1D_nearest_fill( dd_value, id_detect ) 
    2276       IMPLICIT NONE 
    2277       ! Argument 
    2278       REAL(dp)   , DIMENSION(:), INTENT(INOUT) :: dd_value  
    2279       INTEGER(i4), DIMENSION(:), INTENT(INOUT) :: id_detect 
    2280  
    2281       ! local variable 
    2282       INTEGER(i4), DIMENSION(1) :: il_shape 
    2283  
    2284       INTEGER(i4) :: il_i1 
    2285       INTEGER(i4) :: il_i2 
    2286  
    2287       INTEGER(i4) :: il_half1 
    2288        
    2289       ! loop indices 
    2290       INTEGER(i4) :: ji 
    2291       !---------------------------------------------------------------- 
    2292  
    2293       il_shape(:)=SHAPE(dd_value) 
    2294  
    2295       il_i1=1 
    2296       il_i2=il_shape(1) 
    2297  
    2298       il_half1=CEILING(il_shape(1)*0.5) 
    2299        
    2300       DO ji=1,il_half1 
    2301           
    2302          ! lower left point 
    2303          IF(id_detect(ji)==1)THEN 
    2304  
    2305             dd_value( ji)=dd_value(il_i1) 
    2306             id_detect(ji)=0 
    2307  
    2308          ENDIF 
    2309  
    2310          ! lower right point 
    2311          IF(id_detect(il_shape(1)-ji+1)==1)THEN 
    2312  
    2313             dd_value( il_shape(1)-ji+1)=dd_value(il_i2) 
    2314             id_detect(il_shape(1)-ji+1)=0 
    2315  
    2316          ENDIF 
    2317  
    2318       ENDDO 
    2319  
    2320    END SUBROUTINE interp__1D_nearest_fill 
    2321    !> @endcode 
    2322958END MODULE interp 
Note: See TracChangeset for help on using the changeset viewer.