w_zonal_module_sjpack

    function AvrLonLat_xy(xy_data)
      !
      ! 2 次元緯度経度格子点データの全領域平均(1 層用).
      !
      ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
      ! 総和を計算し, x_X_Weight*y_Y_Weight の総和で割ることで平均している. 
      !
      real(8), intent(in)   :: xy_data(0:im-1,1:jm)
      !(in) 2 次元経度緯度格子点データ

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

      AvrLonLat_xy = AvrLon_x(x_AvrLat_xy(xy_data))

    end function AvrLonLat_xy

    function AvrLon_x(x_data)
      !
      ! 1 次元(X)格子点データの経度(X)方向平均(1 層用).
      !
      ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算し, 
      ! x_X_Weight の総和で割ることで平均している. 
      !
      real(8), intent(in) :: x_data(0:im-1)
      !(in)  1 次元経度(X)格子点データ
      real(8)             :: AvrLon_x       
      !(out) 平均値

      AvrLon_x = IntLon_x(x_data)/sum(x_Lon_weight)

    end function AvrLon_x

    function IntLonLat_xy(xy_data)
      !
      ! 2 次元緯度経度格子点データの全領域積分(1 層用). 
      !
      ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), intent(in)   :: xy_data(0:im-1,1:jm)         
      !(in) 2 次元経度緯度格子点データ(0:im-1,1:jm)

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

      IntLonLat_xy = IntLon_x(x_IntLat_xy(xy_data))

    end function IntLonLat_xy

    function IntLon_x(x_data)
      !
      ! 1 次元経度(X)格子点データの X 方向積分(1 層用).
      !
      ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している.
      !
      real(8), intent(in) :: x_data(0:im-1)   !(in)  1 次元経度(X)格子点データ
      real(8)             :: IntLon_x         !(out) 積分値

      IntLon_x = sum(x_data * x_Lon_weight)

    end function IntLon_x

    function Interpolate_w(w_data,alon,alat)
      real(8), intent(IN) :: w_data(nm+1)   ! スペクトルデータ
      real(8), intent(IN) :: alon           ! 補間する位置(経度)
      real(8), intent(IN) :: alat           ! 補間する位置(緯度)
      real(8)             :: Interpolate_w  ! 補間した値
      
      Interpolate_w = Interpolate_alat_w(w_data,alat)

    end function Interpolate_w

    function l_nm_array00(n,m)
      !
      ! 全波数(n)と東西波数(m,ダミー)からそのスペクトルデータの格納位置を返す.
      ! 
      ! 引数 n,m がともに整数値の場合, 整数値を返す. 
      !
      integer               :: l_nm_array00   
      !(out) スペクトルデータの格納位置 

      integer, intent(in)   :: n     !(in) 全波数
      integer, intent(in)   :: m     !(in) 帯状波数           

      if ( m /= 0 ) then
         call MessageNotify('E','l_nm_array00', 'zonal wavenumber should be zero in w_zonal_module')
      end if

      l_nm_array00 = n+1

    end function l_nm_array00

    function l_nm_array01(n,marray)           ! スペクトルデータの格納位置 
      !
      ! 全波数(n)と東西波数(m, ダミー)からそのスペクトルデータの格納位置を返す.
      ! 
      ! 第 1 引数 n が整数, 第 2 引数 marray が整数 1 次元配列の場合, 
      ! marray と同じ大きさの 1 次元整数配列を返す. 
      !
      integer, intent(in)  :: n               !(in) 全波数
      integer, intent(in)  :: marray(:)       !(in) 帯状波数
      integer              :: l_nm_array01(size(marray))
      !(out) スペクトルデータ位置

      integer              :: i 

      do i=1, size(marray)
         l_nm_array01(i) = l_nm_array00(n,marray(i))
      enddo
    end function l_nm_array01

    function l_nm_array10(narray,m)
      !
      ! 全波数(n)と東西波数(m,ダミー)からそのスペクトルデータの格納位置を返す.
      ! 
      ! 第 1 引数 narray が整数 1 次元配列, 第 2 引数  m が整数の場合, 
      ! narray と同じ大きさの 1 次元整数配列を返す. 
      !
      integer, intent(in)  :: narray(:)           !(in) 全波数  
      integer, intent(in)  :: m                   !(in) 帯状波数
      integer              :: l_nm_array10(size(narray))
      !(out) スペクトルデータ位置

      integer              :: i 

      do i=1, size(narray)
         l_nm_array10(i) = l_nm_array00(narray(i),m)
      enddo
    end function l_nm_array10

    function l_nm_array11(narray,marray)
      !
      ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
      ! 
      ! 第 1,2 引数 narray, marray がともに整数 1 次元配列の場合, 
      ! narray, marray と同じ大きさの 1 次元整数配列を返す. 
      ! narray, marray は同じ大きさでなければならない. 
      !
      integer, intent(in)  :: narray(:)          !(in) 全波数  
      integer, intent(in)  :: marray(:)          !(in) 帯状波数
      integer              :: l_nm_array11(size(narray))
      !(out) スペクトルデータ位置

      integer              :: i 

      if ( size(narray) .ne. size(marray) ) then
         call MessageNotify('E','l_nm_array11', 'dimensions of input arrays  n and m are different.')
      endif

      do i=1, size(narray)
         l_nm_array11(i) = l_nm_array00(narray(i),marray(i))
      enddo
    end function l_nm_array11

    function n_EnergyFromStreamfunc_w(w_Strfunc)
      !
      ! 流線関数のスペクトルデータから各全波数のエネルギー成分(スペクトル)を
      ! 計算する(1 層用).
      !
      !  * 全波数 n の流線関数のスペクトル成分ψ(n,m) から
      !    エネルギースペクトルはΣ[m=-nm]^nm(1/2)n(n+1)ψ(n,m)^2 
      !    と計算される.
      !
      !  * 全てのエネルギースペクトル成分の和に 4πをかけたものが
      !    球面上での全エネルギーに等しい.
      !

      real(8), intent(in)      :: w_Strfunc(:)
      !(in) 流線関数(スペクトルデータ)

      real(8), dimension(0:nm) :: n_EnergyFromStreamfunc_w
      !(out) エネルギースペクトル (水平全波数 n 空間) 

      integer :: n

      do n=0,nm
         n_EnergyFromStreamfunc_w(n) = 0.5 * n*(n+1) * w_StrFunc(l_nm(n,0))**2
      enddo

    end function n_EnergyFromStreamfunc_w

    function n_EnstrophyFromStreamfunc_w(w_Strfunc)
      !
      ! 流線関数のスペクトルデータから各全波数のエネルギー成分(スペクトル)を
      ! 計算する(1 層用)
      !
      ! * 全波数 n の流線関数のスペクトル成分ψ(n,m) からエンストロフィー
      !   スペクトルはΣ[m=-nm]^nm(1/2)n^2(n+1)^2ψ(n,m)^2 と計算される.
      !    
      ! * 全てのエネルギースペクトル成分の和に 4π/R^2 をかけたものが
      !   球面上での全エンストフィーに等しい.
      !
      real(8), intent(in)      :: w_Strfunc(:)
      !(in) 流線関数(スペクトルデータ)

      real(8), dimension(0:nm) :: n_EnstrophyFromStreamfunc_w  
      !(out) エンストロフィースペクトル(水平全波数 n 空間)

      integer :: n

      do n=0,nm
         n_EnstrophyFromStreamfunc_w(n) = 0.5 * n**2 * (n+1)**2 * w_StrFunc(l_nm(n,0))**2
      enddo
    end function n_EnstrophyFromStreamfunc_w

    function nm_EnergyFromStreamfunc_w(w_Strfunc)
      ! 
      ! 流線関数のスペクトルデータからエネルギーの球面調和函数成分
      ! (スペクトル)を計算する(1 層用).
      !
      !  * 全波数 n, 帯状波数 m の流線関数のスペクトル成分ψ(n,m) から
      !    エネルギースペクトルは (1/2)n(n+1)ψ(n,m)^2 と計算される.
      !
      !  * 全てのエネルギースペクトル成分の和に4πをかけたものが球面上での
      !    全エネルギーに等しい.
      !
      !  * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される.
      !    欠損値の値はモジュール変数 w_spectrum_VMiss によって設定できる
      !    (初期値は -999.0)
      !
      real(8), intent(in)   :: w_Strfunc(:)
      !(in) 流線関数(スペクトルデータ)

      real(8), dimension(0:nm,-nm:nm) :: nm_EnergyFromStreamfunc_w
      !(out) エネルギースペクトル(水平全波数 n, 帯状波数 m 空間)

      integer :: n, m

      nm_EnergyFromStreamfunc_w = w_spectrum_VMiss

      do n=0,nm
         do m=1,n
            nm_EnergyFromStreamfunc_w(n,m)  = 0.0D0
            nm_EnergyFromStreamfunc_w(n,-m) = 0.0D0
         enddo
         nm_EnergyFromStreamfunc_w(n,0) = 0.5 * n*(n+1) * w_Strfunc(l_nm(n,0))**2
      enddo
    end function nm_EnergyFromStreamfunc_w

    function nm_EnstrophyFromStreamfunc_w(w_Strfunc)
      !
      ! 流線関数のスペクトルデータからエンストロフィーの球面調和函数成分
      ! (スペクトル)を計算する(1 層用). 
      !
      ! * 全波数 n, 帯状波数 m の流線関数のスペクトル成分ψ(n,m) から
      !    エンストロフィースペクトルは (1/2)n^2(n+1)^2ψ(n,m)^2 と計算される.
      !
      ! * 全てのエンストロフィースペクトル成分の和に4π/R^2をかけたものが
      !   球面上での全エンストロフィーに等しい. ここで R は球面の半径である.
      !
      ! * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される.
      !   欠損値の値はモジュール変数 w_spectrum_VMiss によって設定できる
      !   (初期値は -999.0)
      !
      real(8), intent(in)   :: w_Strfunc(:)
      !(in) 流線関数(スペクトルデータ)

      real(8), dimension(0:nm,-nm:nm) :: nm_EnstrophyFromStreamfunc_w
      ! エンストロフィースペクトル (水平全波数 n, 帯状波数 m 空間)

      integer :: n, m

      nm_EnstrophyFromStreamfunc_w = w_spectrum_VMiss

      do n=0,nm
         do m=1,n
            nm_EnstrophyFromStreamfunc_w(n,m) = 0.0
            nm_EnstrophyFromStreamfunc_w(n,-m) = 0.0
         enddo
         nm_EnstrophyFromStreamfunc_w(n,0) = 0.5 * n**2 * (n+1)**2 * w_Strfunc(l_nm(n,0))**2
      enddo

    end function nm_EnstrophyFromStreamfunc_w

    function nm_l_int(l)
      ! 
      ! スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す.
      !
      ! 引数 l が整数値の場合, 対応する全波数と帯状波数を
      ! 長さ 2 の 1 次元整数値を返す. 
      ! nm_l(1) が全波数, nm_l(2) が帯状波数である. 
      !
      integer               :: nm_l_int(2)  !(out) 全波数, 帯状波数
      integer, intent(in)   :: l            !(in) スペクトルデータの格納位置
      
      nm_l_int(1) = l-1
      nm_l_int(2) = 0
    end function nm_l_int

    function nm_l_array(larray)
      ! 
      ! スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す.
      !
      ! 引数 larray が整数 1 次元配列の場合, 
      ! larray に対応する n, m を格納した 2 次元整数配列を返す. 
      ! nm_l_array(:,1) が全波数, nm_l_array(:,2) が帯状波数である. 
      !
      integer, intent(in)  :: larray(:)
      !(out) 全波数, 帯状波数

      integer              :: nm_l_array(size(larray),2)
      !(in) スペクトルデータの格納位置

      integer              :: i

      do i=1, size(larray)
         nm_l_array(i,:) = nm_l_int(larray(i))
      enddo
    end function nm_l_array

    function w_DLon_w(w_data)
      !
      ! スペクトルデータに経度微分 ∂/∂λ を作用させる(1 層用).
      !
      ! スペクトルデータの経度微分とは, 対応する格子点データに
      ! 経度微分∂/∂λを作用させたデータのスペクトル変換のことである.
      ! 
      real(8)              :: w_DLon_w(nm+1)
      !(out) スペクトルデータの経度微分

      real(8), intent(in)  :: w_data(nm+1)
      !(in) 入力スペクトルデータ

      w_DLon_w = 0.0D0

    end function w_DLon_w

    function w_DivLambda_xy(xy_data)
      !
      ! 格子点データに発散型経度微分 1/(1-μ^2)・∂/∂λ (μ=sinφ) 
      ! を作用させてスペクトルデータに変換して返す(1 層用).
      !
      real(8)              :: w_DivLambda_xy(nm+1)
      !(out) 格子点データを発散型経度微分したスペクトルデータ

      real(8), intent(in)  :: xy_data(0:im-1,1:jm)
      !(in) 入力格子点データ

      w_DivLambda_xy = w_xy(xy_data,ipow=2,iflag=-1)

    end function w_DivLambda_xy

    function w_DivLat_xy(xy_data)
      !
      ! 格子点データに発散型緯度微分 1/cosφ・∂(f cosφ)/∂φ を作用させて
      ! スペクトルデータに変換して返す(1 層用).
      !
      real(8)              :: w_DivLat_xy(nm+1)
      !(out) 格子点データを発散型緯度微分したスペクトルデータ

      real(8), intent(in)  :: xy_data(0:im-1,1:jm)
      !(in) 入力格子点データ

      w_DivLat_xy = w_xy(xy_data,ipow=1,iflag=1)

    end function w_DivLat_xy

    function w_DivLon_xy(xy_data)
      !
      ! 格子点データに発散型経度微分 1/cosφ・∂/∂λ を作用させて
      ! スペクトルデータに変換して返す(1 層用).
      !
      real(8)              :: w_DivLon_xy(nm+1)
      !(out) 格子点データを発散型経度微分したスペクトルデータ
      real(8), intent(in)  :: xy_data(0:im-1,1:jm)
      !(in) 入力格子点データ

      w_DivLon_xy = w_xy(xy_data,ipow=1,iflag=-1)

    end function w_DivLon_xy

    function w_DivMu_xy(xy_data)
      !
      ! 格子点データに発散型緯度微分 ∂/∂μ (μ=sinφ)を作用させて
      ! スペクトルデータに変換して返す(1 層用).
      !
      real(8)              :: w_DivMu_xy(nm+1)
      !(out) 格子点データを発散型緯度微分したスペクトルデータ

      real(8), intent(in)  :: xy_data(0:im-1,1:jm)
      !(in) 入力格子点データ

      w_DivMu_xy = w_xy(xy_data,ipow=2,iflag=1)

    end function w_DivMu_xy

    function w_Div_xy_xy(xy_u,xy_v)
      !
      ! 2 つの入力格子点データをベクトル成分とする発散を計算し, 
      ! スペクトルデータとして返す(1 層用).
      !
      real(8)              :: w_Div_xy_xy(nm+1)
      !(out) 2 つの入力格子点データをベクトル成分とする発散のスペクトルデータ

      real(8), intent(in)  :: xy_u(0:im-1,1:jm)
      !(in) ベクトル経度成分の格子点データ

      real(8), intent(in)  :: xy_v(0:im-1,1:jm)
      !(in) ベクトル緯度成分の格子点データ

      w_Div_xy_xy = w_Divlon_xy(xy_u) + w_Divlat_xy(xy_v)

    end function w_Div_xy_xy

    subroutine w_Finalize
      !
      ! モジュールの終了処理(割り付け配列の解放)をおこなう. 
      !
      ! 解像度を変更する際にはこのサブルーチンを呼んで終了処理を
      ! おこなったのちに再度 w_Initial で初期設定しなければ
      ! ならない. 
      !
      if ( .not. w_initialize ) then
         call MessageNotify('W','w_finalize', 'w_zonal_module_sjpack not initialized yet')
         return
      endif

      deallocate(x_Lon)           ! 格子点座標格納配列(経度)
      deallocate(x_Lon_Weight)
      deallocate(xy_Lon)
      deallocate(xy_Lat)          ! 格子点座標格納配列

      deallocate(rn,irm)

      call l_Finalize

      w_initialize = .false.

      call MessageNotify('M','w_finalize', 'w_zonal_module_sjpack (2013/02/24) is finalized')

    end subroutine w_Finalize

    function w_Jacobian_w_w(w_a,w_b)
      ! 2 つのスペクトルデータにヤコビアン
      !
      !   J(f,g) = ∂f/∂λ・∂g/∂μ - ∂g/∂λ・∂f/∂μ
      !          = ∂f/∂λ・1/cosφ・∂g/∂φ
      !             - ∂g/∂λ・1/cosφ・∂f/∂φ
      !
      ! を作用させる(1 層用).

      real(8)             :: w_Jacobian_w_w(nm+1)
      !(out) 2 つのスペクトルデータのヤコビアン

      real(8), intent(in) :: w_a((nm+1)*(nm+1))
      !(in) 1つ目の入力スペクトルデータ
      
      real(8), intent(in) :: w_b((nm+1)*(nm+1))
      !(in) 2つ目の入力スペクトルデータ

      w_Jacobian_w_w = 0.0D0      

    end function w_Jacobian_w_w

    subroutine w_StreamPotential2Vector(w_Psi, w_Chi, xy_U, xy_V)
      !
      ! 流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に
      ! (逆)変換する(1 層用)
      !
      !   u cosφ =      ∂χ/∂λ - cosφ∂ψ/∂φ, 
      !   v cosφ = cosφ∂χ/∂φ +      ∂ψ/∂λ 
      !
      real(8), intent(in)   :: w_Psi(nm+1)
      !(in) 流線関数
      real(8), intent(in)   :: w_Chi(nm+1)
      !(in) 速度ポテンシャル

      real(8), intent(out)   :: xy_U(0:im-1,1:jm)
      !(out) 速度経度成分
      real(8), intent(out)   :: xy_V(0:im-1,1:jm)
      !(out) 速度緯度成分
      !
      !   u cosφ = - cosφ∂ψ/∂φ の計算
      !
      xy_U = - xy_GradLat_w(w_Psi)
      !
      !   v cosφ = cosφ∂χ/∂φ の計算
      !
      xy_V = xy_GradLat_w(w_Chi)

    end subroutine w_StreamPotential2Vector

    subroutine w_Vector2VorDiv(xy_U, xy_V, w_Vor, w_Div)
      !
      ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に
      ! (正)変換する(1 層用)
      ! 
      !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ 
      !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
      !
      real(8), intent(in)   :: xy_U(0:im-1,1:jm)
      !(in) 速度経度成分
      real(8), intent(in)   :: xy_V(0:im-1,1:jm)
      !(in) 速度緯度成分

      real(8), intent(out)   :: w_Vor(nm+1)
      !(out) 流線関数
      real(8), intent(out)   :: w_Div(nm+1)
      !(out) 速度ポテンシャル
      !
      !   ζ =  - 1/cosφ ∂(u cosφ)/∂φ の計算
      !
      w_Vor = - w_DivLat_xy(xy_U)
      !
      !    D =    1/cosφ ∂(v cosφ)/∂φ の計算
      !
      w_Div = w_DivLat_xy(xy_V)

    end subroutine w_Vector2VorDiv

    subroutine w_VectorCosLat2VorDiv(xy_UCosLat, xy_VCosLat, w_Vor, w_Div)
      !
      ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に
      ! (正)変換する(1 層用)
      ! 
      !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ 
      !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
      !
      real(8), intent(in)   :: xy_UCosLat(0:im-1,1:jm)
      !(in) 速度経度成分 * cos(lat)
      real(8), intent(in)   :: xy_VCosLat(0:im-1,1:jm)
      !(in) 速度緯度成分 * cos(lat)

      real(8), intent(out)   :: w_Vor(nm+1)
      !(out) 流線関数
      real(8), intent(out)   :: w_Div(nm+1)
      !(out) 速度ポテンシャル
      !
      !   ζ =  - 1/cosφ ∂(u cosφ)/∂φ の計算
      !
      w_Vor = - w_xy(xy_UCosLat,ipow=2,iflag=1)
      !
      !    D =    1/cosφ ∂(v cosφ)/∂φ の計算
      !
      w_Div = w_xy(xy_VCosLat,ipow=2,iflag=1)

    end subroutine w_VectorCosLat2VorDiv

    subroutine w_initial(n_in,i_in,j_in,np_in)
      !
      ! スペクトル変換の格子点数, 波数および OPENMP 使用時の
      ! 最大スレッド数を設定する.
      !
      ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を
      ! しなければならない. 
      !
      ! np_in は w_module との互換性のためにおいているだけである. 
      ! OPENMP 計算は実装されていない. 
      !
      integer,intent(in) :: i_in              !(in) 格子点数(東西)
      integer,intent(in) :: j_in              !(in) 格子点数(南北)
      integer,intent(in) :: n_in              !(in) 切断波数の設定
      integer,intent(in), optional :: np_in   !(in) OPENMP での最大スレッド数

      integer :: i, j, n

      if ( present (np_in) )then
         call MessageNotify('W','w_initial','OPENMP calculation not supported')
      endif

      if ( i_in /= 1  )then
         call MessageNotify('W','w_initial','Longitudinal dimension should be 1')
      endif

      im = i_in  
       jm = j_in  
       nm = n_in

      allocate(x_Lon(0:im-1))                ! 格子点座標格納配列(経度)
      allocate(x_Lon_Weight(0:im-1))
      allocate(xy_Lon(0:im-1,1:jm))
      allocate(xy_Lat(0:im-1,1:jm))          ! 格子点座標格納配列

      allocate(rn(0:nm,2),irm(nm+1,2))

      call l_initial(n_in,j_in)

      do i=0,im-1
         x_Lon(i)  = 2*pi/im*i               ! 経度座標
         x_Lon_Weight(i) = 2*pi/im           ! 経度座標重み
      enddo

      do j=1,jm
         xy_Lon(:,j) = x_Lon
      enddo

      do i=0,im-1
         xy_Lat(i,:) = y_Lat
      enddo

      do n=0,nm
         rn(n,1) = -n*(n+1)
      enddo

      rn(0,2) = 1.0D0
      do n=1,nm
         rn(n,2) = -1/(n*(n+1))
      enddo

      do n=1,nm+1
         irm(n,1) = n
         irm(n,2) = 0
      enddo

      w_initialize = .true.

      call MessageNotify('M','w_initial', 'w_zonal_module_sjpack (2013/02/24) is initialized')

    end subroutine w_initial

    function w_xy(xy_data,ipow,iflag)
      !
      ! 格子データからスペクトルデータへ(正)変換する(1 層用).
      !
      real(8)               :: w_xy(nm+1)
      !(out) スペクトルデータ

      real(8), intent(in)   :: xy_data(0:im-1,1:jm)
      !(in) 格子点データ

      integer, intent(in), optional  :: ipow
      !(in) 変換時に同時に作用させる 1/cosφ の次数. 省略時は 0.

      integer, intent(in), optional  :: iflag
      ! 変換の種類
      !    0 : 通常の正変換
      !    1 : 緯度微分 1/cosφ・∂(f cos^2φ)/∂φ を作用させた正変換
      !   -1 : 緯度微分を作用させた正変換
      !    2 : sinφを作用させた正変換
      !  省略時は 0.


      integer, parameter  :: ipow_default  = 0    ! スイッチデフォルト値
      integer, parameter  :: iflag_default = 0    ! スイッチデフォルト値

      integer ipval, ifval

      real(8)             :: xy_work(0:im-1,1:jm) ! 格子点データ作業配列

      if (present(ipow)) then
         ipval = ipow
      else
         ipval = ipow_default
      endif

      if (present(iflag)) then
         ifval = iflag
      else
         ifval = iflag_default
      endif

      xy_work = xy_data/cos(xy_Lat)**ipval

      if ( ifval == -1 ) then                         ! 経度微分は 0 
         w_xy = 0.0D0        
      else if ( ifval == 1 ) then                     ! 緯度微分
         xy_work = xy_work * cos(xy_Lat)
         w_xy = w_DivLat_y(xy_work(0,:))
      else if ( ifval == 2 ) then                     ! sinφ をかけた正変換
         xy_work = xy_work * sin(xy_Lat)
         w_xy = w_y(xy_work(0,:))
      else
         w_xy = w_y(xy_work(0,:))
      endif
      
    end function w_xy

    function x_AvrLat_xy(xy_data)
      !
      ! 2 次元緯度経度格子点データの緯度(Y)方向平均(1 層用).
      !
      ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算し, 
      ! y_Y_Weight の総和で割ることで平均している. 
      !
      real(8), intent(in) :: xy_data(0:im-1,1:jm)
      !(in) 2 次元経度緯度格子点データ(0:im-1,1:jm)
      real(8)             :: x_AvrLat_xy(im)
      !(out) 平均された 1 次元経度(X)格子点データ

      x_AvrLat_xy = x_IntLat_xy(xy_data)/sum(y_Lat_weight)

    end function x_AvrLat_xy

    function x_IntLat_xy(xy_data)
      !
      ! 2 次元緯度経度格子点データの緯度(Y)方向積分(1 層用).
      !
      ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算している. 
      !
      real(8), intent(in) :: xy_data(0:im-1,1:jm)           
      !(in) 2 次元経度緯度格子点データ(0:im-1,1:jm)

      real(8)             :: x_IntLat_xy(0:im-1) 
      !(out) 積分された 1 次元経度(X)格子点データ

      integer :: j

      x_IntLat_xy = 0.0D0
      do j=1,jm
         x_IntLat_xy = x_IntLat_xy + xy_data(:,j) * y_Lat_weight(j)
      enddo

    end function x_IntLat_xy

    function xy_GradLambda_w(w_data)
      !
      ! スペクトルデータに勾配型経度微分 ∂/∂λ を作用する(1 層用).
      !
      real(8)              :: xy_GradLambda_w(0:im-1,1:jm)
      !(out) スペクトルデータを勾配型経度微分した格子点データ

      real(8), intent(in)  :: w_data(nm+1)
      !(in) 入力スペクトルデータ
      
      xy_GradLambda_w = xy_w(w_data,ipow=0,iflag=-1)

    end function xy_GradLambda_w

    function xy_GradLat_w(w_data)
      !
      ! スペクトルデータに勾配型緯度微分 ∂/∂φ を作用させて
      ! 格子点データに変換して返す(1 層用).
      !
      real(8)              :: xy_GradLat_w(0:im-1,1:jm)
      !(out) スペクトルデータを勾配型緯度微分した格子点データ

      real(8), intent(in)  :: w_data(nm+1)
      !(in) 入力スペクトルデータ

      xy_GradLat_w = xy_w(w_data,ipow=1,iflag=1)

    end function xy_GradLat_w

    function xy_GradLon_w(w_data)
      !
      ! スペクトルデータに勾配型経度微分 1/cosφ・∂/∂λ を
      ! 作用させた格子点データを返す(1 層用).
      !
      real(8)              :: xy_GradLon_w(0:im-1,1:jm)
      !(out) スペクトルデータを勾配型経度微分した格子点データ

      real(8), intent(in)  :: w_data(nm+1)
      !(in) 入力スペクトルデータ

      xy_GradLon_w = xy_w(w_data,ipow=1,iflag=-1)

    end function xy_GradLon_w

    function xy_GradMu_w(w_data)
      !
      ! スペクトルデータに勾配型緯度微分 (1-μ^2)∂/∂μ  (μ=sinφ)
      ! を作用させて格子点データに変換して返す(1 層用).
      !
      real(8)              :: xy_GradMu_w(0:im-1,1:jm)
      !(out) スペクトルデータを勾配型緯度微分した格子点データ

      real(8), intent(in)  :: w_data(nm+1)
      !(in) 入力スペクトルデータ

      xy_GradMu_w = xy_w(w_data,ipow=0,iflag=1)

    end function xy_GradMu_w

    function xy_w(w_data,ipow,iflag)
      !
      ! スペクトルデータから格子データへ変換する(1 層用).
      !
      real(8)               :: xy_w(0:im-1,1:jm)
      !(out) 格子点データ

      real(8), intent(in)   :: w_data(nm+1)
      !(in) スペクトルデータ

      integer, intent(in), optional  :: ipow      
      !(in) 作用させる 1/cosφ の次数. 省略時は 0. 

      integer, intent(in), optional  :: iflag
      !(in) 変換の種類
      !    0 : 通常の正変換
      !    1 : 緯度微分 cosφ・∂/∂φ を作用させた逆変換
      !   -1 : 経度微分を作用させた逆変換
      !    2 : sinφを作用させた逆変換(非実装)
      !    省略時は 0.
      !
      integer, parameter  :: ipow_default  = 0
      integer, parameter  :: iflag_default = 0

      integer ipval, ifval, i

      if (present(ipow)) then
         ipval = ipow
      else
         ipval = ipow_default
      endif

      if (present(iflag)) then
         ifval = iflag
      else
         ifval = iflag_default
      endif

      if ( ifval == -1 ) then                          ! 経度微分は 0 
         xy_w = 0.0D0        
      else if ( ifval == 1 ) then                    ! 緯度微分
         do i=0,im-1
            xy_w(i,:) = y_GradLat_w(w_data) * cos(y_Lat)
         enddo
      else if ( ifval == 2 ) then                     ! sinφ をかけた逆変換
         do i=0,im-1
            xy_w(i,:) = y_w(w_data)
         enddo
         xy_w = xy_w * sin(xy_Lat)
      else
         do i=0,im-1
            xy_w(i,:) = y_w(w_data)
         enddo
      endif

      if ( ipval /= 0 ) then
         xy_w = xy_w/cos(xy_Lat)**ipval
      end if

    end function xy_w

    function y_AvrLon_xy(xy_data)
      !
      ! 2 次元緯度経度格子点データの経度(X)方向平均(1 層用).
      !
      ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算し, 
      ! x_X_Weight の総和で割ることで平均している. 
      !
      real(8), intent(in) :: xy_data(0:im-1,1:jm)
      !(in) 2 次元経度緯度格子点データ(0:im-1,1:jm)
      real(8)             :: y_AvrLon_xy(1:jm)
      !(out) 平均された 1 次元緯度(Y)格子点

      y_AvrLon_xy = y_IntLon_xy(xy_data)/sum(x_Lon_weight)

    end function y_AvrLon_xy

    function y_IntLon_xy(xy_data)
      !
      ! 2 次元緯度経度格子点データの経度(X)方向積分(1 層用).
      !
      ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している. 
      !
      real(8), intent(in) :: xy_data(0:im-1,1:jm)
      !(in) 2 次元経度緯度格子点データ(0:im-1,1:jm)

      real(8)             :: y_IntLon_xy(1:jm)
      !(out) 積分された 1 次元緯度(Y)格子点データ

      integer :: i

      y_IntLon_xy = 0.0D0
      do i=0,im-1
         y_IntLon_xy = y_IntLon_xy + xy_data(i,:) * x_Lon_weight(i)
      enddo

    end function y_IntLon_xy