!--
!----------------------------------------------------------------------
!     Copyright (c) 2020 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!表題  wb_deriv_mpi_module_sypack
!
!  spml/wb_deriv_mpi_module_sypack モジュールは球面上での流体運動を
!  球面調和函数を用いたスペクトル法によって数値計算するための
!  モジュール wb_mpi_module_sypack の下部モジュールであり, スペクトル法の
!  微分計算のための Fortran90 関数を提供する.
!
!  球面上の 1 層モデル用 w_deriv_mpi_module_sypack モジュールを多層モデル用に
!  拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
!  対する変換が行える.
!
!  内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる.
!  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!  ついては ISPACK3/SYPACK のマニュアルを参照されたい.
!
!履歴  2020/08/01  竹広真一
!
!      制限
!         ・変換する格子点データ, スペクトルデータの配列の大きさは決めうち
!
module wb_mpi_module_deriv
  !
  != wb_mpi_module_deriv
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  !  spml/wb_mpi_module_deriv モジュールは球面上での流体運動を
  !  球面調和函数を用いたスペクトル法によって数値計算するための
  !  モジュール wb_mpi_module の下部モジュールであり, スペクトル法の
  !  微分計算のための Fortran90 関数を提供する.
  !
  !  球面上の 1 層モデル用 w_deriv_mpi_module モジュールを多層モデル用に
  !  拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
  !  対する変換が行える.
  !
  !  内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる.
  !  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  !  ついては ISPACK3/SYPACK のマニュアルを参照されたい.
  !
  use dc_types, only: DP
  use dc_message, only : MessageNotify
  !use w_mpi_module_base, only : im, jm, nm=>nn, jc=>jl
  use w_mpi_module_deriv
  use wb_mpi_module_base, only : &
    & pva_wb, wb_pva, pva_xvb, xvb_pva, wb_xvb, xvb_wb, pv_w, w_pv, &
    & ic, is, ie, jc, nc, kc, km

  implicit none

  ! private im, jm, nm, nc, kc                  ! for Intel Fortran
  private

  public wb_deriv_mpi_Initial
  public wb_deriv_mpi_Finalize
  public wb_Lapla_wb, wb_LaplaInv_wb          ! ラプラシアンと逆演算
  public wb_DLon_wb                           ! 経度微分
  public pva_GradLon_wb, pva_GradLat_wb       ! 勾配型微分
  public wb_DivLon_pva, wb_DivLat_pva         ! 発散型微分
  public wb_Div_pva_pva                       ! 発散型微分
  public wb_JacobianMPI_wb_wb                 ! ヤコビアン
  public pva_GradLambda_wb, pva_GradMu_wb     ! 勾配型微分(λ,μ座標)
  public wb_DivLambda_pva, wb_DivMu_pva       ! 発散型微分(λ,μ座標)
  public w_Lapla_w, w_LaplaInv_w              ! ラプラシアンと逆演算
  public w_DLon_w                             ! 経度微分
  public pv_GradLon_w, pv_GradLat_w           ! 勾配型微分
  public w_DivLon_pv, w_DivLat_pv             ! 発散型微分
  public w_Div_pv_pv                          ! 発散型微分
  public w_JacobianMPI_w_w                    ! ヤコビアン
  public pv_GradLambda_w, pv_GradMu_w         ! 勾配型微分(λ,μ座標)
  public w_DivLambda_pv, w_DivMu_pv           ! 発散型微分(λ,μ座標)

  public rn                                   ! ラプラシアン演算用配列

contains

    !--------------- 初期化 -----------------
    subroutine wb_deriv_mpi_initial
      !
      ! スペクトル微分計算に必要となる作業領域を設定する.
      !
      ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで
      ! 初期設定をしなければならない.
      !
      ! 上位サブルーチン wb_mpi_Initial を使用すること.
      !
!!$      wb_mpi_deriv_initialize=.true.

      call w_deriv_mpi_initial

      call MessageNotify('M','wb_deriv_mpi_initial', &
           'wb_mpi_module_deriv (2020/08/01)is initialized')

    end subroutine wb_deriv_mpi_initial

    !--------------- 微分計算(多層) -----------------
    function wb_Lapla_wb(wb_data)
      !
      ! 入力スペクトルデータにラプラシアン
      !
      !    ▽^2 = 1/cos^2φ・∂^2/∂λ^2 + 1/cosφ・∂/∂φ(cosφ∂/∂φ)
      !
      ! を作用する(多層用).
      !
      ! スペクトルデータのラプラシアンとは, 対応する格子点データに
      ! ラプラシアンを作用させたデータのスペクトル変換のことである.
      !
      real(DP), intent(in)  :: wb_data(:,:)
      !(in) 入力スペクトルデータ
      real(DP)              :: wb_lapla_wb(nc,size(wb_data,2))
      !(out) 入力スペクトルデータのラプラシアン
      integer :: k

      do k=1,size(wb_data,2)
         wb_Lapla_wb(:,k) = rn(:,1)*wb_data(:,k)
      enddo

    end function wb_Lapla_wb

    function wb_LaplaInv_wb(wb_data)
      !
      ! 入力スペクトルデータに逆ラプラシアン
      !
      !    ▽^{-2}
      !      =[1/cos^2φ・∂^2/∂λ^2 + 1/cosφ・∂/∂φ(cosφ∂/∂φ)]^{-1}
      !
      ! を作用する(多層用).
      !
      ! スペクトルデータの逆ラプラシアンとは, 対応する格子点データに
      ! 逆ラプラシアンを作用させたデータのスペクトル変換のことである.
      !
      real(DP), intent(in)  :: wb_data(:,:)
      !(in) 入力スペクトルデータ
      real(DP)              :: wb_LaplaInv_wb(nc,size(wb_data,2))
      !(out) スペクトルデータの逆ラプラシアン
      integer :: k

      do k=1,size(wb_data,2)
         wb_LaplaInv_wb(:,k) = rn(:,2)*wb_data(:,k)
      enddo

    end function wb_LaplaInv_wb

    function wb_DLon_wb(wb_data)
      !
      ! スペクトルデータに経度微分 ∂/∂λ を作用させる(多層用).
      !
      ! スペクトルデータの経度微分とは, 対応する格子点データに
      ! 経度微分∂/∂λを作用させたデータのスペクトル変換のことである.
      !
      real(DP), intent(in)  :: wb_data(:,:)
      !(in) 入力スペクトルデータ
      real(DP)              :: wb_DLon_wb(nc,size(wb_data,2))
      !(out) スペクトルデータの経度微分
      integer :: k

      do k=1,kc
         wb_DLon_wb(:,k) = w_DLon_w(wb_data(:,k))
      enddo

    end function wb_DLon_wb

    function pva_GradLon_wb(wb_data)
      !
      ! スペクトルデータに勾配型経度微分 1/cosφ・∂/∂λ を
      ! 作用させた格子点データを返す(多層用).
      !
      real(DP), intent(in)  :: wb_data(nc,kc)
      !(in) 入力スペクトルデータ
      real(DP)              :: pva_GradLon_wb(ic,1:jc,km)
      !(out) スペクトルデータを勾配型経度微分した格子点データ

      pva_GradLon_wb = pva_wb(wb_data,ipow=1,iflag=-1)

    end function pva_GradLon_wb

    function pva_GradLat_wb(wb_data)
      !
      ! スペクトルデータに勾配型緯度微分 ∂/∂φ を作用させて
      ! 格子点データに変換して返す(多層用).
      !
      real(DP), intent(in)  :: wb_data(nc,kc)
      !(in) 入力スペクトルデータ
      real(DP)              :: pva_GradLat_wb(ic,1:jc,km)
      !(out) スペクトルデータを勾配型緯度微分した格子点データ

      pva_GradLat_wb = pva_wb(wb_data,ipow=1,iflag=1)

    end function pva_GradLat_wb

    function wb_DivLon_pva(pva_data)
      !
      ! 格子点データに発散型経度微分 1/cosφ・∂/∂λ を作用させて
      ! スペクトルデータに変換して返す(多層用).
      !
      real(DP), intent(in)  :: pva_data(ic,jc,km)
      !(in) 入力格子点データ
      real(DP)              :: wb_DivLon_pva(nc,kc)
      !(out) 格子点データを発散型経度微分したスペクトルデータ

      wb_DivLon_pva = wb_pva(pva_data,ipow=1,iflag=-1)

    end function wb_DivLon_pva

    function wb_DivLat_pva(pva_data)
      !
      ! 格子点データに発散型緯度微分 1/cosφ・∂(f cosφ)/∂φ を作用させて
      ! スペクトルデータに変換して返す(多層用).
      !
      real(DP), intent(in)  :: pva_data(ic,jc,km)
      !(in) 入力格子点データ
      real(DP)              :: wb_DivLat_pva(nc,kc)
      !(out) 格子点データを発散型緯度微分したスペクトルデータ

      wb_DivLat_pva = wb_pva(pva_data,ipow=1,iflag=1)

    end function wb_DivLat_pva

    function wb_Div_pva_pva(pva_u,pva_v)
      !
      ! 2 つの入力格子点データをベクトル成分とする発散を計算し,
      ! スペクトルデータとして返す(多層用).
      !
      real(DP), intent(in)  :: pva_u(ic,jc,km)
      !(in) ベクトル経度成分の格子点データ
      real(DP), intent(in)  :: pva_v(ic,jc,km)
      !(in) ベクトル緯度成分の格子点データ
      real(DP)              :: wb_Div_pva_pva(nc,kc)
      !(out) 2 つの入力格子点データをベクトル成分とする発散のスペクトルデータ

      wb_Div_pva_pva = wb_DivLon_pva(pva_u) + wb_DivLat_pva(pva_v)

    end function wb_Div_pva_pva

    function wb_JacobianMPI_wb_wb(wb_a,wb_b)
      ! 2 つのスペクトルデータにヤコビアン
      !
      !   J(f,g) = ∂f/∂λ・∂g/∂μ - ∂g/∂λ・∂f/∂μ
      !          = ∂f/∂λ・1/cosφ・∂g/∂φ
      !             - ∂g/∂λ・1/cosφ・∂f/∂φ
      !
      ! を作用させる(多層用).
      !
      real(DP), intent(in) :: wb_a(nc,kc)
      !(in) 1つ目の入力スペクトルデータ
      real(DP), intent(in) :: wb_b(nc,kc)
      !(in) 2つ目の入力スペクトルデータ
      real(DP)             :: wb_JacobianMPI_wb_wb(nc,kc)
      !(out) 2 つのスペクトルデータのヤコビアン
      integer :: k

      do k=1,kc
         wb_JacobianMPI_wb_wb(:,k) = w_JacobianMPI_w_w(wb_a(:,k),wb_b(:,k))
      end do

    end function wb_JacobianMPI_wb_wb

    !--------------- 微分計算(一層) -----------------
!!$    function w_Lapla_w(w_data)
!!$      !
!!$      ! 入力スペクトルデータにラプラシアン
!!$      !
!!$      !    ▽^2 = 1/cos^2φ・∂^2/∂λ^2 + 1/cosφ・∂/∂φ(cosφ∂/∂φ)
!!$      !
!!$      ! を作用する(多層用).
!!$      !
!!$      ! スペクトルデータのラプラシアンとは, 対応する格子点データに
!!$      ! ラプラシアンを作用させたデータのスペクトル変換のことである.
!!$      !
!!$      real(DP), intent(in)  :: w_data(:)
!!$      !(in) 入力スペクトルデータ
!!$      real(DP)              :: w_lapla_w(nc)
!!$      !(out) 入力スペクトルデータのラプラシアン
!!$
!!$      w_Lapla_w = ru(:,1)*w_data
!!$
!!$    end function w_Lapla_w
!!$
!!$    function w_LaplaInv_w(w_data)
!!$      !
!!$      ! 入力スペクトルデータに逆ラプラシアン
!!$      !
!!$      !    ▽^{-2}
!!$      !      =[1/cos^2φ・∂^2/∂λ^2 + 1/cosφ・∂/∂φ(cosφ∂/∂φ)]^{-1}
!!$      !
!!$      ! を作用する(多層用).
!!$      !
!!$      ! スペクトルデータの逆ラプラシアンとは, 対応する格子点データに
!!$      ! 逆ラプラシアンを作用させたデータのスペクトル変換のことである.
!!$      !
!!$      real(DP), intent(in)  :: w_data(:)
!!$      !(in) 入力スペクトルデータ
!!$      real(DP)              :: w_LaplaInv_w(nc)
!!$      !(out) スペクトルデータの逆ラプラシアン
!!$
!!$      w_LaplaInv_w = ru(:,2)*w_data
!!$
!!$    end function w_LaplaInv_w
!!$
!!$    function w_DLon_w(w_data)
!!$      !
!!$      ! スペクトルデータに経度微分 ∂/∂λ を作用させる(多層用).
!!$      !
!!$      ! スペクトルデータの経度微分とは, 対応する格子点データに
!!$      ! 経度微分∂/∂λを作用させたデータのスペクトル変換のことである.
!!$      !
!!$      real(DP), intent(in)  :: w_data(:)
!!$      !(in) 入力スペクトルデータ
!!$      real(DP)              :: w_DLon_w(nc)
!!$      !(out) スペクトルデータの経度微分
!!$
!!$      real(DP) :: wb_data(nc,1)
!!$
!!$      wb_data(:,1) = w_data
!!$      wb_data = wb_DLon_wb(wb_data)
!!$      w_DLon_w = wb_data(:,1)
!!$
!!$    end function w_DLon_w

    function pv_GradLon_w(w_data)
      !
      ! スペクトルデータに勾配型経度微分 1/cosφ・∂/∂λ を
      ! 作用させた格子点データを返す(多層用).
      !
      real(DP), intent(in)  :: w_data(nc)
      !(in) 入力スペクトルデータ
      real(DP)              :: pv_GradLon_w(ic,1:jc)
      !(out) スペクトルデータを勾配型経度微分した格子点データ

      pv_GradLon_w = pv_w(w_data,ipow=1,iflag=-1)

    end function pv_GradLon_w

    function pv_GradLat_w(w_data)
      !
      ! スペクトルデータに勾配型緯度微分 ∂/∂φ を作用させて
      ! 格子点データに変換して返す(多層用).
      !
      real(DP), intent(in)  :: w_data(nc)
      !(in) 入力スペクトルデータ
      real(DP)              :: pv_GradLat_w(ic,1:jc)
      !(out) スペクトルデータを勾配型緯度微分した格子点データ

      pv_GradLat_w = pv_w(w_data,ipow=1,iflag=1)

    end function pv_GradLat_w

    function w_DivLon_pv(pv_data)
      !
      ! 格子点データに発散型経度微分 1/cosφ・∂/∂λ を作用させて
      ! スペクトルデータに変換して返す(多層用).
      !
      real(DP), intent(in)  :: pv_data(ic,jc)
      !(in) 入力格子点データ
      real(DP)              :: w_DivLon_pv(nc)
      !(out) 格子点データを発散型経度微分したスペクトルデータ

      w_DivLon_pv = w_pv(pv_data,ipow=1,iflag=-1)

    end function w_DivLon_pv

    function w_DivLat_pv(pv_data)
      !
      ! 格子点データに発散型緯度微分 1/cosφ・∂(f cosφ)/∂φ を作用させて
      ! スペクトルデータに変換して返す(多層用).
      !
      real(DP), intent(in)  :: pv_data(ic,jc)
      !(in) 入力格子点データ
      real(DP)              :: w_DivLat_pv(nc)
      !(out) 格子点データを発散型緯度微分したスペクトルデータ

      w_DivLat_pv = w_pv(pv_data,ipow=1,iflag=1)

    end function w_DivLat_pv

    function w_Div_pv_pv(pv_u,pv_v)
      !
      ! 2 つの入力格子点データをベクトル成分とする発散を計算し,
      ! スペクトルデータとして返す(多層用).
      !
      real(DP), intent(in)  :: pv_u(ic,jc)
      !(in) ベクトル経度成分の格子点データ
      real(DP), intent(in)  :: pv_v(ic,jc)
      !(in) ベクトル緯度成分の格子点データ
      real(DP)              :: w_Div_pv_pv(nc)
      !(out) 2 つの入力格子点データをベクトル成分とする発散のスペクトルデータ

      w_Div_pv_pv = w_DivLon_pv(pv_u) + w_DivLat_pv(pv_v)

    end function w_Div_pv_pv

    !--------------- 微分計算 (λ,μ座標系用, 多層) -----------------
    function pva_GradLambda_wb(wb_data)
      !
      ! スペクトルデータに勾配型経度微分 ∂/∂λ を作用する(多層用).
      !
      real(DP), intent(in)  :: wb_data(nc,kc)
      !(in) 入力スペクトルデータ
      real(DP)              :: pva_GradLambda_wb(ic,jc,km)
      !(out) スペクトルデータを勾配型経度微分した格子点データ

      pva_GradLambda_wb = pva_wb(wb_data,ipow=0,iflag=-1)

    end function pva_GradLambda_wb

    function pva_GradMu_wb(wb_data)
      !
      ! スペクトルデータに勾配型緯度微分 (1-μ^2)∂/∂μ  (μ=sinφ)
      ! を作用させて格子点データに変換して返す(多層用).
      !
      real(DP), intent(in)  :: wb_data(nc,kc)
      !(in) 入力スペクトルデータ
      real(DP)              :: pva_GradMu_wb(ic,1:jc,km)
      !(out) スペクトルデータを勾配型緯度微分した格子点データ

      pva_GradMu_wb = pva_wb(wb_data,ipow=0,iflag=1)

    end function pva_GradMu_wb

    function wb_DivLambda_pva(pva_data)
      !
      ! 格子点データに発散型経度微分 1/(1-μ^2)・∂/∂λ (μ=sinφ)
      ! を作用させてスペクトルデータに変換して返す(多層用).
      !
      real(DP), intent(in)  :: pva_data(ic,1:jc,km)
      !(in) 入力格子点データ
      real(DP)              :: wb_DivLambda_pva(nc,kc)
      !(out) 格子点データを発散型経度微分したスペクトルデータ

      wb_DivLambda_pva = wb_pva(pva_data,ipow=2,iflag=-1)

    end function wb_DivLambda_pva

    function wb_DivMu_pva(pva_data)
      !
      ! 格子点データに発散型緯度微分 ∂/∂μ (μ=sinφ)を作用させて
      ! スペクトルデータに変換して返す(多層用).
      !
      real(DP), intent(in)  :: pva_data(ic,1:jc,km)
      !(in) 入力格子点データ
      real(DP)              :: wb_DivMu_pva(nc,kc)
      !(out) 格子点データを発散型緯度微分したスペクトルデータ

      wb_DivMu_pva = wb_pva(pva_data,ipow=2,iflag=1)

    end function wb_DivMu_pva

    !--------------- 微分計算 (λ,μ座標系用, 一層) -----------------
    function pv_GradLambda_w(w_data)
      !
      ! スペクトルデータに勾配型経度微分 ∂/∂λ を作用する(多層用).
      !
      real(DP), intent(in)  :: w_data(nc)
      !(in) 入力スペクトルデータ
      real(DP)              :: pv_GradLambda_w(ic,1:jc)
      !(out) スペクトルデータを勾配型経度微分した格子点データ

      pv_GradLambda_w = pv_w(w_data,ipow=0,iflag=-1)

    end function pv_GradLambda_w

    function pv_GradMu_w(w_data)
      !
      ! スペクトルデータに勾配型緯度微分 (1-μ^2)∂/∂μ  (μ=sinφ)
      ! を作用させて格子点データに変換して返す(多層用).
      !
      real(DP), intent(in)  :: w_data(nc)
      !(in) 入力スペクトルデータ
      real(DP)              :: pv_GradMu_w(ic,1:jc)
      !(out) スペクトルデータを勾配型緯度微分した格子点データ

      pv_GradMu_w = pv_w(w_data,ipow=0,iflag=1)

    end function pv_GradMu_w

    function w_DivLambda_pv(pv_data)
      !
      ! 格子点データに発散型経度微分 1/(1-μ^2)・∂/∂λ (μ=sinφ)
      ! を作用させてスペクトルデータに変換して返す(多層用).
      !
      real(DP), intent(in)  :: pv_data(ic,1:jc)
      !(in) 入力格子点データ
      real(DP)              :: w_DivLambda_pv(nc)
      !(out) 格子点データを発散型経度微分したスペクトルデータ

      w_DivLambda_pv = w_pv(pv_data,ipow=2,iflag=-1)

    end function w_DivLambda_pv

    function w_DivMu_pv(pv_data)
      !
      ! 格子点データに発散型緯度微分 ∂/∂μ (μ=sinφ)を作用させて
      ! スペクトルデータに変換して返す(多層用).
      !
      real(DP), intent(in)  :: pv_data(ic,1:jc)
      !(in) 入力格子点データ
      real(DP)              :: w_DivMu_pv(nc)
      !(out) 格子点データを発散型緯度微分したスペクトルデータ

      w_DivMu_pv = w_pv(pv_data,ipow=2,iflag=1)

    end function w_DivMu_pv

  !--------------- 終了処理 -----------------
    subroutine wb_deriv_mpi_Finalize
      !
      ! モジュールの終了処理(割り付け配列の解放)をおこなう.
      !
      call w_deriv_mpi_Finalize

      call MessageNotify('M','wb_deriv_mpi_finalize',&
           'wb_mpi_module_deriv is finalized (2020/01/25) ')

    end subroutine wb_deriv_mpi_Finalize

end module wb_mpi_module_deriv
