grid2nv.f90

Go to the documentation of this file.

subroutine grid2nv(tin, tout, jpts, acc, igrd, ll)
!---------------------------------------------------------------------------
!  The function tin(tout), given on the grid tout = 2 * pi * (j - 1) / jpts,
!  is inverted to give tout(tin) on the grid tin = 2 * pi * (i - 1) / jpts.
!  This is done by determining the zeros of the function
!     y(t) = t + sum(gf(m) * sin(m * t)) - 2 * pi * (i - 1) / jpts,
!  where gf(m) are the fourier coefficients of g(t) = tin(t) - t.
!  Diagnostic information is printed if l /= 0.
!  The argument ll has the double function of communicating the
!  output unit iout = ll / 10 and the print switch l = ll - iout * 10.
!  iout = 0 (l=ll): file is "output".
!-----------------------------------------------------------------------  

  use itm_types

  implicit none

  integer(itm_i4), intent(in) :: jpts
  real(r8), dimension(jpts), intent(inout) :: tin
  real(r8), dimension(jpts), intent(out) :: tout
  real(r8), intent(in) :: acc
  integer(itm_i4), intent(in) :: igrd, ll

  real(r8), dimension(jpts + 1) :: t, g  
  real(r8), dimension(jpts / 2 - 1) :: gfcos, gfsin

  real(r8), parameter :: pi = 3.1415926535898_r8
  integer(itm_i4), parameter :: ninv = 100 
  real(r8) :: t1, y1, t0, y0, t2, y2  
  real(r8) :: sum1, sum2, gfnul
  integer(itm_i4) :: mharm
  integer(itm_i4) :: iout, iout_pr
  integer(itm_i4) :: l, i, j, jj, j1, n
  integer(itm_i4) :: igrdnv, ifirst, icirc

  mharm = jpts / 2 - 1

  do jj = 2, jpts
    if (tin(jj - 1) > tin(jj))  tin(jj) = tin(jj) + 2._r8 * pi
  end do
  iout = ll / 10
  l = ll - iout * 10
  iout_pr = iout * 10 + l
  if (iout == 0) iout = 6
  if (l /= 0) write(iout, 3)

  do i = 1, jpts
    g(i) = tin(i) - 2._r8 * pi * (i - 1) / dble(jpts)
  end do
  call rft(g, gfnul, gfcos, gfsin, jpts, mharm)
  if (l /= 0) then
    call prarr1('tin(j) : ', tin, jpts, iout_pr)
    call prarr1('g(i):', t, jpts, iout_pr)
    write(iout, 56) gfnul
    call prarr1('gfcos(m):', gfcos, mharm, iout_pr)
    call prarr1('gfsin(m):', gfsin, mharm, iout_pr)
  end if
  do i = 1, jpts + 1
    t(i) = 2._r8 * pi * (i - 1) / dble(jpts)
  end do
  j1 = 1
  igrdnv = 1
  ifirst = 0
  icirc = -(int(tin(1) / (2._r8 * pi) + 10000) -  9999)
  if (abs(tin(1)) < 1.e-12_r8) icirc = 0

  loop_i: do i = 1, jpts
    j = j1
    t1 = t(j) + icirc * 2._r8 * pi
    call fsum2(sum1, t1, gfnul, gfcos, gfsin, mharm)
    y1 = t1 + sum1 - t(i)
    do
      t0 = t1
      y0 = y1
      if (abs(y0) <= acc) then
        tout(i) = t0
        cycle loop_i
      end if
      if (j == jpts + 1) then
        if (ifirst == 0)  then
          j = 1
          icirc = icirc + 1
          ifirst = 1
        else
          write(iout, 70)
          return
        end if
      else
        j = j + 1
      end if
      t1 = t(j) + icirc * 2._r8 * pi
      call fsum2(sum1, t1, gfnul, gfcos, gfsin, mharm)
      y1 = t1 + sum1 - t(i)
      if (sign(1._r8, y0) /= sign(1._r8, y1)) exit
    end do

    j1 = j - 1

    do n = 1, ninv
      t2 = t0 - (t1 - t0) * y0 / (y1 - y0)
      call fsum2(sum2, t2, gfnul, gfcos, gfsin, mharm)
      y2 = t2 + sum2 - t(i)
      if (l /= 0) write(iout, 25) n, t0, t1, y0, y1, t2, y2, j
      if (abs(y2) <= acc) exit
      if (sign(1._r8, y2) /= sign(1._r8, y1)) then
        t0 = t2
        y0 = y2
      else
        t1 = t2
        y1 = y2
      end if
    end do

    tout(i) = t2
    if (l /= 0) write(iout, 55) i, n
    if (n > igrdnv) igrdnv = n

  end do loop_i

  return

    3 format(///1x, 'subroutine gridinv')
   25 format(1x, 'n=', i3, ' t0=', f10.5, ' t1=', f10.5, &
       ' y0=', f10.5, ' y1=', f10.5, ' t2=', f10.5, ' y2=', f10.5, ' j=',i3)
   55 format(1x, 'i=', i3, 5x, 'n=', i3)
   56 format(/1x, 'gfnul = ', 1pe12.4)
   70 format(/1x, '***subroutine gridinv: no zero found ')

end subroutine grid2nv