SUBROUTINE interp_cubic(Vh, n, ndim, bounds, val_flag, X, m, Y, istat) !************************************************************* ! ! Cubic interpolation routine ! ! based upon Keys, R.G. (1981), "Cubic convolution interpolation for digital ! image processing", IEEE, Vol. ASSP-29, pp. 1153-1160. ! ! INPUTS: ! ndim is the number of dimensions of the array Vh ! Bounds(4,ndim) gives the dimensions of the array Vh, where: ! Bounds(1,:) are the lower bounds ! Bounds(2,:) are the upper bounds ! Bounds(3,:) are the number of points on the dimension ! Bounds(4,:) is 1 / step size ! n is the size of array Vh; n = prod(Bounds(3,:)) ! Vh(n) is the data array ! m the number of points to be interpolated ! X(m,ndim) the data points to be interpolated ! OUTPUTS: ! Y(m) the results of the interpolation ! ! code developed for the SIDD microsimulation model ! ! Author: J. van de Ven ! Date: 02 Aug 2016 ! !************************************************************* USE global_ParamStore USE global_CommonAll IMPLICIT none integer(4), intent(in) :: val_flag Integer(4) :: n, m, ndim, istat Real(8) :: Vh(n), Bounds(4,ndim), X(m,ndim), Y(m) !local variables integer(4) :: hh, ii, jj, kk, ll, ndx, ndx_store, dir2, lin_flag integer(4) :: dims(ndim), offset(ndim), mm(ndim), dd0(ndim), dd(ndim), dd_a(ndim), dd1(4**ndim, ndim), nn(ndim), dir(ndim) real(8), parameter :: tol = epsilon(1.0_8) real(8) :: ss(ndim), weight(4**ndim), val(4**ndim), tmp real(8) :: min_chk, max_chk, min_chk2, max_chk2, diff_chk !************************************************************* ! begin code !************************************************************* Y = 0.0 istat = 0 val = -999.0 do ii = 1, ndim dims(ii) = int(Bounds(3,ii)+0.1) if (ii.eq.1_4) then offset(ii) = 1_4 else offset(ii) = product(dims(1:ii-1)) end if end do ii = 1 do while (ii.le.m) ! loop over each iteration point ! identify reference points ! mm is the reference coordinate for the interpolation ! ss is the off-set from the reference point, specified as a fraction of a dimensional step ! dd is a working vector indicating difference off-sets from the reference point ! nn = mm + dd is the co-ordinate of the current test point do jj = 1, ndim ! loop over each interpolation point ss(jj) = 1.0 + (X(ii,jj) - Bounds(1,jj)) * Bounds(4,jj) mm(jj) = floor(ss(jj) + tol) if ( mm(jj).eq.dims(jj) ) then ! upper bound - off-set one step backward dir(jj) = 1 mm(jj) = mm(jj) - 1 dd0(jj) = -1 elseif ( mm(jj).eq.1 ) then ! lower bound - loop backward to facilitate computation of required elements dir(jj) = -1 dd0(jj) = 2 else ! internal - loop forward dir(jj) = 1 dd0(jj) = -1 end if ss(jj) = ss(jj) - real(mm(jj)) end do max_chk = -1.0E+10 min_chk = +1.0E+10 max_chk2 = -1.0E+10 min_chk2 = +1.0E+10 dd = dd0 dd(1) = dd(1) - dir(1) ! here we identify index values, and associated function values ! dd1 is a vector that stores index deviations from mm ! val is a vector that stores function values ! the vectors are organised so that the first element corresponds to the lower ! ordinate of the n-dimensional cube, and the last element to the upper ordinate jj = 1 do while (jj.le.(4**ndim)) ! loop over each test point dd(1) = dd(1) + dir(1) kk = 1 do while ( abs(dd(kk)-dd0(kk)).gt.3 ) dd(kk) = dd0(kk) kk = kk + 1 dd(kk) = dd(kk) + dir(kk) end do ! calculate indices ndx_store = 1 do kk = 1, ndim nn(kk) = mm(kk) + dd(kk) ndx_store = ndx_store + (4**(kk-1)) * (dd(kk)+1) end do dd1(ndx_store,:) = dd ! calculate reference index if ( any(nn.lt.1).or.any(nn.gt.dims) ) then ! outside boundary - need to approximate on basis of internal points kk = 1 dd_a = dd do while ( (nn(kk).ge.1).and.(nn(kk).le.dims(kk)) ) kk = kk + 1 end do if ( nn(kk).lt.1 ) then dir2 = 1 else dir2 = -1 end if do ll = 1, 3 ! loop over each point required to evaluate approximation dd_a(kk) = dd_a(kk) + dir2 ndx = 1 do hh = 1, ndim ndx = ndx + (4**(hh-1)) * (dd_a(hh)+1) end do tmp = val(ndx) if ( abs(tmp + 999.0).lt.1.0E-6 ) then !$OMP CRITICAL(error_1) print *, 'problem in approximating boundary condition in cubic interpolation' pause 'press enter to exit' continue stop continue !$OMP END CRITICAL(error_1) end if if ( ll.eq.1 ) then val(ndx_store) = 3.0 * tmp elseif ( ll.eq.2 ) then val(ndx_store) = val(ndx_store) - 3.0 * tmp else val(ndx_store) = val(ndx_store) + tmp end if end do else ! internal to grid ndx = 1 do kk = 1, ndim ndx = ndx + (nn(kk) - 1) * offset(kk) end do val(ndx_store) = Vh(ndx) if ( abs(val(ndx_store) + 999.0).lt.1.0E-6 ) then istat = 1 jj = 4**ndim + 1 ii = m + 1 else if ( any(dd.gt.1).or.any(dd.lt.0) ) then if ( val(ndx_store).gt.max_chk2 ) then max_chk2 = val(ndx_store) end if if ( val(ndx_store).lt.min_chk2 ) then min_chk2 = val(ndx_store) end if else if ( val(ndx_store).gt.max_chk ) then max_chk = val(ndx_store) end if if ( val(ndx_store).lt.min_chk ) then min_chk = val(ndx_store) end if end if end if end if jj = jj + 1 end do if ( istat.eq.0 ) then diff_chk = max_chk - min_chk lin_flag = 0 if ( diff_chk.lt.(tol*3.0) ) then ! near points define a flat surface - no need to interpolate Y(ii) = max_chk elseif ( (max(abs(min_chk-min_chk2),abs(max_chk2-max_chk))/diff_chk).gt.5.0 ) then ! need to adopt linear interpolation lin_flag = 1 else ! conduct cubic interpolation weight = 1.0 do jj = 1, (4**ndim) dd = dd1(jj,:) do kk = 1, ndim if ( dd(kk).eq.-1 ) then weight(jj) = weight(jj) * ( - ss(kk)**3.0 + 2.0 * ss(kk)**2.0 - ss(kk) ) / 2.0 elseif ( dd(kk).eq.0 ) then weight(jj) = weight(jj) * ( 3.0 * ss(kk)**3.0 - 5.0 * ss(kk)**2.0 + 2.0 ) / 2.0 elseif ( dd(kk).eq.1 ) then weight(jj) = weight(jj) * (-3.0 * ss(kk)**3.0 + 4.0 * ss(kk)**2.0 + ss(kk) ) / 2.0 else weight(jj) = weight(jj) * ( ss(kk)**3.0 - ss(kk)**2.0) / 2.0 end if end do Y(ii) = Y(ii) + val(jj) * weight(jj) end do end if if ( ( (val_flag.eq.1).and.(Y(ii).lt.0.0) ).or.( lin_flag.eq.1 ) ) then ! undertake linear interpolation weight = 0.0 Y(ii) = 0.0 do jj = 1, (4**ndim) if ( any(dd1(jj,:).gt.1).or.any(dd1(jj,:).lt.0) ) then ! skip continue else ! point is within inner call dimension ! calculate weights and indices weight(jj) = 1.0 do kk = 1, ndim weight(jj) = weight(jj) * (1.0 - abs(real(dd1(jj,kk))-ss(kk))) end do Y(ii) = Y(ii) + val(jj) * weight(jj) end if end do tmp = sum(weight) if ( (tmp.gt.1.0E-5).and.(tmp.lt.(1.0+1.0E-5)) ) then Y(ii) = Y(ii) / tmp else !$OMP CRITICAL(error_2) print *, 'problem in weighting values for linear interpolation' pause 'press enter to exit' continue stop continue !$OMP END CRITICAL(error_2) end if end if end if ii = ii + 1 end do END SUBROUTINE interp_cubic