!---------------------------------------------------------------------- ! NEMO system team, System and Interface for oceanic RElocable Nesting !---------------------------------------------------------------------- ! ! MODULE: extrap ! ! DESCRIPTION: !> @brief !> This module !> !> @details !> !> @author !> J.Paul ! REVISION HISTORY: !> @date Nov, 2013 - Initial Version !> !> @note WARNING: FillValue must not be zero (use var_chg_FillValue) !> !> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !> @todo !> - something wrong when computing point to be extralopated !> - take care of ew value in variable structure !---------------------------------------------------------------------- MODULE extrap USE netcdf ! nf90 library USE kind USE phycst USE global USE fct USE logger USE dim USE att USE var IMPLICIT NONE PRIVATE ! NOTE_avoid_public_variables_if_possible ! type and variable ! function and subroutine PUBLIC :: extrap_detect !< detected point to be extrapolated PUBLIC :: extrap_fill_value !< extrapolate value over detected point PUBLIC :: extrap_add_extrabands !< PUBLIC :: extrap_del_extrabands !< PUBLIC :: extrap_deriv_1D !< PUBLIC :: extrap_deriv_2D !< PRIVATE :: extrap__detect_wrapper !< detected point to be extrapolated PRIVATE :: extrap__detect !< detected point to be extrapolated PRIVATE :: extrap__fill_value_wrapper !< extrapolate value over detected point PRIVATE :: extrap__fill_value !< extrapolate value over detected point PRIVATE :: extrap__3D PRIVATE :: extrap_deriv_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. !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] td_var : variable to extrapolate !> @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 !> @return !------------------------------------------------------------------- !> @code FUNCTION extrap__detect( td_var0, td_level1, & & id_offset, id_rho, id_ext ) IMPLICIT NONE ! Argument TYPE(TVAR) , INTENT(INOUT) :: 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_varid INTEGER(i4) , DIMENSION(:,:,:), ALLOCATABLE :: il_detect 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 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) IF( PRESENT(id_offset) )THEN il_offset(1:SIZE(id_offset(:,:),DIM=1),& & 1:SIZE(id_offset(:,:),DIM=2) )= id_offset(:,:) 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 WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=1 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_varid=var_get_id(td_level1(:),TRIM(cl_level)) IF( il_varid == 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 print *,'read ',trim(cl_level) tl_var1=td_level1(il_varid) ENDIF 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(ig_ndim,2) ) ! coarsening fine grid level il_extra(jp_I,1)=CEILING(REAL(il_rho(jp_I)-1,dp)/2._dp) il_extra(jp_I,2)=FLOOR(REAL(il_rho(jp_I)-1,dp)/2._dp) il_extra(jp_J,1)=CEILING(REAL(il_rho(jp_J)-1,dp)/2._dp) il_extra(jp_J,2)=FLOOR(REAL(il_rho(jp_J)-1,dp)/2._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 !il_level1_G0(:,:)=0 !DO jj1=1,tl_var1%t_dim(2)%i_len ! jj0=INT(REAL((jj1+il_offset(jp_J,1)-1)-1,dp)/REAL(il_rho(jp_J),dp)) +1 ! DO ji1=1,tl_var1%t_dim(1)%i_len ! ji0=INT(REAL((ji1+il_offset(jp_I,1)-1)-1,dp)/REAL(il_rho(jp_I),dp)) +1 ! il_level1_G0(ji0,jj0)=MAX(il_level1_G0(ji0,jj0),il_level1(ji1,jj1)) ! ! ENDDO !ENDDO ! clean DEALLOCATE( il_extra ) DEALLOCATE( il_level1 ) ENDIF ! look for sea point !il_detect(:,:,1)=0 DO jk0=1,td_var0%t_dim(3)%i_len WHERE( il_level1_G0(:,:) >= jk0) il_detect(:,:,jk0)=1 END WHERE !il_detect(:,:,jk0)=il_level1_G0(:,:) !WHERE( td_var0%d_value(:,:,jk0,1) /= td_var0%d_fill ) ! il_detect(:,:,1)=jk0-1 !END WHERE ENDDO ! clean DEALLOCATE( il_level1_G0 ) ENDIF ! clean DEALLOCATE( il_offset ) !! select extra point depending on interpolation method !! compute point near grid point already inform !FORALL( ji0=1:il_dim0(1), & !& jj0=1:il_dim0(2), & !& jk0=1:il_dim0(3), & !& il_detect(ji0,jj0,jk0) == 1 ) ! 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 !END FORALL ! do not compute grid point already inform WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=0 ! save result extrap__detect(:,:,:)=il_detect(:,:,:) ! clean DEALLOCATE( il_dim0 ) DEALLOCATE( il_ext ) DEALLOCATE( il_detect ) END FUNCTION extrap__detect !> @endcode !------------------------------------------------------------------- !> @brief !> This function detected point to be extrapolated. !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] td_var : variable to extrapolate !> @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 !> @return !------------------------------------------------------------------- !> @code FUNCTION extrap__detect_wrapper( td_var, td_level, & & id_offset, id_rho, id_ext ) 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), 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 interpolated 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 interpolated 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 interpolated 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 !> @endcode ! !------------------------------------------------------------------- ! !> @brief ! !> This function detected point to be extrapolated. ! !> ! !> @details ! !> ! !> @author J.Paul ! !> - Nov, 2013- Initial Version ! ! ! !> @param[in] td_var : variable to extrapolate ! !> @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 ! !> @return ! ! ! !> @todo ! !------------------------------------------------------------------- ! !> @code ! FUNCTION extrap__detect(td_var, id_iext, id_jext, id_kext) ! IMPLICIT NONE ! ! Argument ! TYPE(TVAR) , INTENT(INOUT) :: td_var ! !INTEGER(i4), DIMENSION(:,:,:), INTENT(OUT ) :: id_detect ! INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext ! INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext ! INTEGER(i4), INTENT(IN ), OPTIONAL :: id_kext ! ! ! 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 ! ! ! local variable ! INTEGER(i4) :: il_iext ! INTEGER(i4) :: il_jext ! INTEGER(i4) :: il_kext ! ! INTEGER(i4), DIMENSION(3) :: il_dim ! ! ! loop indices ! INTEGER(i4) :: ji ! INTEGER(i4) :: jj ! INTEGER(i4) :: jk ! !---------------------------------------------------------------- ! ! ! optional argument ! 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=im_minext ! IF( PRESENT(id_kext) ) il_kext=id_kext ! ! ! init ! extrap__detect(:,:,:)=0 ! ! il_dim(:)=td_var%t_dim(1:3)%i_len ! ! ! compute point near grid point already inform ! FORALL( ji=1:il_dim(1), & ! & jj=1:il_dim(2), & ! & jk=1:il_dim(3), & ! & td_var%d_value(ji,jj,jk,1) /= td_var%d_fill ) ! ! extrap__detect(MAX(1,ji-il_iext):MIN(ji+il_iext,il_dim(1)),& ! & MAX(1,jj-il_jext):MIN(jj+il_jext,il_dim(2)),& ! & MAX(1,jk-il_kext):MIN(jk+il_kext,il_dim(3)) )=1 ! ! ! END FORALL ! ! ! do not compute grid point already inform ! WHERE( td_var%d_value(:,:,:,1) /= td_var%d_fill ) extrap__detect(:,:,:)=0 ! ! END FUNCTION extrap__detect ! !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[inout] td_var : variable structure !> @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_extend : radius of the box used to compute extrapolation !> @param[in] id_maxiter : maximum nuber of iteration !------------------------------------------------------------------- !> @code 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 table of 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) .OR. & & PRESENT(td_level) )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 !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 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 nuber of iteration !------------------------------------------------------------------- !> @code 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 INTEGER(i4) :: il_radius INTEGER(i4) :: il_iter TYPE(TATT) :: tl_att ! loop indices INTEGER(i4) :: jl !---------------------------------------------------------------- !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 logger_warn(" EXTRAP FILL: "//& & TRIM(fct_str(SUM(il_detect(:,:,:))))//& & " point(s) to extrapolate " ) !3- extrapolate DO jl=1,td_var%t_dim(4)%i_len ! td_var%d_value(:,:,:,jl)=il_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) CALL logger_debug(" EXTRAP FILL VALUE: "//& & TRIM(fct_str(SUM(il_detect(:,:,:))))//& & " points to extrapolate " ) CALL extrap__3D(td_var%d_value(:,:,:,jl), td_var%d_fill, & & il_detect(:,:,:), & & cd_method, il_radius ) il_iter=il_iter+1 ENDDO ENDDO IF( SUM(il_detect(:,:,:)) /= 0 )THEN CALL logger_warn(" EXTRAP FILL: still "//& & TRIM(fct_str(SUM(il_detect(:,:,:))))//& & " point(s) to extrapolate " ) ENDIF DEALLOCATE(il_detect) END SUBROUTINE extrap__fill_value !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[inout] dd_value : 3D table of variable to be extrapolated !> @param[in] dd_fill : FillValue of variable !> @param[inout] id_detect : table of point to extrapolate !> @param[in] id_ext : number of point use to compute box !> @param[in] cd_method : extrapolation method !------------------------------------------------------------------- !> @code SUBROUTINE extrap__3D( dd_value, dd_fill, id_detect,& & cd_method, id_ext ) 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_ext ! 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 INTEGER(i4), DIMENSION(3) :: il_dim 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 !---------------------------------------------------------------- il_shape(:)=SHAPE(dd_value) SELECT CASE(TRIM(cd_method)) CASE('min_error') 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(:,:,:), dd_fill, 'I' ) ENDIF ! compute derivative in i-direction dl_dfdy(:,:,:)=dd_fill IF( il_shape(2) > 1 )THEN dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'J' ) ENDIF ! compute derivative in i-direction dl_dfdz(:,:,:)=dd_fill IF( il_shape(3) > 1 )THEN dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'K' ) ENDIF il_dim(1)=2*id_ext+1 IF( il_shape(1) < 2*id_ext+1 ) il_dim(1)=1 il_dim(2)=2*id_ext+1 IF( il_shape(2) < 2*id_ext+1 ) il_dim(2)=1 il_dim(3)=2*id_ext+1 IF( il_shape(3) < 2*id_ext+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))) DO jk=1,il_shape(3) IF( ALL(id_detect(:,:,jk) == 0) ) CYCLE DO jj=1,il_shape(2) IF( ALL(id_detect(:,jj,jk) == 0) ) CYCLE DO ji=1,il_shape(1) IF( id_detect(ji,jj,jk) == 1 )THEN il_imin=MAX(ji-id_ext,1) il_imax=MIN(ji+id_ext,il_shape(1)) IF( il_dim(1) == 1 )THEN il_imin=ji il_imax=ji ENDIF il_jmin=MAX(jj-id_ext,1) il_jmax=MIN(jj+id_ext,il_shape(2)) IF( il_dim(2) == 1 )THEN il_jmin=jj il_jmax=jj ENDIF il_kmin=MAX(jk-id_ext,1) il_kmax=MIN(jk+id_ext,il_shape(3)) IF( il_dim(3) == 1 )THEN il_kmin=jk il_kmax=jk ENDIF dd_value(ji,jj,jk)=extrap__3D_min_error_fill( & & dd_value( il_imin:il_imax, & & il_jmin:il_jmax, & & il_kmin:il_kmax ), dd_fill, id_ext, & & 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) /= dd_fill )THEN id_detect(ji,jj,jk)= 0 ENDIF ENDIF ENDDO ENDDO ENDDO IF( ALLOCATED(dl_dfdx) ) DEALLOCATE( dl_dfdx ) IF( ALLOCATED(dl_dfdy) ) DEALLOCATE( dl_dfdy ) IF( ALLOCATED(dl_dfdz) ) DEALLOCATE( dl_dfdz ) IF( ALLOCATED(dl_coef) ) DEALLOCATE( dl_coef ) CASE DEFAULT ! 'dist_weight' il_dim(1)=2*id_ext+1 IF( il_shape(1) < 2*id_ext+1 ) il_dim(1)=1 il_dim(2)=2*id_ext+1 IF( il_shape(2) < 2*id_ext+1 ) il_dim(2)=1 il_dim(3)=2*id_ext+1 IF( il_shape(3) < 2*id_ext+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)) ) DO jk=1,il_shape(3) IF( ALL(id_detect(:,:,jk) == 0) ) CYCLE DO jj=1,il_shape(2) IF( ALL(id_detect(:,jj,jk) == 0) ) CYCLE DO ji=1,il_shape(1) IF( id_detect(ji,jj,jk) == 1 )THEN il_imin=MAX(ji-id_ext,1) il_imax=MIN(ji+id_ext,il_shape(1)) IF( il_dim(1) == 1 )THEN il_imin=ji il_imax=ji ENDIF il_jmin=MAX(jj-id_ext,1) il_jmax=MIN(jj+id_ext,il_shape(2)) IF( il_dim(2) == 1 )THEN il_jmin=jj il_jmax=jj ENDIF il_kmin=MAX(jk-id_ext,1) il_kmax=MIN(jk+id_ext,il_shape(3)) IF( il_dim(3) == 1 )THEN il_kmin=jk il_kmax=jk ENDIF dd_value(ji,jj,jk)=extrap__3D_dist_weight_fill( & & dd_value( il_imin:il_imax, & & il_jmin:il_jmax, & & il_kmin:il_kmax ), dd_fill, id_ext, & & dl_coef(:,:,:) ) IF( dd_value(ji,jj,jk) /= dd_fill )THEN id_detect(ji,jj,jk)= 0 ENDIF ENDIF ENDDO ENDDO ENDDO IF( ALLOCATED(dl_coef) ) DEALLOCATE( dl_coef ) END SELECT END SUBROUTINE extrap__3D !> @endcode !------------------------------------------------------------------- !> @brief !> This function compute derivative of a function in i- and !> j-direction !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] dd_value : 1D table 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 !------------------------------------------------------------------- !> @code 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 !> @endcode !------------------------------------------------------------------- !> @brief !> This function compute derivative of a function in i- and !> j-direction !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] dd_value : 2D table 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 !------------------------------------------------------------------- !> @code 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 !> @endcode !------------------------------------------------------------------- !> @brief !> This function compute derivative of a function in i- and !> j-direction !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[inout] dd_value : 3D table 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 !------------------------------------------------------------------- !> @code 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 !> @endcode !------------------------------------------------------------------- !> @brief !> This function compute extrapolatd values by calculated minimum error using !> taylor expansion !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] dd_value : 3D table of variable to be extrapolated !> @param[in] dd_fill : FillValue of variable !> @param[in] dd_ideriv : derivative of function in i-direction !> @param[in] dd_jderiv : derivative of function in j-direction !> @param[in] dd_kderiv : derivative of function in k-direction !> @param[in] id_ji : i-direction indices table !> @param[in] id_jj : j-direction indices table !> @param[in] id_ii : i-direction indices of the point to extrapolate !> @param[in] id_ij : j-direction indices of the point to extrapolate !------------------------------------------------------------------- !> @code 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 !> @endcode !------------------------------------------------------------------- !> @brief !> This function compute extrapolatd values by calculated minimum error using !> taylor expansion !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] dd_value : 3D table of variable to be extrapolated !> @param[in] dd_fill : FillValue of variable !> @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 : !------------------------------------------------------------------- !> @code PURE FUNCTION extrap__3D_min_error_fill( dd_value, dd_fill, id_ext, & & 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_ext 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_ext*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 !> @endcode !------------------------------------------------------------------- !> @brief !> This function compute extrapolatd values by calculated minimum error using !> taylor expansion !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] dd_value : 3D table of variable to be extrapolated !------------------------------------------------------------------- !> @code 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 !> @endcode !------------------------------------------------------------------- !> @brief !> This function compute extrapolatd values by calculated minimum error using !> taylor expansion !> !> @details !> !> @author J.Paul !> - Nov, 2013- Initial Version ! !> @param[in] dd_value : 3D table of variable to be extrapolated !> @param[in] dd_fill : FillValue of variable !> @param[in] dd_coef : !------------------------------------------------------------------- !> @code FUNCTION extrap__3D_dist_weight_fill( dd_value, dd_fill, id_ext, & & dd_coef ) IMPLICIT NONE ! Argument REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value REAL(dp) , INTENT(IN) :: dd_fill INTEGER(i4), INTENT(IN) :: id_ext 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_ext*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 FORALL( ji=1:il_shape(1), & & jj=1:il_shape(2), & & jk=1:il_shape(3), & & dd_value(ji,jj,jk) /= dd_fill ) ! 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) END FORALL ! 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 !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine add to the variable (to be extrapolated) an !> extraband of N points at north,south,east and west boundaries. !> !> @author J.Paul !> - 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 !------------------------------------------------------------------- !> @code 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 !> @endcode !------------------------------------------------------------------- !> @brief !> This subroutine remove of the variable an extraband !> of N points at north,south,east and west boundaries. !> !> @author J.Paul !> - 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) !------------------------------------------------------------------- !> @code 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 !> @endcode ! !------------------------------------------------------------------- ! !> @brief ! !> This function ! !> ! !> @details ! !> ! !> @author J.Paul ! !> - Nov, 2013- Initial Version ! ! ! !> @param[in] ! !> @param[out] ! !------------------------------------------------------------------- ! !> @code ! FUNCTION extrap_( ) ! IMPLICIT NONE ! ! Argument ! ! ! local variable ! ! ! loop indices ! !---------------------------------------------------------------- ! END FUNCTION extrap_ ! !> @endcode ! !------------------------------------------------------------------- ! !> @brief ! !> This subroutine ! !> ! !> @details ! !> ! !> @author J.Paul ! !> - Nov, 2013- Initial Version ! ! ! !> @param[in] ! !> @param[out] ! !------------------------------------------------------------------- ! !> @code ! SUBROUTINE extrap_( ) ! IMPLICIT NONE ! ! Argument ! ! ! local variable ! ! ! loop indices ! !---------------------------------------------------------------- ! END SUBROUTINE extrap_ ! !> @endcode END MODULE extrap