!c Description:
!c   収束条件.
!c   ニュートン法の収束条件は, 以下の本を参考にした. 
!c     川上一郎 著, 数値計算の基礎, 日大計算物理教室 (1997)
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Histry: 
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2003-03-03    作成
!c   2.0        2004-07-13    相対誤差, 絶対誤差を利用するように修正
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2004, All rights reserved


module interface_ChkErr
  interface ChkErr
     
     subroutine ChkErr_1d(pre, nxt, ErrR, ErrA, cnvsw)
       real(8), intent(in)  :: pre
       real(8), intent(in)  :: nxt
       real(8), intent(in)  :: ErrR
       real(8), intent(in)  :: ErrA
       logical, intent(out) :: cnvsw
     end subroutine ChkErr_1d
     
     subroutine ChkErr_2d(pre, nxt, ErrR, ErrA, cnvsw)
       real(8), intent(in)  :: pre(:)
       real(8), intent(in)  :: nxt(:)
       real(8), intent(in)  :: ErrR
       real(8), intent(in)  :: ErrA
       logical, intent(out) :: cnvsw
     end subroutine ChkErr_2d
     
  end interface
end module interface_ChkErr




subroutine ChkErr_1d(pre, nxt, ErrR, ErrA, cnvsw)
  implicit none
  real(8), intent(in)  :: pre
  real(8), intent(in)  :: nxt
  real(8), intent(in)  :: ErrR
  real(8), intent(in)  :: ErrA
  logical, intent(out) :: cnvsw

  real(8)              :: chg
  real(8)              :: cnv


  !変数の初期化
  cnvsw = .false. 

  !許容誤差の範囲をチェック. 
  !許容誤差は許容絶対誤差と許容相対誤差から見積もる
  chg = abs(nxt - pre)
  cnv = ErrA + ErrR * (abs(pre) + abs(nxt)) / 2.0d0
  
  !許容誤差の範囲になったら sw を入れる
  if (chg < cnv) cnvsw = .true. 

end subroutine ChkErr_1d


subroutine ChkErr_2d(pre, nxt, ErrR, ErrA, cnvsw)
  use interface_display 

  implicit none
  real(8), intent(in)  :: pre(:)
  real(8), intent(in)  :: nxt(:)
  real(8), intent(in)  :: ErrR
  real(8), intent(in)  :: ErrA
  logical, intent(out) :: cnvsw
  
  real(8), allocatable :: chg(:)
  real(8), allocatable :: cnv(:)
  integer              :: i, n, m
  
  !変数の初期化
  n = size(pre, 1)
  allocate(chg(n), cnv(n))
  cnvsw = .false. 
  m = 0
  
  !許容誤差の範囲をチェック. 
  !許容誤差は許容絶対誤差と許容相対誤差から見積もる
  chg = abs(nxt - pre)
  cnv = ErrA + ErrR * (abs(pre) + abs(nxt))
  
  !許容誤差に収まっている数を数える
  do i = 1, n
     if (chg(i) <= cnv(i)) m = m + 1
  end do

  !全ての値が許容範囲に収まったらスイッチを入れる
  if (m == n) cnvsw = .true. 

!  call display(cnv - chg, "収束チェック(正なら収束)")

end subroutine ChkErr_2d
