!--
!----------------------------------------------------------------------
!     Copyright (c) 2016 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!表題  ua_deriv_mpi_module_svpack
!
!  spml/ua_deriv_mpi_module_svpack モジュールは球面上での流体運動を
!  球面調和函数を用いたスペクトル法によって数値計算するための 
!  モジュール wa_module_svpack の下部モジュールであり, スペクトル法の
!  微分計算のための Fortran90 関数を提供する. 
!
!  球面上の 1 層モデル用 w_deriv_module_svpack モジュールを多層モデル用に
!  拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
!  対する変換が行える.
!
!  内部で ISPACK2 の SVPACK の Fortran77 サブルーチンを呼んでいる. 
!  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!  ついては ISPACK2/SVPACK のマニュアルを参照されたい.
!
!履歴  2016/07/30  竹広真一 wa_deriv_module_sjpack より多層 MPI 化
!
!      制限
!         ・変換する格子点データ, スペクトルデータの配列の大きさは決めうち
!
module ua_deriv_mpi_module_svpack
  !
  != ua_deriv_mpi_module_svpack
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  !  spml/ua_deriv_mpi_module_svpack モジュールは球面上での流体運動を
  !  球面調和函数を用いたスペクトル法によって数値計算するための 
  !  モジュール wa_module_svpack の下部モジュールであり, スペクトル法の
  !  微分計算のための Fortran90 関数を提供する. 
  !
  !  球面上の 1 層モデル用 w_deriv_module_svpack モジュールを多層モデル用に
  !  拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
  !  対する変換が行える.
  !
  !  内部で ISPACK2 の SVPACK の Fortran77 サブルーチンを呼んでいる. 
  !  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  !  ついては ISPACK2/SVPACK のマニュアルを参照されたい.
  !
  use dc_message, only : MessageNotify
  use w_base_module_svpack, only : im, jm, nm=>nn
  use ua_base_mpi_module_svpack, only : xyb_ua, ua_xyb, nc, kc, km, ua_wb, wb_ua
  use w_module_svpack, only : w_Jacobian_w_w
!!$  use w_module_svpack, only : w_xy, xy_w, w_Lapla_w, w_LaplaInv_w, &
!!$                              w_DLon_w, w_Jacobian_w_w

  implicit none
 
  real(8), allocatable  :: D(:)
  ! ラプラシアン演算用配列
  !
  ! スペクトルデータのラプラシアンを計算するための係数
  ! 配列のサイズは(nc*2)
  !

  real(8), allocatable  :: rn(:,:)
  ! ラプラシアン演算用配列(w_module と互換性を保つため)
  !
  ! スペクトルデータのラプラシアンを計算するための係数
  ! 配列のサイズは(nc, 2)
  !
  ! r(L,1) には L 番目の格納位置のスペクトルに対するラプラシアン計算の
  ! 係数 -n(n+1) の値が格納されている.
  !

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

  public ua_deriv_mpi_Initial
  public ua_deriv_mpi_Finalize
  public ua_Lapla_ua, ua_LaplaInv_ua          ! ラプラシアンと逆演算
  public ua_DLon_ua                           ! 経度微分
  public xyb_GradLon_ua, xyb_GradLat_ua       ! 勾配型微分
  public ua_DivLon_xyb, ua_DivLat_xyb         ! 発散型微分
  public ua_Div_xyb_xyb                       ! 発散型微分
  public ua_JacobianMPI_ua_ua                 ! ヤコビアン
  public xyb_GradLambda_ua, xyb_GradMu_ua     ! 勾配型微分(λ,μ座標)
  public ua_DivLambda_xyb, ua_DivMu_xyb       ! 発散型微分(λ,μ座標)

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

  save D, rn

contains

    !--------------- 初期化 -----------------
    subroutine ua_deriv_mpi_initial
      !
      ! スペクトル微分計算に必要となる作業領域を設定する.
      !
      ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで
      ! 初期設定をしなければならない.
      !
      ! 上位サブルーチン ua_mpi_Initial を使用すること.
      !
      allocate(D(nc*2))                      ! ラプラシアン演算用配列
      allocate(rn(nc,2))                     ! ラプラシアン演算用配列

      call suinid(nm,D)

      rn = reshape(D,(/nc,2/))

!!$      ua_mpi_deriv_initialize=.true.

      call MessageNotify('M','ua_deriv_mpi_initial', &
           'ua_deriv_mpi_module_svpack (2016/07/30)is initialized')

    end subroutine ua_deriv_mpi_initial

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

      do k=1,size(ua_data,2)
         call suclap(nm,ua_data(1,k),ua_Lapla_ua(1,k),D,1)
      enddo
      
    end function ua_Lapla_ua

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

      do k=1,size(ua_data,2)
         call suclap(nm,ua_data(1,k),ua_LaplaInv_ua(1,k),D,2)
      enddo

    end function ua_LaplaInv_ua

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

      do k=1,size(ua_data,2)
         call sucs2x(nm,ua_data(1,k),ua_DLon_ua(1,k))
      enddo

    end function ua_DLon_ua

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

      xyb_GradLon_ua = xyb_ua(ua_data,ipow=1,iflag=-1)

    end function xyb_GradLon_ua

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

      xyb_GradLat_ua = xyb_ua(ua_data,ipow=1,iflag=1)

    end function xyb_GradLat_ua

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

      ua_DivLon_xyb = ua_xyb(xyb_data,ipow=1,iflag=-1)

    end function ua_DivLon_xyb

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

      ua_DivLat_xyb = ua_xyb(xyb_data,ipow=1,iflag=1)

    end function ua_DivLat_xyb

    function ua_Div_xyb_xyb(xyb_u,xyb_v)
      !
      ! 2 つの入力格子点データをベクトル成分とする発散を計算し, 
      ! スペクトルデータとして返す(多層用).
      !
      real(8), intent(in)  :: xyb_u(0:im-1,jm,kc)
      !(in) ベクトル経度成分の格子点データ
      real(8), intent(in)  :: xyb_v(0:im-1,jm,kc)
      !(in) ベクトル緯度成分の格子点データ
      real(8)              :: ua_Div_xyb_xyb(nc,km)
      !(out) 2 つの入力格子点データをベクトル成分とする発散のスペクトルデータ

      ua_Div_xyb_xyb = ua_DivLon_xyb(xyb_u) + ua_DivLat_xyb(xyb_v)

    end function ua_Div_xyb_xyb

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

      real(8) :: wb_a((nm+1)**2,kc)
      real(8) :: wb_b((nm+1)**2,kc)
      real(8) :: wb_Jacobian((nm+1)**2,kc)

      wb_a = wb_ua(ua_a)
      wb_b = wb_ua(ua_b)      

      do k=1,kc
         wb_Jacobian(:,k) = w_Jacobian_w_w(wb_a(:,k),wb_b(:,k))
      end do

      ua_JacobianMPI_ua_ua = ua_wb(wb_Jacobian)
      
    end function ua_JacobianMPI_ua_ua

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

      xyb_GradLambda_ua = xyb_ua(ua_data,ipow=0,iflag=-1)

    end function xyb_GradLambda_ua

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

      xyb_GradMu_ua = xyb_ua(ua_data,ipow=0,iflag=1)

    end function xyb_GradMu_ua

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

      ua_DivLambda_xyb = ua_xyb(xyb_data,ipow=2,iflag=-1)

    end function ua_DivLambda_xyb

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

      ua_DivMu_xyb = ua_xyb(xyb_data,ipow=2,iflag=1)

    end function ua_DivMu_xyb

  !--------------- 終了処理 -----------------
    subroutine ua_deriv_mpi_Finalize
      ! 
      !
      ! このサブルーチンを単独で用いるのでなく, 
      ! 上位サブルーチン wa_Finalize を使用すること.
      !
      deallocate(D)
      deallocate(rn)
      
      call MessageNotify('M','ua_deriv_mpi_finalize',&
           'ua_deriv_mpi_module_svpack is finalized (2016/08/15) ')

    end subroutine ua_deriv_mpi_Finalize
    
end module ua_deriv_mpi_module_svpack

