SUBROUTINE interp_lin(Vh, n, ndim, bounds, X, m, Y, istat) !************************************************************* ! ! Linear spline interpolation routine ! ! 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 :: n,m,ndim,istat Real(8) :: Vh(n),Bounds(4,ndim),X(m,ndim),Y(m) !local variables integer(4) :: ii, jj, kk, ndx integer(4) :: dims(ndim), offset(ndim), mm(ndim), nn(ndim), dd(ndim) real(8), parameter :: tol = epsilon(1.0_8) real(8) :: ss(ndim), weight(2**ndim), tmp !************************************************************* ! begin code !************************************************************* Y = 0.0 istat = 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) ! identify reference points do jj = 1, ndim 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 ! step mm one backward mm(jj) = mm(jj) - 1 end if ss(jj) = ss(jj) - real(mm(jj)) end do dd = 0 dd(1) = -1 jj = 1 do while (jj.le.2**ndim) ! loop over each test point dd(1) = dd(1) + 1 kk = 1 do while ( dd(kk).gt.1 ) dd(kk) = 0 kk = kk + 1 dd(kk) = dd(kk) + 1 end do ! calculate weights and indices weight(jj) = 1.0 do kk = 1, ndim nn(kk) = mm(kk) + dd(kk) weight(jj) = weight(jj) * (1.0 - abs(real(dd(kk))-ss(kk))) end do if ( weight(jj).gt.(1.0/(2.0**real(ndim))*1.0E-3) ) then ! take point into consideration ! calculate reference index ndx = 1 do kk = 1, ndim ndx = ndx + (nn(kk) - 1) * offset(kk) end do tmp = Vh(ndx) if ( abs(tmp + 999.0).lt.1.0E-6 ) then istat = 1 ii = m + 1 jj = 2**ndim + 1 else Y(ii) = Y(ii) + tmp * weight(jj) end if end if jj = jj + 1 end do if ( istat.eq.0 ) then 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 ii = ii + 1 end do END SUBROUTINE interp_lin