Class aq_module
In: libsrc/aq_module/aq_module.f90

aq_module

Authors:Shin-ichi Takehiro, Youhei SASAKI
Version:$Id: aq_module.f90 590 2013-08-19 08:48:21Z uwabami $
Copyright&License:See COPYRIGHT

概要

spml/aq_module モジュールは 1 次元有限領域の下での流体運動を Matsushima and Marcus (1994) で提唱された多項式を用いたスペクトル法によって 数値計算するための Fortran90 関数を提供する.

このルーチンは離散化にチェビシェフ—ガウス—ラダウ格子点を適用して おり, 主に 2 次元極座標, 円筒座標, 球座標の原点の特異性を回避しな がらスペクトル計算を行うために用いることを念頭においている.

2 次元データの 1 次元に関して同時にスペクトル計算を実行するための 関数も提供しており, 2, 3 次元領域での計算のベースも提供する.

Matsushima and Marcus (1994) の多項式に関する説明は 動径座標のスペクトル法(spectral_radial.pdf) を参照のこと.

関数・変数の名前と型について

命名法

  • 関数名の先頭 (u_, g_, aq_, ag_) は, 返す値の形を示している. 複数並んだ 2 次元データの第 1 次元の大きさは aq_Initial で 設定する重みの指数の配列 nd の大きさと同じでなければならない.
    q_ :スペクトルデータ
    g_ :1 次元格子点データ
    aq_ :1 次元スペクトルデータが複数並んだ 2 次元データ
    ag_ :1 次元格子点データが複数並んだ 2 次元データ.
  • 関数名の間の文字列(Dr)は, その関数の作用を表している.
  • 関数名の最後 (_e, _aq, _g, _ag) は, 入力変数の形がスペクトルデータ および格子点データであることを示している.
    _q :スペクトルデータ
    _g :1 次元格子点データ
    _aq :1 次元スペクトルデータが複数並んだ 2 次元データ
    _ag :1 次元格子点データが複数並んだ 2 次元データ

各データの種類の説明

  • g : 1 次元格子点データ.
    • 変数の種類と次元は real(8), dimension(im).
    • im は R 座標の格子点数であり, サブルーチン aq_Initial にて あらかじめ設定しておく.
  • q : スペクトルデータ.
    • 変数の種類と次元は real(8), dimension(0:km).
    • km は R 方向の最大波数であり, サブルーチン aq_Initial にて あらかじめ設定しておく. スペクトルデータの格納のされ方については…
  • ag : 1 次元(R)格子点データの並んだ 2 次元データ.
    • 変数の種類と次元は real(8), dimension(size(nd),im). 第 2 次元が R 方向を表す.
  • aq : 1 次元スペクトルデータの並んだ 2 次元データ.
    • 変数の種類と次元は real(8), dimension(size(nd),0:km). 第 2 次元がスペクトルを表す.
  • g_ で始まる関数が返す値は 1 次元格子点データに同じ.
  • q_ で始まる関数が返す値はスペクトルデータに同じ.
  • ag_ で始まる関数が返す値は 1 次元格子点データの並んだ 2 次元データに同じ.
  • aq_ で始まる関数が返す値は 1 次元スペクトルデータの並んだ 2 次元データに同じ.
  • スペクトルデータに対する微分等の作用とは, 対応する格子点データに 微分などを作用させたデータをスペクトル変換したもののことである.

変数・手続き群の要約

初期化

aq_Initial :スペクトル変換の格子点数, 波数, 領域の大きさの設定

座標変数

g_R :格子点座標(R)を格納した 1 次元配列
g_R_Weight :重み座標を格納した 1 次元配列

基本変換

g_q, ag_aq :スペクトルデータから格子データへの変換
q_g, aq_ag :格子データからスペクトルデータへの変換

微分

q_rDr_q, aq_rDr_aq :スペクトルデータに r(d/dR) 微分を作用させる

積

q_r2_q, aq_r2_aq :スペクトルデータに r^2 をかける
q_r2Inv_q, aq_r2Inv_aq :スペクトルデータに r^-2 をかける

積分・平均

a_Int_ag, a_Avr_ag :1 次元格子点データの並んだ 2 次元配列の積分および平均
Int_g, Avr_g :1 次元格子点データの積分および平均

境界値問題

aq_Boundary_D, aq_Boundary_N :外側ディリクレ条件, ノイマン条件
aq_BoundaryTau_D, aq_BoundaryTau_N :外側ディリクレ条件, ノイマン条件(タウ法)
ag_BoundaryGrid_D, ag_BoundaryGrid_N :外側ディリクレ条件, ノイマン条件(選点法)

Methods

Avr_g   Int_g   a_Avr_ag   a_Int_ag   ag_BoundaryGrid_D   ag_BoundaryGrid_D   ag_BoundaryGrid_N   ag_BoundaryGrid_N   ag_aq   aq_BoundaryTau_D   aq_BoundaryTau_D   aq_BoundaryTau_N   aq_BoundaryTau_N   aq_Boundary_D   aq_Boundary_D   aq_Boundary_N   aq_Boundary_N   aq_Initial   aq_ag   aq_r2Inv_aq   aq_r2_aq   aq_rDr_aq   g_R   g_R_Weight   g_q   q_g   q_r2Inv_q   q_r2_q   q_rDr_q  

Included Modules

dc_message lumatrix

Public Instance methods

Avr_g( g ) result(Avr_g)
Function :
Avr_g :real(8)
: (out) 積分値
g :real(8), dimension(im), intent(in)
: (in) 格子点データ

1 次元格子点データの平均 

     \int_0^a f(R) W(R) dR/\int_0^a W(R) dR,
         W(R) = R^beta/(a^2-R^2)^(1-alpha)

を計算する.

[Source]

  function Avr_g(g)
    !
    ! 1 次元格子点データの平均
    !
    !      \int_0^a f(R) W(R) dR/\int_0^a W(R) dR,
    !          W(R) = R^beta/(a^2-R^2)^(1-alpha)
    !
    ! を計算する.
    !
    real(8), dimension(im), intent(in)   :: g
    !(in) 格子点データ

    real(8)                              :: Avr_g
    !(out) 積分値

    Avr_g = Int_g(g)/sum(g_R_Weight)
  end function Avr_g
Int_g( g ) result(Int_g)
Function :
Int_g :real(8)
: (out) 積分値
g :real(8), dimension(im), intent(in)
: (in) 格子点データ

1 次元格子点データの積分 

     \int_0^a f(R) W(R) dR, W(R) = R^beta/(a^2-R^2)^(1-alpha)

を行う

[Source]

  function Int_g(g)
    !
    ! 1 次元格子点データの積分
    !
    !      \int_0^a f(R) W(R) dR, W(R) = R^beta/(a^2-R^2)^(1-alpha)
    !
    ! を行う
    !
    real(8), dimension(im), intent(in)   :: g
    !(in) 格子点データ

    real(8)                              :: Int_g
    !(out) 積分値

    Int_g = sum(g*g_R_Weight)

  end function Int_g
a_Avr_ag( ag ) result(a_Avr_ag)
Function :
a_Avr_ag :real(8), dimension(size(ag,1))
: (out) 平均したデータ
ag :real(8), dimension(:,0:), intent(in)
: (in)入力格子点データ

1 次元格子点データが並んだ 2 次元配列の平均 

     \int_0^a f(R) W(R) dR/\int_0^a W(R) dR,
         W(R) = R^beta/(a^2-R^2)^(1-alpha)

を計算する.

[Source]

  function a_Avr_ag(ag)
    !
    ! 1 次元格子点データが並んだ 2 次元配列の平均
    !
    !      \int_0^a f(R) W(R) dR/\int_0^a W(R) dR,
    !          W(R) = R^beta/(a^2-R^2)^(1-alpha)
    !
    ! を計算する.
    !
    real(8), dimension(:,0:), intent(in)   :: ag
    !(in)入力格子点データ

    real(8), dimension(size(ag,1))         :: a_Avr_ag
    !(out) 平均したデータ

    a_Avr_ag = a_Int_ag(ag)/sum(g_R_Weight)

  end function a_Avr_ag
a_Int_ag( ag ) result(a_Int_ag)
Function :
a_Int_ag :real(8), dimension(size(ag,1))
: (out) 積分したデータ
ag :real(8), dimension(:,:), intent(in)
: (in)入力格子点データ

1 次元格子点データが並んだ 2 次元配列の積分 

     \int_0^a f(R) W(R) dR, W(R) = R^beta/(a^2-R^2)^(1-alpha)

を計算する.

[Source]

  function a_Int_ag(ag)
    !
    ! 1 次元格子点データが並んだ 2 次元配列の積分
    !
    !      \int_0^a f(R) W(R) dR, W(R) = R^beta/(a^2-R^2)^(1-alpha)
    !
    ! を計算する.
    !
    real(8), dimension(:,:), intent(in)      :: ag
    !(in)入力格子点データ

    real(8), dimension(size(ag,1))           :: a_Int_ag
    !(out) 積分したデータ
    integer :: i

    if ( size(ag,2) < im ) then
       call MessageNotify('E','a_Int_ag', 'The Grid points of input data too small.')
    elseif ( size(ag,2) > im ) then
       call MessageNotify('W','a_Int_ag', 'The Grid points of input data too large.')
    endif

    a_Int_ag = 0.0d0
    do i=1,im
       a_Int_ag(:) = a_Int_ag(:) + ag(:,i)*g_R_Weight(i)
    enddo
  end function a_Int_ag
ag_BoundaryGrid_D( ag_data, [value] )
Subroutine :
ag_data :real(8), dimension(:,:),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(jmax,im)
value :real(8), dimension(:), intent(in), optional
: (in) 境界値(jmax)

Dirichlet 型境界条件の適用(選点法, 2 次元配列用)

  • 外側境界(i=im)での値を与える.

[Source]

  subroutine ag_BoundaryGrid_D_2d(ag_data,value)
    !
    ! Dirichlet 型境界条件の適用(選点法, 2 次元配列用)
    ! * 外側境界(i=im)での値を与える.
    !
    real(8), dimension(:,:),intent(inout)         :: ag_data
    !(inout) 境界条件を適用するチェビシェフデータ(jmax,im)

    real(8), dimension(:), intent(in), optional  :: value
    !(in) 境界値(jmax)

    real(8), dimension(size(ag_data,1))       :: value0

    if (.not. present(value)) then
       value0=0.0d0
    else
       value0 = value
    endif
    ag_data(:,im) = value0

  end subroutine ag_BoundaryGrid_D_2d
ag_BoundaryGrid_D( g_data, [value] )
Subroutine :
g_data :real(8), dimension(im),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(im)
value :real(8), intent(in), optional
: (in) 境界値

Dirichlet 型境界条件の適用(タウ法, 1 次元配列用)

  • 両境界での値を与える.

[Source]

  subroutine ag_BoundaryGrid_D_1d(g_data,value)
    !
    ! Dirichlet 型境界条件の適用(タウ法, 1 次元配列用)
    ! * 両境界での値を与える.
    !
    real(8), dimension(im),intent(inout)       :: g_data
    !(inout) 境界条件を適用するチェビシェフデータ(im)

    real(8), intent(in), optional              :: value
    !(in) 境界値

    real(8), dimension(1,im)                   :: ag_work
    real(8), dimension(1)                      :: vwork

    if (.not. present(value)) then
       vwork(1)=0.0d0
    else
       vwork(1) = value
    endif

    ag_work(1,:)=g_data
    call ag_BoundaryGrid_D_2d(ag_work,vwork)
    g_data=ag_work(1,:)

  end subroutine ag_BoundaryGrid_D_1d
ag_BoundaryGrid_N( ag_data, [value] )
Subroutine :
ag_data :real(8), dimension(:,:),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(m,im)
value :real(8), dimension(:), intent(in), optional
: (in) 境界値(m)

外側境界 Neumann 型境界条件の適用(選点法, 2 次元配列用)

  • i=im で勾配の値を与える.

[Source]

  subroutine ag_BoundaryGrid_N_2d(ag_data,value)
    !
    ! 外側境界 Neumann 型境界条件の適用(選点法, 2 次元配列用)
    ! * i=im で勾配の値を与える.
    !
    real(8), dimension(:,:),intent(inout)         :: ag_data
    !(inout) 境界条件を適用するチェビシェフデータ(m,im)

    real(8), dimension(:), intent(in), optional    :: value
    !(in) 境界値(m)

    real(8), dimension(:,:,:), allocatable    :: alu
    integer, dimension(:,:), allocatable      :: kp
    real(8), dimension(size(ag_data,1),im)    :: ag_work
    real(8), dimension(size(ag_data,1))       :: value0           ! 境界値

    integer :: i
    save    :: alu, kp

    if ( size(ag_data,2) < im ) then
       call MessageNotify('E','aq_BoundaryGrid_N', 'The dimension of input data too small.')
    elseif ( size(ag_data,2) > im ) then
       call MessageNotify('W','aq_BoundaryGrid_N', 'The dimension of input data too large.')
    endif

    if (.not. present(value)) then
       value0=0.0D0
    else
       value0 = value
    endif

    if ( first_Grid_N ) then
       first_Grid_N = .false.
       allocate(alu(jmax,im,im),kp(jmax,im))

       alu=0.0D0
       do i=1,im
          alu(:,i,i) = 1.0D0
       enddo

       do i=1,im
          ag_work = 0.0D0
          ag_work(:,i) = 1.0D0
          ag_work = ag_aq(aq_rDr_aq(aq_ag(ag_work)))/spread(g_R,1,jmax)
          alu(:,im,i) = ag_work(:,im)
       enddo

       call ludecomp(alu,kp)
    endif

    ag_data(:,im) = value0

    ag_data = lusolve(alu,kp,ag_data)

  end subroutine ag_BoundaryGrid_N_2d
ag_BoundaryGrid_N( g_data, [value] )
Subroutine :
g_data :real(8), dimension(im),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(0:km)
value :real(8), intent(in), optional
: (in) 境界値

Dirichlet/Neumann 型境界条件の適用(選点法, 1 次元配列用)

  • i=0 で勾配の値を与える.

[Source]

  subroutine ag_BoundaryGrid_N_1d(g_data,value)
    !
    ! Dirichlet/Neumann 型境界条件の適用(選点法, 1 次元配列用)
    ! * i=0 で勾配の値を与える.
    !
    real(8), dimension(im),intent(inout)       :: g_data
    !(inout) 境界条件を適用するチェビシェフデータ(0:km)

    real(8), intent(in), optional                :: value
    !(in) 境界値

    real(8), dimension(1,im)                   :: ag_work
    real(8), dimension(1)                        :: vwork           ! 境界値

    if (.not. present(value)) then
       vwork(1)=0.0D0
    else
       vwork(1) = value
    endif

    ag_work(1,:)=g_data
    call ag_BoundaryGrid_N_2d(ag_work,vwork)
    g_data=ag_work(1,:)

  end subroutine ag_BoundaryGrid_N_1d
ag_aq( aq_data ) result(ag_aq)
Function :
ag_aq :real(8), dimension(size(aq_data,1),im)
: (out) 格子点データ
aq_data :real(8), dimension(:,0:), intent(in)
: (in) スペクトルデータ

スペクトルデータから格子データへ変換する(2 次元配列用).

[Source]

  function ag_aq(aq_data)
    !
    ! スペクトルデータから格子データへ変換する(2 次元配列用).
    !
    real(8), dimension(:,0:), intent(in)    :: aq_data
    !(in) スペクトルデータ

    real(8), dimension(size(aq_data,1),im)  :: ag_aq
    !(out) 格子点データ

    integer :: i, j

    if ( size(aq_data,1) /= jmax ) then
       call MessageNotify('E','ag_aq', '1st dim. of the spectral data should be same as dim. of MD.')
    end if

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','ag_aq', 'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','ag_aq', 'The spectral dimension of input data too large.')
    endif

    ag_aq = 0.0D0
    do i=1,im
!$omp parallel do
       do j=1,jmax
             ag_aq(j,i) = sum(CB(j,i,md(j):km:2)*aq_data(j,md(j):km:2))
       enddo
!$omp end parallel do
    enddo

  end function ag_aq
aq_BoundaryTau_D( aq_data, [value] )
Subroutine :
aq_data :real(8), dimension(:,0:),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(jmax,0:km)
value :real(8), dimension(:), intent(in), optional
: (in) 境界値(jmax)

Dirichlet 型境界条件の適用(タウ法, 2 次元配列用)

  • 外側境界(i=im)での値を与える.

[Source]

  subroutine aq_BoundaryTau_D_2d(aq_data,value)
    !
    ! Dirichlet 型境界条件の適用(タウ法, 2 次元配列用)
    ! * 外側境界(i=im)での値を与える.
    !
    real(8), dimension(:,0:),intent(inout)         :: aq_data
    !(inout) 境界条件を適用するチェビシェフデータ(jmax,0:km)

    real(8), dimension(:), intent(in), optional  :: value
    !(in) 境界値(jmax)

    real(8), dimension(:,:,:), allocatable  :: alu
    integer, dimension(:,:), allocatable    :: kp
    real(8), dimension(size(aq_data,1),0:km):: aq_work
    real(8), dimension(size(aq_data,1),im)  :: ag_work
    real(8), dimension(size(aq_data,1))     :: value0           ! 境界値

    integer :: k, j, kstr, kend
    save    :: alu, kp

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_BoundaryTau_D', 'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_BoundaryTau_D', 'The spectral dimension of input data too large.')
    endif

    if (.not. present(value)) then
       value0=0.0D0
    else
       value0 = value
    endif

    if ( first_Tau_D ) then
       first_Tau_D = .false.

       allocate(alu(size(aq_data,1),0:km,0:km),kp(size(aq_data,1),0:km))

       alu=0.0D0
       do k=0,km
          alu(:,k,k) = 1.0D0
       enddo

       do j=1,jmax
          if ( mod(md(j),2) .eq. mod(km,2) ) then
             kstr = md(j) 
              kend = km
          else
             kstr = md(j) 
              kend = km-1
          endif

          do k=kstr,kend,2
             aq_work = 0.0D0
             aq_work(j,k) = 1.0D0
             ag_work = ag_aq(aq_work)
             alu(j,kend,k) = ag_work(j,im)
          enddo
       enddo

       call ludecomp(alu,kp)
    endif

    do j=1,jmax
       if ( mod(md(j),2) .eq. mod(km,2) ) then
          aq_data(j,km)   = value0(j)
       else
          aq_data(j,km-1)   = value0(j)
       endif
    enddo

    aq_data = lusolve(alu,kp,aq_data)

  end subroutine aq_BoundaryTau_D_2d
aq_BoundaryTau_D( q_data, [value] )
Subroutine :
q_data :real(8), dimension(0:km),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(0:km)
value :real(8), intent(in), optional
: (in) 境界値

Dirichlet 型境界条件の適用(タウ法, 1 次元配列用)

  • 両境界での値を与える.

[Source]

  subroutine aq_BoundaryTau_D_1d(q_data,value)
    !
    ! Dirichlet 型境界条件の適用(タウ法, 1 次元配列用)
    ! * 両境界での値を与える.
    !
    real(8), dimension(0:km),intent(inout)       :: q_data
    !(inout) 境界条件を適用するチェビシェフデータ(0:km)

    real(8), intent(in), optional                :: value
    !(in) 境界値

    real(8), dimension(1,0:km)                   :: aq_work
    real(8), dimension(1)                        :: vwork  ! 境界値

    if (.not. present(value)) then
       vwork(1)=0.0D0
    else
       vwork(1) = value
    endif

    aq_work(1,:)=q_data
    call aq_BoundaryTau_D_2d(aq_work,vwork)
    q_data=aq_work(1,:)

  end subroutine aq_BoundaryTau_D_1d
aq_BoundaryTau_N( aq_data, [value] )
Subroutine :
aq_data :real(8), dimension(:,0:),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(m,0:km)
value :real(8), dimension(:), intent(in), optional
: (in) 境界値(m)

外側境界 Neumann 型境界条件の適用(タウ法, 2 次元配列用)

  • i=im で勾配の値を与える.

[Source]

  subroutine aq_BoundaryTau_N_2d(aq_data,value)
    !
    ! 外側境界 Neumann 型境界条件の適用(タウ法, 2 次元配列用)
    ! * i=im で勾配の値を与える.
    !
    real(8), dimension(:,0:),intent(inout)         :: aq_data
    !(inout) 境界条件を適用するチェビシェフデータ(m,0:km)

    real(8), dimension(:), intent(in), optional    :: value
    !(in) 境界値(m)

    real(8), dimension(:,:,:), allocatable    :: alu
    integer, dimension(:,:), allocatable      :: kp
    real(8), dimension(size(aq_data,1),0:km)  :: aq_work
    real(8), dimension(size(aq_data,1),im)    :: ag_work
    real(8), dimension(size(aq_data,1))       :: value0           ! 境界値

    integer :: j, k, kstr, kend

    save    :: alu, kp

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_BoundaryTau_N', 'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_BoundaryTau_N', 'The spectral dimension of input data too large.')
    endif

    if (.not. present(value)) then
       value0=0.0D0
    else
       value0 = value
    endif

    if ( first_Tau_N ) then
       first_Tau_N = .false.
       allocate(alu(jmax,0:km,0:km),kp(jmax,0:km))

       alu=0.0D0
       do k=0,km
          alu(:,k,k) = 1.0D0
       enddo

       do j=1,jmax
          if ( mod(md(j),2) .eq. mod(km,2) ) then
             kstr=md(j) 
              kend = km
          else
             kstr=md(j) 
              kend = km-1
          endif

          do k=kstr,kend,2
             aq_work = 0.0D0
             aq_work(j,k) = 1.0D0
             ag_work = ag_aq(aq_rDr_aq(aq_work))/spread(g_R,1,jmax)
             alu(j,kend,k) = ag_work(j,im)
          enddo
       enddo

       call ludecomp(alu,kp)
    endif

    do j=1,jmax
       if ( mod(md(j),2) .eq. mod(km,2) ) then
          aq_data(j,km)   = value0(j)
       else
          aq_data(j,km-1)   = value0(j)
       endif
    enddo

    aq_data = lusolve(alu,kp,aq_data)

  end subroutine aq_BoundaryTau_N_2d
aq_BoundaryTau_N( q_data, [value] )
Subroutine :
q_data :real(8), dimension(0:km),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(0:km)
value :real(8), intent(in), optional
: (in) 境界値

Dirichlet/Neumann 型境界条件の適用(タウ法, 1 次元配列用)

  • i=0 で勾配の値を与える.

[Source]

  subroutine aq_BoundaryTau_N_1d(q_data,value)
    !
    ! Dirichlet/Neumann 型境界条件の適用(タウ法, 1 次元配列用)
    ! * i=0 で勾配の値を与える.
    !
    real(8), dimension(0:km),intent(inout)       :: q_data
    !(inout) 境界条件を適用するチェビシェフデータ(0:km)

    real(8), intent(in), optional                :: value
    !(in) 境界値

    real(8), dimension(1,0:km)                   :: aq_work
    real(8), dimension(1)                        :: vwork           ! 境界値

    if (.not. present(value)) then
       vwork(1)=0.0D0
    else
       vwork(1) = value
    endif

    aq_work(1,:)=q_data
    call aq_BoundaryTau_N_2d(aq_work,vwork)
    q_data=aq_work(1,:)

  end subroutine aq_BoundaryTau_N_1d
aq_Boundary_D( aq_data, [value] )
Subroutine :
aq_data :real(8), dimension(:,0:),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(jmax,0:km)
value :real(8), dimension(:), intent(in), optional
: (in) 境界値(jmax)

Dirichlet 型境界条件の適用(タウ法, 2 次元配列用)

  • 外側境界(i=im)での値を与える.

[Source]

  subroutine aq_BoundaryTau_D_2d(aq_data,value)
    !
    ! Dirichlet 型境界条件の適用(タウ法, 2 次元配列用)
    ! * 外側境界(i=im)での値を与える.
    !
    real(8), dimension(:,0:),intent(inout)         :: aq_data
    !(inout) 境界条件を適用するチェビシェフデータ(jmax,0:km)

    real(8), dimension(:), intent(in), optional  :: value
    !(in) 境界値(jmax)

    real(8), dimension(:,:,:), allocatable  :: alu
    integer, dimension(:,:), allocatable    :: kp
    real(8), dimension(size(aq_data,1),0:km):: aq_work
    real(8), dimension(size(aq_data,1),im)  :: ag_work
    real(8), dimension(size(aq_data,1))     :: value0           ! 境界値

    integer :: k, j, kstr, kend
    save    :: alu, kp

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_BoundaryTau_D', 'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_BoundaryTau_D', 'The spectral dimension of input data too large.')
    endif

    if (.not. present(value)) then
       value0=0.0D0
    else
       value0 = value
    endif

    if ( first_Tau_D ) then
       first_Tau_D = .false.

       allocate(alu(size(aq_data,1),0:km,0:km),kp(size(aq_data,1),0:km))

       alu=0.0D0
       do k=0,km
          alu(:,k,k) = 1.0D0
       enddo

       do j=1,jmax
          if ( mod(md(j),2) .eq. mod(km,2) ) then
             kstr = md(j) 
              kend = km
          else
             kstr = md(j) 
              kend = km-1
          endif

          do k=kstr,kend,2
             aq_work = 0.0D0
             aq_work(j,k) = 1.0D0
             ag_work = ag_aq(aq_work)
             alu(j,kend,k) = ag_work(j,im)
          enddo
       enddo

       call ludecomp(alu,kp)
    endif

    do j=1,jmax
       if ( mod(md(j),2) .eq. mod(km,2) ) then
          aq_data(j,km)   = value0(j)
       else
          aq_data(j,km-1)   = value0(j)
       endif
    enddo

    aq_data = lusolve(alu,kp,aq_data)

  end subroutine aq_BoundaryTau_D_2d
aq_Boundary_D( q_data, [value] )
Subroutine :
q_data :real(8), dimension(0:km),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(0:km)
value :real(8), intent(in), optional
: (in) 境界値

Dirichlet 型境界条件の適用(タウ法, 1 次元配列用)

  • 両境界での値を与える.

[Source]

  subroutine aq_BoundaryTau_D_1d(q_data,value)
    !
    ! Dirichlet 型境界条件の適用(タウ法, 1 次元配列用)
    ! * 両境界での値を与える.
    !
    real(8), dimension(0:km),intent(inout)       :: q_data
    !(inout) 境界条件を適用するチェビシェフデータ(0:km)

    real(8), intent(in), optional                :: value
    !(in) 境界値

    real(8), dimension(1,0:km)                   :: aq_work
    real(8), dimension(1)                        :: vwork  ! 境界値

    if (.not. present(value)) then
       vwork(1)=0.0D0
    else
       vwork(1) = value
    endif

    aq_work(1,:)=q_data
    call aq_BoundaryTau_D_2d(aq_work,vwork)
    q_data=aq_work(1,:)

  end subroutine aq_BoundaryTau_D_1d
aq_Boundary_N( aq_data, [value] )
Subroutine :
aq_data :real(8), dimension(:,0:),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(m,0:km)
value :real(8), dimension(:), intent(in), optional
: (in) 境界値(m)

外側境界 Neumann 型境界条件の適用(タウ法, 2 次元配列用)

  • i=im で勾配の値を与える.

[Source]

  subroutine aq_BoundaryTau_N_2d(aq_data,value)
    !
    ! 外側境界 Neumann 型境界条件の適用(タウ法, 2 次元配列用)
    ! * i=im で勾配の値を与える.
    !
    real(8), dimension(:,0:),intent(inout)         :: aq_data
    !(inout) 境界条件を適用するチェビシェフデータ(m,0:km)

    real(8), dimension(:), intent(in), optional    :: value
    !(in) 境界値(m)

    real(8), dimension(:,:,:), allocatable    :: alu
    integer, dimension(:,:), allocatable      :: kp
    real(8), dimension(size(aq_data,1),0:km)  :: aq_work
    real(8), dimension(size(aq_data,1),im)    :: ag_work
    real(8), dimension(size(aq_data,1))       :: value0           ! 境界値

    integer :: j, k, kstr, kend

    save    :: alu, kp

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_BoundaryTau_N', 'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_BoundaryTau_N', 'The spectral dimension of input data too large.')
    endif

    if (.not. present(value)) then
       value0=0.0D0
    else
       value0 = value
    endif

    if ( first_Tau_N ) then
       first_Tau_N = .false.
       allocate(alu(jmax,0:km,0:km),kp(jmax,0:km))

       alu=0.0D0
       do k=0,km
          alu(:,k,k) = 1.0D0
       enddo

       do j=1,jmax
          if ( mod(md(j),2) .eq. mod(km,2) ) then
             kstr=md(j) 
              kend = km
          else
             kstr=md(j) 
              kend = km-1
          endif

          do k=kstr,kend,2
             aq_work = 0.0D0
             aq_work(j,k) = 1.0D0
             ag_work = ag_aq(aq_rDr_aq(aq_work))/spread(g_R,1,jmax)
             alu(j,kend,k) = ag_work(j,im)
          enddo
       enddo

       call ludecomp(alu,kp)
    endif

    do j=1,jmax
       if ( mod(md(j),2) .eq. mod(km,2) ) then
          aq_data(j,km)   = value0(j)
       else
          aq_data(j,km-1)   = value0(j)
       endif
    enddo

    aq_data = lusolve(alu,kp,aq_data)

  end subroutine aq_BoundaryTau_N_2d
aq_Boundary_N( q_data, [value] )
Subroutine :
q_data :real(8), dimension(0:km),intent(inout)
: (inout) 境界条件を適用するチェビシェフデータ(0:km)
value :real(8), intent(in), optional
: (in) 境界値

Dirichlet/Neumann 型境界条件の適用(タウ法, 1 次元配列用)

  • i=0 で勾配の値を与える.

[Source]

  subroutine aq_BoundaryTau_N_1d(q_data,value)
    !
    ! Dirichlet/Neumann 型境界条件の適用(タウ法, 1 次元配列用)
    ! * i=0 で勾配の値を与える.
    !
    real(8), dimension(0:km),intent(inout)       :: q_data
    !(inout) 境界条件を適用するチェビシェフデータ(0:km)

    real(8), intent(in), optional                :: value
    !(in) 境界値

    real(8), dimension(1,0:km)                   :: aq_work
    real(8), dimension(1)                        :: vwork           ! 境界値

    if (.not. present(value)) then
       vwork(1)=0.0D0
    else
       vwork(1) = value
    endif

    aq_work(1,:)=q_data
    call aq_BoundaryTau_N_2d(aq_work,vwork)
    q_data=aq_work(1,:)

  end subroutine aq_BoundaryTau_N_1d
aq_Initial( i_in, k_in, r_in, alpha_in, beta_in, md_in )
Subroutine :
i_in :integer,intent(in)
: (in) 格子点数
k_in :integer,intent(in)
: (in) 切断波数
r_in :real(8),intent(in)
: (in) 外側境界の座標(半径)
alpha_in :real(8),intent(in)
: (in) 展開多項式パラメター
beta_in :real(8),intent(in)
: (in) 展開多項式パラメター
md_in(:) :integer,intent(in)
: (in) 展開多項式上付次数のならび

スペクトル変換の格子点数, 波数, 領域の大きさ, 重みを設定する.

他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで 初期設定をしなければならない.

[Source]

  subroutine aq_Initial(i_in,k_in,r_in,alpha_in,beta_in,md_in)
    !
    ! スペクトル変換の格子点数, 波数, 領域の大きさ, 重みを設定する.
    !
    ! 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで
    ! 初期設定をしなければならない.
    !
    integer,intent(in) :: i_in              !(in) 格子点数
    integer,intent(in) :: k_in              !(in) 切断波数
    real(8),intent(in) :: r_in              !(in) 外側境界の座標(半径)
    real(8),intent(in) :: alpha_in          !(in) 展開多項式パラメター
    real(8),intent(in) :: beta_in           !(in) 展開多項式パラメター
    integer,intent(in) :: md_in(:)          !(in) 展開多項式上付次数のならび

    integer i, j, n

   !--- パラメターのチェック・記憶
    im=i_in 
     km=k_in 
     ra = r_in
     alpha=alpha_in
     beta=beta_in

    if ( km .ge. 2*im ) then
       call MessageNotify('E','aq_initial','KM shoud be less than 2*IM')
    endif

    if ( alpha .le. 0.0D0 ) then
       call MessageNotify('E','aq_initial','alpha must be larger than 0')
    endif
    if ( alpha .gt. 1.0D0 ) then
       call MessageNotify('E','aq_initial','alpha must be smaller equal to 1')
    endif
    if ( beta .le. 0.0D0 ) then
       call MessageNotify('E','aq_initial','beta must be larger than 0')
    endif

    gamma = 2.0D0 * alpha + beta

   !--- 次数情報の計算
    jmax = size(md_in)
    if ( allocated(md) ) deallocate(md)
    allocate(md(jmax))

    md = md_in

   !--- 格子点, 重みの設定
    if ( allocated(g_R) ) deallocate(g_R)
    if ( allocated(g_R_Weight) ) deallocate(g_R_Weight)
    allocate(g_R(im),g_R_Weight(im))

    call gauss_radau(2*im, g_R, g_R_Weight)
    g_R = ra * g_R
    g_R_Weight = ra**(gamma-1) * g_R_Weight

   !--- 変換行列の設定
    if ( allocated(CF) ) deallocate(CF)
    if ( allocated(CB) ) deallocate(CB)
    allocate(CF(jmax,0:km,im),CB(jmax,im,0:km))

    CF = 0.0D0
    do j=1,jmax
       do n=md(j),km,2
          do i=1,im
             CF(j,n,i) = Phi(g_R(i)/ra,n,md(j)) * g_R_Weight(i)/ra**(gamma-1)
          enddo
       enddo
    enddo

    CB = 0.0D0
    do j=1,jmax
       do i=1,im
          do n=md(j),km,2
             CB(j,i,n) = Phi(g_R(i)/ra,n,md(j))
          enddo
       enddo
    enddo

   !--- 各ルーチン変換行列の初期化スイッチ

    first_r2inv  = .true.
    first_Tau_D  = .true.
    first_Tau_N  = .true.
    first_Grid_N = .true.

    call MessageNotify( 'M','aq_initial','aq_module (2009/07/31) is initialized')

  end subroutine aq_Initial
aq_ag( ag_data ) result(aq_ag)
Function :
aq_ag :real(8), dimension(size(ag_data,1),0:km)
: (out) スペクトルデータ
ag_data :real(8), dimension(:,:), intent(in)
: (in) 格子点データ

格子データからスペクトルデータへ変換する(2 次元配列用).

[Source]

  function aq_ag(ag_data)
    !
    ! 格子データからスペクトルデータへ変換する(2 次元配列用).
    !
    real(8), dimension(:,:), intent(in)      :: ag_data
    !(in) 格子点データ

    real(8), dimension(size(ag_data,1),0:km) :: aq_ag
    !(out) スペクトルデータ

    integer :: j, n

    if ( size(ag_data,1) /= jmax ) then
       call MessageNotify('E','aq_ag', '1st dim. of the grid data should be same as dim. of MD.')
    end if

    if ( size(ag_data,2) < im ) then
       call MessageNotify('E','aq_ag', 'The Grid points of input data too small.')
    elseif ( size(ag_data,2) > im ) then
       call MessageNotify('W','aq_ag', 'The Grid points of input data too large.')
    endif

    aq_ag = 0.0D0
!$omp parallel do private(n)
    do j=1,jmax
       do n=md(j),km,2
          aq_ag(j,n) = sum(CF(j,n,:)*ag_data(j,:))
       enddo
    enddo
!$omp end parallel do

  end function aq_ag
aq_r2Inv_aq( aq_data ) result(aq_r2Inv_aq)
Function :
aq_r2Inv_aq :real(8), dimension(size(aq_data,1),0:size(aq_data,2)-1)
: (out) 出力スペクトルデータ
aq_data :real(8), dimension(:,0:), intent(in)
: (in) 入力スペクトルデータ

入力スペクトルデータに対して積 r^-2 のスペクトル係数 を計算する(2 次元配列用). 

  a_n^m = aq_r2Inv_aq/sqrt(Inm),  b_n^m = aq_data/sqrt(Inm)

  b_n^m = (n-|m|)(n+|m|+beta-1)/((2n+gamma-5)(2n+gamma-3)) a_n-2^m
      + (2n(n+gamma-1) + 2|m|(|m|+beta-1)+(gamma-3)(beta+1)
           /((2n+gamma+1)(2n+gamma-3)) a_n^m
      + (n-|m|+gamma-beta)(n+|m|+gamma-1)
           /((2n+gamma+3)(2n+gamma+1)) a_n+2^m

[Source]

  function aq_r2Inv_aq(aq_data)
    !
    ! 入力スペクトルデータに対して積 r^-2 のスペクトル係数
    ! を計算する(2 次元配列用).
    !
    !   a_n^m = aq_r2Inv_aq/sqrt(Inm),  b_n^m = aq_data/sqrt(Inm)
    !
    !   b_n^m = (n-|m|)(n+|m|+beta-1)/((2n+gamma-5)(2n+gamma-3)) a_n-2^m
    !       + (2n(n+gamma-1) + 2|m|(|m|+beta-1)+(gamma-3)(beta+1)
    !            /((2n+gamma+1)(2n+gamma-3)) a_n^m
    !       + (n-|m|+gamma-beta)(n+|m|+gamma-1)
    !            /((2n+gamma+3)(2n+gamma+1)) a_n+2^m
    !
    real(8), dimension(:,0:), intent(in)     :: aq_data
    !(in) 入力スペクトルデータ

    real(8), dimension(size(aq_data,1),0:size(aq_data,2)-1) :: aq_r2Inv_aq
    !(out) 出力スペクトルデータ

    real(8), dimension(:,:,:), allocatable    :: R2MTX
    integer, dimension(:,:), allocatable      :: kp

    integer :: j, n, m
    integer :: nstr, nend
    real(8) :: sqrInp2m, sqrInm, sqrInm2m

    save R2MTX, kp

    if ( size(aq_data,1) /= jmax ) then
       call MessageNotify('E','aq_r2Inv_aq', '1st dim. of the spectral data should be same as dim. of MD.')
    end if

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_r2Inv_aq', 'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_r2Inv_aq', 'The spectral dimension of input data too large.')
    endif

    if ( first_r2inv ) then
       first_r2inv = .false.
       if ( allocated(R2MTX) ) deallocate(R2MTX)
       if ( allocated(kp) ) deallocate(kp)
       allocate(R2MTX(jmax,0:km,0:km),kp(jmax,0:km))

       R2MTX(:,:,:) = 0.0D0
       do n=0,km
          R2MTX(:,n,n) = 1.0D0
       enddo

!$omp parallel do private(nstr,nend,sqrInm,sqrInm2m,sqrInp2m,n,m)
       do j=1,jmax
          if ( mod(md(j),2) .eq. mod(km,2) ) then
             nstr=md(j) 
              nend=km
          else
             nstr=md(j) 
              nend=km-1
          endif

          m = abs(md(j))

          n = nstr

          sqrInm   = sqrt(Inm(n,m))
          sqrInp2m = sqrt(Inm(n+2,m))

          R2MTX(j,n,n)   = (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) /((2*n+gamma+1)*(2*n+gamma-3))

          R2MTX(j,n,n+2) = (n-m+gamma-beta)*(n+m+gamma-1) /((2*n+gamma+3)*(2*n+gamma+1)) * (sqrInm/sqrInp2m)

          do n=nstr+2,nend-2,2
             sqrInm2m = sqrInm
             sqrInm   = sqrInp2m
             sqrInp2m = sqrt(Inm(n+2,m))

             R2MTX(j,n,n-2) = (n-m)*(n+m+beta-1)/((2*n+gamma-5)*(2*n+gamma-3)) *(sqrInm/sqrInm2m)

             R2MTX(j,n,n)   = (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) /((2*n+gamma+1)*(2*n+gamma-3))

             R2MTX(j,n,n+2) = (n-m+gamma-beta)*(n+m+gamma-1) /((2*n+gamma+3)*(2*n+gamma+1)) * (sqrInm/sqrInp2m)
          enddo

          n = nend
          sqrInm2m = sqrInm
          sqrInm   = sqrInp2m
             R2MTX(j,n,n-2) = (n-m)*(n+m+beta-1)/((2*n+gamma-5)*(2*n+gamma-3)) *(sqrInm/sqrInm2m)
             R2MTX(j,n,n)   = (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) /((2*n+gamma+1)*(2*n+gamma-3))
       enddo
!$omp end parallel do

       R2MTX = R2MTX * ra**2

       call LuDecomp(R2MTX,kp)

    end if

    aq_r2Inv_aq = LuSolve(R2MTX,kp,aq_data)

  end function aq_r2Inv_aq
aq_r2_aq( aq_data ) result(aq_r2_aq)
Function :
aq_r2_aq :real(8), dimension(size(aq_data,1),0:size(aq_data,2)-1)
: (out) 出力スペクトルデータ
aq_data :real(8), dimension(:,0:), intent(in)
: (in) 入力スペクトルデータ

入力スペクトルデータに対して積 r^2 のスペクトル係数 を計算する(2 次元配列用). 

  a_n^m = aq_data/sqrt(Inm),  b_n^m = aq_rDr_aq/sqrt(Inm)

  b_n^m = (n-|m|)(n+|m|+beta-1)/((2n+gamma-5)(2n+gamma-3)) a_n-2^m
      + (2n(n+gamma-1) + 2|m|(|m|+beta-1)+(gamma-3)(beta+1)
           /((2n+gamma+1)(2n+gamma-3)) a_n^m
      + (n-|m|+gamma-beta)(n+|m|+gamma-1)
           /((2n+gamma+3)(2n+gamma+1)) a_n+2^m

[Source]

  function aq_r2_aq(aq_data)
    !
    ! 入力スペクトルデータに対して積 r^2 のスペクトル係数
    ! を計算する(2 次元配列用).
    !
    !   a_n^m = aq_data/sqrt(Inm),  b_n^m = aq_rDr_aq/sqrt(Inm)
    !
    !   b_n^m = (n-|m|)(n+|m|+beta-1)/((2n+gamma-5)(2n+gamma-3)) a_n-2^m
    !       + (2n(n+gamma-1) + 2|m|(|m|+beta-1)+(gamma-3)(beta+1)
    !            /((2n+gamma+1)(2n+gamma-3)) a_n^m
    !       + (n-|m|+gamma-beta)(n+|m|+gamma-1)
    !            /((2n+gamma+3)(2n+gamma+1)) a_n+2^m
    !
    real(8), dimension(:,0:), intent(in)     :: aq_data
    !(in) 入力スペクトルデータ

    real(8), dimension(size(aq_data,1),0:size(aq_data,2)-1) :: aq_r2_aq
    !(out) 出力スペクトルデータ

    integer :: j, n, m
    integer :: nstr, nend
    real(8) :: sqrInp2m, sqrInm, sqrInm2m

    if ( size(aq_data,1) /= jmax ) then
       call MessageNotify('E','aq_r2_aq', '1st dim. of the spectral data should be same as dim. of MD.')
    end if

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_r2_aq', 'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_r2_aq', 'The spectral dimension of input data too large.')
    endif

    aq_r2_aq = 0.0D0

!$omp parallel do private(nstr,nend,sqrInm,sqrInm2m,sqrInp2m,n,m)
    do j=1,jmax
       if ( mod(md(j),2) .eq. mod(km,2) ) then
          nstr=md(j) 
           nend=km
       else
          nstr=md(j) 
           nend=km-1
       endif

       m = abs(md(j))

       n = nstr
       sqrInm   = sqrt(Inm(n,m))
       sqrInp2m = sqrt(Inm(n+2,m))
       aq_r2_aq(j,n) = (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) /((2*n+gamma+1)*(2*n+gamma-3)) * aq_data(j,n) + (n-m+gamma-beta)*(n+m+gamma-1) /((2*n+gamma+3)*(2*n+gamma+1)) * aq_data(j,n+2) * (sqrInm/sqrInp2m)

       do n=nstr+2,nend-2,2
          sqrInm2m = sqrInm
          sqrInm   = sqrInp2m
          sqrInp2m = sqrt(Inm(n+2,m))

          aq_r2_aq(j,n) = (n-m)*(n+m+beta-1)/((2*n+gamma-5)*(2*n+gamma-3)) * aq_data(j,n-2)* (sqrInm/sqrInm2m) + (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) /((2*n+gamma+1)*(2*n+gamma-3)) * aq_data(j,n) + (n-m+gamma-beta)*(n+m+gamma-1) /((2*n+gamma+3)*(2*n+gamma+1)) * aq_data(j,n+2) * (sqrInm/sqrInp2m)
       enddo

       n = nend
       sqrInm2m = sqrInm
       sqrInm   = sqrInp2m
       aq_r2_aq(j,n) = (n-m)*(n+m+beta-1)/((2*n+gamma-5)*(2*n+gamma-3)) * aq_data(j,n-2)* (sqrInm/sqrInm2m) + (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) /((2*n+gamma+1)*(2*n+gamma-3)) * aq_data(j,n)
    enddo
!$omp end parallel do

    aq_r2_aq = aq_r2_aq * ra**2

  end function aq_r2_aq
aq_rDr_aq( aq_data ) result(aq_rDr_aq)
Function :
aq_rDr_aq :real(8), dimension(size(aq_data,1),0:size(aq_data,2)-1)
: (out) 出力スペクトルデータ
aq_data :real(8), dimension(:,0:), intent(in)
: (in) 入力スペクトルデータ

入力スペクトルデータに対して微分 r(d/dr) のスペクトル係数 を計算する(2 次元配列用). 

  a_n = aq_data/sqrt(Inm),  b_n = aq_rDr_aq/sqrt(Inm)

  b_n = (2n+gamma-1)/(2n+gamma+3)b_n+2
      + (2n+gamma-1)(n+gamma+1)/(2n+gamma+3)a_n+2 + n a_n

[Source]

  function aq_rDr_aq(aq_data)
    !
    ! 入力スペクトルデータに対して微分 r(d/dr) のスペクトル係数
    ! を計算する(2 次元配列用).
    !
    !   a_n = aq_data/sqrt(Inm),  b_n = aq_rDr_aq/sqrt(Inm)
    !
    !   b_n = (2n+gamma-1)/(2n+gamma+3)b_n+2
    !       + (2n+gamma-1)(n+gamma+1)/(2n+gamma+3)a_n+2 + n a_n
    !
    real(8), dimension(:,0:), intent(in)     :: aq_data
    !(in) 入力スペクトルデータ

    real(8), dimension(size(aq_data,1),0:size(aq_data,2)-1) :: aq_rDr_aq
    !(out) 出力スペクトルデータ

    integer :: j, n
    integer :: nstr, nend
    real(8) :: sqInm, sqInm2

    if ( size(aq_data,1) /= jmax ) then
       call MessageNotify('E','aq_rDr_aq', '1st dim. of the spectral data should be same as dim. of MD.')
    end if

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_rDr_aq', 'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_rDr_aq', 'The spectral dimension of input data too large.')
    endif

    aq_rDr_aq = 0.0D0
!$omp parallel do private(nstr,nend,sqInm,sqInm2,n)
    do j=1,jmax
       if ( mod(md(j),2) .eq. mod(km,2) ) then
          nstr=md(j) 
           nend=km
       else
          nstr=md(j) 
           nend=km-1
       endif

       sqInm = sqrt(Inm(nend,md(j)))
       aq_rDr_aq(j,nend) = nend * aq_data(j,nend)

       do n=nend-2,nstr,-2
          sqInm2 = sqInm
          sqInm = sqrt(Inm(n,md(j)))

          aq_rDr_aq(j,n) = sqInm/sqInm2 * ( (2*n+gamma-1)/(2*n+gamma+3) * aq_rDr_aq(j,n+2) + (2*n+gamma-1)*(n+gamma+1)/(2*n+gamma+3) * aq_data(j,n+2) ) + n * aq_data(j,n)
       enddo
    enddo
!$omp end parallel do

  end function aq_rDr_aq
g_R
Variable :
g_R(:) :real(8), allocatable
: ガウス—ラダウ格子点
g_R_Weight
Variable :
g_R_Weight(:) :real(8), allocatable
: ガウス重み
g_q( q_data ) result(g_q)
Function :
g_q :real(8), dimension(im)
: (out) 格子点データ
q_data :real(8), dimension(:), intent(in)
: (in) スペクトルデータ

スペクトルデータから格子データへ変換する(1 次元配列用).

[Source]

  function g_q(q_data)
    !
    ! スペクトルデータから格子データへ変換する(1 次元配列用).
    !
    real(8), dimension(:), intent(in)  :: q_data
    !(in) スペクトルデータ

    real(8), dimension(im)             :: g_q
    !(out) 格子点データ

    real(8), dimension(1,size(q_data)) :: q_work
    ! 作業用配列
    real(8), dimension(1,im)         :: g_work
    ! 作業用配列

    q_work(1,:) = q_data
    g_work = ag_aq(q_work)
    g_q = g_work(1,:)

  end function g_q
q_g( g_data ) result(q_g)
Function :
q_g :real(8), dimension(0:km)
: (out) スペクトルデータ
g_data :real(8), dimension(:), intent(in)
: (in) 格子点データ

格子データからスペクトルデータへ変換する(1 次元配列用).

[Source]

  function q_g(g_data)
    !
    ! 格子データからスペクトルデータへ変換する(1 次元配列用).
    !
    real(8), dimension(:), intent(in)     :: g_data
    !(in) 格子点データ

    real(8), dimension(0:km)              :: q_g
    !(out) スペクトルデータ

    real(8), dimension(1,size(g_data)) :: ag_work
    real(8), dimension(1,0:km)         :: aq_work

    ag_work(1,:) = g_data
    aq_work = aq_ag(ag_work)
    q_g = aq_work(1,:)

  end function q_g
q_r2Inv_q( q_data ) result(q_r2Inv_q)
Function :
q_r2Inv_q :real(8), dimension(0:km)
: (out) チェビシェフデータの R 微分
q_data :real(8), dimension(:), intent(in)
: (in) 入力チェビシェフデータ

入力スペクトルデータに対して積 r^2 のスペクトル係数 を計算する(1 次元配列用).

[Source]

  function q_r2Inv_q(q_data)
    !
    ! 入力スペクトルデータに対して積 r^2 のスペクトル係数
    ! を計算する(1 次元配列用).
    !
    real(8), dimension(:), intent(in)   :: q_data
    !(in) 入力チェビシェフデータ

    real(8), dimension(0:km)            :: q_r2Inv_q
    !(out) チェビシェフデータの R 微分

    real(8), dimension(1,size(q_data))  :: aq_work
    ! 作業用配列

    aq_work(1,:) = q_data
    aq_work = aq_r2Inv_aq(aq_work)
    q_r2Inv_q = aq_work(1,:)

  end function q_r2Inv_q
q_r2_q( q_data ) result(q_r2_q)
Function :
q_r2_q :real(8), dimension(0:km)
: (out) チェビシェフデータの R 微分
q_data :real(8), dimension(:), intent(in)
: (in) 入力チェビシェフデータ

入力スペクトルデータに対して積 r^2 のスペクトル係数 を計算する(1 次元配列用).

[Source]

  function q_r2_q(q_data)
    !
    ! 入力スペクトルデータに対して積 r^2 のスペクトル係数
    ! を計算する(1 次元配列用).
    !
    real(8), dimension(:), intent(in)   :: q_data
    !(in) 入力チェビシェフデータ

    real(8), dimension(0:km)            :: q_r2_q
    !(out) チェビシェフデータの R 微分

    real(8), dimension(1,size(q_data))  :: aq_work
    ! 作業用配列

    aq_work(1,:) = q_data
    aq_work = aq_r2_aq(aq_work)
    q_r2_q = aq_work(1,:)

  end function q_r2_q
q_rDr_q( q_data ) result(q_rDr_q)
Function :
q_rDr_q :real(8), dimension(0:km)
: (out) チェビシェフデータの R 微分
q_data :real(8), dimension(:), intent(in)
: (in) 入力チェビシェフデータ

入力スペクトルデータに r(d/dR) 微分を作用する(1 次元配列用).

スペクトルデータの r(d/dR) 微分とは, 対応する格子点データに R 微分を 作用させたデータのスペクトル変換のことである.

[Source]

  function q_rDr_q(q_data)
    !
    ! 入力スペクトルデータに r(d/dR) 微分を作用する(1 次元配列用).
    !
    ! スペクトルデータの r(d/dR) 微分とは, 対応する格子点データに R 微分を
    ! 作用させたデータのスペクトル変換のことである.
    !
    !
    real(8), dimension(:), intent(in)   :: q_data
    !(in) 入力チェビシェフデータ

    real(8), dimension(0:km)            :: q_rDr_q
    !(out) チェビシェフデータの R 微分

    real(8), dimension(1,size(q_data))  :: aq_work
    ! 作業用配列

    aq_work(1,:) = q_data
    aq_work = aq_rDr_aq(aq_work)
    q_rDr_q = aq_work(1,:)

  end function q_rDr_q