mnewtax.f90

Go to the documentation of this file.

subroutine mnewtax(ntrial, naxis, x, n, tolx, tolf, errx, errf)
!------------------------------------------------------------------------- 
! routine to solve two nonlinear equations using newtons method from
! numerical recipes.
! lu decomposition replaced by explicit solution of 2x2 matrix. 
!-------------------------------------------------------------------------

  use mod_corners
  use mod_mesh
  use itm_types
  use mod_interpolation

  implicit none

  real(r8), intent(in) :: tolx, tolf
  real(r8), intent(out) :: errx, errf
  integer(itm_i4), intent(in) :: ntrial, naxis, n
  real(r8), intent(inout) :: x(2)
  
  real(r8) :: fvec(2), p(2), fjac(2,2)
  real(r8) :: r, s
  real(r8) :: zpsi, zpsir, zpsis, zpsirs, zpsirr, zpsiss
  real(r8) :: temp, dis
  integer(itm_i4) :: i, k
  integer(itm_i4) :: n1, n2, n3, n4

  do k = 1, ntrial
!----------------------------- usrfun iserted here -----------------------
    do i = 1, n
      if (abs(x(i)) > 2.0_r8) then
        errx = 2._r8
        errf = 2._r8
        return
      end if
    end do

    r = x(1)
    s = x(2)
    n1 = nodeno(naxis, 1)
    n2 = nodeno(naxis, 2)
    n3 = nodeno(naxis, 3)
    n4 = nodeno(naxis, 4)
    call interpolation(2, psi(4 * (n1 - 1) + 1 : 4 * n1), psi(4 * (n2 - 1) &
     + 1 : 4 * n2), psi(4 * (n3 - 1) + 1 : 4 * n3), psi(4 * (n4 - 1) + 1 &
     : 4 * n4), r, s, zpsi, zpsir, zpsis, zpsirs, zpsirr, zpsiss)
    fvec(1) = zpsir
    fvec(2) = zpsis
    fjac(1, 1) = zpsirr
    fjac(1, 2) = zpsirs
    fjac(2, 1) = zpsirs
    fjac(2, 2) = zpsiss
!-------------------------------------------------------------------------
    errf = 0._r8
    do i = 1, n
      errf = errf + abs(fvec(i))
    end do
    if (errf <= tolf) return
    do i = 1, n
      p(i) = -fvec(i)
    end do
    temp = p(1)
    dis = fjac(2, 2) * fjac(1, 1) - fjac(1, 2) * fjac(2, 1)
    p(1) = (fjac(2, 2) * p(1) - fjac(1, 2) * p(2)) / dis
    p(2) = (fjac(1, 1) * p(2) - fjac(2, 1) * temp) / dis      
    
    p = min(p, 0.5_r8)
    p = max(p, -0.5_r8)

    errx = 0._r8
    do i = 1, n
      errx = errx + abs(p(i))
      x(i) = x(i) + p(i)
    end do
    if (errx <= tolx) return
  end do

  return
end subroutine mnewtax

!  (c) copr. 1986-92 numerical recipes software *n*1v45_lt+v'.