!---------------------------------------------------------------------- ! NEMO system team, System and Interface for oceanic RElocable Nesting !---------------------------------------------------------------------- ! ! MODULE: interp ! ! DESCRIPTION: !> @brief !> This module manage interpolation on regular grid. !> @note It is used to work on ORCA grid, as we work only with grid indices. !> !> @warning due to the use of second derivative when using cubic interpolation !> you should add 2 extrabands !> !> @details !> @author !> J.Paul ! REVISION HISTORY: !> @date Nov, 2013 - Initial Version ! !> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !> @todo !> - interp 3D !> - see issue when fill value is zero for check cubic_fill.. !---------------------------------------------------------------------- MODULE interp USE netcdf ! nf90 library USE global ! global variable USE kind ! F90 kind parameter USE logger ! log file manager USE fct ! basic useful function ! USE date ! date manager USE att ! attribute manager USE dim ! dimension manager USE var ! variable manager ! USE file ! file manager ! USE iom ! I/O manager ! USE dom ! domain manager USE grid ! grid manager USE extrap ! extrapolation manager ! USE interp ! interpolation manager ! USE filter ! filter manager ! USE mpp ! MPP manager ! USE iom_mpp ! MPP I/O manager IMPLICIT NONE PRIVATE ! NOTE_avoid_public_variables_if_possible ! type and variable ! function and subroutine PUBLIC :: interp_detect !< detected point to be interpolated PUBLIC :: interp_fill_value !< interpolate value over detectected point PUBLIC :: interp_create_mixed_grid !< create mixed grid PUBLIC :: interp_clean_mixed_grid !< clean mixed grid PRIVATE :: interp__detect !< detected point to be interpolated PRIVATE :: interp__detect_wrapper !< detected point to be interpolated PRIVATE :: interp__fill_value_wrapper !< interpolate value over detectected point PRIVATE :: interp__fill_value !< interpolate value over detectected point PRIVATE :: interp__clean_even_grid !< clean even mixed grid PRIVATE :: interp__del_offset !< remove offset from interpolated grid PRIVATE :: interp__check_method !< check if interpolation method available ! PRIVATE :: interp__3D !< interpolate 3D grid PRIVATE :: interp__2D !< interpolate 2D grid PRIVATE :: interp__1D !< interpolate 1D grid PRIVATE :: interp__2D_cubic_coef !< compute coefficient for bicubic interpolation PRIVATE :: interp__2D_cubic_fill !< compute bicubic interpolation PRIVATE :: interp__2D_linear_coef !< compute coefficient for bilinear interpolation PRIVATE :: interp__2D_linear_fill !< compute bilinear interpolation PRIVATE :: interp__2D_nearest_fill !< compute nearest interpolation PRIVATE :: interp__1D_cubic_coef !< compute coefficient for cubic interpolation PRIVATE :: interp__1D_cubic_fill !< compute cubic interpolation PRIVATE :: interp__1D_linear_coef !< compute coefficient for linear interpolation PRIVATE :: interp__1D_linear_fill !< compute linear interpolation PRIVATE :: interp__1D_nearest_fill !< compute nearest interpolation ! PRIVATE :: interp__longitude TYPE TINTERP !CHARACTER(LEN=lc) :: c_name = 'unknown' !< interpolation method name !CHARACTER(LEN=lc) :: c_factor = 'unknown' !< interpolation factor !CHARACTER(LEN=lc) :: c_divisor = 'unknown' !< interpolation divisor CHARACTER(LEN=lc) :: c_name = '' !< interpolation method name CHARACTER(LEN=lc) :: c_factor = '' !< interpolation factor CHARACTER(LEN=lc) :: c_divisor= '' !< interpolation divisor END TYPE TINTERP INTERFACE interp_detect MODULE PROCEDURE interp__detect_wrapper !< detected point to be interpolated END INTERFACE interp_detect INTERFACE interp_fill_value MODULE PROCEDURE interp__fill_value_wrapper !< detected point to be interpolated END INTERFACE interp_fill_value CONTAINS !------------------------------------------------------------------- !> @brief !> This function check if interpolation method available. !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] cd_method : interpolation method !> @return !> @todo see extrap_detect !------------------------------------------------------------------- !> @code FUNCTION interp__check_method( cd_method ) IMPLICIT NONE ! Argument CHARACTER(LEN=lc) :: cd_method ! function LOGICAL :: interp__check_method ! local variable CHARACTER(LEN=lc) :: cl_interp CHARACTER(LEN=lc) :: cl_method ! loop indices INTEGER(I4) :: ji !---------------------------------------------------------------- cl_method=fct_lower(cd_method) interp__check_method=.FALSE. DO ji=1,ig_ninterp cl_interp=fct_lower(cg_interp_list(ji)) IF( TRIM(cl_interp) == TRIM(cl_method) )THEN interp__check_method=.TRUE. EXIT ENDIF ENDDO END FUNCTION interp__check_method !> @endcode !------------------------------------------------------------------- !> @brief !> This function detected point to be interpolated. !> !> @details !> Actually it checks, the number of dimension used for this variable !> and launch interp__detect which detected point to be interpolated. !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] td_mix : mixed grid variable (to interpolate) !> @param[in] id_rho : table of refinement factor !> @return table of detected point to be interpolated !------------------------------------------------------------------- !> @code FUNCTION interp__detect_wrapper( td_mix, id_rho ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(IN) :: td_mix INTEGER(I4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho ! function INTEGER(i4), DIMENSION(td_mix%t_dim(1)%i_len,& & td_mix%t_dim(2)%i_len,& & td_mix%t_dim(3)%i_len ) :: interp__detect_wrapper ! local variable ! loop indices !---------------------------------------------------------------- IF( .NOT. ANY(td_mix%t_dim(1:3)%l_use) )THEN ! no dimension I-J-K used CALL logger_debug(" INTERP DETECT: nothing done for variable"//& & TRIM(td_mix%c_name) ) interp__detect_wrapper(:,:,:)=0 ELSE IF( ALL(td_mix%t_dim(1:3)%l_use) )THEN ! detect point to be interpolated on I-J-K CALL logger_debug(" INTERP DETECT: detect point "//& & TRIM(td_mix%c_point)//" for variable "//& & TRIM(td_mix%c_name) ) interp__detect_wrapper(:,:,:)=interp__detect( td_mix, id_rho(:) ) ELSE IF( ALL(td_mix%t_dim(1:2)%l_use) )THEN ! detect point to be interpolated on I-J CALL logger_debug(" INTERP DETECT: detect point "//& & TRIM(td_mix%c_point)//" for variable "//& & TRIM(td_mix%c_name) ) interp__detect_wrapper(:,:,1:1)=interp__detect( td_mix, id_rho(:)) ELSE IF( td_mix%t_dim(3)%l_use )THEN ! detect point to be interpolated on K CALL logger_debug(" INTERP DETECT: detect vertical point "//& & " for variable "//TRIM(td_mix%c_name) ) interp__detect_wrapper(1:1,1:1,:)=interp__detect( td_mix, id_rho(:) ) ENDIF END FUNCTION interp__detect_wrapper !> @endcode !------------------------------------------------------------------- !> @brief !> This function detected point to be interpolated. !> !> @details !> A special case is done for even refinement on ARAKAWA-C grid. !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] td_mix : mixed grid variable (to interpolate) !> @param[in] id_rho : table of refinement factor !> @return table of detected point to be interpolated !------------------------------------------------------------------- !> @code FUNCTION interp__detect( td_mix, id_rho ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(IN) :: td_mix INTEGER(I4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho ! function INTEGER(i4), DIMENSION(td_mix%t_dim(1)%i_len,& & td_mix%t_dim(2)%i_len,& & td_mix%t_dim(3)%i_len ) :: interp__detect ! local variable INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho INTEGER(I4) :: il_xextra INTEGER(I4) :: il_yextra INTEGER(I4) :: il_zextra INTEGER(i4), DIMENSION(3) :: il_dim LOGICAL , DIMENSION(3) :: ll_even ! loop indices INTEGER(I4) :: ji INTEGER(I4) :: jj INTEGER(I4) :: jk !---------------------------------------------------------------- ALLOCATE( il_rho(ip_maxdim) ) il_rho(:)=1 IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:) ! special case for even refinement on ARAKAWA-C grid ll_even(:)=.FALSE. IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE. IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE. IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE. SELECT CASE(TRIM(td_mix%c_point)) CASE('U') ll_even(1)=.FALSE. CASE('V') ll_even(2)=.FALSE. CASE('F') ll_even(:)=.FALSE. END SELECT IF( ll_even(1) ) il_rho(jp_I)=il_rho(jp_I)+1 IF( ll_even(2) ) il_rho(jp_J)=il_rho(jp_J)+1 IF( ll_even(3) ) il_rho(jp_K)=il_rho(jp_K)+1 ! special case for cubic interpolation il_xextra=0 il_yextra=0 il_zextra=0 SELECT CASE(TRIM(td_mix%c_interp(1))) CASE('cubic') ! those points can not be compute cause cubic interpolation ! need second derivative. IF( il_rho(jp_I) /= 1 ) il_xextra=3*il_rho(jp_I) IF( il_rho(jp_J) /= 1 ) il_yextra=3*il_rho(jp_J) IF( il_rho(jp_K) /= 1 ) il_zextra=3*il_rho(jp_K) END SELECT il_dim(:)=td_mix%t_dim(1:3)%i_len interp__detect(:,:,:)=1 ! do not compute coarse grid point interp__detect(1:td_mix%t_dim(1)%i_len:il_rho(jp_I), & & 1:td_mix%t_dim(2)%i_len:il_rho(jp_J), & & 1:td_mix%t_dim(3)%i_len:il_rho(jp_K) ) = 0 ! do not compute point near fill value FORALL( ji=1:il_dim(1):il_rho(jp_I), & & jj=1:il_dim(2):il_rho(jp_J), & & jk=1:il_dim(3):il_rho(jp_K), & & td_mix%d_value(ji,jj,jk,1) == td_mix%d_fill ) ! i-direction interp__detect(MAX(1,ji-il_xextra):MIN(ji+il_xextra,il_dim(1)),& & MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),& & MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0 ! j-direction interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),& & MAX(1,jj-il_yextra):MIN(jj+il_yextra,il_dim(2)),& & MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0 ! k-direction interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),& & MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),& & MAX(1,jk-il_zextra):MIN(jk+il_zextra,il_dim(3)) )=0 END FORALL END FUNCTION interp__detect !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine create mixed grid. !> !> @details !> Created grid is fine resolution grid. !> First and last point are coasre grid point. !> !> A special case is done for even refinement on ARAKAWA-C grid. !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] td_var : coarse grid variable (should be extrapolated) !> @param[out] td_mix : mixed grid variable !> @param[in] id_rho : table of refinment factor !> !> @todo !------------------------------------------------------------------- !> @code SUBROUTINE interp_create_mixed_grid( td_var, td_mix, id_rho ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(IN ) :: td_var TYPE(TVAR) , INTENT( OUT) :: td_mix INTEGER(I4), DIMENSION(:), INTENT(IN ), OPTIONAL :: id_rho ! local variable INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho INTEGER(i4) :: il_xextra INTEGER(i4) :: il_yextra INTEGER(i4) :: il_zextra LOGICAL, DIMENSION(3) :: ll_even ! loop indices !---------------------------------------------------------------- ALLOCATE(il_rho(ip_maxdim)) il_rho(:)=1 IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:) ! special case for even refinement on ARAKAWA-C grid ll_even(:)=.FALSE. IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE. IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE. IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE. SELECT CASE(TRIM(td_var%c_point)) CASE('U') ll_even(1)=.FALSE. CASE('V') ll_even(2)=.FALSE. CASE('F') ll_even(:)=.FALSE. END SELECT IF( ll_even(1) ) il_rho(jp_I)=il_rho(jp_I)+1 IF( ll_even(2) ) il_rho(jp_J)=il_rho(jp_J)+1 IF( ll_even(3) ) il_rho(jp_K)=il_rho(jp_K)+1 ! copy variable td_mix=td_var ! compute new dimension length il_xextra=il_rho(jp_I)-1 td_mix%t_dim(1)%i_len=td_mix%t_dim(1)%i_len*il_rho(jp_I)-il_xextra il_yextra=il_rho(jp_J)-1 td_mix%t_dim(2)%i_len=td_mix%t_dim(2)%i_len*il_rho(jp_J)-il_yextra il_zextra=il_rho(jp_K)-1 td_mix%t_dim(3)%i_len=td_mix%t_dim(3)%i_len*il_rho(jp_K)-il_zextra IF( ASSOCIATED(td_mix%d_value) ) DEALLOCATE( td_mix%d_value ) ALLOCATE( td_mix%d_value( td_mix%t_dim(1)%i_len, & & td_mix%t_dim(2)%i_len, & & td_mix%t_dim(3)%i_len, & & td_mix%t_dim(4)%i_len) ) ! initialise to FillValue td_mix%d_value(:,:,:,:)=td_mix%d_fill ! quid qd coord ou bathy fourni par user ?? (offset ??) td_mix%d_value(1:td_mix%t_dim(1)%i_len:il_rho(jp_I), & & 1:td_mix%t_dim(2)%i_len:il_rho(jp_J), & & 1:td_mix%t_dim(3)%i_len:il_rho(jp_K), :) = & & td_var%d_value(:,:,:,:) END SUBROUTINE interp_create_mixed_grid !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine remove points added to mixed grid to compute !> interpolation in the special case of even refinement on ARAKAWA-C grid. !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] td_mix : mixed grid variable !> @param[in] id_rho : table of refinment factor !> !> @todo !------------------------------------------------------------------- !> @code SUBROUTINE interp__clean_even_grid( td_mix, id_rho ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(INOUT) :: td_mix INTEGER(I4), DIMENSION(:), INTENT(IN ), OPTIONAL :: id_rho ! local variable INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho INTEGER(i4) :: il_xextra INTEGER(i4) :: il_yextra INTEGER(i4) :: il_zextra LOGICAL, DIMENSION(3) :: ll_even LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ll_mask REAL(dp), DIMENSION(:) , ALLOCATABLE :: dl_vect TYPE(TVAR) :: tl_mix ! loop indices !---------------------------------------------------------------- ALLOCATE(il_rho(ip_maxdim)) il_rho(:)=1 IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:) ! special case for even refinement on ARAKAWA-C grid ll_even(:)=.FALSE. IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE. IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE. IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE. SELECT CASE(TRIM(td_mix%c_point)) CASE('U') ll_even(1)=.FALSE. CASE('V') ll_even(2)=.FALSE. CASE('F') ll_even(:)=.FALSE. END SELECT ! copy variable tl_mix=td_mix ! remove some point only if refinement in some direction is even IF( ANY(ll_even(:)) )THEN ALLOCATE( ll_mask( tl_mix%t_dim(1)%i_len, & & tl_mix%t_dim(2)%i_len, & & tl_mix%t_dim(3)%i_len, & & tl_mix%t_dim(4)%i_len) ) ll_mask(:,:,:,:)=.TRUE. IF( tl_mix%t_dim(1)%l_use .AND. ll_even(1) )THEN il_rho(jp_I)=il_rho(jp_I)+1 ! locate wrong point on mixed grid ll_mask(1:td_mix%t_dim(1)%i_len:il_rho(jp_I),:,:,:)=.FALSE. ! compute coasre grid dimension length il_xextra=il_rho(jp_I)-1 td_mix%t_dim(1)%i_len=(tl_mix%t_dim(1)%i_len+il_xextra)/il_rho(jp_I) il_rho(jp_I)=il_rho(jp_I)-1 ! compute right fine grid dimension length td_mix%t_dim(1)%i_len=td_mix%t_dim(1)%i_len*il_rho(jp_I)-il_xextra ENDIF IF( tl_mix%t_dim(2)%l_use .AND. ll_even(2) )THEN il_rho(jp_J)=il_rho(jp_J)+1 ! locate wrong point on mixed grid ll_mask(:,1:tl_mix%t_dim(2)%i_len:il_rho(jp_J),:,:)=.FALSE. ! compute coasre grid dimension length il_yextra=il_rho(jp_J)-1 td_mix%t_dim(2)%i_len=(tl_mix%t_dim(2)%i_len+il_yextra)/il_rho(jp_J) il_rho(jp_J)=il_rho(jp_J)-1 ! compute right fine grid dimension length td_mix%t_dim(2)%i_len=td_mix%t_dim(2)%i_len*il_rho(jp_J)-il_yextra ENDIF IF( tl_mix%t_dim(3)%l_use .AND. ll_even(3) )THEN il_rho(jp_K)=il_rho(jp_K)+1 ! locate wrong point on mixed grid ll_mask(:,:,1:tl_mix%t_dim(3)%i_len:il_rho(jp_K),:)=.FALSE. ! compute coasre grid dimension length il_zextra=il_rho(jp_K)-1 td_mix%t_dim(3)%i_len=(tl_mix%t_dim(3)%i_len+il_zextra)/il_rho(jp_K) il_rho(jp_K)=il_rho(jp_K)-1 ! compute right fine grid dimension length td_mix%t_dim(3)%i_len=td_mix%t_dim(3)%i_len*il_rho(jp_K)-il_zextra ENDIF IF( ASSOCIATED(td_mix%d_value) ) DEALLOCATE( td_mix%d_value ) ALLOCATE( td_mix%d_value( td_mix%t_dim(1)%i_len, & & td_mix%t_dim(2)%i_len, & & td_mix%t_dim(3)%i_len, & & td_mix%t_dim(4)%i_len) ) ! initialise to FillValue td_mix%d_value(:,:,:,:)=td_mix%d_fill IF( COUNT(ll_mask(:,:,:,:)) /= SIZE(td_mix%d_value(:,:,:,:)) )THEN CALL logger_error("INTERP CLEAN EVEN GRID: output value size "//& & " and mask count differ ") ELSE ALLOCATE( dl_vect(COUNT(ll_mask(:,:,:,:))) ) dl_vect(:)= PACK( tl_mix%d_value(:,:,:,:), & & MASK=ll_mask(:,:,:,:) ) td_mix%d_value(:,:,:,:)=RESHAPE( dl_vect(:), & & SHAPE=td_mix%t_dim(:)%i_len ) DEALLOCATE( dl_vect ) ENDIF DEALLOCATE( ll_mask ) ENDIF CALL var_clean(tl_mix) END SUBROUTINE interp__clean_even_grid !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine save interpolated value over domain. !> And so remove points added on mixed grid !> to compute interpolation !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] td_var : table of mixed grid variable (to interpolate) !> @param[in] id_rho : table of refinement factor !> @return !------------------------------------------------------------------- !> @code SUBROUTINE interp_clean_mixed_grid( td_mix, td_var, & & id_rho, & & id_offset ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(IN ) :: td_mix TYPE(TVAR) , INTENT( OUT) :: td_var INTEGER(I4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho INTEGER(I4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset ! local variable INTEGER(i4) :: il_imin0 INTEGER(i4) :: il_imax0 INTEGER(i4) :: il_jmin0 INTEGER(i4) :: il_jmax0 INTEGER(i4) :: il_imin1 INTEGER(i4) :: il_jmin1 INTEGER(i4) :: il_imax1 INTEGER(i4) :: il_jmax1 INTEGER(i4), DIMENSION(2,2) :: il_offset REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value TYPE(TVAR) :: tl_mix ! loop indices !---------------------------------------------------------------- il_offset(:,:)=0 IF( PRESENT(id_offset) )THEN IF( ANY( SHAPE(id_offset(:,:)) /= SHAPE(il_offset(:,:)) ) )THEN CALL logger_error("INTERP CLEAN MIXED GRID: invalid dimension of"//& & " offset table") ELSE il_offset(:,:)=id_offset(:,:) ENDIF ENDIF ! copy mixed variable in temporary structure tl_mix=td_mix ! remove unusefull points over mixed grid for even refinement CALL interp__clean_even_grid(tl_mix, id_rho(:)) ! copy cleaned mixed variable td_var=tl_mix IF( ALL(td_var%t_dim(1:2)%l_use) )THEN ! delete table of value CALL var_del_value(td_var) ! compute domain indices il_imin0=1 ; il_imax0=td_var%t_dim(1)%i_len il_jmin0=1 ; il_jmax0=td_var%t_dim(2)%i_len il_imin1=il_imin0+il_offset(1,1) il_jmin1=il_jmin0+il_offset(2,1) il_imax1=il_imax0-il_offset(1,2) il_jmax1=il_jmax0-il_offset(2,2) SELECT CASE(TRIM(td_var%c_point)) CASE('U') il_imin1=il_imin0+(il_offset(1,1)-1) il_imax1=il_imax0-(il_offset(1,2)+1) CASE('V') il_jmin1=il_jmin0+(il_offset(2,1)-1) il_jmax1=il_jmax0-(il_offset(2,2)+1) CASE('F') il_imin1=il_imin0+(il_offset(1,1)-1) il_imax1=il_imax0-(il_offset(1,2)+1) il_jmin1=il_jmin0+(il_offset(2,1)-1) il_jmax1=il_jmax0-(il_offset(2,2)+1) END SELECT ! compute new dimension td_var%t_dim(1)%i_len=il_imax1-il_imin1+1 td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1 ALLOCATE(dl_value(td_var%t_dim(1)%i_len, & & td_var%t_dim(2)%i_len, & & td_var%t_dim(3)%i_len, & & td_var%t_dim(4)%i_len) ) dl_value( 1:td_var%t_dim(1)%i_len, & & :,:,:) = tl_mix%d_value( il_imin1:il_imax1, & & il_jmin1:il_jmax1, & & :, : ) ! add variable value CALL var_add_value(td_var,dl_value(:,:,:,:)) ! save variable type td_var%i_type=tl_mix%i_type DEALLOCATE(dl_value) ENDIF CALL var_clean(tl_mix) END SUBROUTINE interp_clean_mixed_grid !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine save interpolated value over domain. !> And so remove points added on mixed grid !> to compute interpolation !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] td_var : table of mixed grid variable (to interpolate) !> @param[in] id_offset : table of offset !> @return !------------------------------------------------------------------- !> @code SUBROUTINE interp__del_offset( td_var, id_offset ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(INOUT) :: td_var INTEGER(i4), DIMENSION(:,:), INTENT(IN ) :: id_offset ! local variable INTEGER(i4) :: il_imin1 INTEGER(i4) :: il_jmin1 INTEGER(i4) :: il_imax1 INTEGER(i4) :: il_jmax1 REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value ! loop indices !---------------------------------------------------------------- IF( ALL(td_var%t_dim(1:2)%l_use) )THEN il_imin1=1+id_offset(1,1) il_jmin1=1+id_offset(2,1) il_imax1=td_var%t_dim(1)%i_len-id_offset(2,1) il_jmax1=td_var%t_dim(2)%i_len-id_offset(2,2) ! compute new dimension td_var%t_dim(1)%i_len=il_imax1-il_imin1+1 td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1 ALLOCATE( dl_value( td_var%t_dim(1)%i_len, & & td_var%t_dim(2)%i_len, & & td_var%t_dim(3)%i_len, & & td_var%t_dim(4)%i_len) ) dl_value(:,:,:,:)=td_var%d_value( il_imin1:il_imax1, & & il_jmin1:il_jmax1, & & :,: ) ! delete table of value CALL var_del_value(td_var) ! add variable value CALL var_add_value(td_var,dl_value(:,:,:,:)) DEALLOCATE(dl_value) ENDIF END SUBROUTINE interp__del_offset !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine interpolate detected point. !> !> @details !> Actually it checks, the number of dimension used for this variable !> and launch interp__fill_value. !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] td_var : variable to be interpolated !> @param[in] id_rho : table of refinement factor !------------------------------------------------------------------- !> @code SUBROUTINE interp__fill_value_wrapper( td_var, & & id_rho, & & id_offset ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(INOUT) :: td_var INTEGER(I4), DIMENSION(:), INTENT(IN ), OPTIONAL :: id_rho INTEGER(I4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset ! local variable CHARACTER(LEN=lc) :: cl_method INTEGER(i4) , DIMENSION(:), ALLOCATABLE :: il_rho INTEGER(i4) , DIMENSION(2,2) :: il_offset ! loop indices !---------------------------------------------------------------- SELECT CASE(TRIM(td_var%c_interp(1))) CASE('cubic','linear','nearest') cl_method=TRIM(td_var%c_interp(1)) CASE DEFAULT CALL logger_warn("INTERP FILL: interpolation method unknown."//& & " use linear interpolation") cl_method='linear' ! update variable structure value td_var%c_interp(1)='linear' END SELECT ALLOCATE( il_rho(ip_maxdim) ) il_rho(:)=1 IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:) IF( ANY(il_rho(:) < 0) )THEN CALL logger_error("INTERP FILL VALUE: invalid "//& & " refinement factor ") ENDIF il_offset(:,:)=0 IF( PRESENT(id_offset) )THEN IF( ANY(SHAPE(id_offset(:,:)) /= (/2,2/)) )THEN CALL logger_error("INTERP FILL VALUE: invalid table of offset") ELSE il_offset(:,:)=id_offset(:,:) ENDIF ENDIF IF( (il_rho(jp_I) /= 1 .AND. td_var%t_dim(1)%l_use) .OR. & & (il_rho(jp_J) /= 1 .AND. td_var%t_dim(2)%l_use) .OR. & & (il_rho(jp_K) /= 1 .AND. td_var%t_dim(3)%l_use) )THEN CALL logger_info("INTERP FILL: interpolate "//TRIM(td_var%c_name)//& & " using "//TRIM(cl_method)//" method." ) CALL logger_info("INTERP FILL: refinement factor "//& & TRIM(fct_str(il_rho(jp_I)))//& & " "//TRIM(fct_str(il_rho(jp_J)))//& & " "//TRIM(fct_str(il_rho(jp_K))) ) CALL interp__fill_value( td_var, cl_method, & & il_rho(:), & & il_offset(:,:) ) SELECT CASE(TRIM(td_var%c_interp(2))) CASE('/rhoi') WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill ) td_var%d_value(:,:,:,:) = & & td_var%d_value(:,:,:,:) / REAL(il_rho(jp_I),dp) END WHERE CASE('/rhoj') WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill ) td_var%d_value(:,:,:,:) = & & td_var%d_value(:,:,:,:) / REAL(il_rho(jp_J),dp) END WHERE CASE('/rhok') WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill ) td_var%d_value(:,:,:,:) = & & td_var%d_value(:,:,:,:) / REAL(il_rho(jp_K),dp) END WHERE CASE('*rhoi') WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill ) td_var%d_value(:,:,:,:) = & & td_var%d_value(:,:,:,:) * REAL(il_rho(jp_I),dp) END WHERE CASE('*rhoj') WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill ) td_var%d_value(:,:,:,:) = & & td_var%d_value(:,:,:,:) * REAL(il_rho(jp_J),dp) END WHERE CASE('*rhok') WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill ) td_var%d_value(:,:,:,:) = & & td_var%d_value(:,:,:,:) * REAL(il_rho(jp_K),dp) END WHERE CASE DEFAULT td_var%c_interp(2)='' END SELECT ENDIF END SUBROUTINE interp__fill_value_wrapper !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine interpolate value over mixed grid. !> !> @details !> !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[inout] td_var : variable !> @param[in] cd_method : interpolation method !> @param[in] id_rho : table of refinment factor !------------------------------------------------------------------- !> @code SUBROUTINE interp__fill_value( td_var, cd_method, & & id_rho, & & id_offset ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(INOUT) :: td_var CHARACTER(LEN=*) , INTENT(IN ) :: cd_method INTEGER(I4) , DIMENSION(:) , INTENT(IN ) :: id_rho INTEGER(I4) , DIMENSION(:,:), INTENT(IN ) :: id_offset ! local variable CHARACTER(LEN=lc) :: cl_interp INTEGER(I4) , DIMENSION(:) , ALLOCATABLE :: il_rho INTEGER(i4) , DIMENSION(:,:,:) , ALLOCATABLE :: il_detect REAL(dp) :: dl_min REAL(dp) :: dl_max LOGICAL , DIMENSION(3) :: ll_even LOGICAL :: ll_discont TYPE(TVAR) :: tl_mix TYPE(TATT) :: tl_att ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj INTEGER(i4) :: jk INTEGER(i4) :: jl !---------------------------------------------------------------- !1- create mixed grid CALL interp_create_mixed_grid( td_var, tl_mix, id_rho(:) ) ! clean variable structure CALL var_clean(td_var) !2- detect point to be interpolated ALLOCATE( il_detect( tl_mix%t_dim(1)%i_len, & & tl_mix%t_dim(2)%i_len, & & tl_mix%t_dim(3)%i_len) ) il_detect(:,:,:)=0 il_detect(:,:,:)=interp_detect(tl_mix, id_rho(:) ) ! add attribute to variable cl_interp=fct_concat(tl_mix%c_interp(:)) tl_att=att_init('interpolation',cl_interp) CALL var_move_att(tl_mix, tl_att) ! special case for even refinement on ARAKAWA-C grid ll_even(:)=.FALSE. IF( MOD(id_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE. IF( MOD(id_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE. IF( MOD(id_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE. SELECT CASE(TRIM(tl_mix%c_point)) CASE('U') ll_even(1)=.FALSE. CASE('V') ll_even(2)=.FALSE. CASE('F') ll_even(:)=.FALSE. END SELECT ALLOCATE(il_rho(ip_maxdim)) il_rho(:)=id_rho(:) IF( ll_even(1) ) il_rho(jp_I)=id_rho(jp_I)+1 IF( ll_even(2) ) il_rho(jp_J)=id_rho(jp_J)+1 IF( ll_even(3) ) il_rho(jp_K)=id_rho(jp_K)+1 ! special case for longitude ll_discont=.FALSE. IF( TRIM(tl_mix%c_units) == 'degrees_east' )THEN dl_min=MINVAL( tl_mix%d_value(:,:,:,:), & & tl_mix%d_value(:,:,:,:)/=tl_mix%d_fill) dl_max=MAXVAL( tl_mix%d_value(:,:,:,:), & & tl_mix%d_value(:,:,:,:)/=tl_mix%d_fill) IF( dl_min < -170_dp .AND. dl_max > 170_dp .OR. & & dl_min < 10_dp .AND. dl_max > 350_dp )THEN ll_discont=.TRUE. ENDIF ENDIF !3- interpolate DO jl=1,tl_mix%t_dim(4)%i_len IF( il_rho(jp_K) /= 1 )THEN ! CALL interp__3D(tl_mix%d_value(:,:,:,jl), tl_mix%d_fill, & ! & il_detect(:,:,:), cd_method, & ! & il_rhoi, il_rhoj, il_rhok, & ! & ll_even(:), ll_discont ) CALL logger_error("INTERP FILL: can not interpolate "//& & "vertically for now ") ENDIF IF( ANY(il_detect(:,:,:)==1) )THEN ! I-J plan DO jk=1,tl_mix%t_dim(3)%i_len CALL interp__2D(tl_mix%d_value(:,:,jk,jl), tl_mix%d_fill, & & il_detect(:,:,jk), cd_method, & & il_rho(jp_I), il_rho(jp_J), & & ll_even(1:2), ll_discont) ENDDO IF( ALL(il_detect(:,:,:)==0) ) CYCLE IF( il_rho(jp_K) /= 1 )THEN ! I-K plan DO jj=1,tl_mix%t_dim(2)%i_len CALL interp__2D(tl_mix%d_value(:,jj,:,jl), tl_mix%d_fill, & & il_detect(:,jj,:), cd_method, & & il_rho(jp_J), il_rho(jp_K), & & ll_even(1:3:2), ll_discont ) ENDDO IF( ALL(il_detect(:,:,:)==0) ) CYCLE ! J-K plan DO ji=1,tl_mix%t_dim(1)%i_len CALL interp__2D(tl_mix%d_value(ji,:,:,jl), tl_mix%d_fill, & & il_detect(ji,:,:), cd_method, & & il_rho(jp_J), il_rho(jp_K), & & ll_even(2:3), ll_discont ) ENDDO ENDIF IF( ANY(il_detect(:,:,:)==1) )THEN ! I direction DO jk=1,tl_mix%t_dim(3)%i_len DO jj=1,tl_mix%t_dim(2)%i_len CALL interp__1D( tl_mix%d_value(:,jj,jk,jl), & & tl_mix%d_fill, & & il_detect(:,jj,jk), cd_method, & & il_rho(jp_I), & & ll_even(1), ll_discont ) ENDDO ENDDO IF( ALL(il_detect(:,:,:)==0) ) CYCLE ! J direction DO jk=1,tl_mix%t_dim(3)%i_len DO ji=1,tl_mix%t_dim(1)%i_len CALL interp__1D( tl_mix%d_value(ji,:,jk,jl), & & tl_mix%d_fill, & & il_detect(ji,:,jk), cd_method, & & il_rho(jp_J), & & ll_even(2), ll_discont ) ENDDO ENDDO IF( il_rho(jp_K) /= 1 )THEN IF( ALL(il_detect(:,:,:)==0) ) CYCLE ! K direction DO jj=1,tl_mix%t_dim(2)%i_len DO ji=1,tl_mix%t_dim(1)%i_len CALL interp__1D( tl_mix%d_value(ji,jj,:,jl), & & tl_mix%d_fill, & & il_detect(ji,jj,:), cd_method, & & il_rho(jp_K), & & ll_even(3), ll_discont ) ENDDO ENDDO ENDIF ENDIF ENDIF ENDDO IF( ANY(il_detect(:,:,:)==1) )THEN CALL logger_warn("INTERP FILL: some points can not be interpolated "//& & "for variable "//TRIM(tl_mix%c_name) ) ENDIF DEALLOCATE(il_detect) !4- save useful domain (remove offset) CALL interp_clean_mixed_grid( tl_mix, td_var, & & id_rho(:), & & id_offset(:,:) ) ! clean variable structure CALL var_clean(tl_mix) END SUBROUTINE interp__fill_value !> @endcode ! !------------------------------------------------------------------- ! !> @brief ! !> This subroutine ! !> ! !> @details ! !> ! !> @author J.Paul ! !> - Nov, 2013- Initial Version ! ! ! !> @param[in] ! !> @param[out] ! !------------------------------------------------------------------- ! !> @code ! SUBROUTINE interp__3D( dd_value, dd_fill, & ! & id_detect, cd_method, & ! & id_rhoi, id_rhoj, id_rhok, & ! & ld_even, ld_discont ) ! IMPLICIT NONE ! ! Argument ! REAL(dp) , DIMENSION(:,:,:), INTENT(INOUT) :: dd_value ! REAL(dp) , INTENT(IN ) :: dd_fill ! INTEGER(I4) , DIMENSION(:,:,:), INTENT(INOUT) :: id_detect ! CHARACTER(LEN=*) , INTENT(IN ) :: cd_method ! INTEGER(I4) , INTENT(IN ) :: id_rhoi ! INTEGER(I4) , INTENT(IN ) :: id_rhoj ! INTEGER(I4) , INTENT(IN ) :: id_rhok ! LOGICAL , DIMENSION(:) , INTENT(IN ) :: ld_even ! LOGICAL , INTENT(IN ), OPTIONAL :: ld_discont ! ! ! local variable ! INTEGER(i4), DIMENSION(3) :: il_shape ! INTEGER(i4), DIMENSION(3) :: il_dim ! ! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coarse ! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx ! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdy ! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_d2fdxy ! ! LOGICAL :: ll_discont ! ! ! loop indices !! INTEGER(i4) :: ji !! INTEGER(i4) :: jj !! INTEGER(i4) :: ii !! INTEGER(i4) :: ij ! !---------------------------------------------------------------- ! ll_discont=.FALSE. ! IF( PRESENT(ld_discont) ) ll_discont=ld_discont ! ! il_shape(:)=SHAPE(dd_value) ! il_dim(:)=(il_shape(:)-1)/2 ! ! ALLOCATE( dl_coarse(il_dim(1),il_dim(2),il_dim(3)) ) ! ! value on coarse grid ! dl_coarse(:,:,:)=dd_value( 1:il_shape(1):id_rhoi, & ! & 1:il_shape(2):id_rhoj, & ! & 1:il_shape(3):id_rhok ) ! SELECT CASE(TRIM(cd_method)) ! ! CASE('cubic') ! ! ALLOCATE( dl_dfdx( il_dim(1),il_dim(2),il_dim(3)), & ! & dl_dfdy( il_dim(1),il_dim(2),il_dim(3)), & ! & dl_d2fdxy(il_dim(1),il_dim(2),il_dim(3)) ) ! !! ! compute derivative on coarse grid !! dl_dfdx(:,:)=extrap_deriv_2D(dl_coarse(:,:,:), dd_fill, 'I') !! dl_dfdy(:,:)=extrap_deriv_2D(dl_coarse(:,:,:), dd_fill, 'J') !! !! ! compute cross derivative on coarse grid !! dl_d2fdxy(:,:)=extrap_deriv_2D(dl_dfdx(:,:,:), dd_fill, 'J') !! !! DO jj=1,il_shape(2)-1,id_rhoj !! ij=((jj-1)/id_rhoj)+1 !! DO ji=1,il_shape(1)-1,id_rhoi !! ii=((ji-1)/id_rhoi)+1 !! !! IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE !! ! compute bicubic coefficient !! dl_coef(:)=interp__2D_cubic_coef(dl_coarse(ii:ii+1,ij:ij+1),& !! & dl_dfdx( ii:ii+1,ij:ij+1),& !! & dl_dfdy( ii:ii+1,ij:ij+1),& !! & dl_d2fdxy(ii:ii+1,ij:ij+1) ) !! !! ! compute value on detetected point !! CALL interp__2D_cubic_fill(dl_coef(:), & !! & dd_value( ji:ji+id_rhoi, & !! & jj:jj+id_rhoj ), & !! & id_detect(ji:ji+id_rhoi, & !! & jj:jj+id_rhoj ) ) !! !! ENDDO !! ENDDO ! ! DEALLOCATE( dl_dfdx, & ! & dl_dfdy, & ! & dl_d2fdxy ) ! ! CASE('nearest') ! ! CASE DEFAULT ! linear ! !! DO jj=1,il_shape(2)-1,id_rhoj !! ij=((jj-1)/id_rhoj)+1 !! DO ji=1,il_shape(1)-1,id_rhoi !! ii=((ji-1)/id_rhoi)+1 !! !! IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE !! ! compute bilinear coefficient !! dl_coef(:)=interp__2D_linear_coef(dl_coarse(ii:ii+1,ij:ij+1)) !! !! ! compute value on detetected point !! CALL interp__2D_linear_fill(dl_coef(:), & !! & dd_value( ji:ji+id_rhoi, & !! & jj:jj+id_rhoj ), & !! & id_detect(ji:ji+id_rhoi, & !! & jj:jj+id_rhoj ) ) !! !! ENDDO !! ENDDO ! ! END SELECT ! ! DEALLOCATE( dl_coarse ) ! ! END SUBROUTINE interp__3D ! !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code SUBROUTINE interp__2D( dd_value, dd_fill, & & id_detect, cd_method, & & id_rhoi, id_rhoj, & & ld_even, ld_discont ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_value REAL(dp) , INTENT(IN ) :: dd_fill INTEGER(I4) , DIMENSION(:,:), INTENT(INOUT) :: id_detect CHARACTER(LEN=*) , INTENT(IN ) :: cd_method INTEGER(I4) , INTENT(IN ) :: id_rhoi INTEGER(I4) , INTENT(IN ) :: id_rhoj LOGICAL , DIMENSION(:) , INTENT(IN ) :: ld_even LOGICAL , INTENT(IN ), OPTIONAL :: ld_discont ! local variable INTEGER(I4) :: il_xextra INTEGER(I4) :: il_yextra INTEGER(i4), DIMENSION(2) :: il_shape INTEGER(i4), DIMENSION(2) :: il_dim REAL(dp) :: dl_min REAL(dp) :: dl_max REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_coef REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_coarse REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_tmp REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_dfdx REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_dfdy REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_d2fdxy LOGICAL :: ll_discont ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj INTEGER(i4) :: ii INTEGER(i4) :: ij !---------------------------------------------------------------- ll_discont=.FALSE. IF( PRESENT(ld_discont) ) ll_discont=ld_discont CALL logger_debug("INTERP 2D: interpolation method "//TRIM(cd_method)//& & " discont "//TRIM(fct_str(ll_discont)) ) il_shape(:)=SHAPE(dd_value) ! compute coarse grid dimension il_xextra=id_rhoi-1 il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi il_yextra=id_rhoj-1 il_dim(2)=(il_shape(2)+il_yextra)/id_rhoj ALLOCATE( dl_coarse(il_dim(1),il_dim(2)) ) ! value on coarse grid dl_coarse(:,:)=dd_value( 1:il_shape(1):id_rhoi, & & 1:il_shape(2):id_rhoj ) SELECT CASE(TRIM(cd_method)) CASE('cubic') ALLOCATE( dl_dfdx( il_dim(1),il_dim(2)), & & dl_dfdy( il_dim(1),il_dim(2)), & & dl_d2fdxy(il_dim(1),il_dim(2)) ) ! compute derivative on coarse grid dl_dfdx(:,:)=extrap_deriv_2D(dl_coarse(:,:), dd_fill, 'I', ll_discont) dl_dfdy(:,:)=extrap_deriv_2D(dl_coarse(:,:), dd_fill, 'J', ll_discont) ! compute cross derivative on coarse grid dl_d2fdxy(:,:)=extrap_deriv_2D(dl_dfdx(:,:), dd_fill, 'J', ll_discont) ALLOCATE( dl_tmp(2,2) ) ALLOCATE( dl_coef(16) ) DO jj=1,il_shape(2)-1,id_rhoj ij=((jj-1)/id_rhoj)+1 DO ji=1,il_shape(1)-1,id_rhoi ii=((ji-1)/id_rhoi)+1 IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill) .OR. & & ANY( dl_dfdx(ii:ii+1,ij:ij+1)==dd_fill) .OR. & & ANY( dl_dfdy(ii:ii+1,ij:ij+1)==dd_fill) .OR. & & ANY(dl_d2fdxy(ii:ii+1,ij:ij+1)==dd_fill) ) CYCLE dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1) IF( ll_discont )THEN dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill ) dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill ) IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN WHERE( dl_tmp(:,:) < 0_dp ) dl_tmp(:,:) = dl_tmp(:,:)+360._dp END WHERE ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN WHERE( dl_tmp(:,:) > 180_dp ) dl_tmp(:,:) = dl_tmp(:,:)-180._dp END WHERE ENDIF ENDIF ! compute bicubic coefficient dl_coef(:)=interp__2D_cubic_coef(dl_tmp(:,:),& & dl_dfdx( ii:ii+1,ij:ij+1),& & dl_dfdy( ii:ii+1,ij:ij+1),& & dl_d2fdxy( ii:ii+1,ij:ij+1),& & dd_fill ) ! compute value on detetected point CALL interp__2D_cubic_fill(dd_value( ji:ji+id_rhoi, & & jj:jj+id_rhoj ), & & id_detect(ji:ji+id_rhoi, & & jj:jj+id_rhoj ), & & dl_coef(:), dd_fill, & & ld_even(:), id_rhoi, id_rhoj ) IF( ll_discont )THEN WHERE( dd_value( ji:ji+id_rhoi, & & jj:jj+id_rhoj ) >= 180._dp .AND. & & dd_value( ji:ji+id_rhoi, & & jj:jj+id_rhoj ) /= dd_fill ) dd_value( ji:ji+id_rhoi, & & jj:jj+id_rhoj ) = & & dd_value( ji:ji+id_rhoi, & & jj:jj+id_rhoj ) - 360._dp END WHERE ENDIF ENDDO ENDDO DEALLOCATE(dl_coef) DEALLOCATE(dl_tmp ) DEALLOCATE(dl_dfdx, & & dl_dfdy, & & dl_d2fdxy ) CASE('nearest') DO jj=1,il_shape(2)-1,id_rhoj ij=((jj-1)/id_rhoj)+1 DO ji=1,il_shape(1)-1,id_rhoi ii=((ji-1)/id_rhoi)+1 ! compute value on detetected point CALL interp__2D_nearest_fill(dd_value( ji:ji+id_rhoi, & & jj:jj+id_rhoj ), & & id_detect(ji:ji+id_rhoi, & & jj:jj+id_rhoj ) ) ENDDO ENDDO CASE DEFAULT ! linear ALLOCATE( dl_coef(4) ) ALLOCATE( dl_tmp(2,2) ) DO jj=1,il_shape(2)-1,id_rhoj ij=((jj-1)/id_rhoj)+1 DO ji=1,il_shape(1)-1,id_rhoi ii=((ji-1)/id_rhoi)+1 IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1) IF( ll_discont )THEN dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill ) dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill ) IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN WHERE( dl_tmp(:,:) < 0_dp ) dl_tmp(:,:) = dl_tmp(:,:)+360._dp END WHERE ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN WHERE( dl_tmp(:,:) > 180_dp ) dl_tmp(:,:) = dl_tmp(:,:)-180._dp END WHERE ENDIF ENDIF ! compute bilinear coefficient dl_coef(:)=interp__2D_linear_coef(dl_tmp(:,:), dd_fill) ! compute value on detetected point CALL interp__2D_linear_fill(dd_value( ji:ji+id_rhoi, & & jj:jj+id_rhoj ), & & id_detect(ji:ji+id_rhoi, & & jj:jj+id_rhoj ), & & dl_coef(:), dd_fill, & & ld_even(:), id_rhoi, id_rhoj ) IF( ll_discont )THEN WHERE( dd_value( ji:ji+id_rhoi, & & jj:jj+id_rhoj ) >= 180._dp .AND. & & dd_value( ji:ji+id_rhoi, & & jj:jj+id_rhoj ) /= dd_fill ) dd_value( ji:ji+id_rhoi, & & jj:jj+id_rhoj ) = & & dd_value( ji:ji+id_rhoi, & & jj:jj+id_rhoj ) - 360._dp END WHERE ENDIF ENDDO ENDDO DEALLOCATE(dl_coef) DEALLOCATE(dl_tmp ) END SELECT DEALLOCATE( dl_coarse ) END SUBROUTINE interp__2D !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code SUBROUTINE interp__1D( dd_value, dd_fill, & & id_detect, cd_method, & & id_rhoi, & & ld_even, ld_discont ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value REAL(dp) , INTENT(IN ) :: dd_fill INTEGER(I4) , DIMENSION(:), INTENT(INOUT) :: id_detect CHARACTER(LEN=*) , INTENT(IN ) :: cd_method INTEGER(I4) , INTENT(IN ) :: id_rhoi LOGICAL , INTENT(IN ) :: ld_even LOGICAL , INTENT(IN ), OPTIONAL :: ld_discont ! local variable INTEGER(i4), DIMENSION(1) :: il_shape INTEGER(i4), DIMENSION(1) :: il_dim INTEGER(I4) :: il_xextra REAL(dp) :: dl_min REAL(dp) :: dl_max REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_coarse REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_tmp REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_dfdx REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_coef LOGICAL :: ll_discont ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: ii !---------------------------------------------------------------- ll_discont=.FALSE. IF( PRESENT(ld_discont) ) ll_discont=ld_discont CALL logger_debug("INTERP 1D: interpolation method "//TRIM(cd_method)//& & " discont "//TRIM(fct_str(ll_discont)) ) il_shape(:)=SHAPE(dd_value) il_xextra=id_rhoi-1 il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi ALLOCATE( dl_coarse(il_dim(1)) ) ! value on coarse grid dl_coarse(:)=dd_value( 1:il_shape(1):id_rhoi ) SELECT CASE(TRIM(cd_method)) CASE('cubic') ALLOCATE( dl_dfdx( il_dim(1)) ) ! compute derivative on coarse grid dl_dfdx(:)=extrap_deriv_1D(dl_coarse(:), dd_fill, ll_discont) ALLOCATE( dl_coef(4)) ALLOCATE( dl_tmp(2) ) DO ji=1,il_shape(1)-1,id_rhoi ii=((ji-1)/id_rhoi)+1 IF( ANY( dl_tmp(:)==dd_fill ) .OR. & & ANY(dl_dfdx(:)==dd_fill ) ) CYCLE dl_tmp(:)=dl_coarse(ii:ii+1) IF( ll_discont )THEN dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill ) dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill ) IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN WHERE( dl_tmp(:) < 0_dp ) dl_tmp(:) = dl_tmp(:)+360._dp END WHERE ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN WHERE( dl_tmp(:) > 180_dp ) dl_tmp(:) = dl_tmp(:)-180._dp END WHERE ENDIF ENDIF ! compute bicubic coefficient dl_coef(:)=interp__1D_cubic_coef(dl_tmp(:),& & dl_dfdx( ii:ii+1), & & dd_fill ) ! compute value on detetected point CALL interp__1D_cubic_fill(dd_value( ji:ji+id_rhoi),& & id_detect(ji:ji+id_rhoi),& & dl_coef(:), dd_fill, & & ld_even, id_rhoi ) IF( ll_discont )THEN WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. & & dd_value( ji:ji+id_rhoi ) /= dd_fill ) dd_value( ji:ji+id_rhoi ) = & & dd_value( ji:ji+id_rhoi ) - 360._dp END WHERE ENDIF ENDDO DEALLOCATE(dl_coef) DEALLOCATE(dl_tmp ) CASE('nearest') DO ji=1,il_shape(1)-1,id_rhoi ii=((ji-1)/id_rhoi)+1 ! compute value on detetected point CALL interp__1D_nearest_fill(dd_value( ji:ji+id_rhoi), & & id_detect(ji:ji+id_rhoi) ) ENDDO CASE DEFAULT ! linear ALLOCATE(dl_coef(2)) ALLOCATE( dl_tmp(2) ) DO ji=1,il_shape(1)-1,id_rhoi ii=((ji-1)/id_rhoi)+1 IF( ANY(dl_coarse(ii:ii+1)==dd_fill ) ) CYCLE dl_tmp(:)=dl_coarse(ii:ii+1) IF( ll_discont )THEN dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill ) dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill ) IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN WHERE( dl_tmp(:) < 0_dp ) dl_tmp(:) = dl_tmp(:)+360._dp END WHERE ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN WHERE( dl_tmp(:) > 180_dp ) dl_tmp(:) = dl_tmp(:)-180._dp END WHERE ENDIF ENDIF ! compute bilinear coefficient dl_coef(:)=interp__1D_linear_coef(dl_tmp(:), dd_fill) ! compute value on detetected point CALL interp__1D_linear_fill( dd_value( ji:ji+id_rhoi),& & id_detect(ji:ji+id_rhoi),& & dl_coef(:), dd_fill, & & ld_even, id_rhoi ) IF( ll_discont )THEN WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. & & dd_value( ji:ji+id_rhoi ) /= dd_fill ) dd_value( ji:ji+id_rhoi ) = & & dd_value( ji:ji+id_rhoi ) - 360._dp END WHERE ENDIF ENDDO DEALLOCATE(dl_coef) END SELECT DEALLOCATE( dl_coarse ) END SUBROUTINE interp__1D !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code FUNCTION interp__2D_cubic_coef( dd_value, & & dd_dfdx, & & dd_dfdy, & & dd_d2fdxy,& & dd_fill ) IMPLICIT NONE ! Argument REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_value REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdx REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdy REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_d2fdxy REAL(dp) , INTENT(IN) :: dd_fill ! function REAL(dp), DIMENSION(16) :: interp__2D_cubic_coef ! local variable REAL(dp), DIMENSION(16,16), PARAMETER :: dp_matrix = RESHAPE( & & (/ 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,& 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,& -3 , 3 , 0 , 0 ,-2 ,-1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,& 2 ,-2 , 0 , 0 , 1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,& 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,& 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 ,& 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,-3 , 3 , 0 , 0 ,-2 ,-1 , 0 , 0 ,& 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 ,-2 , 0 , 0 , 1 , 1 , 0 , 0 ,& -3 , 0 , 3 , 0 , 0 , 0 , 0 , 0 ,-2 , 0 ,-1 , 0 , 0 , 0 , 0 , 0 ,& 0 , 0 , 0 , 0 ,-3 , 0 , 3 , 0 , 0 , 0 , 0 , 0 ,-2 , 0 ,-1 , 0 ,& 9 ,-9 ,-9 , 9 , 6 , 3 ,-6 ,-3 , 6 ,-6 , 3 ,-3 , 4 , 2 , 2 , 1 ,& -6 , 6 , 6 ,-6 ,-3 ,-3 , 3 , 3 ,-4 , 4 ,-2 , 2 ,-2 ,-2 ,-1 ,-1 ,& 2 , 0 ,-2 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 ,& 0 , 0 , 0 , 0 , 2 , 0 ,-2 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 ,& -6 , 6 , 6 ,-6 ,-4 ,-2 , 4 , 2 ,-3 , 3 ,-3 , 3 ,-2 ,-1 ,-2 ,-1 ,& 4 ,-4 ,-4 , 4 , 2 , 2 ,-2 ,-2 , 2 ,-2 , 2 ,-2 , 1 , 1 , 1 , 1 /), & & (/ 16, 16 /) ) REAL(dp), DIMENSION(16) :: dl_vect !---------------------------------------------------------------- ! init interp__2D_cubic_coef(:)=dd_fill IF( ANY(SHAPE( dd_value(:,:))/= 2) .OR. & & ANY(SHAPE( dd_dfdx(:,:))/= 2) .OR. & & ANY(SHAPE( dd_dfdy(:,:))/= 2) .OR. & & ANY(SHAPE(dd_d2fdxy(:,:))/= 2) )THEN CALL logger_error("INTERP CUBIC COEF: invalid dimension of "//& & "input tables. shape should be (/2,2/)") ELSEIF( ANY( dd_value(:,:) == dd_fill) .OR. & & ANY( dd_dfdx(:,:) == dd_fill) .OR. & & ANY( dd_dfdy(:,:) == dd_fill) .OR. & & ANY(dd_d2fdxy(:,:) == dd_fill) )THEN CALL logger_warn("INTERP CUBIC COEF: fill value detected. "//& & "can not compute coefficient ") ELSE dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. ) dl_vect( 5: 8)=PACK(dd_dfdx(:,:),.TRUE. ) dl_vect( 9:12)=PACK(dd_dfdy(:,:),.TRUE. ) dl_vect(13:16)=PACK(dd_d2fdxy(:,:),.TRUE. ) interp__2D_cubic_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:)) ENDIF END FUNCTION interp__2D_cubic_coef !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code SUBROUTINE interp__2D_cubic_fill( dd_value, id_detect, dd_coef, dd_fill, & & ld_even, id_rhoi, id_rhoj ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_value INTEGER(i4) , DIMENSION(:,:), INTENT(INOUT) :: id_detect REAL(dp) , DIMENSION(:) , INTENT(IN ) :: dd_coef REAL(dp) , INTENT(IN ) :: dd_fill LOGICAL , DIMENSION(:), INTENT(IN ) :: ld_even INTEGER(I4) , INTENT(IN ) :: id_rhoi INTEGER(I4) , INTENT(IN ) :: id_rhoj ! local variable INTEGER(i4), DIMENSION(2) :: il_shape REAL(dp) , DIMENSION(16) :: dl_vect REAL(dp) :: dl_dx REAL(dp) :: dl_dy REAL(dp) :: dl_x REAL(dp) :: dl_x2 REAL(dp) :: dl_x3 REAL(dp) :: dl_y REAL(dp) :: dl_y2 REAL(dp) :: dl_y3 ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj !---------------------------------------------------------------- IF( SIZE(dd_coef(:)) /= 16 )THEN CALL logger_error("INTERP CUBIC FILL: invalid dimension of "//& & "coef table. shape should be (/16/)") ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN CALL logger_error("INTERP CUBIC FILL: fill value detected in coef . "//& & "can not compute interpolation.") ELSE il_shape(:)=SHAPE(dd_value(:,:)) dl_dx=1./REAL(id_rhoi) dl_dy=1./REAL(id_rhoj) DO jj=1,il_shape(2) IF( ld_even(2) )THEN dl_y=(jj-1)*dl_dy - dl_dy*0.5 ELSE ! odd refinement dl_y=(jj-1)*dl_dy ENDIF dl_y2=dl_y*dl_y dl_y3=dl_y2*dl_y DO ji=1,il_shape(1) IF(id_detect(ji,jj)==1)THEN IF( ld_even(1) )THEN dl_x=(ji-1)*dl_dx - dl_dx*0.5 ELSE ! odd refinement dl_x=(ji-1)*dl_dx ENDIF dl_x2=dl_x*dl_x dl_x3=dl_x2*dl_x dl_vect(:)=(/1._dp, dl_x , dl_x2 , dl_x3 , & & dl_y , dl_x*dl_y , dl_x2*dl_y , dl_x3*dl_y , & & dl_y2, dl_x*dl_y2, dl_x2*dl_y2, dl_x3*dl_y2, & & dl_y3, dl_x*dl_y3, dl_x2*dl_y3, dl_x3*dl_y3 /) dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dl_vect(:)) id_detect(ji,jj)=0 ENDIF ENDDO ENDDO ENDIF END SUBROUTINE interp__2D_cubic_fill !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code FUNCTION interp__2D_linear_coef( dd_value, dd_fill ) IMPLICIT NONE ! Argument REAL(dp), DIMENSION(:,:) , INTENT(IN) :: dd_value REAL(dp) , INTENT(IN) :: dd_fill ! function REAL(dp), DIMENSION(4) :: interp__2D_linear_coef ! local variable REAL(dp), DIMENSION(4,4), PARAMETER :: dp_matrix = RESHAPE( & & (/ 1 , 0 , 0 , 0 ,& -1 , 1 , 0 , 0 ,& -1 , 0 , 1 , 0 ,& 1 ,-1 ,-1 , 1 /), & & (/ 4, 4 /) ) REAL(dp), DIMENSION(4) :: dl_vect !---------------------------------------------------------------- IF( ANY(SHAPE(dd_value(:,:))/= 2) )THEN CALL logger_error("INTERP LINEAR COEF: invalid dimension of "//& & "input tables. shape should be (/2,2/)") ELSEIF( ANY(dd_value(:,:)==dd_fill) )THEN CALL logger_error("INTERP LINEAR COEF: fill value detected. "//& & "can not compute coefficient.") ELSE dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. ) interp__2D_linear_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:)) ENDIF END FUNCTION interp__2D_linear_coef !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code SUBROUTINE interp__2D_linear_fill( dd_value, id_detect, dd_coef, dd_fill, & & ld_even, id_rhoi, id_rhoj ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_value INTEGER(i4) , DIMENSION(:,:), INTENT(INOUT) :: id_detect REAL(dp) , DIMENSION(:) , INTENT(IN ) :: dd_coef REAL(dp) , INTENT(IN ) :: dd_fill LOGICAL , DIMENSION(:), INTENT(IN ) :: ld_even INTEGER(I4) , INTENT(IN ) :: id_rhoi INTEGER(I4) , INTENT(IN ) :: id_rhoj ! local variable INTEGER(i4), DIMENSION(2) :: il_shape REAL(dp) , DIMENSION(4) :: dl_vect REAL(dp) :: dl_dx REAL(dp) :: dl_dy REAL(dp) :: dl_x REAL(dp) :: dl_y ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj !---------------------------------------------------------------- IF( SIZE(dd_coef(:)) /= 4 )THEN CALL logger_error("INTERP LINEAR FILL: invalid dimension of "//& & "coef table. shape should be (/4/)") ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN CALL logger_error("INTERP LINEAR FILL: fill value detected in coef. "//& & "can not compute interpolation.") ELSE il_shape(:)=SHAPE(dd_value(:,:)) dl_dx=1./REAL(id_rhoi) dl_dy=1./REAL(id_rhoj) DO jj=1,il_shape(2) IF( ld_even(2) )THEN dl_y=(jj-1)*dl_dy - dl_dy*0.5 ELSE ! odd refinement dl_y=(jj-1)*dl_dy ENDIF DO ji=1,il_shape(1) IF(id_detect(ji,jj)==1)THEN IF( ld_even(1) )THEN dl_x=(ji-1)*dl_dx - dl_dx*0.5 ELSE ! odd refinement dl_x=(ji-1)*dl_dx ENDIF dl_vect(:)=(/1._dp, dl_x, dl_y, dl_x*dl_y /) dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dl_vect(:)) id_detect(ji,jj)=0 ENDIF ENDDO ENDDO ENDIF END SUBROUTINE interp__2D_linear_fill !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code SUBROUTINE interp__2D_nearest_fill( dd_value, id_detect ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:,:) , INTENT(INOUT) :: dd_value INTEGER(i4), DIMENSION(:,:) , INTENT(INOUT) :: id_detect ! local variable INTEGER(i4), DIMENSION(2) :: il_shape INTEGER(i4) :: il_i1 INTEGER(i4) :: il_i2 INTEGER(i4) :: il_j1 INTEGER(i4) :: il_j2 INTEGER(i4) :: il_half1 INTEGER(i4) :: il_half2 ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj !---------------------------------------------------------------- il_shape(:)=SHAPE(dd_value(:,:)) il_i1=1 il_i2=il_shape(1) il_j1=1 il_j2=il_shape(2) il_half1=CEILING(il_shape(1)*0.5) il_half2=CEILING(il_shape(2)*0.5) DO jj=1,il_half2 DO ji=1,il_half1 ! lower left point IF(id_detect(ji,jj)==1)THEN dd_value( ji,jj)=dd_value(il_i1,il_j1) id_detect(ji,jj)=0 ENDIF ! lower right point IF(id_detect(il_shape(1)-ji+1,jj)==1)THEN dd_value( il_shape(1)-ji+1,jj)=dd_value(il_i2,il_j1) id_detect(il_shape(1)-ji+1,jj)=0 ENDIF ! upper left point IF(id_detect(ji,il_shape(2)-jj+1)==1)THEN dd_value( ji,il_shape(2)-jj+1)=dd_value(il_i1,il_j2) id_detect(ji,il_shape(2)-jj+1)=0 ENDIF ! upper right point IF(id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)==1)THEN dd_value( il_shape(1)-ji+1,il_shape(2)-jj+1)=dd_value(il_i2,il_j2) id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)=0 ENDIF ENDDO ENDDO END SUBROUTINE interp__2D_nearest_fill !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code FUNCTION interp__1D_cubic_coef( dd_value, & & dd_dfdx, & & dd_fill ) IMPLICIT NONE ! Argument REAL(dp), DIMENSION(:) , INTENT(IN) :: dd_value REAL(dp), DIMENSION(:) , INTENT(IN) :: dd_dfdx REAL(dp) , INTENT(IN) :: dd_fill ! function REAL(dp), DIMENSION(4) :: interp__1D_cubic_coef ! local variable REAL(dp), DIMENSION(4,4), PARAMETER :: dp_matrix = RESHAPE( & & (/ 1 , 0 , 0 , 0 ,& -1 , 1 , 0 , 0 ,& -3 , 3 ,-2 ,-1 ,& 2 ,-2 , 1 , 1 /), & & (/ 4, 4 /) ) REAL(dp), DIMENSION(4) :: dl_vect !---------------------------------------------------------------- IF( SIZE(dd_value(:))/= 2 .OR. & & SIZE(dd_dfdx(:) )/= 2 )THEN CALL logger_error("INTERP CUBIC COEF: invalid dimension of "//& & "input tables. shape should be (/2,2/)") ELSEIF( ANY(dd_value(:)==dd_fill) .OR. & & ANY( dd_dfdx(:)==dd_fill) )THEN CALL logger_error("INTERP CUBIC COEF: fill value detected. "//& & "can not compute coefficient.") ELSE dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. ) dl_vect( 3: 4)=PACK(dd_dfdx(:),.TRUE. ) interp__1D_cubic_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:)) ENDIF END FUNCTION interp__1D_cubic_coef !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code SUBROUTINE interp__1D_cubic_fill( dd_value, id_detect, dd_coef, dd_fill, & & ld_even, id_rhoi ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value INTEGER(i4) , DIMENSION(:), INTENT(INOUT) :: id_detect REAL(dp) , DIMENSION(:), INTENT(IN ) :: dd_coef REAL(dp) , INTENT(IN ) :: dd_fill LOGICAL , INTENT(IN ) :: ld_even INTEGER(I4) , INTENT(IN ) :: id_rhoi ! local variable INTEGER(i4), DIMENSION(1) :: il_shape REAL(dp) , DIMENSION(4) :: dl_vect REAL(dp) :: dl_dx REAL(dp) :: dl_x REAL(dp) :: dl_x2 REAL(dp) :: dl_x3 ! loop indices INTEGER(i4) :: ji !---------------------------------------------------------------- IF( SIZE(dd_coef(:)) /= 4 )THEN CALL logger_error("INTERP CUBIC FILL: invalid dimension of "//& & "coef table. shape should be (/4/)") !ELSEIF( ANY(dd_value(:)==dd_fill .AND. id_detect(:)==0 ) .OR. & ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN CALL logger_error("INTERP CUBIC FILL: fill value detected. "//& & "can not compute interpolation") ELSE il_shape(:)=SHAPE(dd_value(:)) dl_dx=1./REAL(id_rhoi) DO ji=1,il_shape(1) IF(id_detect(ji)==1)THEN IF( ld_even )THEN dl_x=(ji-1)*dl_dx - dl_dx*0.5 ELSE ! odd refinement dl_x=(ji-1)*dl_dx ENDIF dl_x2=dl_x*dl_x dl_x3=dl_x2*dl_x dl_vect(:)=(/1._dp, dl_x, dl_x2, dl_x3 /) dd_value(ji)=DOT_PRODUCT(dd_coef(:),dl_vect(:)) id_detect(ji)=0 ENDIF ENDDO ENDIF END SUBROUTINE interp__1D_cubic_fill !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code FUNCTION interp__1D_linear_coef( dd_value, dd_fill ) IMPLICIT NONE ! Argument REAL(dp), DIMENSION(:) , INTENT(IN) :: dd_value REAL(dp) , INTENT(IN) :: dd_fill ! function REAL(dp), DIMENSION(2) :: interp__1D_linear_coef ! local variable REAL(dp), DIMENSION(2,2), PARAMETER :: dp_matrix = RESHAPE( & & (/ 1 , 0 ,& -1 , 1 /), & & (/ 2, 2 /) ) REAL(dp), DIMENSION(2) :: dl_vect !---------------------------------------------------------------- IF( ANY(SHAPE(dd_value(:))/= 2) )THEN CALL logger_error("INTERP LINEAR COEF: invalid dimension of "//& & "input tables. shape should be (/2/)") ELSEIF( ANY(dd_value(:)==dd_fill) )THEN CALL logger_error("INTERP LINEAR COEF: fill value detected. "//& & "can not compute coefficient.") ELSE dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. ) interp__1D_linear_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:)) ENDIF END FUNCTION interp__1D_linear_coef !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code SUBROUTINE interp__1D_linear_fill( dd_value, id_detect, dd_coef, dd_fill, & & ld_even, id_rhoi ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value INTEGER(i4) , DIMENSION(:), INTENT(INOUT) :: id_detect REAL(dp) , DIMENSION(:), INTENT(IN ) :: dd_coef REAL(dp) , INTENT(IN ) :: dd_fill LOGICAL , INTENT(IN ) :: ld_even INTEGER(I4) , INTENT(IN ) :: id_rhoi ! local variable INTEGER(i4), DIMENSION(1) :: il_shape REAL(dp) , DIMENSION(2) :: dl_vect REAL(dp) :: dl_dx REAL(dp) :: dl_x ! loop indices INTEGER(i4) :: ji !---------------------------------------------------------------- IF( SIZE(dd_coef(:)) /= 2 )THEN CALL logger_error("INTERP LINEAR FILL: invalid dimension of "//& & "coef table. shape should be (/2/)") !ELSEIF( ANY(dd_value(:)==dd_fill .AND. id_detect(:)==0 ) .OR. & ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN CALL logger_error("INTERP LINEAR FILL: fill value detected. "//& & "can not compute interpolation") ELSE il_shape(:)=SHAPE(dd_value) dl_dx=1./REAL(id_rhoi) DO ji=1,il_shape(1) IF(id_detect(ji)==1)THEN IF( ld_even )THEN dl_x=(ji-1)*dl_dx - dl_dx*0.5 ELSE ! odd refinement dl_x=(ji-1)*dl_dx ENDIF dl_vect(:)=(/1._dp, dl_x /) dd_value(ji)=DOT_PRODUCT(dd_coef(:),dl_vect(:)) id_detect(ji)=0 ENDIF ENDDO ENDIF END SUBROUTINE interp__1D_linear_fill !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] !> @param[out] !------------------------------------------------------------------- !> @code SUBROUTINE interp__1D_nearest_fill( dd_value, id_detect ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value INTEGER(i4), DIMENSION(:), INTENT(INOUT) :: id_detect ! local variable INTEGER(i4), DIMENSION(1) :: il_shape INTEGER(i4) :: il_i1 INTEGER(i4) :: il_i2 INTEGER(i4) :: il_half1 ! loop indices INTEGER(i4) :: ji !---------------------------------------------------------------- il_shape(:)=SHAPE(dd_value) il_i1=1 il_i2=il_shape(1) il_half1=CEILING(il_shape(1)*0.5) DO ji=1,il_half1 ! lower left point IF(id_detect(ji)==1)THEN dd_value( ji)=dd_value(il_i1) id_detect(ji)=0 ENDIF ! lower right point IF(id_detect(il_shape(1)-ji+1)==1)THEN dd_value( il_shape(1)-ji+1)=dd_value(il_i2) id_detect(il_shape(1)-ji+1)=0 ENDIF ENDDO END SUBROUTINE interp__1D_nearest_fill !> @endcode END MODULE interp