!= ガウス重み, 分点の計算
!
!= Calculate Gauss node and Gaussian weight
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: gauss_quad.f90,v 1.3 2011-06-19 11:05:23 yot 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.

  ! 物理・数学定数設定
  ! 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.3 2011-06-19 11:05:23 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

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

  subroutine GauLeg( x1, x2, n, a_x, a_w )

    real(DP), intent(in ) :: x1,x2
    integer , intent(in ) :: n
    real(DP), intent(out) :: a_x(n)
    real(DP), intent(out) :: a_w(n)



    call GAUSS( n, a_x, a_w )

    a_w = a_w * 2.0_DP

    ! Change integration domain from [-1,1] to [x1,x2]
    a_x = ( x2 - x1 ) * 0.5_DP * a_x + ( x1 + x2 ) * 0.5_DP
    a_w = a_w * ( x2 - x1 ) * 0.5_DP


  end subroutine GauLeg

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

  ! Subroutine below is provided by Ishioka-san.

  !**********************************************************************
  !     CALCULATE GAUSSIAN LATITUDES AND WEIGHTS               
  !**********************************************************************
  !     X(J): sin(\phi_j)
  !     W(J): w_j/2
  !**********************************************************************

  SUBROUTINE GAUSS(JM,X,W)

!!$    IMPLICIT REAL*8(A-H,O-Z)
!!$    PARAMETER(PI=3.1415926535897932385D0)
    integer, parameter :: NB=64
    integer, intent(in ) :: JM
!!$    DIMENSION X(JM),W(JM),E(NB)
    real(DP), intent(out) :: X(JM)
    real(DP), intent(out) :: W(JM)

    real(DP) :: E(NB)
    real(DP) :: EPS
    real(DP) :: Z
    real(DP) :: P0
    real(DP) :: P1
    real(DP) :: DPTMP
    real(DP) :: DZ
    integer  :: JH
    integer  :: IFLAG
    integer  :: I
    integer  :: J
    integer  :: N


    JH=JM/2

    EPS=1
    DO I=1,NB
      EPS=EPS/2
      E(I)=EPS+1
    END DO

    I=0
    EPS=1
10  CONTINUE
    I=I+1
    EPS=EPS/2
    IF(E(I).GT.1) GOTO 10

    EPS=EPS*16

    IF(MOD(JM,2).EQ.0) THEN
      DO J=1,JH
        Z=SIN(PI*(2*J-1)/(2*JM+1))
        IFLAG=0
20      CONTINUE
        P0=0
        P1=1
        DO N=1,JM-1,2
          P0=((2*N-1)*Z*P1-(N-1)*P0)/N
          P1=((2*N+1)*Z*P0-N*P1)/(N+1)
        END DO
        DPTMP=JM*(P0-Z*P1)/(1-Z*Z)
        DZ=P1/DPTMP
        Z=Z-DZ
        IF(IFLAG.EQ.0) THEN
          IF(ABS(DZ/Z).LE.EPS) THEN
            IFLAG=1
            X(JM-JH+J)=Z
          END IF
          GOTO 20
        END IF
        W(JM-JH+J)=1/(DPTMP*DPTMP)/(1-X(JM-JH+J)*X(JM-JH+J))
        W(JH+1-J)=W(JM-JH+J)
        X(JH+1-J)=-X(JM-JH+J)
      END DO
    ELSE
      DO J=1,JH
        Z=SIN(PI*2*J/(2*JM+1))
        IFLAG=0
30      CONTINUE
        P0=1
        P1=Z
        DO N=2,JM-1,2
          P0=((2*N-1)*Z*P1-(N-1)*P0)/N
          P1=((2*N+1)*Z*P0-N*P1)/(N+1)
        END DO
        DPTMP=JM*(P0-Z*P1)/(1-Z*Z)
        DZ=P1/DPTMP
        Z=Z-DZ
        IF(IFLAG.EQ.0) THEN
          IF(ABS(DZ/Z).LE.EPS) THEN
            IFLAG=1
            X(JM-JH+J)=Z
          END IF
          GOTO 30
        END IF
        W(JM-JH+J)=1/(DPTMP*DPTMP)/(1-X(JM-JH+J)*X(JM-JH+J))
        W(JH+1-J)=W(JM-JH+J)
        X(JH+1-J)=-X(JM-JH+J)
      END DO
      P0=1
      DO N=2,JM-1,2
        P0=-(N-1)*P0/N
      END DO
      DPTMP=JM*P0
      W(JH+1)=1/(DPTMP*DPTMP)
      X(JH+1)=0
    END IF

  END SUBROUTINE GAUSS

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

end module gauss_quad