!= ガウス重み, 分点の計算
!
!= Calculate Gauss node and Gaussian weight
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: gauss_quad.f90,v 1.4 2014/05/07 09:39:20 murashin Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module gauss_quad
  !
  != ガウス重み, 分点の計算
  !
  != Calculate Gauss node and Gaussian weight
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 
  !
  ! 
  !
  !== References
  !
  !== Procedures List
  !
!!$  ! RadiationFluxDennouAGCM :: 放射フラックスの計算
!!$  ! RadiationDTempDt        :: 放射フラックスによる温度変化の計算
!!$  ! RadiationFluxOutput     :: 放射フラックスの出力
!!$  ! RadiationFinalize       :: 終了処理 (モジュール内部の変数の割り付け解除)
!!$  ! ------------            :: ------------
!!$  ! RadiationFluxDennouAGCM :: Calculate radiation flux
!!$  ! RadiationDTempDt        :: Calculate temperature tendency with radiation flux
!!$  ! RadiationFluxOutput     :: Output radiation fluxes
!!$  ! RadiationFinalize       :: Termination (deallocate variables in this module)
  !
  !== NAMELIST
  !
!!$  ! NAMELIST#radiation_DennouAGCM_nml
  !

  ! USE statements
  !

  ! 
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! Double precision.
    &                 STRING, &  ! Strings.
    &                 TOKEN      ! Keywords.

  ! 物理・数学定数設定
  ! Physical and mathematical constants settings
  !
  use constants0, only: &
    & PI                    ! $ \pi $ .
                            ! Circular constant

  ! Declaration statements
  !
  implicit none
  private

  ! 
  ! Public procedure
  !
  public:: GauLeg


  character(*), parameter:: module_name = 'gauss_quad'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: gauss_quad.f90,v 1.4 2014/05/07 09:39:20 murashin Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

  !--------------------------------------------------------------------------------------
  ! This code is the same as sh_gauleg subrouting in sh_module.f90.
  ! This code was obtained from Press et al. [1993], p141.
  !--------------------------------------------------------------------------------------

  subroutine GauLeg( x1, x2, n, x, w )

    real(DP), intent(in ) :: x1,x2
    integer , intent(in ) :: n
    real(DP), intent(out) :: x(n),w(n)


    real(DP):: eps
    real(DP):: z1, z, xm, xl, pp, p3, p2, p1
    integer :: i, j, m


    !
    ! changed at 2005/09/14
    !
!!$      eps = 1.0e-11
    eps = 1.0d-15

    m  = ( n+1 ) / 2
    xm = ( x2+x1 ) / 2.0d0
    xl = ( x2-x1 ) / 2.0d0

    do i = 1, m
      z = cos( pi*dble(i-0.25d0)/(dble(n+0.5d0)) )

100   continue
      p1 = 1.0d0
      p2 = 0.0d0
!!$         do j = 1, n
      do j = 1, n
        p3 = p2
        p2 = p1
        p1 = ( (2.0d0*dble(j)-1.0d0)*z*p2-(dble(j)-1.0d0)*p3 ) / dble(j)
      enddo

      pp = dble(n) * (z*p1-p2) / (z*z-1.0d0)
      z1 = z
      z  = z1 - p1 / pp

      !
      ! changed at 2005/09/14
      !
!!$         if( abs( z-z1 ) .gt. eps ) go to 100
!!$         if( abs( z-z1 ) / z1 .gt. eps ) go to 100
      !
      ! changed before 2008/08/03
      !
      if( abs( z-z1 ) / abs( z1 + 1.0d-200 ) .gt. eps ) go to 100

      x( i     ) = xm-xl*z
      x( n+1-i ) = xm+xl*z
      w( i     ) = 2.0d0*xl/((1.0d0-z*z)*pp*pp)
      w( n+1-i ) = w(i)
    end do


  end subroutine GauLeg

  !--------------------------------------------------------------------------------------

end module gauss_quad
