ガウス重み, 分点の計算

Calculate Gauss node and Gaussian weight

Note that Japanese and English are described in parallel.

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

Methods

GAUSS GauLeg module_name version

Included Modules

dc_types constants0

Public Instance methods

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

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

[Source]

  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

Private Instance methods

GAUSS( JM, X, W )

Subroutine :

JM :

integer, intent(in )

:	!$ DIMENSION X(JM),W(JM),E(NB)

X(JM) : real(DP), intent(out)

W(JM) : real(DP), intent(out)

[Source]

  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

module_name

Constant :

module_name = ‘gauss_quad‘ :

character(*), parameter

:	モジュールの名称. Module name

version

Constant :

version = ’$Name: $’ // ’$Id: gauss_quad.f90,v 1.3 2011-06-19 11:05:23 yot Exp $’ :

character(*), parameter

:	モジュールのバージョン Module version