!--
!----------------------------------------------------------------------
! Copyright(c) 2017 SPMDODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!表題  scee_module_fftj
!
!    spml/scee_module_fftj モジュールは平行平板間での 3 次元流体運動を
!    スペクトル法によって数値計算するための Fortran90 関数を提供する
!    ものである.
!
!    水平方向にフーリエ数変換および上下の境界壁を扱うための
!    sin/cos 変換を用いる場合のスペクトル計算のためのさまざまな
!    関数を提供する.
!
!    内部で ee_module_fftj, asc_module を用いている. 最下部ではフーリエ
!    およびチェビシェフ変換のエンジンとして ISPACK の Fortran77
!    サブルーチンを用いている.
!
!
!履歴  2017/05/26  竹広真一  作成
!
!凡例
!      データ種類と index
!        x : 水平(X)        y : 水平(Y)        z : 鉛直
!        e : フーリエ変換スペクトル
!        l : フーリエ関数スペクトル(X 方向波数)
!        m : フーリエ関数スペクトル(Y 方向波数)
!        s : sin 関数スペクトル
!        c : cos 関数スペクトル
!        a : 任意の次元
!
!        zyx : 3 次元格子点データ
!        yx  : 水平 2 次元格子点データ
!        zy  : 鉛直 2 次元格子点データ
!        zx  : 鉛直 2 次元格子点データ
!
!        zee : 水平スペクトル鉛直格子点データ
!        see : スペクトルデータ鉛直 sin 変換格子
!        cee : スペクトルデータ鉛直 cos 変換格子
!
!++
module scee_module_fftj
  !
  != scee_module_fftj
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  !    spml/scee_module_fftj モジュールは平行平板間での 3 次元流体運動を
  !    スペクトル法によって数値計算するための Fortran90 関数を提供する
  !    ものである.
  !
  !    水平方向にフーリエ数変換および上下の境界壁を扱うための
  !    sin/cos 変換を用いる場合のスペクトル計算のためのさまざまな
  !    関数を提供する.
  !
  !    内部で ee_module_fftj, asc_module を用いている. 最下部ではフーリエ
  !    およびチェビシェフ変換のエンジンとして ISPACK の Fortran77
  !    サブルーチンを用いている.
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !
  ! * 関数名の先頭 (see_, cee_, zyx_, zee_, ee_, yx_, x_, y_, z_, a_) は,
  !   返す値の形を示している.
  !   see_  :: スペクトルデータ(2 重フーリエ・sin変換)
  !   cee_  :: スペクトルデータ(2 重フーリエ・cos変換)
  !   zyx_ :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・)
  !   zee_ :: 水平スペクトル, 鉛直格子点データ
  !   e2a_ :: 1 次元化した水平スペクトル, 任意座標データ
  !   aee_ :: 任意座標データ, 水平スペクトルデータ
  !
  ! * 関数名の間の文字列(Dx, Dy, Dz, Lapla,..)
  !   は, その関数の作用を表している.
  !
  ! * 関数名の最後 (see_, cee_, zyx_, zee_, ee_, yx_, x_, y_, z_, a_) は, 
  !   入力変数の形がスペクトルデータおよび格子点データであることを示している.
  !   _see :: スペクトルデータ(2 重フーリエ・sin変換)
  !   _cee :: スペクトルデータ(2 重フーリエ・cos変換)  
  !   _zyx :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・)
  !   _zee :: 水平スペクトル, 鉛直格子点データ
  !   _e2a :: 1 次元化した水平スペクトル, 任意座標データ
  !   _aee :: 任意座標データ, 水平スペクトルデータ
  !
  !=== 各データの種類の説明
  !
  ! * zyx : 3 次元格子点データ(鉛直, 水平 2 次元)
  !   * 変数の種類と次元は real(8), dimension(0:km,0:jm-1,0:im-1).
  !   * im, jm, km はそれぞれ水平 X, Y, 鉛直 Z 座標の格子点数であり,
  !     サブルーチン tee_Initial にてあらかじめ設定しておく.
  !
  ! * see : スペクトルデータ
  !   * 変数の種類と次元は real(8), dimension(nm,-mm:mm,-lm:lm).
  !   * lm, mm は X,Y 方向最大波数, nm は sin 級数の最大次数
  !     であり, サブルーチン tee_Initial にてあらかじめ設定しておく.
  !
  ! * cee : スペクトルデータ
  !   * 変数の種類と次元は real(8), dimension(0:nm,-mm:mm,-lm:lm).
  !   * lm, mm は X,Y 方向最大波数, nm は cos 級数の最大次数
  !     であり, サブルーチン tee_Initial にてあらかじめ設定しておく.
  !
  ! * zee : 水平スペクトル, 鉛直格子点データ.
  !   * 変数の種類と次元は real(8), dimension(0:km,-mm:mm,-lm:lm).
  !
  ! * e2a : 1 次元化した水平スペクトル, 任意座標データ
  !   * 変数の種類と次元は real(8), dimension((2*mm+1)*(2*lm+1),:)
  !
  ! * aee :   任意座標データ, 水平スペクトルデータ
  !   * 変数の種類と次元は real(8), dimension(:,-mm:mm,-lm:lm).
  !
  ! * see_, cee_ で始まる関数が返す値はスペクトルデータに同じ.
  !
  ! * zyx_ で始まる関数が返す値は 3 次元格子点データに同じ.
  !
  ! * zee_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したものことである.
  !
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化
  !
  ! scee_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
  !
  !==== 座標変数
  !
  ! x_X, y_Y, z_Z                :: 格子点座標(水平 X,Y, 鉛直 Z 座標)を
  !                                 格納した1 次元配列
  ! x_X_Weight, y_X_Weight, z_Z_Weight :: 重み座標を格納した 1 次元配列
  ! zyx_X, zyx_Y, zyx_Z          :: 格子点データの水平鉛直座標(X,Y,Z)
  !                                 (格子点データ型 3 次元配列)
  ! yx_X, yx_Y                   :: 格子点データの水平座標(X,Y)
  ! zy_Z, zy_Y                   :: 格子点データの鉛直水平座標(Z,Y)
  ! zx_Z, zx_X                   :: 格子点データの鉛直水平座標(Z,X)
  !                                 (格子点データ型 2 次元配列)
  !
  !==== 基本変換
  !
  ! zyx_see, see_zyx :: スペクトルデータと 3 次元格子データの間の変換
  !                     (2 重フーリエ, sin 変換)
  ! zyx_cee, cee_zyx :: スペクトルデータと 3 次元格子データの間の変換
  !                     (2 重フーリエ, cos 変換)
  ! zyx_zee, zee_zyx :: 3 次元格子データと水平スペクトル・鉛直格子データとの間
  !                     の変換 (2 重フーリエ変換)
  ! zee_see, see_zee :: スペクトルデータと水平スペクトル・鉛直格子データとの間
  !                     の変換 (sin変換)
  ! zee_cee, cee_zee :: スペクトルデータと水平スペクトル・鉛直格子データとの間
  !                     の変換 (cos変換)
  ! ee_yx, yx_ee     :: スペクトルデータと 2 次元水平格子データの間の変換
  !                     (2 重フーリエ変換)
  ! az_as, as_az     :: 格子データと sin 変換データの間の変換を同時に複数個行う
  ! az_ac, ac_az     :: 格子データと cos 変換データの間の変換を同時に複数個行う
  ! e2a_aee, aee_e2a :: 水平スペクトル軸を 1 次元化し転置, 転置 ２次元化する.
  !
  !==== 微分
  !
  ! see_Dx_see          :: スペクトルデータに動径微分∂/∂x を作用させる
  ! cee_Dx_cee          :: スペクトルデータに動径微分∂/∂x を作用させる
  ! see_Dy_see          :: スペクトルデータに動径微分∂/∂y を作用させる
  ! cee_Dy_cee          :: スペクトルデータに動径微分∂/∂y を作用させる
  ! see_Dz_cee          :: スペクトルデータに動径微分∂/∂z を作用させる
  ! cee_Dz_see          :: スペクトルデータに動径微分∂/∂z を作用させる
  !
  ! see_Lapla_see       :: スペクトルデータにラプラシアンを作用させる
  ! cee_Lapla_cee       :: スペクトルデータにラプラシアンを作用させる
  ! see_LaplaInv_see    :: スペクトルデータに逆ラプラシアンを作用させる
  ! cee_LaplaInv_cee    :: スペクトルデータに逆ラプラシアンを作用させる
  ! see_LaplaH_see      :: スペクトルデータに水平ラプラシアンを作用させる
  ! cee_LaplaH_cee      :: スペクトルデータに水平ラプラシアンを作用させる
  ! see_LaplaHInv_see   :: スペクトルデータに逆水平ラプラシアンを作用させる
  ! cee_LaplaHInv_cee   :: スペクトルデータに逆水平ラプラシアンを作用させる
  !
  ! see_Div_zyx_zyx_zyx :: ベクトル成分である 3 つの格子データに
  !                        発散を作用させる
  ! cee_Div_zyx_zyx_zyx :: ベクトル成分である 3 つの格子データに
  !                        発散を作用させる
  !
  !==== トロイダルポロイダル計算用微分
  !
  ! see_ZRot_zyx_zyx         :: ベクトル v の渦度と動径ベクトル r の内積
  !                             z・(▽×v) を計算する
  ! cee_ZRot_zyx_zyx         :: ベクトル v の渦度と動径ベクトル r の内積
  !                             z・(▽×v) を計算する
  ! see_ZRotRot_zyx_zyx_zyx  :: ベクトルの v の z・(▽×▽×v) を計算する
  ! cee_ZRotRot_zyx_zyx_zyx  :: ベクトルの v の z・(▽×▽×v) を計算する
  ! scee_Potential2Vector_see_cee :: トロイダルポロイダルポテンシャルから
  !                                  ベクトル場を計算する
  ! scee_Potential2Vector_cee_see :: トロイダルポロイダルポテンシャルから
  !                                  ベクトル場を計算する
  ! scee_Potential2Rotation_see_cee :: トロイダルポロイダルポテンシャルで表される
  !                                    非発散ベクトル場の回転の各成分を計算する
  ! scee_Potential2Rotation_cee_see :: トロイダルポロイダルポテンシャルで表される
  !                                    非発散ベクトル場の回転の各成分を計算する
  !
  !==== ポロイダル/トロイダルモデル用スペクトル解析
  !
  ! zee_ToroidalEnergySpectrum_cee, zee_ToroidalEnergySpectrum_see ::
  !     トロイダルポテンシャルからエネルギーのフーリエ各成分を計算する
  ! zee_PoloidalEnergySpectrum_tee, zee_PoloidalEnergySpectrum_cee ::
  !     ポロイダルポテンシャルからエネルギーのフーリエ各成分を計算する
  !
  !==== 積分・平均(3 次元データ)
  !
  ! IntZYX_zyx, AvrZYX_zyx     :: 3 次元格子点データの全領域積分および平均
  ! z_IntYX_zyx, z_AvrYX_zyx   :: 3 次元格子点データの水平積分および平均
  ! y_IntZX_zyx, y_AvrZX_zyx   :: 3 次元格子点データのZX積分および平均
  ! z_IntZY_zyx, z_AvrZY_zyx   :: 3 次元格子点データのZY積分および平均
  ! zy_IntX_zyx, zy_AvrX_zyx   :: 3 次元格子点データの水平X方向積分および平均
  ! zx_IntY_zyx, zx_AvrY_zyx   :: 3 次元格子点データの水平Y方向積分および平均
  ! zx_IntZ_zyx, zx_AvrZ_zyx   :: 3 次元格子点データの鉛直方向積分および平均
  !
  !==== 積分・平均(2 次元データ)
  !
  ! IntYX_yx, AvrYX_yx :: 2 次元格子点データの水平積分および平均
  ! IntZX_zx, AvrZX_zx :: 2 次元(ZX)格子点データのZX積分および平均
  ! IntZY_zy, AvrZY_zy :: 2 次元(ZY)格子点データのZY積分および平均
  ! y_IntX_yx, y_AvrX_yx   :: 水平 2 次元格子点データのX方向積分および平均
  ! x_IntY_yx, x_AvrY_yx   :: 水平2 次元格子点データのY方向積分および平均
  ! z_IntX_zx, z_AvrX_zx   :: 2 次元(ZX)格子点データのX方向積分および平均
  ! x_IntZ_zx, x_AvrZ_zx   :: 2 次元(ZX)格子点データのZ方向積分および平均
  ! z_IntY_zy, z_AvrY_zy   :: 2 次元(ZY)格子点データのY方向積分および平均
  ! y_IntZ_zy, y_AvrZ_zy   :: 2 次元(ZY)格子点データのZ方向積分および平均
  !
  !==== 積分・平均(1 次元データ)
  !
  ! IntX_x, AvrX_x  :: 1 次元(X)格子点データのX方向積分および平均
  ! IntY_y, AvrY_y  :: 1 次元(Y)格子点データのY方向積分および平均
  ! IntZ_z, AvrZ_z  :: 1 次元(Z)格子点データのZ方向積分および平均
  !
  !==== 補間計算
  !
  ! Interpolate_cee :: スペクトルデータから任意の点の値を補間する.
  ! Interpolate_see :: スペクトルデータから任意の点の値を補間する.  
  !
  use dc_message
  use lumatrix
  use ee_module_fftj
  use asc_module, z_Z => g_X, z_Z_WEIGHT => g_X_WEIGHT, &
                  as_az => as_ag, az_as => ag_as, &
                  ac_az => ac_ag, az_ac => ag_ac, &
                  s_z => s_g, z_s => g_s, &
                  c_z => c_g, z_c => g_c, &                 
                  c_Dz_s => c_Dx_s, ac_Dz_as => ac_Dx_as, &
                  s_Dz_c => s_Dx_c, as_Dz_ac => as_Dx_ac
  implicit none
  private

  public scee_Initial

  public x_X, x_X_Weight
  public y_Y, y_Y_Weight
  public z_Z, z_Z_Weight
  public yx_X, yx_Y, zy_Z, zy_Y, zx_X, zx_Z
  public zyx_X, zyx_Y, zyx_Z
  public zee_Z
  public scee_VMiss

  public ee_yx, yx_ee
  public ac_Dz_as, as_Dz_ac, c_Dz_s, s_Dz_c
  public az_ac, az_as, ac_az, as_az, z_c, z_s, c_z, s_z
  public zyx_see, see_zyx, zee_see, see_zee
  public zyx_cee, cee_zyx, zee_cee, cee_zee
  public zyx_zee, zee_zyx 
  public e2z_e2s, e2s_e2z
  public e2z_e2c, e2c_e2z   
  public aee_e2a, e2a_aee, ee_e2, e2_ee
  public cee_Dx_cee, see_Dx_see, cee_Dy_cee, see_Dy_see, see_Dz_cee, cee_Dz_see
  public cee_Lapla_cee, see_Lapla_see
  public cee_LaplaInv_cee, see_LaplaInv_see
  public cee_LaplaH_cee, see_LaplaH_see, cee_LaplaHInv_cee, see_LaplaHInv_see
  public cee_Div_zyx_zyx_zyx, see_Div_zyx_zyx_zyx

  public zy_IntX_zyx, zx_IntY_zyx, yx_IntZ_zyx
  public x_IntZY_zyx, y_IntZX_zyx, z_IntYX_zyx
  public IntZYX_zyx

  public x_IntY_yx, y_IntX_yx, IntYX_yx
  public z_IntY_zy, y_IntZ_zy, IntZY_zy
  public z_IntX_zx, x_IntZ_zx, IntZX_zx
  public IntX_x, IntY_y, IntZ_z

  public zy_AvrX_zyx, zx_AvrY_zyx, yx_AvrZ_zyx
  public x_AvrZY_zyx, y_AvrZX_zyx, z_AvrYX_zyx
  public AvrZYX_zyx

  public x_AvrY_yx, y_AvrX_yx, AvrYX_yx
  public z_AvrY_zy, y_AvrZ_zy, AvrZY_zy
  public z_AvrX_zx, x_AvrZ_zx, AvrZX_zx
  public AvrX_x, AvrY_y, AvrZ_z

  public cee_ZRot_zyx_zyx, see_ZRot_zyx_zyx
  public cee_ZRotRot_zyx_zyx_zyx, see_ZRotRot_zyx_zyx_zyx
  public scee_Potential2vector_cee_see, scee_Potential2vector_see_cee
  public scee_Potential2Rotation_cee_see, scee_Potential2Rotation_see_cee

  public Interpolate_cee, Interpolate_see

  public zee_ToroidalEnergySpectrum_cee, zee_ToroidalEnergySpectrum_see ! nz_ToroidalEnergySpectrum_wt
  public zee_PoloidalEnergySpectrum_cee, zee_PoloidalEnergySpectrum_see ! nz_PoloidalEnergySpectrum_wt

  integer            :: im=32, jm=32, km=16  ! 格子点の設定(X,Y,Z)
  integer            :: lm=10, mm=10, nm=16  ! 切断波数の設定(水平X,Y, 鉛直Z)
  real(8)            :: xl, yl, zl           ! 領域の大きさ(水平X,Y, 鉛直Z)

  real(8), parameter :: pi=3.1415926535897932385D0

  real(8), dimension(:,:,:), allocatable :: zyx_X, zyx_Y, zyx_Z    ! 座標
  real(8), dimension(:,:),   allocatable :: zy_Z, zy_Y             ! 座標
  real(8), dimension(:,:),   allocatable :: zx_Z, zx_X             ! 座標
  real(8), dimension(:,:,:), allocatable :: zee_Z                  ! 座標

  real(8) :: scee_VMiss = -999.0        ! 欠損値

  save im, jm, km, lm, mm, nm, xl, yl, zl

contains
  !--------------- 初期化 -----------------
  subroutine scee_Initial(i,j,k,l,m,n,xmin,xmax,ymin,ymax,zmin,zmax)
     !
     ! スペクトル変換の格子点数, 波数, 各座標の範囲を設定する.
     !
     ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を
     ! しなければならない.
     !
     integer,intent(in) :: i              ! 格子点数(水平X)
     integer,intent(in) :: j              ! 格子点数(水平Y)
     integer,intent(in) :: k              ! 格子点数(鉛直Z)
     integer,intent(in) :: l              ! 切断波数(水平X波数)
     integer,intent(in) :: m              ! 切断波数(水平Y波数)
     integer,intent(in) :: n              ! 切断波数(鉛直Z波数)

     real(8),intent(in) :: xmin, xmax     ! X 方向領域
     real(8),intent(in) :: ymin, ymax     ! Y 方向領域
     real(8),intent(in) :: zmin, zmax     ! Z 方向領域

     integer            :: id

     im = i  ; jm = j ; km = k
     lm = l  ; mm = m ; nm = n
     xl = xmax - xmin
     yl = ymax - ymin
     zl = zmax - zmin

     call ee_Initial(im,jm,lm,mm,xmin,xmax,ymin,ymax,id)

     call asc_Initial(km,nm,zmin,zmax)

     allocate(zyx_X(0:km,0:jm-1,0:im-1))
     allocate(zyx_Y(0:km,0:jm-1,0:im-1))
     allocate(zyx_Z(0:km,0:jm-1,0:im-1))

     allocate(zy_Z(0:km,0:jm-1))
     allocate(zy_Y(0:km,0:jm-1))
     allocate(zx_Z(0:km,0:im-1))
     allocate(zx_X(0:km,0:im-1))

     allocate(zee_Z(0:km,-mm:mm,-lm:lm))

     zyx_X = spread(yx_X,1,km+1)
     zyx_Y = spread(yx_Y,1,km+1)
     zyx_Z = spread(spread(z_Z,2,jm),3,im)

     zy_Z = spread(z_Z,2,jm)
     zy_Y = spread(y_Y,1,km+1)

     zx_Z = spread(z_Z,2,im)
     zx_X = spread(x_X,1,km+1)

     zee_Z = spread(spread(z_Z,2,2*mm+1),3,2*lm+1)

     call MessageNotify('M','tee_initial', &
          'scee_module_fftj (2017/05/28) is initialized')

   end subroutine scee_Initial

   !--------------- 基本変換 -----------------
   
   function zyx_see(see)
     !
     ! スペクトルデータから 3 次元格子点データへ(逆)変換する.
     !
     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see
     !(in) 2 重フーリエ sin スペクトルデータ
     real(8), dimension(0:km,0:jm-1,0:im-1)           :: zyx_see
     !(out) 3 次元格子点データ

     zyx_see = zyx_zee(zee_see(see))

   end function zyx_see

   function see_zyx(zyx)
     !
     ! 3 次元格子点データからスペクトルデータへ(正)変換する.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元格子点データ
     real(8), dimension(nm,-mm:mm,-lm:lm)               :: see_zyx
     !(out) 2 重フーリエ sin スペクトルデータ

     see_zyx = see_zee(zee_zyx(zyx))

   end function see_zyx

   function zyx_cee(cee)
     !
     ! スペクトルデータから 3 次元格子点データへ(逆)変換する.
     !
     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee
     !(in) 2 重フーリエ cos スペクトルデータ
     real(8), dimension(0:km,0:jm-1,0:im-1)             :: zyx_cee
     !(out) 3 次元格子点データ

     zyx_cee = zyx_zee(zee_cee(cee))

   end function zyx_cee

   function cee_zyx(zyx)
     !
     ! 3 次元格子点データからスペクトルデータへ(正)変換する.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元格子点データ
     real(8), dimension(0:nm,-mm:mm,-lm:lm)             :: cee_zyx
     !(out) 2 重フーリエ cos スペクトルデータ

     cee_zyx = cee_zee(zee_zyx(zyx))

   end function cee_zyx

   function zyx_zee(zee)
     !
     ! 水平スペクトル・鉛直格子点データから 3 次元格子点データへ(逆)変換する.
     !
     real(8), dimension(0:km,-mm:mm,-lm:lm), intent(in) :: zee
     !(in) 2 次元水平スペクトル・鉛直格子点データ
     real(8), dimension(0:km,0:jm-1,0:im-1)             :: zyx_zee
     !(out) 3 次元格子点データ
     real(8), dimension(-mm:mm,-lm:lm)                  :: ee
     real(8), dimension(0:jm-1,0:im-1)                  :: yx

     integer :: k

     do k = 0, km
        !zyx_zee(k,:,:) = yx_ee(zee(k,:,:))
        !
        ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
        !
        ee = zee(k,-mm:mm,-lm:lm)
        yx = yx_ee(ee)
        zyx_zee(k,:,:) = yx
     enddo

   end function zyx_zee

   function zee_zyx(zyx)
     !
     ! 3 次元格子データから水平スペクトル・鉛直格子点データへ(正)変換する.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元格子点データ
     real(8), dimension(0:km,-mm:mm,-lm:lm)             :: zee_zyx
     !(out) 2 次元スペクトル・鉛直格子点データ
     real(8), dimension(-mm:mm,-lm:lm)                  :: ee
     real(8), dimension(0:jm-1,0:im-1)                  :: yx

     integer :: k

     do k = 0, km
        !zee_zyx(k,:,:) = ee_yx(zyx(k,:,:))
        !
        ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
        !
        yx = zyx(k,0:jm-1,0:im-1)
        ee = ee_yx(yx)
        zee_zyx(k,:,:) = ee
     enddo

   end function zee_zyx

   function zee_see(see)
     !
     ! スペクトルデータから水平スペクトル・鉛直格子点データへ(正)変換する.
     !
     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see
     !(in) 2 次元水平スペクトル鉛直 sin スペクトルデータ
     real(8), dimension(0:km,-mm:mm,-lm:lm)           :: zee_see
     !(out) 2 次元水平スペクトル・鉛直格子点データ

     zee_see = aee_e2a(e2z_e2s(e2a_aee(see)))

   end function zee_see

   function see_zee(zee)
     !
     ! 水平スペクトル・鉛直格子点データからスペクトルデータへ(正)変換する.
     !
     real(8), dimension(0:km,-mm:mm,-lm:lm), intent(in) :: zee
     !(in) 2 次元水平スペクトル・鉛直格子点データ
     real(8), dimension(nm,-mm:mm,-lm:lm)               :: see_zee
     !(out) 2 次元水平鉛直 sin スペクトルデータ

     see_zee = aee_e2a(e2s_e2z(e2a_aee(zee)))

   end function see_zee

   function zee_cee(cee)
     !
     ! スペクトルデータから水平スペクトル・鉛直格子点データへ(正)変換する.
     !
     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee
     !(in) 2 次元水平スペクトル鉛直 cos スペクトルデータ
     real(8), dimension(0:km,-mm:mm,-lm:lm)             :: zee_cee
     !(out) 2 次元水平スペクトル・鉛直格子点データ

     zee_cee = aee_e2a(e2z_e2c(e2a_aee(cee)))

   end function zee_cee

   function cee_zee(zee)
     !
     ! 水平スペクトル・鉛直格子点データからスペクトルデータへ(正)変換する.
     !
     real(8), dimension(0:km,-mm:mm,-lm:lm), intent(in) :: zee
     !(in) 2 次元水平スペクトル・鉛直格子点データ
     real(8), dimension(0:nm,-mm:mm,-lm:lm)             :: cee_zee
     !(out) 2 次元水平鉛直 cos スペクトルデータ

     cee_zee = aee_e2a(e2c_e2z(e2a_aee(zee)))

   end function cee_zee

   function e2z_e2s(e2s)
     !
     ! 鉛直スペクトルを格子点に変換する
     !
     real(8), dimension((2*lm+1)*(2*mm+1),nm)    :: e2s
     !(in)  1 次元化された水平スペクトル・鉛直スペクトル
     real(8), dimension((2*lm+1)*(2*mm+1),0:km)  :: e2z_e2s
     !(out) 1 次元化された水平スペクトル・鉛直格子座標

     e2z_e2s = az_as(e2s)

   end function e2z_e2s

   function e2s_e2z(e2z)
     !
     ! 鉛直格子点をスペクトルに変換する
     !
     real(8), dimension((2*lm+1)*(2*mm+1),0:km)  :: e2z
     !(in)  1 次元化された水平スペクトル・鉛直格子座標
     real(8), dimension((2*lm+1)*(2*mm+1),nm)    :: e2s_e2z
     !(out) 1 次元化された水平スペクトル・鉛直スペクトル

     e2s_e2z = as_az(e2z)

   end function e2s_e2z

   function e2z_e2c(e2c)
     !
     ! 鉛直スペクトルを格子点に変換する
     !
     real(8), dimension((2*lm+1)*(2*mm+1),0:nm)    :: e2c
     !(in)  1 次元化された水平スペクトル・鉛直スペクトル
     real(8), dimension((2*lm+1)*(2*mm+1),0:km)    :: e2z_e2c
     !(out) 1 次元化された水平スペクトル・鉛直格子座標

     e2z_e2c = az_ac(e2c)

   end function e2z_e2c

   function e2c_e2z(e2z)
     !
     ! 鉛直格子点をスペクトルに変換する
     !
     real(8), dimension((2*lm+1)*(2*mm+1),0:km)  :: e2z
     !(in)  1 次元化された水平スペクトル・鉛直格子座標
     real(8), dimension((2*lm+1)*(2*mm+1),0:nm)  :: e2c_e2z
     !(out) 1 次元化された水平スペクトル・鉛直スペクトル

     e2c_e2z = ac_az(e2z)

   end function e2c_e2z
   
   function e2a_aee(aee)
     !
     ! 水平スペクトルを 1 次元化転置する.
     !
     real(8), dimension(:,-mm:,-lm:), intent(in)        :: aee
     !(in) 任意座標・水平 2 次元スペクトルデータ dimension(:,-mm:mm,-lm:lm)
     real(8), dimension((2*lm+1)*(2*mm+1),size(aee,1))  :: e2a_aee
     !(out) 1 次元化された水平スペクトル・任意座標データ

     integer :: j,k,l,m
     !
     !e2a_aee = transpose(reshape(aee,(/size(aee,1),(2*lm+1)*(2*mm+1)/)))
     !
     ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
     !

     do k=1,size(aee,1)
        do l=-lm,lm
           do m=-mm,mm
              j = (m+mm+1)+ (2*mm+1)*(l+lm)
              e2a_aee(j,k) = aee(k,m,l)
           enddo
        enddo
     enddo

   end function e2a_aee

   function aee_e2a(e2a)
     !
     ! 水平スペクトルを転置展開する.
     !
     real(8), dimension(:,:),intent(in)                 :: e2a
     !(in) 1 次元化された水平スペクトル・任意座標データ
     !     dimmension((2*mm*1)*(2*lm*1),:)
     real(8), dimension(size(e2a,2),-mm:mm,-lm:lm)      :: aee_e2a
     !(out) 任意座標・水平 2 次元スペクトルデータ

     integer :: j,k,l,m
     !
     !aee_e2a = reshape(transpose(e2a),(/size(e2a,2),2*mm+1,2*lm+1/))
     !
     ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
     !

     do k=1,size(e2a,2)
        do l=-lm,lm
           do m=-mm,mm
              j = (m+mm+1) + (2*mm+1)*(l+lm)
              aee_e2a(k,m,l) = e2a(j,k)
           enddo
        enddo
     enddo

   end function aee_e2a

   function e2_ee(ee)
     !
     ! 水平スペクトルを 1 次元化転置する.
     !
     real(8), dimension(-mm:,-lm:), intent(in)        :: ee
     !(in) 任意座標・水平 2 次元スペクトルデータ dimension(-mm:mm,-lm:lm)
     real(8), dimension((2*lm+1)*(2*mm+1))            :: e2_ee
     !(out) 1 次元化された水平スペクトル・任意座標データ

     integer :: j,l,m

     !e2_ee = reshape(e2a_aee(reshape(ee,(/1,2*mm+1,2*lm+1/))), &
     !                (/(2*lm+1)*(2*mm+1)/))
     !
     ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
     !

     do l=-lm,lm
        do m=-mm,mm
           j = (m+mm+1)+ (2*mm+1)*(l+lm)
           e2_ee(j) = ee(m,l)
        enddo
     enddo

   end function e2_ee

   function ee_e2(e2)
     !
     ! 水平スペクトルを転置展開する.
     !
     real(8), dimension(:),intent(in)                   :: e2
     !(in) 1 次元化された水平スペクトル・任意座標データ
     !     dimmension((2*mm*1)*(2*lm*1),:)
     real(8), dimension(-mm:mm,-lm:lm)                  :: ee_e2
     !(out) 任意座標・水平 2 次元スペクトルデータ

     integer :: j,l,m

     !ee_e2 = reshape(aee_e2a(reshape(e2,(/(2*mm+1)*(2*lm+1),1/))), &
     !                (/2*mm+1,2*lm+1/))
     !
     ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
     !

     do l=-lm,lm
        do m=-mm,mm
           j = (m+mm+1) + (2*mm+1)*(l+lm)
           ee_e2(m,l) = e2(j)
        enddo
     enddo

   end function ee_e2

   !--------------- 微分計算 -----------------
   function cee_Dx_cee(cee)
     !
     ! 入力スペクトルデータに水平微分 ∂/∂x を作用する.
     !
     ! スペクトルデータのX微分とは, 対応する格子点データにX微分を
     ! 作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee
     !(in) 入力スペクトルデータ

     real(8), dimension(0:nm,-mm:mm,-lm:lm)             :: cee_Dx_cee
     !(in) X微分されたスペクトルデータ

     integer :: n

     do n=0,nm
        cee_Dx_cee(n,:,:) = ee_Dx_ee(cee(n,:,:))
     enddo

   end function cee_Dx_cee

   function see_Dx_see(see)
     !
     ! 入力スペクトルデータに水平微分 ∂/∂x を作用する.
     !
     ! スペクトルデータのX微分とは, 対応する格子点データにX微分を
     ! 作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see
     !(in) 入力スペクトルデータ

     real(8), dimension(nm,-mm:mm,-lm:lm)             :: see_Dx_see
     !(in) X微分されたスペクトルデータ

     integer :: n

     do n=1,nm
        see_Dx_see(n,:,:) = ee_Dx_ee(see(n,:,:))
     enddo

   end function see_Dx_see

   function cee_Dy_cee(cee)
     !
     ! 入力スペクトルデータに水平微分 ∂/∂y を作用する.
     !
     ! スペクトルデータのY微分とは, 対応する格子点データにY微分を
     ! 作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee
     !(in) 入力スペクトルデータ

     real(8), dimension(0:nm,-mm:mm,-lm:lm)             :: cee_Dy_cee
     !(in) X微分されたスペクトルデータ

     integer :: n

     do n=0,nm
        cee_Dy_cee(n,:,:) = ee_Dy_ee(cee(n,:,:))
     enddo

   end function cee_Dy_cee

   function see_Dy_see(see)
     !
     ! 入力スペクトルデータに水平微分 ∂/∂y を作用する.
     !
     ! スペクトルデータのY微分とは, 対応する格子点データにY微分を
     ! 作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see
     !(in) 入力スペクトルデータ

     real(8), dimension(nm,-mm:mm,-lm:lm)             :: see_Dy_see
     !(in) X微分されたスペクトルデータ

     integer :: n

     do n=1,nm
        see_Dy_see(n,:,:) = ee_Dy_ee(see(n,:,:))
     enddo

   end function see_Dy_see

   function see_Dz_cee(cee)
     !
     ! 入力スペクトルデータに鉛直微分 ∂/∂z を作用する.
     !
     ! スペクトルデータの鉛直微分とは, 対応する格子点データに鉛直微分を
     ! 作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee
     !(in) 入力スペクトルデータ

     real(8), dimension(nm,-mm:mm,-lm:lm)               :: see_Dz_cee
     !(in) 鉛直微分されたスペクトルデータ

     see_Dz_cee = aee_e2a(e2s_Dz_e2c(e2a_aee(cee)))

   end function see_Dz_cee

   function cee_Dz_see(see)
     !
     ! 入力スペクトルデータに鉛直微分 ∂/∂z を作用する.
     !
     ! スペクトルデータの鉛直微分とは, 対応する格子点データに鉛直微分を
     ! 作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see
     !(in) 入力スペクトルデータ

     real(8), dimension(0:nm,-mm:mm,-lm:lm)           :: cee_Dz_see
     !(in) 鉛直微分されたスペクトルデータ

     cee_Dz_see = aee_e2a(e2c_Dz_e2s(e2a_aee(see)))     

   end function cee_Dz_see

   function e2s_Dz_e2c(e2c)
     !
     ! 入力スペクトルデータに鉛直微分 ∂/∂z を作用する.
     !
     ! スペクトルデータの鉛直微分とは, 対応する格子点データに鉛直微分を
     ! 作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension((2*lm+1)*(2*mm+1),0:nm), intent(IN) :: e2c
     !(in) 入力スペクトルデータ

     real(8), dimension((2*lm+1)*(2*mm+1),nm)           :: e2s_Dz_e2c
     !(in) 鉛直微分されたスペクトルデータ
     
     e2s_Dz_e2c = as_Dz_ac(e2c)
   end function e2s_Dz_e2c
     
   function e2c_Dz_e2s(e2s)
     !
     ! 入力スペクトルデータに鉛直微分 ∂/∂z を作用する.
     !
     ! スペクトルデータの鉛直微分とは, 対応する格子点データに鉛直微分を
     ! 作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension((2*lm+1)*(2*mm+1),nm), intent(IN) :: e2s
     !(in) 入力スペクトルデータ

     real(8), dimension((2*lm+1)*(2*mm+1),0:nm)           :: e2c_Dz_e2s 
     !(in) 鉛直微分されたスペクトルデータ
     
     e2c_Dz_e2s = ac_Dz_as(e2s)
   end function e2c_Dz_e2s
     
   function cee_Lapla_cee(cee)
     !
     ! 入力スペクトルデータにラプラシアン
     !
     !     ▽^2 = ∂^2/∂X^2 + ∂^2/∂Y^2 + ∂^2/∂Z^2
     !
     ! を作用する.
     !
     ! スペクトルデータのラプラシアンとは, 対応する格子点データに
     ! ラプラシアンを作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee
     !(in) 2 次元球面調和函数 cos スペクトルデータ

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

     cee_Lapla_cee = cee_LaplaH_cee(cee) + cee_Dz_see(see_Dz_cee(cee))

   end function cee_Lapla_cee

   function see_Lapla_see(see)
     !
     ! 入力スペクトルデータにラプラシアン
     !
     !     ▽^2 = ∂^2/∂X^2 + ∂^2/∂Y^2 + ∂^2/∂Z^2
     !
     ! を作用する.
     !
     ! スペクトルデータのラプラシアンとは, 対応する格子点データに
     ! ラプラシアンを作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see
     !(in) 2 次元球面調和函数 sin スペクトルデータ

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

     see_Lapla_see = see_LaplaH_see(see) + see_Dz_cee(cee_Dz_see(see))

   end function see_Lapla_see

   function cee_LaplaInv_cee(cee)
     !
     ! スペクトルデータに逆ラプラシアンを作用させる. 
     !
     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee
     !(in) 2 次元球面調和函数 cos スペクトルデータ

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

     real(8), allocatable, dimension(:,:,:) :: cee_LaplaInvFactor
     logical :: first=.true.

     save first, cee_LaplaInvFactor

     if ( first ) then
        first = .false.
        allocate(cee_LaplaInvFactor(0:nm,-mm:mm,-lm:lm))
        cee_LaplaInvFactor = 1.0D0
        cee_LaplaInvFactor = 1.0D0/cee_Lapla_cee(cee_LaplaInvFactor)
        cee_LaplaInvFactor(0,0,0) = 0.0D0
     end if

     cee_LaplaInv_cee = cee_LaplaInvFactor * cee

   end function cee_LaplaInv_cee
   
   function see_LaplaInv_see(see)
     !
     ! スペクトルデータに逆ラプラシアンを作用させる. 
     !
     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see
     !(in) 2 次元球面調和函数 cos スペクトルデータ

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

     real(8), allocatable, dimension(:,:,:) :: see_LaplaInvFactor
     logical :: first=.true.

     save first, see_LaplaInvFactor

     if ( first ) then
        first = .false.
        allocate(see_LaplaInvFactor(nm,-mm:mm,-lm:lm))
        see_LaplaInvFactor = 1.0D0
        see_LaplaInvFactor = 1.0D0/see_Lapla_see(see_LaplaInvFactor)
     end if

     see_LaplaInv_see = see_LaplaInvFactor * see

   end function see_LaplaInv_see
   
   function cee_LaplaH_cee(cee)
     !
     ! 入力スペクトルデータに水平ラプラシアン
     !
     !     ▽^2_H = ∂^2/∂X^2 + ∂^2/∂Y^2
     !
     ! を作用する.
     !
     ! スペクトルデータの水平ラプラシアンとは, 対応する格子点データに
     ! 水平ラプラシアンを作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee
     !(in) 入力スペクトルデータ

     real(8), dimension(0:nm,-mm:mm,-lm:lm)             :: cee_LaplaH_cee
     !(out) 水平ラプラシアンを作用された 2 次元スペクトルデータ

     integer :: n

     do n=0,nm
        cee_LaplaH_cee(n,:,:) = ee_Lapla_ee(cee(n,:,:))
     enddo

   end function cee_LaplaH_cee

   function see_LaplaH_see(see)
     !
     ! 入力スペクトルデータに水平ラプラシアン
     !
     !     ▽^2_H = ∂^2/∂X^2 + ∂^2/∂Y^2
     !
     ! を作用する.
     !
     ! スペクトルデータの水平ラプラシアンとは, 対応する格子点データに
     ! 水平ラプラシアンを作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see
     !(in) 入力スペクトルデータ

     real(8), dimension(nm,-mm:mm,-lm:lm)             :: see_LaplaH_see
     !(out) 水平ラプラシアンを作用された 2 次元スペクトルデータ

     integer :: n

     do n=1,nm
        see_LaplaH_see(n,:,:) = ee_Lapla_ee(see(n,:,:))
     enddo

   end function see_LaplaH_see

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

     real(8), dimension(0:nm,-mm:mm,-lm:lm)             :: cee_LaplaHInv_cee
     !(out) 逆水平ラプラシアンを作用された 2 次元スペクトルデータ

     integer :: n

     do n=0,nm
        cee_LaplaHInv_cee(n,:,:) = ee_LaplaInv_ee(cee(n,:,:))
     enddo

   end function cee_LaplaHInv_cee

   function see_LaplaHInv_see(see)
     !
     ! 入力スペクトルデータに逆水平ラプラシアン
     !
     !     ▽^-2_H = (∂^2/∂X^2 + ∂^2/∂Y^2)^-1
     !
     ! を作用する.
     !
     ! スペクトルデータの逆水平ラプラシアンとは, 対応する格子点データに
     ! 逆水平ラプラシアンを作用させたデータのスペクトル変換のことである.
     !
     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see
     !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

     real(8), dimension(nm,-mm:mm,-lm:lm)             :: see_LaplaHInv_see
     !(out) 逆水平ラプラシアンを作用された 2 次元スペクトルデータ

     integer :: n

     do n=1,nm
        see_LaplaHInv_see(n,:,:) = ee_LaplaInv_ee(see(n,:,:))
     enddo

   end function see_LaplaHInv_see

   function see_Div_zyx_zyx_zyx(zyx_VX,zyx_VY,zyx_VZ)
     !
     ! べクトル成分である 3 つの格子データに発散を作用させた
     ! スペクトルデータを返す.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX
     !(in) ベクトル場のX成分
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY
     !(in) ベクトル場のY成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VZ
     !(in) ベクトル場のZ成分

     real(8), dimension(nm,-mm:mm,-lm:lm)     :: see_Div_zyx_zyx_zyx
     !(out) ベクトル場の発散

     see_Div_zyx_zyx_zyx =   see_Dx_see(see_zyx(zyx_VX)) &
          + see_Dy_see(see_zyx(zyx_VY)) &
          + see_Dz_cee(cee_zyx(zyx_VZ))

   end function see_Div_zyx_zyx_zyx

   function cee_Div_zyx_zyx_zyx(zyx_VX,zyx_VY,zyx_VZ)
     !
     ! べクトル成分である 3 つの格子データに発散を作用させた
     ! スペクトルデータを返す.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX
     !(in) ベクトル場のX成分
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY
     !(in) ベクトル場のY成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VZ
     !(in) ベクトル場のZ成分

     real(8), dimension(0:nm,-mm:mm,-lm:lm)     :: cee_Div_zyx_zyx_zyx
     !(out) ベクトル場の発散

     cee_Div_zyx_zyx_zyx =   cee_Dx_cee(cee_zyx(zyx_VX)) &
          + cee_Dy_cee(cee_zyx(zyx_VY)) &
          + cee_Dz_see(see_zyx(zyx_VZ))

   end function cee_Div_zyx_zyx_zyx

   !--------------- 積分計算 -----------------
   !----(入力データ zyx)---
   function zy_IntX_zyx(zyx)  ! X積分
     !
     ! 3 次元格子点データのX方向積分.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元格子点データ

     real(8), dimension(0:km,0:jm-1)  :: zy_IntX_zyx
     !(out) X方向積分された 2 次元ZY格子点データ

     integer :: i

     zy_IntX_zyx = 0.0d0
     do i=0,im-1
        zy_IntX_zyx(:,:) = zy_IntX_zyx(:,:) &
             + zyx(:,:,i) * x_X_Weight(i)
     enddo
   end function zy_IntX_zyx

   function zx_IntY_zyx(zyx)
     !
     ! 3 次元格子点データのY方向域積分.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元格子点データ

     real(8), dimension(0:km,0:im-1)  :: zx_IntY_zyx
     !(out) Y積分された 2 次元ZY格子点データ.

     integer :: j

     zx_IntY_zyx = 0.0d0
     do j=0,jm-1
        zx_IntY_zyx(:,:) = zx_IntY_zyx(:,:) &
             + zyx(:,j,:) * y_Y_Weight(j)
     enddo
   end function zx_IntY_zyx

   function yx_IntZ_zyx(zyx)  ! Z積分
     !
     ! 3 次元格子点データのZ方向積分.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元格子点データ

     real(8), dimension(0:jm-1,0:im-1)  :: yx_IntZ_zyx
     !(out) Z積分された 2 次元YX(水平)格子点データ

     integer :: k

     yx_IntZ_zyx = 0.0d0
     do k=0,km
        yx_IntZ_zyx(:,:) = yx_IntZ_zyx(:,:) &
             + zyx(k,:,:) * z_Z_Weight(k)
     enddo
   end function yx_IntZ_zyx

   function x_IntZY_zyx(zyx)
     !
     ! 3 次元格子点データのZY積分
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元ZYX格子点データ

     real(8), dimension(0:im-1)     :: x_IntZY_zyx
     !(out) ZY(子午面)積分された 1 次元X格子点データ

     integer :: j, k

     x_IntZY_zyx = 0.0D0
     do j=0,jm-1
        do k=0,km
           x_IntZY_zyx = x_IntZY_zyx &
                + zyx(k,j,:) * y_Y_Weight(j) * z_Z_Weight(k)
        enddo
     enddo
   end function x_IntZY_zyx

   function y_IntZX_zyx(zyx)
     !
     ! 3 次元格子点データのZX積分.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元ZYX格子点データ

     real(8), dimension(0:jm-1)       :: y_IntZX_zyx
     !(out) ZX積分された 1 次元Y格子点データ

     integer :: i, k

     y_IntZX_zyx = 0
     do i=0,im-1
        do k=0,km
           y_IntZX_zyx = y_IntZX_zyx &
                + zyx(k,:,i) * x_X_Weight(i) * z_Z_Weight(k)
        enddo
     enddo
   end function y_IntZX_zyx

   function z_IntYX_zyx(zyx)  ! YX(水平)積分
     !
     ! 3 次元格子点データのYX(水平)積分
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元ZYX格子点データ

     real(8), dimension(0:km)     :: z_IntYX_zyx
     !(out) YX(水平)積分された 1 次元Z格子点データ

     integer :: i, j

     z_IntYX_zyx = 0
     do j=0,jm-1
        do i=0,im-1
           z_IntYX_zyx = z_IntYX_zyx &
                + zyx(:,j,i) * x_X_Weight(i) * y_Y_Weight(j)
        enddo
     enddo
   end function z_IntYX_zyx

   function IntZYX_zyx(zyx) ! ZYX(全領域)積分
     !
     ! 3 次元格子点データのZYX(全領域)積分
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元格子点データ

     real(8)                     :: IntZYX_zyx
     !(out) 全領域積分値

     integer :: i, j, k

     IntZYX_zyx = 0
     do i=0,im-1
        do j=0,jm-1
           do k=0,km
              IntZYX_zyx = IntZYX_zyx &
                   + zyx(k,j,i) * x_X_Weight(i) &
                   * y_Y_Weight(j) * z_Z_Weight(k)
           enddo
        enddo
     enddo
   end function IntZYX_zyx

   !----(入力データ zy)---
   function z_IntY_zy(zy)  ! Y積分
     !
     ! 2 次元(ZY)格子点データのY方向積分.
     !
     real(8), dimension(0:km,0:jm-1), intent(in) :: zy
     !(in) 2 次元ZY(子午面)格子点データ

     real(8), dimension(0:km)  :: z_IntY_zy
     !(out) Y積分された 1 次元Z格子点データ

     integer :: j

     z_IntY_zy = 0.0d0
     do j=0,jm-1
        z_IntY_zy(:) = z_IntY_zy(:) + zy(:,j) * y_Y_Weight(j)
     enddo
   end function z_IntY_zy

   function y_IntZ_zy(zy)  ! Z積分
     !
     ! 2 次元(ZY)格子点データのZ方向域積分.
     !
     real(8), dimension(0:km,0:jm-1), intent(in) :: zy
     !(in) 2 次元ZY格子点データ

     real(8), dimension(0:jm-1)  :: y_IntZ_zy
     !(out) Z積分された 1 次元Y格子点データ

     integer :: k

     y_IntZ_zy = 0.0d0
     do k=0,km
        y_IntZ_zy(:) = y_IntZ_zy(:) &
             + zy(k,:) * z_Z_Weight(k)
     enddo
   end function y_IntZ_zy

   function IntZY_zy(zy)
     !
     ! 2 次元(ZY)格子点データのZY積分
     !
     real(8), dimension(0:km,0:jm-1), intent(in) :: zy
     !(in) 2 次元ZY(子午面)格子点データ

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

     IntZY_zy = 0
     do j=0,jm-1
        do k=0,km
           IntZY_zy = IntZY_zy &
                + zy(k,j) * y_Y_Weight(j) * z_Z_Weight(k)
        enddo
     enddo
   end function IntZY_zy

   !----(入力データ zx)---
   function z_IntX_zx(zx)
     !
     ! 2 次元(ZX)格子点データのX方向積分.
     !
     real(8), dimension(0:km,0:im-1), intent(in) :: zx
     !(in) 2 次元ZY格子点データ

     real(8), dimension(0:km)  :: z_IntX_zx
     !(out) X積分された 1 次元Z格子点データ

     integer :: i

     z_IntX_zx = 0.0d0
     do i=0,im-1
        z_IntX_zx(:) = z_IntX_zx(:) + zx(:,i) * x_X_Weight(i)
     enddo

   end function z_IntX_zx

   function x_IntZ_zx(zx)
     !
     ! 2 次元(ZX)格子点データのZ方向積分.
     !
     real(8), dimension(0:km,0:im-1), intent(in) :: zx
     !(in) 2 次元ZY格子点データ

     real(8), dimension(0:im-1)  :: x_IntZ_zx
     !(out) Z積分された 1 次元X格子点データ

     integer :: k

     x_IntZ_zx = 0.0d0
     do k=0,km
        x_IntZ_zx(:) = x_IntZ_zx(:) &
             + zx(k,:) * z_Z_Weight(k)
     enddo

   end function x_IntZ_zx

   function IntZX_zx(zx)  ! ZX積分
     !
     ! 2 次元(ZX)格子点データのZX積分
     !
     real(8), dimension(0:km,0:im-1), intent(in) :: zx
     !(in) 2 次元ZY格子点データ

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

     integer :: i, k

     IntZX_zx = 0
     do i=0,im-1
        do k=0,km
           IntZX_zx = IntZX_zx &
                + zx(k,i) * x_X_Weight(i) * z_Z_Weight(k)
        enddo
     enddo
   end function IntZX_zx

   !----(入力データ z)---
   function IntZ_z(z)  ! Z積分
     !
     ! 1 次元(Z)格子点データのZ方向積分.
     !
     real(8), dimension(0:km), intent(in) :: z
     !(in) 1 次元Z格子点データ

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

     integer :: k

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

   !--------------- 平均計算 -----------------
   !----(入力データ zyx)---
   function zy_AvrX_zyx(zyx)  ! X積分
     !
     ! 3 次元格子点データのX方向平均.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元ZYX格子点データ

     real(8), dimension(0:km,0:jm-1)  :: zy_AvrX_zyx
     !(out) X方向平均された 2 次元子午面格子点データ

     zy_AvrX_zyx = zy_IntX_zyx(zyx)/sum(x_X_Weight)

   end function zy_AvrX_zyx

   function zx_AvrY_zyx(zyx)  ! Y平均
     !
     ! 3 次元格子点データのY方向域平均.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元格子点データ

     real(8), dimension(0:km,0:im-1)  :: zx_AvrY_zyx
     !(out) Y平均された 2 次元ZX格子点データ

     zx_AvrY_zyx = zx_IntY_zyx(zyx)/sum(y_Y_Weight)

   end function zx_AvrY_zyx

   function yx_AvrZ_zyx(zyx)
     !
     ! 3 次元格子点データのZ方向平均.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元格子点データ

     real(8), dimension(0:jm-1,0:im-1)  :: yx_AvrZ_zyx
     !(out) Z平均された 2 次元YX(水平)格子点データ

     yx_AvrZ_zyx = yx_IntZ_zyx(zyx)/sum(z_Z_Weight)

   end function yx_AvrZ_zyx

   function x_AvrZY_zyx(zyx)  ! ZY積分
     !
     ! 3 次元格子点データのZY平均
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元ZYX格子点データ

     real(8), dimension(0:im-1)     :: x_AvrZY_zyx
     !(out) ZY平均された 1 次元X格子点データ

     x_AvrZY_zyx = x_IntZY_zyx(zyx) &
          /( sum(y_Y_Weight)*sum(z_Z_Weight) )

   end function x_AvrZY_zyx

   function y_AvrZX_zyx(zyx)  ! ZX積分
     !
     ! 3 次元格子点データのZX平均.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元格子点データ

     real(8), dimension(0:jm-1)       :: y_AvrZX_zyx
     !(out) ZX平均された 1 次元Y格子点データ

     y_AvrZX_zyx = y_IntZX_zyx(zyx) &
          /(sum(x_X_Weight)*sum(z_Z_Weight))

   end function y_AvrZX_zyx

   function z_AvrYX_zyx(zyx)  ! YX(水平)積分
     !
     ! 3 次元格子点データのYX(水平)積分
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元格子点データ

     real(8), dimension(0:km)     :: z_AvrYX_zyx
     !(out) YX(水平)平均された 1 次元Z格子点データ

     z_AvrYX_zyx = z_IntYX_zyx(zyx) &
          /(sum(x_X_Weight)*sum(y_Y_Weight))

   end function z_AvrYX_zyx

   function AvrZYX_zyx(zyx) ! ZYX(全領域)積分
     !
     ! 3 次元格子点データのZYX(全領域)積分
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
     !(in) 3 次元ZYX格子点データ

     real(8)                     :: AvrZYX_zyx
     !(out) 全領域平均値

     AvrZYX_zyx = IntZYX_zyx(zyx) &
          /(sum(x_X_Weight)*sum(y_Y_Weight) * sum(z_Z_Weight))

   end function AvrZYX_zyx

   !----(入力データ zy)---
   function z_AvrY_zy(zy)
     !
     ! 2 次元(ZY)格子点データのY方向平均.
     !
     real(8), dimension(0:km,0:jm-1), intent(in) :: zy
     !(in) 2 次元ZY格子点データ

     real(8), dimension(0:km)  :: z_AvrY_zy
     !(out) Y平均された 1 次元Z格子点データ

     z_AvrY_zy = z_IntY_zy(zy)/sum(y_Y_Weight)

   end function z_AvrY_zy

   function y_AvrZ_zy(zy)
     !
     ! 2 次元(ZY)格子点データのZ方向平均.
     !
     real(8), dimension(0:km,0:jm-1), intent(in) :: zy
     !(in) 2 次元ZY格子点データ

     real(8), dimension(0:jm-1)  :: y_AvrZ_zy
     !(out) Z平均された 1 次元Y格子点データ

     y_AvrZ_zy = y_IntZ_zy(zy)/sum(z_Z_Weight)

   end function y_AvrZ_zy

   function AvrZY_zy(zy)  ! ZY平均
     !
     ! 2 次元(ZY)格子点データのZY平均
     !
     real(8), dimension(0:km,0:jm-1), intent(in) :: zy
     !(in) 2 次元ZY(子午面)格子点データ

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

     AvrZY_zy = IntZY_zy(zy)/(sum(y_Y_Weight)*sum(z_Z_Weight))

   end function AvrZY_zy

   !----(入力データ zx)---
   function z_AvrX_zx(zx)  ! X積分
     !
     ! 2 次元(ZX)格子点データのX方向平均.
     !
     real(8), dimension(0:km,0:im-1), intent(in) :: zx
     !(in) 2 次元ZY格子点データ

     real(8), dimension(0:km)  :: z_AvrX_zx
     !(out) X平均された 1 次元Z格子点データ

     z_AvrX_zx = z_IntX_zx(zx)/sum(x_X_Weight)

   end function z_AvrX_zx

   function x_AvrZ_zx(zx)  ! Z平均
     !
     ! 2 次元(ZX)格子点データのZ方向平均.
     !
     real(8), dimension(0:km,0:im-1), intent(in) :: zx
     !(in) 2 次元ZY格子点データ

     real(8), dimension(0:im-1)  :: x_AvrZ_zx
     !(out) Z平均された 1 次元X格子点データ

     x_AvrZ_zx = x_IntZ_zx(zx)/sum(z_Z_Weight)

   end function x_AvrZ_zx

   function AvrZX_zx(zx)  ! ZX積分
     !
     ! 2 次元(ZX)格子点データのZX平均
     !
     real(8), dimension(0:km,0:im-1), intent(in) :: zx
     ! (in) 2 次元格子点データ
     real(8)                                 :: AvrZX_zx
     ! 平均値

     AvrZX_zx = IntZX_zx(zx)/(sum(x_X_Weight)*sum(z_Z_Weight))

   end function AvrZX_zx

   !----(入力データ z)---
   function AvrZ_z(z)
     !
     ! 1 次元(Z)格子点データのZ方向平均.
     !
     real(8), dimension(0:km), intent(in) :: z
     !(in) 1 次元Z格子点データ
     real(8)                              :: AvrZ_z
     !(out) 平均値

     AvrZ_z = IntZ_z(z)/sum(z_Z_Weight)

   end function AvrZ_z


   !--------------- ポロイダル/トロイダルモデル用微分 -----------------
!!$   function wt_KxRGrad_wt(wt)
!!$     !
!!$     ! 入力スペクトルデータにX微分 k×r・▽ = ∂/∂λを作用する.
!!$     !
!!$     real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
!!$     !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$     real(8), dimension(-lm:lm,-mm,mm,0:nm)             :: wt_KxRGrad_wt
!!$     !(out) X微分を作用された 2 次元スペクトルデータ
!!$
!!$     wt_KxRGrad_wt =  wa_Dlon_wa(wt)
!!$
!!$   end function wt_KxRGrad_wt
!!$
!!$   function zyx_KGrad_wt(wt)    ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r
!!$     !
!!$     ! 入力スペクトルデータに対応する格子データに軸方向微分
!!$     !
!!$     !    k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r
!!$     !
!!$     ! を作用させた格子データが返される.
!!$     ! ここでベクトル k は球の中心から北極向きの単位ベクトルである.
!!$     !
!!$     real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
!!$     !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$     real(8), dimension(0:km,0:jm-1,0:im-1)                     :: yzx_KGrad_wt
!!$     !(out) 軸方向微分を作用された 2 次元スペクトルデータ
!!$
!!$     xzy_KGrad_wt =  cos(xyz_Y)*xyz_GradY_wt(wt) &
!!$          + sin(xyz_Y)*xyz_wt(wt_Drad_wt(wt))
!!$
!!$   end function zyx_KGrad_wt
!!$
!!$   function wt_L2_wt(wt)
!!$     !
!!$     ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する.
!!$     !
!!$     ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる.
!!$     !  入力スペクトルデ ータに対応する格子点データに演算子
!!$     !
!!$     !     L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ)
!!$     !
!!$     ! を作用させたデータのスペクトル変換が返される.
!!$     !
!!$     real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
!!$     !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$     real(8), dimension(-lm:lm,-mm,mm,0:nm)             :: wt_L2_wt
!!$     !(out) L^2 演算子を作用された 2 次元スペクトルデータ
!!$
!!$     wt_L2_wt = -wa_Lapla_wa(wt)
!!$
!!$   end function wt_L2_wt
!!$
!!$   function wt_L2Inv_wt(wt)
!!$     !
!!$     ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を
!!$     ! 作用する.
!!$     !
!!$     ! スペクトルデータに L^2 演算子を作用させる関数 wt_L2_wt の逆計算を
!!$     ! 行う関数である.
!!$     !
!!$     real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
!!$     !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$     real(8), dimension(-lm:lm,-mm,mm,0:nm)             :: wt_L2Inv_wt
!!$     !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ
!!$
!!$     wt_L2Inv_wt = -wa_LaplaInv_wa(wt)
!!$
!!$   end function wt_L2Inv_wt
!!$
!!$   function wt_QOperator_wt(wt)
!!$     !
!!$     ! 入力スペクトルデータに対応する格子点データに演算子
!!$     !
!!$     !    Q=(k・▽-1/2(L2 k・▽+ k・▽L2))
!!$     !
!!$     ! を作用させたデータのスペクトル変換が返される.
!!$     !
!!$     real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
!!$     !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$     real(8), dimension(-lm:lm,-mm,mm,0:nm)             :: wt_QOperator_wt
!!$     !(out) Q 演算子を作用された 2 次元スペクトルデータ
!!$
!!$     wt_QOperator_wt = &
!!$          wt_xyz(xyz_KGrad_wt(wt) - xyz_KGrad_wt(wt_L2_wt(wt))/2) &
!!$          - wt_L2_wt(wt_xyz(xyz_KGrad_wt(wt)))/2
!!$
!!$   end function wt_QOperator_wt

   function cee_ZRot_zyx_zyx(zyx_VX,zyx_VY)  ! z・(▽×v)
     !
     ! ベクトルの渦度とZベクトルの内積 z・(▽×v) を計算する.
     !
     ! 第 1, 2 引数(v[x], v[y])がそれぞれベクトルのX成分, Y成分を表す.
     !
     !    z・(▽×v) = ∂v[y]/∂x - ∂(v[x])/∂y
     !
     ! のスペクトル データが返される.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX
     !(in) ベクトルのX成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY
     !(in) ベクトルのY成分

     real(8), dimension(0:nm,-mm:mm,-lm:lm)     :: cee_ZRot_zyx_zyx
     !(out) ベクトルの渦度とZベクトルの内積

     cee_ZRot_zyx_zyx = cee_DX_cee(cee_zyx(zyx_VY)) &
          - cee_DY_cee(cee_zyx(zyx_VX))

   end function cee_ZRot_zyx_zyx

   function see_ZRot_zyx_zyx(zyx_VX,zyx_VY)  ! z・(▽×v)
     !
     ! ベクトルの渦度とZベクトルの内積 z・(▽×v) を計算する.
     !
     ! 第 1, 2 引数(v[x], v[y])がそれぞれベクトルのX成分, Y成分を表す.
     !
     !    z・(▽×v) = ∂v[y]/∂x - ∂(v[x])/∂y
     !
     ! のスペクトル データが返される.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX
     !(in) ベクトルのX成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY
     !(in) ベクトルのY成分

     real(8), dimension(nm,-mm:mm,-lm:lm)     :: see_ZRot_zyx_zyx
     !(out) ベクトルの渦度とZベクトルの内積

     see_ZRot_zyx_zyx = see_DX_see(see_zyx(zyx_VY)) &
          - see_DY_see(see_zyx(zyx_VX))

   end function see_ZRot_zyx_zyx

   function cee_ZRotRot_zyx_zyx_zyx(zyx_VX,zyx_VY,zyx_VZ)
     !
     ! ベクトル v に対して z・(▽×▽×v) を計算する.
     !
     ! 第 1, 2, 3 引数(v[x], v[y], v[z])がそれぞれベクトルのX成分,
     ! Y成分, Z成分を表す.
     !
     !    z・(▽×▽×v)  = ∂/∂z (∂v[x]/∂x + ∂(v[y])/∂y ) )
     !                    - ▽_H^2 v[z]
     !
     ! のスペクトルデータが返される.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX
     !(in) ベクトルのX成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY
     !(in) ベクトルのY成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VZ
     !(in) ベクトルのZ成分

     real(8), dimension(0:nm,-mm:mm,-lm:lm)     :: cee_ZRotRot_zyx_zyx_zyx
     !(out) ベクトル v の z・(▽×▽×v)

     cee_ZRotRot_zyx_zyx_zyx = &
          cee_DZ_see(   see_DX_see(see_zyx(zyx_VX))   &
                      + see_DY_see(see_zyx(zyx_VY)) ) &
        - cee_LaplaH_cee(cee_zyx(zyx_VZ))

   end function cee_ZRotRot_zyx_zyx_zyx

   function see_ZRotRot_zyx_zyx_zyx(zyx_VX,zyx_VY,zyx_VZ)
     !
     ! ベクトル v に対して z・(▽×▽×v) を計算する.
     !
     ! 第 1, 2, 3 引数(v[x], v[y], v[z])がそれぞれベクトルのX成分,
     ! Y成分, Z成分を表す.
     !
     !    z・(▽×▽×v)  = ∂/∂z (∂v[x]/∂x + ∂(v[y])/∂y ) )
     !                    - ▽_H^2 v[z]
     !
     ! のスペクトルデータが返される.
     !
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX
     !(in) ベクトルのX成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY
     !(in) ベクトルのY成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VZ
     !(in) ベクトルのZ成分

     real(8), dimension(nm,-mm:mm,-lm:lm)     :: see_ZRotRot_zyx_zyx_zyx
     !(out) ベクトル v の z・(▽×▽×v)

     see_ZRotRot_zyx_zyx_zyx = &
          see_DZ_cee(   cee_DX_cee(cee_zyx(zyx_VX))   &
                      + cee_DY_cee(cee_zyx(zyx_VY)) ) &
        - see_LaplaH_see(see_zyx(zyx_VZ))

   end function see_ZRotRot_zyx_zyx_zyx

   subroutine scee_Potential2Vector_cee_see(&
        zyx_VX,zyx_VY,zyx_VZ,cee_Torvel,see_Polvel)
     !
     ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
     !
     !     v = ▽x(Ψz) + ▽x▽x(Φz)
     !
     ! の各成分を計算する
     !
     real(8), dimension(0:km,0:jm-1,0:im-1)     :: zyx_VX
     !(out) ベクトル場のX成分

     real(8), dimension(0:km,0:jm-1,0:im-1)     :: zyx_VY
     !(out) ベクトル場のY成分

     real(8), dimension(0:km,0:jm-1,0:im-1)     :: zyx_VZ
     !(out) ベクトル場のZ成分

     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee_Torvel
     !(in) トロイダルポテンシャル

     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in)   :: see_Polvel
     !(in) ポロイダルポテンシャル

     zyx_VX = zyx_cee(   cee_DY_cee(cee_Torvel) &
                       + cee_DX_cee(cee_DZ_see(see_Polvel))  )
     zyx_VY = zyx_cee( - cee_DX_cee(cee_Torvel) &
                       + cee_DY_cee(cee_DZ_see(see_Polvel))  )
     zyx_VZ = -zyx_see(see_LaplaH_see(see_Polvel))

   end subroutine scee_Potential2Vector_cee_see

   subroutine scee_Potential2Vector_see_cee(&
        zyx_VX,zyx_VY,zyx_VZ,see_Torvel,cee_Polvel)
     !
     ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
     !
     !     v = ▽x(Ψz) + ▽x▽x(Φz)
     !
     ! の各成分を計算する
     !
     real(8), dimension(0:km,0:jm-1,0:im-1)     :: zyx_VX
     !(out) ベクトル場のX成分

     real(8), dimension(0:km,0:jm-1,0:im-1)     :: zyx_VY
     !(out) ベクトル場のY成分

     real(8), dimension(0:km,0:jm-1,0:im-1)     :: zyx_VZ
     !(out) ベクトル場のZ成分

     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in)   :: see_Torvel
     !(in) トロイダルポテンシャル

     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee_Polvel
     !(in) ポロイダルポテンシャル

     zyx_VX = zyx_see(   see_DY_see(see_Torvel) &
                       + see_DX_see(see_DZ_cee(cee_Polvel))  )
     zyx_VY = zyx_see( - see_DX_see(see_Torvel) &
                       + see_DY_see(see_DZ_cee(cee_Polvel))  )
     zyx_VZ = -zyx_cee(cee_LaplaH_cee(cee_Polvel))

   end subroutine scee_Potential2Vector_see_cee

   subroutine scee_Potential2Rotation_cee_see(&
        zyx_RotVX,zyx_RotVY,zyx_RotVZ,cee_Torvel,see_Polvel)
     !
     ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
     !
     !     v = ▽x(Ψz) + ▽x▽x(Φz)
     !
     ! に対して, その回転
     !
     !     ▽xv = ▽x▽x(Ψz) + ▽x▽x▽x(Φz) = ▽x▽x(Ψz) - ▽x((▽^2Φ)z)
     !
     ! を計算する.

     ! ベクトル場の回転
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVX
     !(out) 回転のX成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVY
     !(out) 回転のY成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVZ
     !(out) 回転のZ成分

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

     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in)   :: see_Polvel
     !(in) ポロイダルポテンシャル


     call scee_Potential2Vector_see_cee(&
          zyx_RotVX,zyx_RotVY,zyx_RotVZ,-see_Lapla_see(see_Polvel),cee_Torvel)
     
   end subroutine scee_Potential2Rotation_cee_see

   subroutine scee_Potential2Rotation_see_cee(&
        zyx_RotVX,zyx_RotVY,zyx_RotVZ,see_Torvel,cee_Polvel)
     !
     ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
     !
     !     v = ▽x(Ψz) + ▽x▽x(Φz)
     !
     ! に対して, その回転
     !
     !     ▽xv = ▽x▽x(Ψz) + ▽x▽x▽x(Φz) = ▽x▽x(Ψz) - ▽x((▽^2Φ)z)
     !
     ! を計算する.

     ! ベクトル場の回転
     real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVX
     !(out) 回転のX成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVY
     !(out) 回転のY成分

     real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVZ
     !(out) 回転のZ成分

     ! 入力ベクトル場を表すポテンシャル
     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in)   :: see_Torvel
     !(in) トロイダルポテンシャル

     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee_Polvel
     !(in) ポロイダルポテンシャル


     call scee_Potential2Vector_cee_see(&
          zyx_RotVX,zyx_RotVY,zyx_RotVZ,-cee_Lapla_cee(cee_Polvel),see_Torvel)
     
   end subroutine scee_Potential2Rotation_see_cee

   !--------------- 補間計算 -----------------
   function Interpolate_cee(cee_data,x,y,z)
     !
     ! (x,y,z) における関数値を
     ! そのスペクトル係数 tee_data から補間計算する
     !
     real(8), intent(IN) :: cee_data(0:nm,-mm:mm,-lm:lm)  ! スペクトルデータ
     real(8), intent(IN) :: x                             ! 補間する位置(X)
     real(8), intent(IN) :: y                             ! 補間する位置(Y)
     real(8), intent(IN) :: z                             ! 補間する位置(Z)
     real(8) :: Interpolate_cee                           ! 補間した値

     Interpolate_cee = &
          Interpolate_ee(ee_e2(a_Interpolate_ac(e2a_aee(cee_data),z)),x,y)

   end function Interpolate_cee

   function Interpolate_see(see_data,x,y,z)
     !
     ! (x,y,z) における関数値を
     ! そのスペクトル係数 tee_data から補間計算する
     !
     real(8), intent(IN) :: see_data(nm,-mm:mm,-lm:lm)  ! スペクトルデータ
     real(8), intent(IN) :: x                           ! 補間する位置(X)
     real(8), intent(IN) :: y                           ! 補間する位置(Y)
     real(8), intent(IN) :: z                           ! 補間する位置(Z)
     real(8) :: Interpolate_see                         ! 補間した値

     Interpolate_see = &
          Interpolate_ee(ee_e2(a_Interpolate_as(e2a_aee(see_data),z)),x,y)

   end function Interpolate_see

   
   !--------------- ポロイダル/トロイダルモデル用スペクトル解析 ----------------
   
   function zee_ToroidalEnergySpectrum_cee(cee_TORPOT)
     !
     ! トロイダルポテンシャルから, トロイダルエネルギーの
     ! 水平 X 波数 l, Y  波数 m の各成分を計算する
     !
     !  * X 波数 l, Y 波数 m のトロイダルポテンシャルのスペクトル成分
     !    ψ(l,m,z)から水平 X 波数 l, Y 波数 m 成分のトロイダルエネルギー
     !    スペクトルは  (1/2)(l**2+m**2)|ψ(l,m,z)|^2  と計算される.
     !
     !  * 全てのエネルギースペクトル成分の和をZ積分したものが領域内での
     !    水平平均エネルギーに等しい.
     !
     real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee_TORPOT
     !(in) トロイダルポテンシャル

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

     real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_Work
     integer :: l, m

     zee_Work = zee_cee(cee_Torpot)
     do l=-lm,lm
        do m=-mm,mm
           zee_ToroidalEnergySpectrum_cee(:,m,l) &
                = 0.5 * ((2*PI*l/xl)**2+(2*PI*m/yl)**2) &
                * ( zee_Work(:,m,l)**2 + zee_Work(:,-m,-l)**2 )
        enddo
     enddo

   end function zee_ToroidalEnergySpectrum_cee

   function zee_ToroidalEnergySpectrum_see(see_TORPOT)
     !
     ! トロイダルポテンシャルから, トロイダルエネルギーの
     ! 水平 X 波数 l, Y  波数 m の各成分を計算する
     !
     !  * X 波数 l, Y 波数 m のトロイダルポテンシャルのスペクトル成分
     !    ψ(l,m,z)から水平 X 波数 l, Y 波数 m 成分のトロイダルエネルギー
     !    スペクトルは  (1/2)(l**2+m**2)|ψ(l,m,z)|^2  と計算される.
     !
     !  * 全てのエネルギースペクトル成分の和をZ積分したものが領域内での
     !    水平平均エネルギーに等しい.
     !
     real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see_TORPOT
     !(in) トロイダルポテンシャル

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

     real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_Work
     integer :: l, m

     zee_Work = zee_see(see_Torpot)
     do l=-lm,lm
        do m=-mm,mm
           zee_ToroidalEnergySpectrum_see(:,m,l) &
                = 0.5 * ((2*PI*l/xl)**2+(2*PI*m/yl)**2) &
                * ( zee_Work(:,m,l)**2 + zee_Work(:,-m,-l)**2 )
        enddo
     enddo

   end function zee_ToroidalEnergySpectrum_see

!!$    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 と計算される.
!!$      !
!!$      ! * 全てのエネルギースペクトル成分の和をZ積分したもの(r^2の重み無し)
!!$      !    が球殻内での全エネルギーに等しい.
!!$      !
!!$      real(8), dimension(0:nm,-mm:mm,-lm: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_Z**2 &
!!$                * sum(wz_Data(l_nm(n,(/(m,m=-n,n)/)),:)**2,1)
!!$      enddo
!!$
!!$    end function nz_ToroidalEnergySpectrum_wt

    function zee_PoloidalEnergySpectrum_cee(cee_POLPOT)
      !
      ! ポロイダルポテンシャルから, ポロイダルエネルギーの
      ! 水平 X 波数 l, Y 波数 m の各成分を計算する.
      !
      !  * X 波数 l, Y 波数 m のポロイダルポテンシャルのスペクトル成分
      !    φ(l,m,z)から X 波数 l, Y 波数 m 成分のポロイダルエネルギー
      !    スペクトルは
      !
      !      (1/2)(l**2+m**2){[dφ(n,m,z)/dz]^2 + (l**2+m**2)φ(n,m,z)^2}
      !
      !    と計算される.
      !
      !  * 全てのエネルギースペクトル成分の和をZ積分したもの
      !    が水平平均エネルギーに等しい.
      !
      real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee_POLPOT
      !(in) ポロイダルポテンシャル

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


      real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_Data   ! 作業領域
      real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_DData  ! 作業領域
      integer :: l, m

      zee_Data = zee_cee(cee_POLPOT)
      zee_DData = zee_see(see_DZ_cee(cee_POLPOT))

      do l=-lm,lm
         do m=-mm,mm
            zee_PoloidalEnergySpectrum_cee(:,m,l) =                   &
                 + 0.5* ((2*pi*l/xl)**2+(2*pi*m/yl)**2)               &
                 *(   zee_DData(:,m,l)**2 + zee_DData(:,-m,-l)**2     &
                    + ((2*pi*l/xl)**2+(2*pi*m/yl)**2)                 &
                         *( zee_Data(:,m,l)**2 + zee_Data(:,-m,-l)**2))
         enddo
      enddo

    end function zee_PoloidalEnergySpectrum_cee

    function zee_PoloidalEnergySpectrum_see(see_POLPOT)
      !
      ! ポロイダルポテンシャルから, ポロイダルエネルギーの
      ! 水平 X 波数 l, Y 波数 m の各成分を計算する.
      !
      !  * X 波数 l, Y 波数 m のポロイダルポテンシャルのスペクトル成分
      !    φ(l,m,z)から X 波数 l, Y 波数 m 成分のポロイダルエネルギー
      !    スペクトルは
      !
      !      (1/2)(l**2+m**2){[dφ(n,m,z)/dz]^2 + (l**2+m**2)φ(n,m,z)^2}
      !
      !    と計算される.
      !
      !  * 全てのエネルギースペクトル成分の和をZ積分したもの
      !    が水平平均エネルギーに等しい.
      !
      real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see_POLPOT
      !(in) ポロイダルポテンシャル

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


      real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_Data   ! 作業領域
      real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_DData  ! 作業領域
      integer :: l, m

      zee_Data = zee_see(see_POLPOT)
      zee_DData = zee_cee(cee_DZ_see(see_POLPOT))

      do l=-lm,lm
         do m=-mm,mm
            zee_PoloidalEnergySpectrum_see(:,m,l) =                   &
                 + 0.5* ((2*pi*l/xl)**2+(2*pi*m/yl)**2)               &
                 *(   zee_DData(:,m,l)**2 + zee_DData(:,-m,-l)**2     &
                    + ((2*pi*l/xl)**2+(2*pi*m/yl)**2)                 &
                         *( zee_Data(:,m,l)**2 + zee_Data(:,-m,-l)**2))
         enddo
      enddo

    end function zee_PoloidalEnergySpectrum_see

!!$    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}
!!$      !
!!$      !    と計算される.
!!$      !
!!$      !  * 全ての全波数に対してのエネルギースペクトル成分の和をZ積分したもの
!!$      !    (r^2の重み無し)が球殻内での全エネルギーに等しい.
!!$      !
!!$      real(8), dimension(0:nm,-mm:mm,-lm: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_Z*wz_wt(wt_DZ_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

end module scee_module_fftj
