!---------------------------------------------------------------------- ! NEMO system team, System and Interface for oceanic RElocable Nesting !---------------------------------------------------------------------- ! ! MODULE: extrap ! ! DESCRIPTION: !> @brief !> This module manage extrapolation. !> !> @details !> Extrapolation method to be used is specify inside variable !> strcuture, as array of string character.
!> - td_var\%c_extrap(1) string character is the interpolation name choose between: !> - 'dist_weight' !> - 'min_error' !> !> @note Extrapolation method could be specify for each variable in namelist _namvar_, !> defining string character _cn\_varinfo_. By default _dist_weight_.
!> Example: !> - cn_varinfo='varname1:dist_weight', 'varname2:min_error' !> !> to detect point to be extrapolated:
!> @code !> il_detect(:,:,:)=extrap_detect(td_var, [td_level], [id_offset,] [id_rho,] [id_ext]) !> @endcode !> - il_detect(:,:,:) is 3D array of point to be extrapolated !> - td_var is coarse grid variable to be extrapolated !> - td_level is fine grid array of level (see vgrid_get_level) [optional] !> - id_offset is array of offset between fine and coarse grid [optional] !> - id_rho is array of refinment factor [optional] !> - id_ext is array of number of points to be extrapolated [optional] !> !> to extrapolate variable:
!> @code !> CALL extrap_fill_value( td_var, [td_level], [id_offset], [id_rho], [id_iext], [id_jext], [id_kext], [id_radius], [id_maxiter]) !> @endcode !> - td_var is coarse grid variable to be extrapolated !> - td_level is fine grid array of level (see vgrid_get_level) [optional] !> - id_offset is array of offset between fine and coarse grid [optional] !> - id_rho is array of refinment factor [optional] !> - id_iext is number of points to be extrapolated in i-direction [optional] !> - id_jext is number of points to be extrapolated in j-direction [optional] !> - id_kext is number of points to be extrapolated in k-direction [optional] !> - id_radius is radius of the halo used to compute extrapolation [optional] !> - id_maxiter is maximum number of iteration [optional] !> !> to add extraband to the variable (to be extrapolated):
!> @code !> CALL extrap_add_extrabands(td_var, [id_isize,] [id_jsize] ) !> @endcode !> - td_var is variable structure !> - id_isize : i-direction size of extra bands [optional] !> - id_jsize : j-direction size of extra bands [optional] !> !> to delete extraband of a variable:
!> @code !> CALL extrap_del_extrabands(td_var, [id_isize,] [id_jsize] ) !> @endcode !> - td_var is variable structure !> - id_isize : i-direction size of extra bands [optional] !> - id_jsize : j-direction size of extra bands [optional] !> !> to compute first derivative of 1D array:
!> @code !> dl_value(:)=extrap_deriv_1D( dd_value(:), dd_fill, [ld_discont] ) !> @endcode !> - dd_value is 1D array of variable !> - dd_fill is FillValue of variable !> - ld_discont is logical to take into account longitudinal East-West discontinuity [optional] !> !> to compute first derivative of 2D array:
!> @code !> dl_value(:,:)=extrap_deriv_2D( dd_value(:,:), dd_fill, cd_dim, [ld_discont] ) !> @endcode !> - dd_value is 2D array of variable !> - dd_fill is FillValue of variable !> - cd_dim is character to compute derivative on first (I) or second (J) dimension !> - ld_discont is logical to take into account longitudinal East-West discontinuity [optional] !> !> to compute first derivative of 3D array:
!> @code !> dl_value(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, cd_dim, [ld_discont] ) !> @endcode !> - dd_value is 3D array of variable !> - dd_fill is FillValue of variable !> - cd_dim is character to compute derivative on first (I), second (J), or third (K) dimension !> - ld_discont is logical to take into account longitudinal East-West discontinuity [optional] !> !> @warning _FillValue must not be zero (use var_chg_FillValue()) !> !> @author !> J.Paul ! REVISION HISTORY: !> @date Nov, 2013 - Initial Version !> @date September, 2014 !> - add header !> !> @todo !> - create module for each extrapolation method !> !> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !---------------------------------------------------------------------- MODULE extrap USE netcdf ! nf90 library USE kind ! F90 kind parameter USE phycst ! physical constant USE global ! global variable USE fct ! basic useful function USE date ! date manager USE logger ! log file manager USE att ! attribute manager USE dim ! dimension manager USE var ! variable manager IMPLICIT NONE ! NOTE_avoid_public_variables_if_possible ! type and variable PRIVATE :: im_maxiter !< default maximum number of iteration PRIVATE :: im_minext !< default minumum number of point to extrapolate PRIVATE :: im_mincubic !< default minumum number of point to extrapolate for cubic interpolation ! function and subroutine PUBLIC :: extrap_detect !< detected point to be extrapolated PUBLIC :: extrap_fill_value !< extrapolate value over detected point PUBLIC :: extrap_add_extrabands !< add extraband to the variable (to be extrapolated) PUBLIC :: extrap_del_extrabands !< delete extraband of the variable PUBLIC :: extrap_deriv_1D !< compute first derivative of 1D array PUBLIC :: extrap_deriv_2D !< compute first derivative of 2D array PUBLIC :: extrap_deriv_3D !< compute first derivative of 3D array PRIVATE :: extrap__detect_wrapper ! detected point to be extrapolated wrapper PRIVATE :: extrap__detect ! detected point to be extrapolated PRIVATE :: extrap__fill_value_wrapper ! extrapolate value over detected point wrapper PRIVATE :: extrap__fill_value ! extrapolate value over detected point PRIVATE :: extrap__3D ! PRIVATE :: extrap__3D_min_error_coef ! PRIVATE :: extrap__3D_min_error_fill ! PRIVATE :: extrap__3D_dist_weight_coef ! PRIVATE :: extrap__3D_dist_weight_fill ! INTEGER(i4), PARAMETER :: im_maxiter = 10 !< default maximum number of iteration INTEGER(i4), PARAMETER :: im_minext = 2 !< default minumum number of point to extrapolate INTEGER(i4), PARAMETER :: im_mincubic= 4 !< default minumum number of point to extrapolate for cubic interpolation INTERFACE extrap_detect MODULE PROCEDURE extrap__detect_wrapper !< detected point to be extrapolated END INTERFACE extrap_detect INTERFACE extrap_fill_value MODULE PROCEDURE extrap__fill_value_wrapper !< detected point to be interpolated END INTERFACE extrap_fill_value CONTAINS !------------------------------------------------------------------- !> @brief !> This function detected point to be extrapolated, given variable structure. !> !> @details !> optionaly, you could sepcify fine grid level, refinment factor (default 1), !> offset between fine and coarse grid (default compute from refinment factor !> as offset=(rho-1)/2), number of point to be extrapolated in each direction !> (default im_minext).
!> !> First coarsening fine grid level, if need be, then select point near !> grid point already inform. !> !> @note point to be extrapolated are selected using FillValue, !> so to avoid mistake FillValue should not be zero (use var_chg_FillValue) !> !> @author J.Paul !> - November, 2013- Initial Version ! !> @param[in] td_var0 coarse grid variable to extrapolate !> @param[in] td_level1 fine grid array of level !> @param[in] id_offset array of offset between fine and coarse grid !> @param[in] id_rho array of refinment factor !> @param[in] id_ext array of number of points to be extrapolated !> @return array of point to be extrapolated !------------------------------------------------------------------- FUNCTION extrap__detect( td_var0, td_level1, & & id_offset, id_rho, id_ext ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(IN ) :: td_var0 TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level1 INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_ext ! function INTEGER(i4), DIMENSION(td_var0%t_dim(1)%i_len,& & td_var0%t_dim(2)%i_len,& & td_var0%t_dim(3)%i_len ) :: extrap__detect ! local variable CHARACTER(LEN=lc) :: cl_level INTEGER(i4) :: il_ind INTEGER(i4) , DIMENSION(:,:,:), ALLOCATABLE :: il_detect INTEGER(i4) , DIMENSION(:,:,:), ALLOCATABLE :: il_tmp INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_offset INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_level1 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_level1_G0 INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_extra INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_ext INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_rho INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_dim0 TYPE(TVAR) :: tl_var1 ! loop indices INTEGER(i4) :: ji0 INTEGER(i4) :: jj0 INTEGER(i4) :: jk0 INTEGER(i4) :: ji1 INTEGER(i4) :: jj1 INTEGER(i4) :: ji1m INTEGER(i4) :: jj1m INTEGER(i4) :: ji1p INTEGER(i4) :: jj1p !---------------------------------------------------------------- ! init extrap__detect(:,:,:)=0 ALLOCATE( il_dim0(3) ) il_dim0(:)=td_var0%t_dim(1:3)%i_len ! optional argument ALLOCATE( il_rho(ip_maxdim) ) il_rho(:)=1 IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:) ALLOCATE( il_offset(ip_maxdim,2) ) il_offset(:,:)=0 IF( PRESENT(id_offset) )THEN il_offset(1:SIZE(id_offset(:,:),DIM=1),& & 1:SIZE(id_offset(:,:),DIM=2) )= id_offset(:,:) ELSE il_offset(jp_I,:)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5) il_offset(jp_J,:)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5) ENDIF ALLOCATE( il_ext(ip_maxdim) ) il_ext(:)=im_minext IF( PRESENT(id_ext) ) il_ext(1:SIZE(id_ext(:)))=id_ext(:) ALLOCATE( il_detect(il_dim0(1),& & il_dim0(2),& & il_dim0(3)) ) il_detect(:,:,:)=0 ! select point already inform DO jk0=1,td_var0%t_dim(3)%i_len DO jj0=1,td_var0%t_dim(2)%i_len DO ji0=1,td_var0%t_dim(1)%i_len IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill ) il_detect(ji0,jj0,jk0)=1 ENDDO ENDDO ENDDO IF( PRESENT(td_level1) )THEN SELECT CASE(TRIM(td_var0%c_point)) CASE DEFAULT !'T' cl_level='tlevel' CASE('U') cl_level='ulevel' CASE('V') cl_level='vlevel' CASE('F') cl_level='flevel' END SELECT il_ind=var_get_index(td_level1(:),TRIM(cl_level)) IF( il_ind == 0 )THEN CALL logger_error("EXTRAP DETECT: can not compute point to be "//& & "extrapolated for variable "//TRIM(td_var0%c_name)//& & ". can not find "//& & "level for variable point "//TRIM(TRIM(td_var0%c_point))) ELSE tl_var1=var_copy(td_level1(il_ind)) ALLOCATE( il_level1_G0( il_dim0(1), il_dim0(2)) ) IF( ALL(tl_var1%t_dim(1:2)%i_len == il_dim0(1:2)) )THEN ! variable to be extrapolated use same resolution than level il_level1_G0(:,:)=INT(tl_var1%d_value(:,:,1,1),i4) ELSE ! variable to be extrapolated do not use same resolution than level ALLOCATE( il_level1(tl_var1%t_dim(1)%i_len, & & tl_var1%t_dim(2)%i_len) ) ! match fine grid vertical level with coarse grid il_level1(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)/il_rho(jp_K) ALLOCATE( il_extra(ip_maxdim,2) ) ! coarsening fine grid level il_extra(jp_I,1)=CEILING(REAL(il_rho(jp_I)-1,dp)*0.5_dp) il_extra(jp_I,2)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5_dp) il_extra(jp_J,1)=CEILING(REAL(il_rho(jp_J)-1,dp)*0.5_dp) il_extra(jp_J,2)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5_dp) DO jj0=1,td_var0%t_dim(2)%i_len jj1=(jj0-1)*il_rho(jp_J)+1-il_offset(jp_J,1) jj1m=MAX( jj1-il_extra(jp_J,1), 1 ) jj1p=MIN( jj1+il_extra(jp_J,2), & & tl_var1%t_dim(2)%i_len-il_offset(jp_J,2) ) DO ji0=1,td_var0%t_dim(1)%i_len ji1=(ji0-1)*il_rho(jp_I)+1-id_offset(jp_I,1) ji1m=MAX( ji1-il_extra(jp_I,1), 1 ) ji1p=MIN( ji1+il_extra(jp_I,2), & & tl_var1%t_dim(1)%i_len-id_offset(jp_I,2) ) il_level1_G0(ji0,jj0)=MAXVAL(il_level1(ji1m:ji1p,jj1m:jj1p)) ENDDO ENDDO ! clean DEALLOCATE( il_extra ) DEALLOCATE( il_level1 ) ENDIF ! look for sea point DO jk0=1,td_var0%t_dim(3)%i_len WHERE( il_level1_G0(:,:) >= jk0) il_detect(:,:,jk0)=1 END WHERE ENDDO ! clean DEALLOCATE( il_level1_G0 ) CALL var_clean(tl_var1) ENDIF ENDIF ! clean DEALLOCATE( il_offset ) ALLOCATE( il_tmp(il_dim0(1),& & il_dim0(2),& & il_dim0(3)) ) il_tmp(:,:,:)=il_detect(:,:,:) ! select extra point depending on interpolation method ! compute point near grid point already inform DO jk0=1,il_dim0(3) DO jj0=1,il_dim0(2) DO ji0=1,il_dim0(1) IF( il_tmp(ji0,jj0,jk0) == 1 )THEN il_detect( & & MAX(1,ji0-il_ext(jp_I)):MIN(ji0+il_ext(jp_I),il_dim0(1)),& & MAX(1,jj0-il_ext(jp_J)):MIN(jj0+il_ext(jp_J),il_dim0(2)),& & MAX(1,jk0-il_ext(jp_K)):MIN(jk0+il_ext(jp_K),il_dim0(3)) & & ) = 1 ENDIF ENDDO ENDDO ENDDO ! clean DEALLOCATE( il_tmp ) ! do not compute grid point already inform DO jk0=1,td_var0%t_dim(3)%i_len DO jj0=1,td_var0%t_dim(2)%i_len DO ji0=1,td_var0%t_dim(1)%i_len IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill ) il_detect(ji0,jj0,jk0)=0 ENDDO ENDDO ENDDO ! save result extrap__detect(:,:,:)=il_detect(:,:,:) ! clean DEALLOCATE( il_dim0 ) DEALLOCATE( il_ext ) DEALLOCATE( il_detect ) DEALLOCATE( il_rho ) END FUNCTION extrap__detect !------------------------------------------------------------------- !> @brief !> This function sort variable to be extrapolated, depending on number of !> dimentsion, then detected point to be extrapolated. !> !> @author J.Paul !> - November, 2013- Initial Version !> !> @param[in] td_var coarse grid variable to extrapolate !> @param[in] td_level fine grid array of level !> @param[in] id_offset array of offset between fine and coarse grid !> @param[in] id_rho array of refinment factor !> @param[in] id_ext array of number of points to be extrapolated !> @return 3D array of point to be extrapolated !------------------------------------------------------------------- FUNCTION extrap__detect_wrapper( td_var, td_level, & & id_offset, id_rho, id_ext ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(IN ) :: td_var TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_ext ! function INTEGER(i4), DIMENSION(td_var%t_dim(1)%i_len,& & td_var%t_dim(2)%i_len,& & td_var%t_dim(3)%i_len ) :: extrap__detect_wrapper ! local variable ! loop indices !---------------------------------------------------------------- ! init extrap__detect_wrapper(:,:,:)=0 IF( .NOT. ANY(td_var%t_dim(1:3)%l_use) )THEN ! no dimension I-J-K used CALL logger_debug(" EXTRAP DETECT: nothing done for variable"//& & TRIM(td_var%c_name) ) ELSE IF( ALL(td_var%t_dim(1:3)%l_use) )THEN ! detect point to be extrapolated on I-J-K CALL logger_debug(" EXTRAP DETECT: detect point "//& & " for variable "//TRIM(td_var%c_name) ) extrap__detect_wrapper(:,:,:)=extrap__detect( td_var, td_level, & & id_offset, & & id_rho, & & id_ext ) ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN ! detect point to be extrapolated on I-J CALL logger_debug(" EXTRAP DETECT: detect horizontal point "//& & " for variable "//TRIM(td_var%c_name) ) extrap__detect_wrapper(:,:,1:1)=extrap__detect( td_var , td_level,& & id_offset, & & id_rho, & & id_ext ) ELSE IF( td_var%t_dim(3)%l_use )THEN ! detect point to be extrapolated on K CALL logger_debug(" EXTRAP DETECT: detect vertical point "//& & " for variable "//TRIM(td_var%c_name) ) extrap__detect_wrapper(1:1,1:1,:)=extrap__detect( td_var , td_level, & & id_offset, & & id_rho, & & id_ext ) ENDIF CALL logger_debug(" EXTRAP DETECT: "//& & TRIM(fct_str(SUM(extrap__detect_wrapper(:,:,:))))//& & " points to be extrapolated" ) END FUNCTION extrap__detect_wrapper !------------------------------------------------------------------- !> @brief !> This subroutine select method to be used for extrapolation. !> If need be, increase number of points to be extrapolated. !> Finally launch extrap__fill_value. !> !> @details !> optionaly, you could specify :
!> - refinment factor (default 1) !> - offset between fine and coarse grid (default compute from refinment factor !> as offset=(rho-1)/2) !> - number of point to be extrapolated in each direction (default im_minext) !> - radius of the halo used to compute extrapolation !> - maximum number of iteration !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[inout] td_var variable structure !> @param[in] td_level fine grid array of level !> @param[in] id_offset array of offset between fine and coarse grid !> @param[in] id_rho array of refinment factor !> @param[in] id_iext number of points to be extrapolated in i-direction !> @param[in] id_jext number of points to be extrapolated in j-direction !> @param[in] id_kext number of points to be extrapolated in k-direction !> @param[in] id_radius radius of the halo used to compute extrapolation !> @param[in] id_maxiter maximum number of iteration !------------------------------------------------------------------- SUBROUTINE extrap__fill_value_wrapper( td_var, td_level, & & id_offset, & & id_rho, & & id_iext, id_jext, id_kext, & & id_radius, id_maxiter ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(INOUT) :: td_var TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext INTEGER(i4), INTENT(IN ), OPTIONAL :: id_kext INTEGER(i4), INTENT(IN ), OPTIONAL :: id_radius INTEGER(i4), INTENT(IN ), OPTIONAL :: id_maxiter ! local variable INTEGER(i4) :: il_iext INTEGER(i4) :: il_jext INTEGER(i4) :: il_kext INTEGER(i4) :: il_radius INTEGER(i4) :: il_maxiter CHARACTER(LEN=lc) :: cl_method ! loop indices !---------------------------------------------------------------- IF( .NOT. ASSOCIATED(td_var%d_value) )THEN CALL logger_error("EXTRAP FILL VALUE: no value "//& & "associted to variable "//TRIM(td_var%c_name) ) ELSE SELECT CASE(TRIM(td_var%c_extrap(1))) CASE('min_error') cl_method='min_error' CASE DEFAULT cl_method='dist_weight' !update variable structure td_var%c_extrap(1)='dist_weight' END SELECT il_iext=im_minext IF( PRESENT(id_iext) ) il_iext=id_iext il_jext=im_minext IF( PRESENT(id_jext) ) il_jext=id_jext il_kext=0 IF( PRESENT(id_kext) ) il_kext=id_kext IF( TRIM(td_var%c_interp(1)) == 'cubic')THEN IF( il_iext > 0 .AND. il_iext < im_mincubic ) il_iext=im_mincubic IF( il_jext > 0 .AND. il_jext < im_mincubic ) il_jext=im_mincubic ENDIF IF( il_iext < 0 )THEN CALL logger_error("EXTRAP FILL VALUE: invalid "//& & " number of points to be extrapolated in i-direction "//& & "("//TRIM(fct_str(il_iext))//")") ENDIF IF( il_jext < 0 )THEN CALL logger_error("EXTRAP FILL VALUE: invalid "//& & " number of points to be extrapolated in j-direction "//& & "("//TRIM(fct_str(il_jext))//")") ENDIF IF( il_kext < 0 )THEN CALL logger_error("EXTRAP FILL VALUE: invalid "//& & " number of points to be extrapolated in k-direction "//& & "("//TRIM(fct_str(il_kext))//")") ENDIF IF( (il_iext /= 0 .AND. td_var%t_dim(1)%l_use) .OR. & & (il_jext /= 0 .AND. td_var%t_dim(2)%l_use) .OR. & & (il_kext /= 0 .AND. td_var%t_dim(3)%l_use) )THEN ! number of point use to compute box il_radius=1 IF( PRESENT(id_radius) ) il_radius=id_radius IF( il_radius < 0 )THEN CALL logger_error("EXTRAP FILL VALUE: invalid "//& & " radius of the box used to compute extrapolation "//& & "("//TRIM(fct_str(il_radius))//")") ENDIF ! maximum number of iteration il_maxiter=im_maxiter IF( PRESENT(id_maxiter) ) il_maxiter=id_maxiter IF( il_maxiter < 0 )THEN CALL logger_error("EXTRAP FILL VALUE: invalid "//& & " maximum nuber of iteration "//& & "("//TRIM(fct_str(il_maxiter))//")") ENDIF CALL logger_info("EXTRAP FILL: extrapolate "//TRIM(td_var%c_name)//& & " using "//TRIM(cl_method)//" method." ) CALL extrap__fill_value( td_var, cl_method, & & il_iext, il_jext, il_kext, & & il_radius, il_maxiter, & & td_level, & & id_offset, id_rho ) ENDIF ENDIF END SUBROUTINE extrap__fill_value_wrapper !------------------------------------------------------------------- !> @brief !> This subroutine compute point to be extrapolated, then extrapolate point. !> !> @details !> optionaly, you could specify :
!> - refinment factor (default 1) !> - offset between fine and coarse grid (default compute from refinment factor !> as offset=(rho-1)/2) !> !> @author J.Paul !> - November, 2013- Initial Version ! !> @param[inout] td_var variable structure !> @param[in] cd_method extrapolation method !> @param[in] id_iext number of points to be extrapolated in i-direction !> @param[in] id_jext number of points to be extrapolated in j-direction !> @param[in] id_kext number of points to be extrapolated in k-direction !> @param[in] id_radius radius of the halo used to compute extrapolation !> @param[in] id_maxiter maximum number of iteration !> @param[in] td_level fine grid array of level !> @param[in] id_offset array of offset between fine and coarse grid !> @param[in] id_rho array of refinment factor !------------------------------------------------------------------- SUBROUTINE extrap__fill_value( td_var, cd_method, & & id_iext, id_jext, id_kext, & & id_radius, id_maxiter, & & td_level, & & id_offset, & & id_rho ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(INOUT) :: td_var CHARACTER(LEN=*), INTENT(IN ) :: cd_method INTEGER(i4) , INTENT(IN ) :: id_iext INTEGER(i4) , INTENT(IN ) :: id_jext INTEGER(i4) , INTENT(IN ) :: id_kext INTEGER(i4) , INTENT(IN ) :: id_radius INTEGER(i4) , INTENT(IN ) :: id_maxiter TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level INTEGER(i4) , DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset INTEGER(i4) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho ! local variable CHARACTER(LEN=lc) :: cl_extrap INTEGER(i4), DIMENSION(:,:,:) , ALLOCATABLE :: il_detect TYPE(TATT) :: tl_att ! loop indices !---------------------------------------------------------------- !1- detect point to be extrapolated ALLOCATE( il_detect( td_var%t_dim(1)%i_len, & & td_var%t_dim(2)%i_len, & & td_var%t_dim(3)%i_len) ) il_detect(:,:,:) = extrap_detect( td_var, td_level, & & id_offset, & & id_rho, & & id_ext=(/id_iext, id_jext, id_kext/) ) !2- add attribute to variable cl_extrap=fct_concat(td_var%c_extrap(:)) tl_att=att_init('extrapolation',cl_extrap) CALL var_move_att(td_var, tl_att) CALL att_clean(tl_att) CALL logger_info(" EXTRAP FILL: "//& & TRIM(fct_str(SUM(il_detect(:,:,:))))//& & " point(s) to extrapolate " ) !3- extrapolate CALL extrap__3D(td_var%d_value(:,:,:,:), td_var%d_fill, & & il_detect(:,:,:), & & cd_method, id_radius, id_maxiter ) DEALLOCATE(il_detect) END SUBROUTINE extrap__fill_value !------------------------------------------------------------------- !> @brief !> This subroutine compute point to be extrapolated in 3D array. !> !> @details !> in case of 'min_error' method:
!> - compute derivative in i-, j- and k- direction !> - compute minimum error coefficient (distance to center of halo) !> - compute extrapolatd values by calculated minimum error using taylor expansion !> in case of 'dist_weight' method:
!> - compute distance weight coefficient (inverse of distance to center of halo) !> - compute extrapolatd values using Inverse Distance Weighting !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[inout] dd_value 3D array of variable to be extrapolated !> @param[in] dd_fill FillValue of variable !> @param[inout] id_detect array of point to extrapolate !> @param[in] cd_method extrapolation method !> @param[in] id_radius radius of the halo used to compute extrapolation !------------------------------------------------------------------- SUBROUTINE extrap__3D( dd_value, dd_fill, id_detect,& & cd_method, id_radius, id_maxiter ) 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_radius INTEGER(i4), INTENT(IN ) :: id_maxiter ! local variable INTEGER(i4) :: il_imin INTEGER(i4) :: il_imax INTEGER(i4) :: il_jmin INTEGER(i4) :: il_jmax INTEGER(i4) :: il_kmin INTEGER(i4) :: il_kmax INTEGER(i4) :: il_iter INTEGER(i4) :: il_radius INTEGER(i4), DIMENSION(4) :: il_shape INTEGER(i4), DIMENSION(3) :: il_dim INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdy REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdz REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj INTEGER(i4) :: jk INTEGER(i4) :: jl !---------------------------------------------------------------- il_shape(:)=SHAPE(dd_value) ALLOCATE( il_detect( il_shape(1), il_shape(2), il_shape(3)) ) SELECT CASE(TRIM(cd_method)) CASE('min_error') DO jl=1,il_shape(4) ! intitialise table of poitn to be extrapolated il_detect(:,:,:)=id_detect(:,:,:) il_iter=1 DO WHILE( ANY(il_detect(:,:,:)==1) ) ! change extend value to minimize number of iteration il_radius=id_radius+(il_iter/id_maxiter) ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) ) ALLOCATE( dl_dfdy(il_shape(1), il_shape(2), il_shape(3)) ) ALLOCATE( dl_dfdz(il_shape(1), il_shape(2), il_shape(3)) ) ! compute derivative in i-direction dl_dfdx(:,:,:)=dd_fill IF( il_shape(1) > 1 )THEN dl_dfdx(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'I' ) ENDIF ! compute derivative in j-direction dl_dfdy(:,:,:)=dd_fill IF( il_shape(2) > 1 )THEN dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'J' ) ENDIF ! compute derivative in k-direction dl_dfdz(:,:,:)=dd_fill IF( il_shape(3) > 1 )THEN dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:,jl), dd_fill, 'K' ) ENDIF il_dim(1)=2*il_radius+1 IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1 il_dim(2)=2*il_radius+1 IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1 il_dim(3)=2*il_radius+1 IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1 ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) dl_coef(:,:,:)=extrap__3D_min_error_coef(dd_value( 1:il_dim(1), & & 1:il_dim(2), & & 1:il_dim(3), & & jl )) DO jk=1,il_shape(3) IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE DO jj=1,il_shape(2) IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE DO ji=1,il_shape(1) IF( il_detect(ji,jj,jk) == 1 )THEN il_imin=MAX(ji-il_radius,1) il_imax=MIN(ji+il_radius,il_shape(1)) IF( il_dim(1) == 1 )THEN il_imin=ji il_imax=ji ENDIF il_jmin=MAX(jj-il_radius,1) il_jmax=MIN(jj+il_radius,il_shape(2)) IF( il_dim(2) == 1 )THEN il_jmin=jj il_jmax=jj ENDIF il_kmin=MAX(jk-il_radius,1) il_kmax=MIN(jk+il_radius,il_shape(3)) IF( il_dim(3) == 1 )THEN il_kmin=jk il_kmax=jk ENDIF dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( & & dd_value( il_imin:il_imax, & & il_jmin:il_jmax, & & il_kmin:il_kmax,jl ), dd_fill, il_radius, & & dl_dfdx( il_imin:il_imax, & & il_jmin:il_jmax, & & il_kmin:il_kmax ), & & dl_dfdy( il_imin:il_imax, & & il_jmin:il_jmax, & & il_kmin:il_kmax ), & & dl_dfdz( il_imin:il_imax, & & il_jmin:il_jmax, & & il_kmin:il_kmax ), & & dl_coef(:,:,:) ) IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN il_detect(ji,jj,jk)= 0 ENDIF ENDIF ENDDO ENDDO ENDDO DEALLOCATE( dl_dfdx ) DEALLOCATE( dl_dfdy ) DEALLOCATE( dl_dfdz ) DEALLOCATE( dl_coef ) il_iter=il_iter+1 ENDDO ENDDO CASE DEFAULT ! 'dist_weight' DO jl=1,il_shape(4) ! intitialise table of poitn to be extrapolated il_detect(:,:,:)=id_detect(:,:,:) il_iter=1 DO WHILE( ANY(il_detect(:,:,:)==1) ) ! change extend value to minimize number of iteration il_radius=id_radius+(il_iter/id_maxiter) il_dim(1)=2*il_radius+1 IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1 il_dim(2)=2*il_radius+1 IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1 il_dim(3)=2*il_radius+1 IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1 ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) ) dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), & & 1:il_dim(2), & & 1:il_dim(3), & & jl ) ) DO jk=1,il_shape(3) IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE DO jj=1,il_shape(2) IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE DO ji=1,il_shape(1) IF( il_detect(ji,jj,jk) == 1 )THEN il_imin=MAX(ji-il_radius,1) il_imax=MIN(ji+il_radius,il_shape(1)) IF( il_dim(1) == 1 )THEN il_imin=ji il_imax=ji ENDIF il_jmin=MAX(jj-il_radius,1) il_jmax=MIN(jj+il_radius,il_shape(2)) IF( il_dim(2) == 1 )THEN il_jmin=jj il_jmax=jj ENDIF il_kmin=MAX(jk-il_radius,1) il_kmax=MIN(jk+il_radius,il_shape(3)) IF( il_dim(3) == 1 )THEN il_kmin=jk il_kmax=jk ENDIF dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( & & dd_value( il_imin:il_imax, & & il_jmin:il_jmax, & & il_kmin:il_kmax, & & jl), dd_fill, il_radius, & & dl_coef(:,:,:) ) IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN il_detect(ji,jj,jk)= 0 ENDIF ENDIF ENDDO ENDDO ENDDO DEALLOCATE( dl_coef ) il_iter=il_iter+1 ENDDO ENDDO END SELECT DEALLOCATE( il_detect ) END SUBROUTINE extrap__3D !------------------------------------------------------------------- !> @brief !> This function compute derivative of 1D array. !> !> @details !> optionaly you could specify to take into account east west discontinuity !> (-180° 180° or 0° 360° for longitude variable) !> !> @author J.Paul !> - November, 2013- Initial Version ! !> @param[in] dd_value 1D array of variable to be extrapolated !> @param[in] dd_fill FillValue of variable !> @param[in] ld_discont logical to take into account east west discontinuity !------------------------------------------------------------------- PURE FUNCTION extrap_deriv_1D( dd_value, dd_fill, ld_discont ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:), INTENT(IN) :: dd_value REAL(dp) , INTENT(IN) :: dd_fill LOGICAL , INTENT(IN), OPTIONAL :: ld_discont ! function REAL(dp), DIMENSION(SIZE(dd_value,DIM=1) ) :: extrap_deriv_1D ! local variable INTEGER(i4) :: il_imin INTEGER(i4) :: il_imax INTEGER(i4), DIMENSION(1) :: il_shape REAL(dp) :: dl_min REAL(dp) :: dl_max REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_value LOGICAL :: ll_discont ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: i1 INTEGER(i4) :: i2 !---------------------------------------------------------------- ! init extrap_deriv_1D(:)=dd_fill ll_discont=.FALSE. IF( PRESENT(ld_discont) ) ll_discont=ld_discont il_shape(:)=SHAPE(dd_value(:)) ALLOCATE( dl_value(3)) ! compute derivative in i-direction DO ji=1,il_shape(1) il_imin=MAX(ji-1,1) il_imax=MIN(ji+1,il_shape(1)) IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN i1=1 ; i2=3 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN i1=1 ; i2=2 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN i1=2 ; i2=3 ENDIF dl_value(i1:i2)=dd_value(il_imin:il_imax) IF( il_imin == 1 )THEN dl_value(:)=EOSHIFT( dl_value(:), & & DIM=1, & & SHIFT=-1, & & BOUNDARY=dl_value(1) ) ENDIF IF( il_imax == il_shape(1) )THEN dl_value(:)=EOSHIFT( dl_value(:), & & DIM=1, & & SHIFT=1, & & BOUNDARY=dl_value(3)) ENDIF IF( ll_discont )THEN dl_min=MINVAL( dl_value(:), dl_value(:)/=dd_fill ) dl_max=MAXVAL( dl_value(:), dl_value(:)/=dd_fill ) IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN WHERE( dl_value(:) < 0._dp ) dl_value(:) = dl_value(:)+360._dp END WHERE ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN WHERE( dl_value(:) > 180._dp ) dl_value(:) = dl_value(:)-180._dp END WHERE ENDIF ENDIF IF( dl_value( 2) /= dd_fill .AND. & ! ji & dl_value( 3) /= dd_fill .AND. & ! ji+1 & dl_value( 1) /= dd_fill )THEN ! ji-1 extrap_deriv_1D(ji)=& & ( dl_value(3) - dl_value(1) ) / & & REAL( il_imax-il_imin ,dp) ENDIF ENDDO DEALLOCATE( dl_value ) END FUNCTION extrap_deriv_1D !------------------------------------------------------------------- !> @brief !> This function compute derivative of 2D array. !> you have to specify in which direction derivative have to be computed: !> first (I) or second (J) dimension. !> !> @details !> optionaly you could specify to take into account east west discontinuity !> (-180° 180° or 0° 360° for longitude variable) !> !> @author J.Paul !> - November, 2013- Initial Version ! !> @param[in] dd_value 2D array of variable to be extrapolated !> @param[in] dd_fill FillValue of variable !> @param[in] cd_dim compute derivative on first (I) or second (J) dimension !> @param[in] ld_discont logical to take into account east west discontinuity !------------------------------------------------------------------- FUNCTION extrap_deriv_2D( dd_value, dd_fill, cd_dim, ld_discont ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_value REAL(dp) , INTENT(IN) :: dd_fill CHARACTER(LEN=*) , INTENT(IN) :: cd_dim LOGICAL , INTENT(IN), OPTIONAL :: ld_discont ! function REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), & & SIZE(dd_value,DIM=2) ) :: extrap_deriv_2D ! local variable INTEGER(i4) :: il_imin INTEGER(i4) :: il_imax INTEGER(i4) :: il_jmin INTEGER(i4) :: il_jmax INTEGER(i4), DIMENSION(2) :: il_shape REAL(dp) :: dl_min REAL(dp) :: dl_max REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_value LOGICAL :: ll_discont ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj INTEGER(i4) :: i1 INTEGER(i4) :: i2 INTEGER(i4) :: j1 INTEGER(i4) :: j2 !---------------------------------------------------------------- ! init extrap_deriv_2D(:,:)=dd_fill ll_discont=.FALSE. IF( PRESENT(ld_discont) ) ll_discont=ld_discont il_shape(:)=SHAPE(dd_value(:,:)) SELECT CASE(TRIM(fct_upper(cd_dim))) CASE('I') ALLOCATE( dl_value(3,il_shape(2)) ) ! compute derivative in i-direction DO ji=1,il_shape(1) ! init dl_value(:,:)=dd_fill il_imin=MAX(ji-1,1) il_imax=MIN(ji+1,il_shape(1)) IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN i1=1 ; i2=3 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN i1=1 ; i2=2 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN i1=2 ; i2=3 ENDIF dl_value(i1:i2,:)=dd_value(il_imin:il_imax,:) IF( il_imin == 1 )THEN dl_value(:,:)=EOSHIFT( dl_value(:,:), & & DIM=1, & & SHIFT=-1, & & BOUNDARY=dl_value(1,:) ) ENDIF IF( il_imax == il_shape(1) )THEN dl_value(:,:)=EOSHIFT( dl_value(:,:), & & DIM=1, & & SHIFT=1, & & BOUNDARY=dl_value(3,:)) ENDIF IF( ll_discont )THEN dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN WHERE( dl_value(:,:) < 0_dp ) dl_value(:,:) = dl_value(:,:)+360._dp END WHERE ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN WHERE( dl_value(:,:) > 180 ) dl_value(:,:) = dl_value(:,:)-180._dp END WHERE ENDIF ENDIF WHERE( dl_value(2,:) /= dd_fill .AND. & ! ji & dl_value(3,:) /= dd_fill .AND. & ! ji+1 & dl_value(1,:) /= dd_fill ) ! ji-1 extrap_deriv_2D(ji,:)=& & ( dl_value(3,:) - dl_value(1,:) ) / & & REAL( il_imax-il_imin,dp) END WHERE ENDDO CASE('J') ALLOCATE( dl_value(il_shape(1),3) ) ! compute derivative in j-direction DO jj=1,il_shape(2) il_jmin=MAX(jj-1,1) il_jmax=MIN(jj+1,il_shape(2)) IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN j1=1 ; j2=3 ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN j1=1 ; j2=2 ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN j1=2 ; j2=3 ENDIF dl_value(:,j1:j2)=dd_value(:,il_jmin:il_jmax) IF( il_jmin == 1 )THEN dl_value(:,:)=EOSHIFT( dl_value(:,:), & & DIM=2, & & SHIFT=-1, & & BOUNDARY=dl_value(:,1)) ENDIF IF( il_jmax == il_shape(2) )THEN dl_value(:,:)=EOSHIFT( dl_value(:,:), & & DIM=2, & & SHIFT=1, & & BOUNDARY=dl_value(:,3)) ENDIF IF( ll_discont )THEN dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill ) IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN WHERE( dl_value(:,:) < 0_dp ) dl_value(:,:) = dl_value(:,:)+360._dp END WHERE ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN WHERE( dl_value(:,:) > 180 ) dl_value(:,:) = dl_value(:,:)-180._dp END WHERE ENDIF ENDIF WHERE( dl_value(:, 2) /= dd_fill .AND. & ! jj & dl_value(:, 3) /= dd_fill .AND. & ! jj+1 & dl_value(:, 1) /= dd_fill ) ! jj-1 extrap_deriv_2D(:,jj)=& & ( dl_value(:,3) - dl_value(:,1) ) / & & REAL(il_jmax-il_jmin,dp) END WHERE ENDDO END SELECT DEALLOCATE( dl_value ) END FUNCTION extrap_deriv_2D !------------------------------------------------------------------- !> @brief !> This function compute derivative of 3D array. !> you have to specify in which direction derivative have to be computed: !> first (I), second (J) or third (K) dimension. !> !> @details !> optionaly you could specify to take into account east west discontinuity !> (-180° 180° or 0° 360° for longitude variable) !> !> @author J.Paul !> - November, 2013- Initial Version ! !> @param[inout] dd_value 3D array of variable to be extrapolated !> @param[in] dd_fill FillValue of variable !> @param[in] cd_dim compute derivative on first (I) second (J) or third (K) dimension !> @param[in] ld_discont logical to take into account east west discontinuity !------------------------------------------------------------------- PURE FUNCTION extrap_deriv_3D( dd_value, dd_fill, cd_dim, ld_discont ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value REAL(dp) , INTENT(IN) :: dd_fill CHARACTER(LEN=*) , INTENT(IN) :: cd_dim LOGICAL , INTENT(IN), OPTIONAL :: ld_discont ! function REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), & & SIZE(dd_value,DIM=2), & & SIZE(dd_value,DIM=3)) :: extrap_deriv_3D ! local variable INTEGER(i4) :: il_imin INTEGER(i4) :: il_imax INTEGER(i4) :: il_jmin INTEGER(i4) :: il_jmax INTEGER(i4) :: il_kmin INTEGER(i4) :: il_kmax INTEGER(i4), DIMENSION(3) :: il_shape REAL(dp) :: dl_min REAL(dp) :: dl_max REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_value LOGICAL :: ll_discont ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj INTEGER(i4) :: jk INTEGER(i4) :: i1 INTEGER(i4) :: i2 INTEGER(i4) :: j1 INTEGER(i4) :: j2 INTEGER(i4) :: k1 INTEGER(i4) :: k2 !---------------------------------------------------------------- ! init extrap_deriv_3D(:,:,:)=dd_fill ll_discont=.FALSE. IF( PRESENT(ld_discont) ) ll_discont=ld_discont il_shape(:)=SHAPE(dd_value(:,:,:)) SELECT CASE(TRIM(fct_upper(cd_dim))) CASE('I') ALLOCATE( dl_value(3,il_shape(2),il_shape(3)) ) ! compute derivative in i-direction DO ji=1,il_shape(1) il_imin=MAX(ji-1,1) il_imax=MIN(ji+1,il_shape(1)) IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN i1=1 ; i2=3 ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN i1=1 ; i2=2 ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN i1=2 ; i2=3 ENDIF dl_value(i1:i2,:,:)=dd_value(il_imin:il_imax,:,:) IF( il_imin == 1 )THEN dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & & DIM=1, & & SHIFT=-1, & & BOUNDARY=dl_value(1,:,:) ) ENDIF IF( il_imax == il_shape(1) )THEN dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & & DIM=1, & & SHIFT=1, & & BOUNDARY=dl_value(3,:,:)) ENDIF IF( ll_discont )THEN dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN WHERE( dl_value(:,:,:) < 0_dp ) dl_value(:,:,:) = dl_value(:,:,:)+360._dp END WHERE ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN WHERE( dl_value(:,:,:) > 180 ) dl_value(:,:,:) = dl_value(:,:,:)-180._dp END WHERE ENDIF ENDIF WHERE( dl_value(2,:,:) /= dd_fill .AND. & ! ji & dl_value(3,:,:) /= dd_fill .AND. & !ji+1 & dl_value(1,:,:) /= dd_fill ) !ji-1 extrap_deriv_3D(ji,:,:)= & & ( dl_value(3,:,:) - dl_value(1,:,:) ) / & & REAL( il_imax-il_imin ,dp) END WHERE ENDDO CASE('J') ALLOCATE( dl_value(il_shape(1),3,il_shape(3)) ) ! compute derivative in j-direction DO jj=1,il_shape(2) il_jmin=MAX(jj-1,1) il_jmax=MIN(jj+1,il_shape(2)) IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN j1=1 ; j2=3 ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN j1=1 ; j2=2 ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN j1=2 ; j2=3 ENDIF dl_value(:,j1:j2,:)=dd_value(:,il_jmin:il_jmax,:) IF( il_jmin == 1 )THEN dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & & DIM=2, & & SHIFT=-1, & & BOUNDARY=dl_value(:,1,:) ) ENDIF IF( il_jmax == il_shape(2) )THEN dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & & DIM=2, & & SHIFT=1, & & BOUNDARY=dl_value(:,3,:)) ENDIF IF( ll_discont )THEN dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN WHERE( dl_value(:,:,:) < 0_dp ) dl_value(:,:,:) = dl_value(:,:,:)+360._dp END WHERE ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN WHERE( dl_value(:,:,:) > 180 ) dl_value(:,:,:) = dl_value(:,:,:)-180._dp END WHERE ENDIF ENDIF WHERE( dl_value(:, 2,:) /= dd_fill .AND. & ! jj & dl_value(:, 3,:) /= dd_fill .AND. & ! jj+1 & dl_value(:, 1,:) /= dd_fill ) ! jj-1 extrap_deriv_3D(:,jj,:)=& & ( dl_value(:,3,:) - dl_value(:,1,:) ) / & & REAL( il_jmax - il_jmin ,dp) END WHERE ENDDO CASE('K') ! compute derivative in k-direction DO jk=1,il_shape(3) il_kmin=MAX(jk-1,1) il_kmax=MIN(jk+1,il_shape(3)) IF( il_kmin==jk-1 .AND. il_kmax==jk+1 )THEN k1=1 ; k2=3 ELSEIF( il_kmin==jk .AND. il_kmax==jk+1 )THEN k1=1 ; k2=2 ELSEIF( il_kmin==jk-1 .AND. il_kmax==jk )THEN k1=2 ; k2=3 ENDIF dl_value(:,:,k1:k2)=dd_value(:,:,il_kmin:il_kmax) IF( il_kmin == 1 )THEN dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & & DIM=3, & & SHIFT=-1, & & BOUNDARY=dl_value(:,:,1) ) ENDIF IF( il_kmax == il_shape(3) )THEN dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), & & DIM=3, & & SHIFT=1, & & BOUNDARY=dl_value(:,:,3)) ENDIF IF( ll_discont )THEN dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill ) IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN WHERE( dl_value(:,:,:) < 0_dp ) dl_value(:,:,:) = dl_value(:,:,:)+360._dp END WHERE ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN WHERE( dl_value(:,:,:) > 180 ) dl_value(:,:,:) = dl_value(:,:,:)-180._dp END WHERE ENDIF ENDIF WHERE( dl_value(:,:, 2) /= dd_fill .AND. & ! jk & dl_value(:,:, 3) /= dd_fill .AND. & ! jk+1 & dl_value(:,:, 1) /= dd_fill ) ! jk-1 extrap_deriv_3D(:,:,jk)=& & ( dl_value(:,:,3) - dl_value(:,:,1) ) / & & REAL( il_kmax-il_kmin,dp) END WHERE ENDDO END SELECT DEALLOCATE( dl_value ) END FUNCTION extrap_deriv_3D !------------------------------------------------------------------- !> @brief !> This function compute coefficient for min_error extrapolation. !> !> @details !> coefficients are "grid distance" to the center of the box choosed to compute !> extrapolation. !> !> @author J.Paul !> - November, 2013- Initial Version ! !> @param[in] dd_value 3D array of variable to be extrapolated !> @return 3D array of coefficient for minimum error extrapolation !------------------------------------------------------------------- PURE FUNCTION extrap__3D_min_error_coef( dd_value ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value ! function REAL(dp), DIMENSION(SIZE(dd_value(:,:,:),DIM=1), & & SIZE(dd_value(:,:,:),DIM=2), & & SIZE(dd_value(:,:,:),DIM=3) ) :: extrap__3D_min_error_coef ! local variable INTEGER(i4), DIMENSION(3) :: il_shape INTEGER(i4) :: il_imid INTEGER(i4) :: il_jmid INTEGER(i4) :: il_kmid REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dist ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj INTEGER(i4) :: jk !---------------------------------------------------------------- ! init extrap__3D_min_error_coef(:,:,:)=0 il_shape(:)=SHAPE(dd_value(:,:,:)) il_imid=INT(REAL(il_shape(1),sp)*0.5+1) il_jmid=INT(REAL(il_shape(2),sp)*0.5+1) il_kmid=INT(REAL(il_shape(3),sp)*0.5+1) ALLOCATE( dl_dist(il_shape(1),il_shape(2),il_shape(3)) ) DO jk=1,il_shape(3) DO jj=1,il_shape(2) DO ji=1,il_shape(1) ! compute distance dl_dist(ji,jj,jk) = (ji-il_imid)**2 + & & (jj-il_jmid)**2 + & & (jk-il_kmid)**2 IF( dl_dist(ji,jj,jk) /= 0 )THEN dl_dist(ji,jj,jk)=SQRT( dl_dist(ji,jj,jk) ) ENDIF ENDDO ENDDO ENDDO WHERE( dl_dist(:,:,:) /= 0 ) extrap__3D_min_error_coef(:,:,:)=dl_dist(:,:,:) END WHERE DEALLOCATE( dl_dist ) END FUNCTION extrap__3D_min_error_coef !------------------------------------------------------------------- !> @brief !> This function compute extrapolatd value by calculated minimum error using !> taylor expansion !> !> @author J.Paul !> - November, 2013- Initial Version !> !> @param[in] dd_value 3D array of variable to be extrapolated !> @param[in] dd_fill FillValue of variable !> @param[in] id_radius radius of the halo used to compute extrapolation !> @param[in] dd_dfdx derivative of function in i-direction !> @param[in] dd_dfdy derivative of function in j-direction !> @param[in] dd_dfdz derivative of function in k-direction !> @param[in] dd_coef array of coefficient for min_error extrapolation !> @return extrapolatd value !------------------------------------------------------------------- PURE FUNCTION extrap__3D_min_error_fill( dd_value, dd_fill, id_radius, & & dd_dfdx, dd_dfdy, dd_dfdz, & & dd_coef ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value REAL(dp) , INTENT(IN) :: dd_fill INTEGER(i4), INTENT(IN) :: id_radius REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdx REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdy REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdz REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_coef ! function REAL(dp) :: extrap__3d_min_error_fill ! local variable INTEGER(i4), DIMENSION(3) :: il_shape INTEGER(i4), DIMENSION(3) :: il_ind INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_mask REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_error INTEGER(i4) :: il_min ! loop indices !---------------------------------------------------------------- ! init extrap__3D_min_error_fill=dd_fill il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2)) IF( COUNT(dd_value(:,:,:) /= dd_fill) >= il_min )THEN il_shape(:)=SHAPE(dd_value(:,:,:)) ALLOCATE( il_mask( il_shape(1),il_shape(2),il_shape(3)) ) ALLOCATE( dl_error(il_shape(1),il_shape(2),il_shape(3)) ) ! compute error dl_error(:,:,:)=0. il_mask(:,:,:)=0 WHERE( dd_dfdx(:,:,:) /= dd_fill ) dl_error(:,:,:)=dd_coef(:,:,:)*dd_dfdx(:,:,:) il_mask(:,:,:)=1 END WHERE WHERE( dd_dfdy(:,:,:) /= dd_fill ) dl_error(:,:,:)=(dl_error(:,:,:)+dd_coef(:,:,:)*dd_dfdy(:,:,:)) il_mask(:,:,:)=1 END WHERE WHERE( dd_dfdz(:,:,:) /= dd_fill ) dl_error(:,:,:)=(dl_error(:,:,:)+dd_coef(:,:,:)*dd_dfdz(:,:,:)) il_mask(:,:,:)=1 END WHERE ! get minimum error indices il_ind(:)=MINLOC(dl_error(:,:,:),il_mask(:,:,:)==1) ! return value IF( ALL(il_ind(:)/=0) )THEN extrap__3D_min_error_fill=dd_value(il_ind(1),il_ind(2),il_ind(3)) ENDIF DEALLOCATE( il_mask ) DEALLOCATE( dl_error ) ENDIF END FUNCTION extrap__3D_min_error_fill !------------------------------------------------------------------- !> @brief !> This function compute coefficient for inverse distance weighted method !> !> @details !> coefficients are inverse "grid distance" to the center of the box choosed to compute !> extrapolation. !> !> @author J.Paul !> - November, 2013- Initial Version ! !> @param[in] dd_value 3D array of variable to be extrapolated !> @return 3D array of coefficient for inverse distance weighted extrapolation !------------------------------------------------------------------- PURE FUNCTION extrap__3D_dist_weight_coef( dd_value ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value ! function REAL(dp), DIMENSION(SIZE(dd_value(:,:,:),DIM=1), & & SIZE(dd_value(:,:,:),DIM=2), & & SIZE(dd_value(:,:,:),DIM=3) ) :: extrap__3D_dist_weight_coef ! local variable INTEGER(i4), DIMENSION(3) :: il_shape INTEGER(i4) :: il_imid INTEGER(i4) :: il_jmid INTEGER(i4) :: il_kmid REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dist ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj INTEGER(i4) :: jk !---------------------------------------------------------------- ! init extrap__3D_dist_weight_coef(:,:,:)=0 il_shape(:)=SHAPE(dd_value(:,:,:)) il_imid=INT(REAL(il_shape(1),sp)*0.5+1,i4) il_jmid=INT(REAL(il_shape(2),sp)*0.5+1,i4) il_kmid=INT(REAL(il_shape(3),sp)*0.5+1,i4) ALLOCATE( dl_dist(il_shape(1),il_shape(2),il_shape(3)) ) DO jk=1,il_shape(3) DO jj=1,il_shape(2) DO ji=1,il_shape(1) ! compute distance dl_dist(ji,jj,jk) = (ji-il_imid)**2 + & & (jj-il_jmid)**2 + & & (jk-il_kmid)**2 IF( dl_dist(ji,jj,jk) /= 0 )THEN dl_dist(ji,jj,jk)=SQRT( dl_dist(ji,jj,jk) ) ENDIF ENDDO ENDDO ENDDO WHERE( dl_dist(:,:,:) /= 0 ) extrap__3D_dist_weight_coef(:,:,:)=1./dl_dist(:,:,:) END WHERE DEALLOCATE( dl_dist ) END FUNCTION extrap__3D_dist_weight_coef !------------------------------------------------------------------- !> @brief !> This function compute extrapolatd value using inverse distance weighted !> method !> !> @details !> !> @author J.Paul !> - November, 2013- Initial Version ! !> @param[in] dd_value 3D array of variable to be extrapolated !> @param[in] dd_fill FillValue of variable !> @param[in] id_radius radius of the halo used to compute extrapolation !> @param[in] dd_coef 3D array of coefficient for inverse distance weighted extrapolation !> @return extrapolatd value !------------------------------------------------------------------- FUNCTION extrap__3D_dist_weight_fill( dd_value, dd_fill, id_radius, & & dd_coef ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value REAL(dp) , INTENT(IN) :: dd_fill INTEGER(i4), INTENT(IN) :: id_radius REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_coef ! function REAL(dp) :: extrap__3D_dist_weight_fill ! local variable INTEGER(i4), DIMENSION(3) :: il_shape REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_value REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef INTEGER(i4) :: il_min ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: jj INTEGER(i4) :: jk !---------------------------------------------------------------- ! init extrap__3D_dist_weight_fill=dd_fill il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2)) IF( COUNT(dd_value(:,:,:)/= dd_fill) >= il_min )THEN il_shape(:)=SHAPE(dd_value(:,:,:)) ALLOCATE( dl_value( il_shape(1),il_shape(2),il_shape(3)) ) ALLOCATE( dl_coef( il_shape(1),il_shape(2),il_shape(3)) ) dl_value(:,:,:)=0 dl_coef(:,:,:)=0 DO jk=1,il_shape(3) DO jj=1,il_shape(2) DO ji=1,il_shape(1) IF( dd_value(ji,jj,jk) /= dd_fill )THEN ! compute factor dl_value(ji,jj,jk)=dd_coef(ji,jj,jk)*dd_value(ji,jj,jk) dl_coef(ji,jj,jk)=dd_coef(ji,jj,jk) ENDIF ENDDO ENDDO ENDDO ! return value IF( SUM( dl_coef(:,:,:) ) /= 0 )THEN extrap__3D_dist_weight_fill = & & SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) ) ENDIF DEALLOCATE( dl_value ) DEALLOCATE( dl_coef ) ENDIF END FUNCTION extrap__3D_dist_weight_fill !------------------------------------------------------------------- !> @brief !> This subroutine add to the variable (to be extrapolated) an !> extraband of N points at north,south,east and west boundaries. !> !> @details !> optionaly you could specify size of extra bands in i- and j-direction !> !> @author J.Paul !> - November, 2013-Initial version ! !> @param[inout] td_var variable !> @param[in] id_isize i-direction size of extra bands (default=im_minext) !> @param[in] id_jsize j-direction size of extra bands (default=im_minext) !> @todo !> - invalid special case for grid with north fold !------------------------------------------------------------------- SUBROUTINE extrap_add_extrabands(td_var, id_isize, id_jsize ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(INOUT) :: td_var INTEGER(i4), INTENT(IN ), OPTIONAL :: id_isize INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jsize ! local variable REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value INTEGER(i4) :: il_isize INTEGER(i4) :: il_jsize INTEGER(i4) :: il_tmp ! loop indices INTEGER(i4) :: ji INTEGER(i4) :: ij !---------------------------------------------------------------- il_isize=im_minext IF(PRESENT(id_isize)) il_isize=id_isize IF( il_isize < im_minext .AND. & & TRIM(td_var%c_interp(1)) == 'cubic' )THEN CALL logger_warn("EXTRAP ADD EXTRABANDS: size of extrabands "//& & "should be at least "//TRIM(fct_str(im_minext))//" for "//& & " cubic interpolation ") ENDIF il_jsize=im_minext IF(PRESENT(id_jsize)) il_jsize=id_jsize IF( il_jsize < im_minext .AND. & & TRIM(td_var%c_interp(1)) == 'cubic' )THEN CALL logger_warn("EXTRAP ADD EXTRABANDS: size of extrabands "//& & "should be at least "//TRIM(fct_str(im_minext))//" for "//& & " cubic interpolation ") ENDIF IF( .NOT. td_var%t_dim(1)%l_use ) il_isize=0 IF( .NOT. td_var%t_dim(2)%l_use ) il_jsize=0 CALL logger_trace( "EXTRAP ADD EXTRABANDS: dimension change "//& & "in variable "//TRIM(td_var%c_name) ) ! add extrabands in variable 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(:,:,:,:) td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len + 2*il_isize td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len + 2*il_jsize DEALLOCATE(td_var%d_value) ALLOCATE( td_var%d_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 ) ) ! intialise td_var%d_value(:,:,:,:)=td_var%d_fill ! fill center td_var%d_value( 1+il_isize:td_var%t_dim(1)%i_len-il_isize, & & 1+il_jsize:td_var%t_dim(2)%i_len-il_jsize, & & :,:) = dl_value(:,:,:,:) ! special case for overlap IF( td_var%i_ew >= 0 .AND. il_isize /= 0 )THEN DO ji=1,il_isize ! from east to west il_tmp=td_var%t_dim(1)%i_len-td_var%i_ew+ji-2*il_isize td_var%d_value(ji,:,:,:) = td_var%d_value(il_tmp,:,:,:) ! from west to east ij=td_var%t_dim(1)%i_len-ji+1 il_tmp=td_var%i_ew-ji+2*il_isize+1 td_var%d_value(ij,:,:,:) = td_var%d_value(il_tmp,:,:,:) ENDDO ENDIF DEALLOCATE( dl_value ) END SUBROUTINE extrap_add_extrabands !------------------------------------------------------------------- !> @brief !> This subroutine remove of the variable an extraband !> of N points at north,south,east and west boundaries. !> !> @details !> optionaly you could specify size of extra bands in i- and j-direction !> !> @author J.Paul !> - November, 2013-Initial version !> !> @param[inout] td_var variable !> @param[in] id_isize i-direction size of extra bands (default=im_minext) !> @param[in] id_jsize j-direction size of extra bands (default=im_minext) !------------------------------------------------------------------- SUBROUTINE extrap_del_extrabands(td_var, id_isize, id_jsize ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(INOUT) :: td_var INTEGER(i4), INTENT(IN ), OPTIONAL :: id_isize INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jsize ! local variable REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value INTEGER(i4) :: il_isize INTEGER(i4) :: il_jsize INTEGER(i4) :: il_imin INTEGER(i4) :: il_imax INTEGER(i4) :: il_jmin INTEGER(i4) :: il_jmax ! loop indices !---------------------------------------------------------------- il_isize=im_minext IF(PRESENT(id_isize)) il_isize=id_isize il_jsize=im_minext IF(PRESENT(id_jsize)) il_jsize=id_jsize IF( .NOT. td_var%t_dim(1)%l_use ) il_isize=0 IF( .NOT. td_var%t_dim(2)%l_use ) il_jsize=0 CALL logger_trace( "EXTRAP DEL EXTRABANDS: dimension change "//& & "in variable "//TRIM(td_var%c_name) ) ! add extrabands in variable 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(:,:,:,:) ! fill center il_imin=1+il_isize il_imax=td_var%t_dim(1)%i_len-il_isize il_jmin=1+il_jsize il_jmax=td_var%t_dim(2)%i_len-il_jsize td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len - 2*il_isize td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len - 2*il_jsize DEALLOCATE(td_var%d_value) ALLOCATE( td_var%d_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 ) ) ! intialise td_var%d_value(:,:,:,:)=td_var%d_fill td_var%d_value(:,:,:,:)=dl_value(il_imin:il_imax,& & il_jmin:il_jmax,& & :,:) DEALLOCATE( dl_value ) END SUBROUTINE extrap_del_extrabands END MODULE extrap