subroutine xint(x, f, n, iflag, lflag, xi, fi, ierr) ! This routine applies a quadratic interpolation procedure ! to f(x) between x(1) and x(n). f(x) is assumed to be ! represented by quadratic polynomials between the points ! x(i). The polynomials are chosen so that they equal f(i) ! at the points x(i), the first derviatives on either ! side of the interior x(i) match at x(i), and the second ! derivative of the approximated function integrated ! over the domain is minimized. ! This version is for interpolating longitude ! Input: x(1),x(2) ... x(n) The x values (must be sequential) ! f(1),f(2) ... f(n) The function values ! n The number of x,f pairs ! iflag Flag for initialization ! =1 for coefficient calculation ! =0 to use previous coefficients ! lflag Flag for linear interpolation ! =0 to perform linear interpolation ! =1 to perform quadratic interpolation ! xi The x value at which to interpolate f ! Output: fi The interpolated function value ! ierr Error flag ! =0 Normal return ! =1 Parameter nmax is too small ! =2 The x values are not sequential ! =3 Coefficient iteration did not ! converge ! =4 Mix-up finding coefficients ! =5 if xi .gt. x(n) or .lt. x(1), ! xi is set to nearest endpoint ! before the interpolation ! Note: fi is set to -99.9 if ! ierr=1,2,3 or 4 ! Modifications: ! 5/13/21 (J. Dostalek) Started this Fortran 90 version !------------------------------------------------------------------------------- implicit none integer, parameter :: nmax = 1000 real, intent(in), dimension(n) :: x real, intent(in), dimension(n) :: f integer, intent(in) :: n integer, intent(in) :: iflag integer, intent(in) :: lflag real, intent(inout) :: xi real, intent(out) :: fi integer, intent(out) :: ierr real, dimension(nmax) :: ax, bx, cx real, dimension(nmax) :: df, dx, gm, ct integer :: idbug integer :: lutest real :: thresh real :: d real :: eps, aa, bb, dev integer :: nit real :: slt, cfsave real :: cf0, rel, dsdc1, cft, den, slo integer :: i, j, k, ii common/xsave/ax, bx, cx ! Specify unit number for debug write statements ! and debug flag idbug = 0 lutest = 6 ! Initialize error flag and output value ierr = 4 fi = -99.9 ! Specify minimum reduction in cost function for convergence thresh = 1.0e-10 ! Check to make sure nmax is large enough, and n is .gt. 1 if (n > nmax .or. n < 2) then ierr = 1 return end if ! Check to make sure x is sequential do i = 1, n - 1 if (x(i) >= x(i + 1)) then ierr = 2 return end if end do if (iflag == 1) then ! Perform the initialization for later interpolation ! Calculate x and f differences do i = 1, n - 1 df(i) = f(i + 1) - f(i) dx(i) = x(i + 1) - x(i) end do ! Calculate domain size d = x(n) - x(1) ! Check for linearity of input points eps = 1.0e-10 bb = (f(2) - f(1))/(x(2) - x(1)) aa = f(1) - bb*x(1) dev = 0.0 do i = 3, n dev = dev + abs(aa + bb*x(i) - f(i)) end do if (dev < eps .or. lflag == 0 .or. n == 2) then ! Linear Interpolation do i = 1, n - 1 cx(i) = 0.0 end do ierr = 0 else ! Quadratic Interpolation ! Iterate to find the c-coefficients cx(1) = 0.0 nit = 100 slt = 0.01 cfsave = 1.0e+10 do k = 1, nit ! Calculate c values do i = 2, n - 1 cx(i) = -cx(i - 1)*dx(i - 1)/dx(i) - df(i - 1)/(dx(i)*dx(i - 1)) + & df(i)/(dx(i)*dx(i)) end do ! Calculate current value of cost function cf0 = 0.0 do i = 1, n - 1 cf0 = cf0 + cx(i)*cx(i)*dx(i) end do cf0 = 0.5*cf0/d if (idbug /= 0) then write (lutest, 101) cf0 101 format(/, ' cf0=', e13.6) end if ! Check for convergence rel = abs(cf0 - cfsave)/abs(cfsave) if (rel < thresh) then if (idbug /= 0) then write (lutest, 104) 104 format(/, ' Iteration converged') end if ierr = 0 exit end if cfsave = cf0 ! Calculate values of Lagrange multipliers gm(n - 1) = cx(n - 1)*dx(n - 1)/d if (n > 3) then do i = n - 2, 2, -1 gm(i) = cx(i)*dx(i)/d - gm(i + 1)*dx(i)/dx(i + 1) end do end if ! Calculate gradient of cost function with respect to c1 dsdc1 = dx(1)*(cx(1)/d - gm(2)/dx(2)) ! Adjust cx(1) using trial step ct(1) = cx(1) - slt*dsdc1 ! Calculate remaining c values at trial step do i = 2, n - 1 ct(i) = -ct(i - 1)*dx(i - 1)/dx(i) - df(i - 1)/(dx(i)*dx(i - 1)) + & df(i)/(dx(i)*dx(i)) end do ! Calculate cost function at trial step cft = 0.0 do i = 1, n - 1 cft = cft + ct(i)*ct(i)*dx(i) end do cft = 0.5*cft/d ! Calculate optimal step length and re-adjust cx(1) den = 2.0*((cft - cf0) + slt*dsdc1*dsdc1) if (den /= 0.0) then slo = dsdc1*dsdc1*slt*slt/den else slo = 0.0 end if ! Adjust slo if desired slo = 1.0*slo cx(1) = cx(1) - slo*dsdc1 if (idbug /= 0) then write (lutest, 100) k, cft, slt, slo 100 format(' Iteration=', i4, ' cf1=', e11.4, ' slt=', e11.4, ' slo=', e11.4) ! do j = 1, n - 1 write (lutest, 102) j, cx(j) 102 format(' i=', i2, ' c=', f8.4) end do end if ! Calculate trial step for next time step slt = 0.5*slo end do ! If iteration did not converge, exit subroutine if (ierr /= 0) then ierr = 3 return end if end if ! Calculate b and a coefficients do i = 1, n - 1 bx(i) = df(i)/dx(i) - cx(i)*(x(i + 1) + x(i)) ax(i) = f(i) - bx(i)*x(i) - cx(i)*x(i)*x(i) end do end if ! Interpolate the function ! Check for xi out of bounds if (xi < x(1)) then xi = x(1) ierr = 5 end if if (xi > x(n)) then xi = x(n) ierr = 5 end if ! Find the interval for the interpolation ii = 1 do i = 2, n if (xi <= x(i)) then ii = i - 1 fi = ax(ii) + bx(ii)*xi + cx(ii)*xi*xi if (ierr /= 5) ierr = 0 exit end if end do return end subroutine xint