wt_module

    function AvrLatRad_yz(yz)  ! 緯度動径(子午面)積分
      !
      ! 2 次元(YZ)格子点データの緯度動径(子午面)平均
      !
      ! 2 次元データ f(φ,r) に対して
      !
      !    ∫f(φ,r) r^2cosφ dφdr /(2(r[o]^3-r[i]^3)/3) 
      !
      ! を計算する.
      !
      real(8), dimension(jm,0:km), intent(in) :: yz
      !(in) 2 次元緯度動径(子午面)格子点データ

      real(8)                   :: AvrLatRad_yz
      !(out) 平均値

      AvrLatRad_yz = IntLatRad_yz(yz)/(sum(y_Lat_Weight)*sum(z_Rad_Weight))

    end function AvrLatRad_yz

    function AvrLonLatRad_xyz(xyz) ! 緯度経度動径(全球)積分
      !
      ! 3 次元格子点データの緯度経度動径(全球)積分
      !
      ! 3 次元データ f(λ,φ,r) に対して
      !
      !    ∫f(λ,φ,r) r^2cosφ dλdφdr /(4π(r[o]^3-r[i]^3)/3) 
      !
      ! を計算する.
      !
      real(8), dimension(im,jm,0:km), intent(in) :: xyz
      !(in) 3 次元経度緯度動径格子点データ

      real(8)                     :: AvrLonLatRad_xyz
      !(out) 全球平均値

      AvrLonLatRad_xyz = IntLonLatRad_xyz(xyz) /(sum(x_Lon_Weight)*sum(y_Lat_Weight) * sum(z_Rad_Weight))

    end function AvrLonLatRad_xyz

    function AvrLonRad_xz(xz)  ! 経度動径(緯度円)積分
      !
      ! 2 次元(XZ)格子点データの経度動径平均
      !
      ! 2 次元データ f(λ,r) に対して 
      ! 
      !    ∫f(λ,r) r^2dλdr /(2π(r[o]^3-r[i]^3)/3)
      !
      ! を計算する.
      !
      real(8), dimension(im,0:km), intent(in) :: xz    ! 2 次元格子点データ
      real(8)                                 :: AvrLonRad_xz      ! 積分値

      AvrLonRad_xz = IntLonRad_xz(xz)/(sum(x_Lon_Weight)*sum(z_Rad_Weight))

    end function AvrLonRad_xz

    function AvrRad_z(z)
      !
      ! 1 次元(Z)格子点データの動径方向域平均.
      !
      ! 1 次元データ f(r) に対して ∫f(r) r^2dr /((r[o]^3-r[i]^3)/3) を
      ! 計算する.
      !
      real(8), dimension(im,0:km), intent(in) :: z
      !(in) 1 次元動径格子点データ

      real(8)                                 :: AvrRad_z
      !(out) 平均値

      AvrRad_z = IntRad_z(z)/sum(z_Rad_Weight)

    end function AvrRad_z

    function IntLatRad_yz(yz)
      !
      ! 2 次元(YZ)格子点データの緯度動径積分(子午面)および平均
      !
      ! 2 次元データ f(φ,r) に対して ∫f(φ,r) r^2cosφ dφdr を計算する.
      !
      real(8), dimension(jm,0:km), intent(in) :: yz
      !(in) 2 次元緯度動径(子午面)格子点データ

      real(8)                   :: IntLatRad_yz
      !(out) 積分値
      integer :: j, k

      IntLatRad_yz = 0
      do k=0,km
         do j=1,jm
            IntLatRad_yz = IntLatRad_yz + yz(j,k) * y_Lat_Weight(j) * z_Rad_Weight(k)
         enddo
      enddo
    end function IntLatRad_yz

    function IntLonLatRad_xyz(xyz) ! 緯度経度動径(全球)積分
      !
      ! 3 次元格子点データの緯度経度動径(全球)積分
      !
      ! 3 次元データ f(λ,φ,r) に対して
      !
      !     ∫f(λ,φ,r) r^2cosφ dλdφdr 
      !
      ! を計算する.
      !
      real(8), dimension(im,jm,0:km), intent(in) :: xyz 
      !(in) 3 次元経度緯度動径格子点データ

      real(8)                     :: IntLonLatRad_xyz 
      !(out) 全球積分値

      integer :: i, j, k

      IntLonLatRad_xyz = 0
      do k=0,km
         do j=1,jm
            do i=1,im
               IntLonLatRad_xyz = IntLonLatRad_xyz + xyz(i,j,k) * x_Lon_Weight(i) * y_Lat_Weight(j) * z_Rad_Weight(k)
            enddo
         enddo
      enddo
    end function IntLonLatRad_xyz

    function IntLonRad_xz(xz)  ! 経度動径(緯度円)積分
      !
      ! 2 次元(XZ)格子点データの経度動径積分
      !
      ! 2 次元データ f(λ,r) に対して∫f(λ,r) r^2dλdr を計算する.
      !
      real(8), dimension(im,0:km), intent(in) :: xz
      !(in) 2 次元緯度動径格子点データ

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

      integer :: i, k

      IntLonRad_xz = 0
      do k=0,km
         do i=1,im
            IntLonRad_xz = IntLonRad_xz + xz(i,k) * x_Lon_Weight(i) * z_Rad_Weight(k)
         enddo
      enddo
    end function IntLonRad_xz

    function IntRad_z(z)  ! 動径積分
      !
      ! 1 次元(Z)格子点データの動径方向域積分.
      !
      ! 1 次元データ f(r) に対して ∫f(r) r^2dr を計算する.
      !
      real(8), dimension(0:km), intent(in) :: z
      !(in) 1 次元動径格子点データ

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

      integer :: k

      IntRad_z = 0.0d0
      do k=0,km
         IntRad_z = IntRad_z + z(k) * z_Rad_Weight(k) 
      enddo
    end function IntRad_z

    function Interpolate_array000_wt(wt_data,alon,alat,arad)
      !
      ! 緯度 alon, 経度 alat 動径 arad における関数値を
      ! その球面調和変換係数 wa_data から補間計算する
      !
      real(8), intent(IN) :: wt_data((nm+1)**2,0:km)  ! スペクトルデータ
      real(8), intent(IN) :: alon                     ! 補間する位置(経度)
      real(8), intent(IN) :: alat                     ! 補間する位置(緯度)
      real(8), intent(IN) :: arad                     ! 補間する位置(動径)
      real(8) :: Interpolate_array000_wt               ! 補間した値
      
      Interpolate_array000_wt = Interpolate_w(a_Interpolate_at(wt_data,arad),alon,alat)

    end function Interpolate_array000_wt

    function Interpolate_array111_wt(wt_data,alon,alat,arad)
      !
      ! 緯度 alon(:), 経度 alat(:) 動径 arad(:) における関数値を
      ! その球面調和変換係数 wa_data から補間計算する. 
      ! alon, alat, arad は同じ大きさでなければならない. 
      ! かえす値の形状は dimension(size(alon))
      !
      real(8), intent(IN) :: wt_data((nm+1)**2,0:km)  ! スペクトルデータ
      real(8), intent(IN) :: alon(:)                  ! 補間する位置(経度)
      real(8), intent(IN) :: alat(:)                  ! 補間する位置(緯度)
      real(8), intent(IN) :: arad(:)                  ! 補間する位置(動径)
      real(8) :: Interpolate_array111_wt(size(alon))  ! 補間した値
      integer :: i

      if ( size(alon) == size(alat) .AND. size(alon) == size(arad))then
         do i=1,size(alon)
            Interpolate_array111_wt(i) = Interpolate_array000_wt(wt_data,alon(i),alat(i),arad(i))
         enddo
      else
         call MessageNotify('E','Interpolate_array111_wt', 'Dimensions of lon,lat,rad miss match.')
      endif
    end function Interpolate_array111_wt

    function nmz_PoloidalEnergySpectrum_wt(wt_POLPOT)
      !
      ! ポロイダルポテンシャルから, ポロイダルエネルギーの
      ! 球面調和函数全波数 n, 帯状波数 m の各成分を計算する.
      !
      !  * 全波数 n, 帯状波数 m のポロイダルポテンシャルのスペクトル成分
      !    φ(n,m,r)から全波数 n, 帯状波数 m 成分のポロイダルエネルギー
      !    スペクトルは 
      !
      !      (1/2)n(n+1)4πr^2{[d(rφ(n,m,r))/dr]^2 + n(n+1)φ(n,m,r)^2} 
      !
      !    と計算される.
      !
      !  * 全てのエネルギースペクトル成分の和を動径積分したもの(r^2の重み無し)
      !    が球殻内での全エネルギーに等しい.
      !
      !  * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される.
      !    欠損値の値はモジュール変数 wt_VMiss によって設定できる
      !    (初期値は -999.0)
      !
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt_POLPOT
      !(in) ポロイダルポテンシャル

      real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_PoloidalEnergySpectrum_wt 
      !(out) エネルギースペクトルポロイダル成分


      real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA1   ! 作業領域
      real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA2   ! 作業領域
      integer :: n, m

      nmz_PoloidalEnergySpectrum_wt = wt_VMiss

      wz_Data1 = wz_wt(wt_POLPOT)
      wz_Data2 = wz_Rad*wz_wt(wt_DRad_wt(wt_POLPOT))     ! d(rφ)/dr + wz_wt(wt_POLPOT)                         ! = rdφ/dr+φ

      do n=0,nm
         do m=-n,n
            nmz_PoloidalEnergySpectrum_wt(n,m,:) = + 0.5* n*(n+1)* (4*pi) *( wz_Data2(l_nm(n,m),:)**2 + n*(n+1)*wz_Data1(l_nm(n,m),:)**2 )
         enddo
      enddo

    end function nmz_PoloidalEnergySpectrum_wt

    function nmz_ToroidalEnergySpectrum_wt(wt_TORPOT)
      !
      ! トロイダルポテンシャルから, トロイダルエネルギーの
      ! 球面調和函数全波数 n, 帯状波数 m の各成分を計算する
      !
      !  * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分
      !    ψ(n,m,r)から全波数 n, 帯状波数 m 成分のトロイダルエネルギー
      !    スペクトルは  (1/2)n(n+1)4πr^2ψ(n,m,r)^2  と計算される.
      !
      !  * 全てのエネルギースペクトル成分の和を動径積分したもの(r^2の重み無し)
      !    が球殻内での全エネルギーに等しい.
      !    
      !  * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される.
      !    wt_VMiss によって設定できる (初期値は -999.0)
      !
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt_TORPOT
      !(in) トロイダルポテンシャル

      real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_ToroidalEnergySpectrum_wt
      !(out) エネルギースペクトルトロイダル成分

      real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA   ! 作業領域
      integer :: n, m

      nmz_ToroidalEnergySpectrum_wt = wt_VMiss

      wz_DATA = wz_wt(wt_TORPOT)
      do n=0,nm
         do m=-n,n
            nmz_ToroidalEnergySpectrum_wt(n,m,:) = 0.5 * n*(n+1)* (4*pi) * z_Rad**2 * wz_DATA(l_nm(n,m),:)**2
         enddo
      enddo

    end function nmz_ToroidalEnergySpectrum_wt

    function nz_PoloidalEnergySpectrum_wt(wt_POLPOT)
      !
      ! ポロイダルポテンシャルから, ポロイダルエネルギーの
      ! 球面調和函数全波数の各成分を計算する
      !
      !  * 全波数 n, 帯状波数 m のポロイダルポテンシャルのスペクトル成分
      !    φ(n,m,r)から全波数 n 成分のポロイダルエネルギースペクトルは
      !
      !      Σ[m=-n]^n ((1/2)n(n+1)4πr^2{[d(rφ(n,m,r))/dr]^2 
      !                 + n(n+1)φ(n,m,r)^2} 
      !
      !    と計算される.
      !
      !  * 全ての全波数に対してのエネルギースペクトル成分の和を動径積分したもの
      !    (r^2の重み無し)が球殻内での全エネルギーに等しい.
      !
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt_POLPOT
      !(in) ポロイダルポテンシャル

      real(8), dimension(0:nm,0:km) :: nz_PoloidalEnergySpectrum_wt
      !(out) エネルギースペクトルポロイダル成分

      real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA1   ! 作業領域
      real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA2   ! 作業領域
      integer :: n, m

      wz_Data1 = wz_wt(wt_POLPOT)
      wz_Data2 = wz_Rad*wz_wt(wt_DRad_wt(wt_POLPOT))     ! d(rφ)/dr + wz_wt(wt_POLPOT)                         ! = rdφ/dr+φ

      do n=0,nm
         nz_PoloidalEnergySpectrum_wt(n,:) = + 0.5* n*(n+1)* (4*pi) *( sum(wz_Data2(l_nm(n,(/(m,m=-n,n)/)),:)**2,1) + n*(n+1)*sum(wz_Data1(l_nm(n,(/(m,m=-n,n)/)),:)**2,1) )
      enddo

    end function nz_PoloidalEnergySpectrum_wt

    function nz_ToroidalEnergySpectrum_wt(wt_TORPOT)
      !
      ! トロイダルポテンシャルから, トロイダルエネルギーの
      ! 球面調和函数全波数の各成分を計算する.
      !
      !  * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分
      !    ψ(n,m,r)から全波数 n 成分のトロイダルエネルギースペクトルは
      !    Σ[m=-n]^n(1/2)n(n+1)4πr^2ψ(n,m,r)^2 と計算される.
      !
      ! * 全てのエネルギースペクトル成分の和を動径積分したもの(r^2の重み無し)
      !    が球殻内での全エネルギーに等しい.
      !
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt_TORPOT
      !(in) トロイダルポテンシャル

      real(8), dimension(0:nm,0:km) :: nz_ToroidalEnergySpectrum_wt 
      !(out) エネルギースペクトルトロイダル成分

      real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA   ! 作業領域
      integer :: n, m

      wz_DATA = wz_wt(wt_TORPOT)
      do n=0,nm
         nz_ToroidalEnergySpectrum_wt(n,:) = 0.5 * n*(n+1)* (4*pi) * z_Rad**2 * sum(wz_Data(l_nm(n,(/(m,m=-n,n)/)),:)**2,1)
      enddo

    end function nz_ToroidalEnergySpectrum_wt

    subroutine wt_BoundariesGrid(wt,values,cond)
      !
      ! スペクトルデータにディリクレ・ノイマン境界条件を適用する
      ! 実空間での境界条件適用
      !
      ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
      ! 条件を課している(選点法). このルーチンを用いるためには 
      ! wt_Initial にて設定するチェビシェフ切断波数(lm)と鉛直格子点数(km)を
      ! 等しくしておく必要がある.
      !
      real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)      :: wt
              !(inout) 境界条件を適用するデータ. 修正された値を返す. 

      real(8), dimension((nm+1)*(nm+1),2), intent(in), optional :: values
              !(in) 境界での 値/勾配 分布を水平スペクトル変換したものを与える. 
              !    省略時は値/勾配 0 となる. 

      character(len=2), intent(in), optional             :: cond
              !(in) 境界条件. 省略時は 'DD'
              !        DD : 両端ディリクレ条件
              !        DN : 上端ディリクレ, 下端ノイマン条件
              !        ND : 上端ノイマン, 下端ディリクレ条件
              !        NN : 両端ノイマン条件

      if (.not. present(cond)) then
         if (present(values)) then
            call at_boundariesGrid_DD(wt,values)
         else
            call at_boundariesGrid_DD(wt)
         endif
         return
      endif

      select case(cond)
      case ('NN')
         if (present(values)) then
            call at_BoundariesGrid_NN(wt,values)
         else
            call at_BoundariesGrid_NN(wt)
         endif
      case ('DN')
         if (present(values)) then
            call at_BoundariesGrid_DN(wt,values)
         else
            call at_BoundariesGrid_DN(wt)
         endif
      case ('ND')
         if (present(values)) then
            call at_BoundariesGrid_ND(wt,values)
         else
            call at_BoundariesGrid_ND(wt)
         endif
      case ('DD')
         if (present(values)) then
            call at_BoundariesGrid_DD(wt,values)
         else
            call at_BoundariesGrid_DD(wt)
         endif
      case default
         call MessageNotify('E','wt_BoundariesGrid','B.C. not supported')
      end select

    end subroutine wt_BoundariesGrid

    subroutine wt_BoundariesTau(wt,values,cond)
      !
      ! スペクトルデータにディリクレ・ノイマン境界条件を適用する
      ! Chebyshev 空間での境界条件適用(タウ法)
      !
      ! チェビシェフ空間において境界条件を満たすべく高次の係数を
      ! 定める方法をとっている(タウ法).
      !
      real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)      :: wt
              !(inout) 境界条件を適用するデータ. 修正された値を返す. 

      real(8), dimension((nm+1)*(nm+1),2), intent(in), optional :: values
              !(in) 境界での 値/勾配 分布を水平スペクトル変換したものを与える. 
              !     省略時は値/勾配 0 となる. 

      character(len=2), intent(in), optional             :: cond
              !(in) 境界条件. 省略時は 'DD'
              !        DD : 両端ディリクレ条件
              !        DN : 上端ディリクレ, 下端ノイマン条件
              !        ND : 上端ノイマン, 下端ディリクレ条件
              !        NN : 両端ノイマン条件

      if (.not. present(cond)) then
         if (present(values)) then
            call at_BoundariesTau_DD(wt,values)
         else
            call at_BoundariesTau_DD(wt)
         endif
         return
      endif

      select case(cond)
      case ('NN')
         if (present(values)) then
            call at_BoundariesTau_NN(wt,values)
         else
            call at_BoundariesTau_NN(wt)
         endif
      case ('DN')
         if (present(values)) then
            call at_BoundariesTau_DN(wt,values)
         else
            call at_BoundariesTau_DN(wt)
         endif
      case ('ND')
         if (present(values)) then
            call at_BoundariesTau_ND(wt,values)
         else
            call at_BoundariesTau_ND(wt)
         endif
      case ('DD')
         if (present(values)) then
            call at_BoundariesTau_DD(wt,values)
         else
            call at_BoundariesTau_DD(wt)
         endif
      case default
         call MessageNotify('E','wt_BoundariesTau','B.C. not supported')
      end select

    end subroutine wt_BoundariesTau

    function wt_DRad_wt(wt)
      !
      ! 入力スペクトルデータに動径微分 ∂/∂r を作用する.
      !
      ! スペクトルデータの動径微分とは, 対応する格子点データに動径微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(8), dimension((nm+1)*(nm+1),0:lm)             :: wt_DRad_wt
      !(in) 動径微分された2 次元球面調和函数チェビシェフスペクトルデータ

      wt_DRad_wt = at_Dr_at(wt)

    end function wt_DRad_wt

    function wt_DivLat_xyz(xyz)
      !
      ! 格子データに発散型緯度微分 1/rcosφ・∂(f cosφ)/∂φ を
      ! 作用させたスペクトルデータを返す.
      !
      real(8), dimension(im,jm,0:km), intent(in)   :: xyz
      !(in) 3 次元経度緯度動径格子点データ

      real(8), dimension((nm+1)*(nm+1),0:lm)       :: wt_DivLat_xyz
      !(out) 発散型緯度微分を作用された 2 次元スペクトルデータ

      wt_DivLat_xyz = wt_wz(wa_divlat_xya(xyz/xyz_Rad))
    end function wt_DivLat_xyz

    function wt_DivLon_xyz(xyz)
      ! 
      ! 格子点データに発散型経度微分 1/rcosφ・∂/∂λ を作用させた
      ! スペクトルデータを返す.
      !
      real(8), dimension(im,jm,0:km), intent(in)   :: xyz
      !(in) 3 次元経度緯度動径格子点データ

      real(8), dimension((nm+1)*(nm+1),0:lm)       :: wt_DivLon_xyz
      !(out) 発散型経度微分を作用された 2 次元スペクトルデータ

      wt_DivLon_xyz = wt_wz(wa_DivLon_xya(xyz/xyz_Rad))
    end function wt_DivLon_xyz

    function wt_DivRad_wt(wt)
      ! 
      ! 入力スペクトルデータに発散型動径微分
      !
      !       1/r^2 ∂/∂r (r^2 .)= ∂/∂r + 2/r
      !
      ! を作用する.
      !
      ! スペクトルデータの発散型動径微分とは, 対応する格子点データに
      ! 発散型動径微分を作用させたデータのスペクトル変換のことである. 
      !
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(8), dimension((nm+1)*(nm+1),0:lm)             :: wt_DivRad_wt
      !(out) 発散型動径微分を作用された 2 次元スペクトルデータ

      wt_DivRad_wt = wt_Drad_wt(wt) + wt_wz(2/wz_rad*wz_wt(wt))


    end function wt_DivRad_wt

    function wt_Div_xyz_xyz_xyz(xyz_Vlon,xyz_Vlat,xyz_Vrad)
      !
      ! べクトル成分である 3 つの格子データに発散を作用させた
      ! スペクトルデータを返す.
      !
      ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, 
      ! 動径成分を表し, 発散は 
      !
      !      1/rcosφ・∂u/∂λ + 1/rcosφ・∂(v cosφ)/∂φ 
      !    + 1/r^2 ∂/∂r (r^2 w)
      !
      ! と計算される.
      !
      real(8), dimension(im,jm,0:km), intent(in) :: xyz_Vlon
      !(in) ベクトル場の経度成分

      real(8), dimension(im,jm,0:km), intent(in) :: xyz_Vlat
      !(in) ベクトル場の緯度成分

      real(8), dimension(im,jm,0:km), intent(in) :: xyz_Vrad
      !(in) ベクトル場の動径成分

      real(8), dimension((nm+1)*(nm+1),0:lm)     :: wt_Div_xyz_xyz_xyz
      !(out) ベクトル場の発散

      wt_Div_xyz_xyz_xyz =   wt_DivLon_xyz(xyz_Vlon) + wt_DivLat_xyz(xyz_Vlat) + wt_DivRad_wt(wt_xyz(xyz_Vrad))

    end function wt_Div_xyz_xyz_xyz

   subroutine wt_Initial(i,j,k,n,l,r_in,r_out,np)
     !
     ! スペクトル変換の格子点数, 波数, 動径座標の範囲を設定する.
     !
     ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を
     ! しなければならない. 
     !
     ! np に 1 より大きな値を指定すれば ISPACK の球面調和函数変換 
     ! OPENMP 並列計算ルーチンが用いられる. 並列計算を実行するには, 
     ! 実行時に環境変数 OMP_NUM_THREADS を np 以下の数字に設定する等の
     ! システムに応じた準備が必要となる. 
     !
     ! np に 1 より大きな値を指定しなければ並列計算ルーチンは呼ばれない.
     !
     !
     integer,intent(in) :: i              ! 格子点数(経度λ)
     integer,intent(in) :: j              ! 格子点数(緯度φ)
     integer,intent(in) :: k              ! 格子点数(動径 r)
     integer,intent(in) :: n              ! 切断波数(水平全波数)
     integer,intent(in) :: l              ! 切断波数(動径波数)

     real(8),intent(in) :: r_in           ! 球殻内半径
     real(8),intent(in) :: r_out          ! 球殻外半径

     integer,intent(in), optional :: np   ! OPENMP での最大スレッド数

     im = i  
      jm = j 
      km = k
     nm = n  
      lm = l
     ri = r_in 
      ro = r_out

     if ( present(np) ) then
        call wa_Initial(nm,im,jm,km+1,np)
     else
        call wa_Initial(nm,im,jm,km+1)
     endif

     call at_Initial(km,lm,r_in,r_out)

     allocate(xyz_Lon(im,jm,0:km))
     allocate(xyz_Lat(im,jm,0:km))
     allocate(xyz_Rad(im,jm,0:km))

     allocate(wz_Rad((nm+1)*(nm+1),0:km))

     xyz_Lon = spread(xy_Lon,3,km+1)
     xyz_Lat = spread(xy_Lat,3,km+1)
     xyz_Rad = spread(spread(z_Rad,1,jm),1,im)

     wz_Rad = spread(z_Rad,1,(nm+1)*(nm+1))

     z_Rad_Weight = z_Rad_Weight * z_Rad**2       ! r^2 dr の積分重み
   end subroutine wt_initial

    function wt_KxRGrad_wt(wt)
      !
      ! 入力スペクトルデータに経度微分 k×r・▽ = ∂/∂λを作用する.
      !
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(8), dimension((nm+1)*(nm+1),0:lm)             :: wt_KxRGrad_wt
      !(out) 経度微分を作用された 2 次元スペクトルデータ

      wt_KxRGrad_wt =  wa_Dlon_wa(wt)
    end function wt_KxRGrad_wt

    function wt_L2Inv_wt(wt)
      !
      ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を
      ! 作用する.
      !
      ! スペクトルデータに L^2 演算子を作用させる関数 wt_L2_wt の逆計算を
      ! 行う関数である.
      !
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(8), dimension((nm+1)*(nm+1),0:lm)             :: wt_L2Inv_wt
      !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ

      wt_L2Inv_wt = -wa_LaplaInv_wa(wt)
    end function wt_L2Inv_wt

    function wt_L2_wt(wt)
      !
      ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する.
      !
      ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる. 
      !  入力スペクトルデ ータに対応する格子点データに演算子 
      !
      !     L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ)
      !
      ! を作用させたデータのスペクトル変換が返される.
      !
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(8), dimension((nm+1)*(nm+1),0:lm)             :: wt_L2_wt
      !(out) L^2 演算子を作用された 2 次元スペクトルデータ

      wt_L2_wt = -wa_Lapla_wa(wt)
    end function wt_L2_wt

    function wt_LaplaPol2PolGrid_wt(wt,cond)
      !
      ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する.
      ! チェビシェフ格子点空間で境界条件を適用している. 
      !
      ! この関数を用いるためには wt_Initial にて設定する
      ! チェビシェフ切断波数(lm)と鉛直格子点数(km)を等しく
      ! しておく必要がある. 
      !
      ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は
      !
      !    ▽^2Φ = f
      !      Φ = const. at boundaries.
      !      ∂Φ/∂r = 0 at boundaries          (粘着条件) 
      !      or ∂^2Φ/∂r^2 = 0 at boundaries   (応力なし条件)
      !
      ! 最初に呼ばれるときの境界条件で以後計算される(要仕様変更). 
      !
      ! 最終的にチェビシェフ係数の解が欲しい場合には, wz_LaplaPol2Pol_wz に
      ! 比べてチェビシェフ -- 格子点変換が 1 回分少なくて済む.
      !
      real(8), dimension((nm+1)*(nm+1),0:lm),intent(in)  :: wt
              !(in) 入力▽^2φ分布

      real(8), dimension((nm+1)*(nm+1),0:lm)             :: wt_LaplaPol2PolGrid_wt
              !(out) 出力ポロイダルポテンシャル分布

      character(len=2), intent(in), optional  :: cond
              !(in) 境界条件スイッチ. 省略時は 'RR'
              !     RR    : 両端粘着条件
              !     RF    : 上端粘着, 下端応力なし条件
              !     FR    : 上端応力なし, 下端粘着条件
              !     FF    : 両端応力なし条件

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

      real(8), dimension((nm+1)*(nm+1),0:km)  :: wz_work
      real(8), dimension((nm+1)*(nm+1),0:lm)  :: wt_work
      real(8), dimension(0:lm,0:lm)           :: tt_I
      real(8), dimension(0:lm,0:km)           :: tz_work
      logical                                 :: rigid1, rigid2   ! 境界条件

      logical :: first = .true.
      integer :: l,n
      save    :: alu, kp, first

      if (.not. present(cond)) then
         rigid1=.TRUE. 
          rigid2=.TRUE.
      else
         select case (cond)
         case ('RR')
            rigid1 = .TRUE.  
             rigid2 = .TRUE.
         case ('RF')
            rigid1 = .TRUE.  
             rigid2 = .FALSE.
         case ('FR')
            rigid1 = .FALSE. 
             rigid2 = .TRUE.
         case ('FF')
            rigid1 = .FALSE. 
             rigid2 = .FALSE.
         case default
            call MessageNotify('E','wt_LaplaPol2PolGrid_wt','B.C. not supported')
         end select
      endif

      if ( first ) then
         first = .false.

         if ( lm /= km ) then
            call MessageNotify('E','wt_LaplaPol2PolGrid_wt', 'Chebyshev truncation and number of grid points should be same.')
         endif

         allocate(alu((nm+1)*(nm+1),0:km,0:lm),kp((nm+1)*(nm+1),0:lm))

         do l=0,lm
            wt_work = 0 
             wt_work(:,l) = 1

            ! 各水平波数に関して独立の式
            alu(:,:,l) = wz_wt(wt_Lapla_wt(wt_work))
         enddo

         ! 運動学的条件. 流線は境界で一定
         tt_I = 0.0
         do l=0,lm
            tt_I(l,l)=1
         enddo

         ! 非電気伝導体
         tz_work = az_at(tt_I)

         do n=1,(nm+1)*(nm+1)
            alu(n,0,:)  = tz_work(:,0)
            alu(n,km,:) = tz_work(:,km)
         enddo

         ! 力学的条件粘着条件 
         if ( rigid1 ) then
            tz_work=az_at(at_Dr_at(tt_I))
         else
            tz_work=az_at(at_Dr_at(at_Dr_at(tt_I)))
         endif
         do n=1,(nm+1)*(nm+1)
            alu(n,1,:) = tz_work(:,0)
         enddo

         ! 力学的条件粘着条件 
         if ( rigid2 ) then
            tz_work=az_at(at_Dr_at(tt_I))
         else
            tz_work=az_at(at_Dr_at(at_Dr_at(tt_I)))
         endif
         do n=1,(nm+1)*(nm+1)
            alu(n,km-1,:) = tz_work(:,km)
         enddo

         call ludecomp(alu,kp)
      endif

      wz_work         = wz_wt(wt)
      wz_work(:,1)    = 0               ! 力学的条件
      wz_work(:,km-1) = 0               ! 力学的条件
      wz_work(:,0)    = 0               ! 運動学的条件
      wz_work(:,km)   = 0               ! 運動学的条件 

      wt_LaplaPol2PolGrid_wt = lusolve(alu,kp,wz_work)

    end function wt_LaplaPol2PolGrid_wt

    function wt_Lapla_wt(wt)
      ! 入力スペクトルデータにラプラシアン
      !
      !     ▽^2 =   1/r^2 cos^2φ・∂^2/∂λ^2 
      !            + 1/r^2 cosφ・∂/∂φ(cosφ∂/∂φ) 
      !            + 1/r^2 ∂/∂r (r^2 ∂/∂r) 
      !
      ! を作用する.
      !
      ! スペクトルデータのラプラシアンとは, 対応する格子点データに
      ! ラプラシアンを作用させたデータのスペクトル変換のことである. 
      !
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(8), dimension((nm+1)*(nm+1),0:lm)             :: wt_Lapla_wt
      !(out) ラプラシアンを作用された 2 次元スペクトルデータ

      wt_Lapla_wt = wt_DivRad_wt(wt_Drad_wt(wt)) + wt_wz(wz_wt(wa_Lapla_wa(wt))/wz_Rad**2)

    end function wt_Lapla_wt

    subroutine wt_PolmagBoundariesGrid(wt_POL)
      !
      ! 磁場ポロイダルポテンシャルに対して境界条件を適用する.
      ! 鉛直実空間での境界条件適用. 
      !
      ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
      ! 条件を課している(選点法). このルーチンを用いるためには 
      ! wt_Initial にて設定するチェビシェフ切断波数(lm)と鉛直格子点数(km)を
      ! 等しくしておく必要がある. 
      !
      ! 現在のところ境界物質が非電気伝導体の場合のみ対応している. 
      ! その場合, 磁場ポロイダルポテンシャルの各水平スペクトル成分 h に
      ! たいして境界条件が与えられ,
      !
      !  * 外側境界 : dh/dr + (n+1)h/r = 0
      !  * 内側境界 : dh/dr - nh/r = 0
      !
      ! である. ここで n は h の水平全波数である. 
      !
      ! 最初に呼ばれるときの境界条件で以後計算される(要仕様変更).
      !
      !
      real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)   :: wt_POL
              !(inout) 境界条件を適用するデータ. 修正された値を返す. 

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

      real(8), dimension(:,:), allocatable    :: tt_I
      real(8), dimension(:,:), allocatable    :: tz_PSI
      real(8), dimension(:,:), allocatable    :: tz_DPSIDR
      real(8), dimension((nm+1)*(nm+1),0:km)  :: wz_POL

      logical :: first = .true.
      integer  :: l, n, nn(2)
      save     :: alu, kp, first

      if ( first ) then
         first = .false.

         if ( lm /= km ) then
            call MessageNotify('E','PolMagBoundariesGrid', 'Chebyshev truncation and number of grid points should be same.')
         endif

         allocate(alu((nm+1)*(nm+1),0:lm,0:lm),kp((nm+1)*(nm+1),0:lm))
         allocate(tt_I(0:lm,0:lm),tz_PSI(0:lm,0:km),tz_DPSIDR(0:lm,0:km))

         tt_I = 0.0
         do l=0,lm
            tt_I(l,l)=1
         enddo
         do n=1,(nm+1)*(nm+1)
            alu(n,:,:) = transpose(az_at(tt_I))  ! 内部領域は値を保存
         enddo

         ! 非電気伝導体
         tz_PSI = az_at(tt_I)
         tz_DPSIDR = az_at(at_dr_at(tt_I))

         do n=1,(nm+1)*(nm+1)
            nn=nm_l(n)
            alu(n,0,:)  = tz_DPSIDR(:,0) + (nn(1)+1) * tz_PSI(:,0)/z_RAD(0)
            alu(n,km,:) = tz_DPSIDR(:,km) - nn(1) * tz_PSI(:,km)/z_RAD(km)
         enddo
         call ludecomp(alu,kp)

         deallocate(tt_I,tz_PSI,tz_DPSIDR)
      endif

      wz_POL       = wz_wt(wt_POL)
      wz_POL(:,0)  = 0.0
      wz_POL(:,km) = 0.0
      wt_POL = lusolve(alu,kp,wz_POL)

    end subroutine wt_PolmagBoundariesGrid

    subroutine wt_PolmagBoundariesTau(wt_POL)
      !
      ! 磁場ポロイダルポテンシャルに対して境界条件を適用する.
      ! Chebyshev 空間での境界条件適用
      !
      ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を
      ! とっている(タウ法). 現在のところ境界物質が非電気伝導体の場合のみ
      ! 対応している. その場合, 磁場ポロイダルポテンシャルの各水平スペクトル
      ! 成分 h にたいして境界条件が与えられ,
      !
      !  * 外側境界 : dh/dr + (n+1)h/r = 0
      !  * 内側境界 : dh/dr - nh/r = 0
      !
      ! である. ここで n は h の水平全波数である. 
      !
      ! 最初に呼ばれるときの境界条件で以後計算される(要仕様変更).
      !
      real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)   :: wt_POL
              !(inout) 境界条件を適用するデータ. 修正された値を返す. 

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

      real(8), dimension(:,:), allocatable    :: tt_I
      real(8), dimension(:,:), allocatable    :: tz_PSI
      real(8), dimension(:,:), allocatable    :: tz_DPSIDR

      logical :: first = .true.
      integer  :: l, n, nn(2)
      save     :: alu, kp, first

      if ( first ) then
         first = .false.

         allocate(alu((nm+1)*(nm+1),0:lm,0:lm),kp((nm+1)*(nm+1),0:lm))
         allocate(tt_I(0:lm,0:lm),tz_PSI(0:lm,0:km),tz_DPSIDR(0:lm,0:km))

         tt_I = 0.0
         do l=0,lm
            tt_I(l,l)=1
         enddo
         do n=1,(nm+1)*(nm+1)
            alu(n,:,:) = tt_I
         enddo

         ! 非電気伝導体
         tz_PSI = az_at(tt_I)
         tz_DPSIDR = az_at(at_dr_at(tt_I))

         do n=1,(nm+1)*(nm+1)
            nn=nm_l(n)
            alu(n,lm-1,:) = tz_DPSIDR(:,0) + (nn(1)+1) * tz_PSI(:,0)/z_RAD(0)
            alu(n,lm,:)   = tz_DPSIDR(:,km) - nn(1) * tz_PSI(:,km)/z_RAD(km)
         enddo
         call ludecomp(alu,kp)

         deallocate(tt_I,tz_PSI,tz_DPSIDR)
      endif

      wt_POL(:,lm-1) = 0.0
      wt_POL(:,lm)   = 0.0
      wt_POL = lusolve(alu,kp,wt_POL)

    end subroutine wt_PolmagBoundariesTau

    subroutine wt_Potential2Rotation( xyz_RotVLON,xyz_RotVLAT,xyz_RotVRAD,wt_TORPOT,wt_POLPOT)
      !
      ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
      !
      !     v = ▽x(Ψr) + ▽x▽x(Φr) 
      !
      ! に対して, その回転
      !
      !     ▽xv = ▽x▽x(Ψr) + ▽x▽x▽x(Φr) = ▽x▽x(Ψr) - ▽x((▽^2Φ)r)
      !
      ! を計算する. 
      
      ! ベクトル場の回転
      real(8), dimension(im,jm,0:km), intent(OUT) :: xyz_RotVLON
      !(out) 回転の経度成分

      real(8), dimension(im,jm,0:km), intent(OUT) :: xyz_RotVLAT
      !(out) 回転の緯度成分

      real(8), dimension(im,jm,0:km), intent(OUT) :: xyz_RotVRAD
      !(out) 回転の動径成分

      ! 入力ベクトル場を表すポテンシャル
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt_TORPOT
      !(in) トロイダルポテンシャル

      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt_POLPOT
      !(in) ポロイダルポテンシャル

      call wt_Potential2Vector( xyz_RotVLON,xyz_RotVLAT,xyz_RotVRAD, -wt_Lapla_wt(wt_POLPOT), wt_TORPOT)

    end subroutine wt_Potential2Rotation