! -*- mode: f90; coding: utf-8 -*-
!-------------------------------------------------------------------------
! Copyright (c) 2002-2019 SPMODEL Development Group. All rights reserved.
!
! 履歴  2019/03/13  竹広真一
!       2019/10/04  佐々木洋平
!-------------------------------------------------------------------------
!> @copyright Copyright (C) SPMODEL Development Group, 2002-2019.
!>            License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @author Shin-ichi Takehiro, Youhei SASAKI
!> @brief
!> @ref wa_module の下部モジュール: 基本変換
!>
!> @warning
!> 本モジュールを直接呼ぶ事は想定されていないことに注意されたい．
!> ユーザは @ref wa_module を呼ぶこと．
!--
!----------------------------------------------------------------------
!     Copyright (c) 2002-2013 SPMDOEL Development Group
!----------------------------------------------------------------------
!表題  wa_base_module_svpack
!
!  spml/wa_base_module_svpack モジュールは球面上での流体運動を
!  球面調和函数を用いたスペクトル法によって数値計算するための
!  モジュール wa_module_svpack の下部モジュールであり, スペクトル計算の
!  基本的な Fortran90 関数を提供する.
!
!  球面上の 1 層モデル用 w_base_module_svpack モジュールを多層モデル用に
!  拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
!  対する変換が行える.
!
!  内部で ISPACK の SVPACK の Fortran77 サブルーチンを呼んでいる.
!  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!  ついては ISPACK/SVPACK のマニュアルを参照されたい.
!
!  このモジュールを使うためには前もって w_base_initial を呼んで
!  切断波数, 格子点数の設定をしておく必要がある.
!
!++
module wa_module_base
  !
  != wa_base_module_svpack
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/wa_base_module_svpack モジュールは球面上での流体運動を
  ! 球面調和函数を用いたスペクトル法によって数値計算するための
  ! モジュール wa_module_svpack の下部モジュールであり, スペクトル計算の
  ! 基本的な Fortran90 関数を提供する.
  !
  ! 球面上の 1 層モデル用 w_base_module_svpack モジュールを多層モデル用に
  ! 拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
  ! 対する変換が行える.
  !
  ! 内部で ISPACK の SVPACK の Fortran77 サブルーチンを呼んでいる.
  ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  ! ついては ISPACK/SVPACK のマニュアルを参照されたい.
  !
  ! このモジュールを使うためには前もって w_base_initial を呼んで
  ! 切断波数, 格子点数の設定をしておく必要がある.
  !
  use dc_types, only: DP
  use dc_message, only: MessageNotify
  use w_module_base, only : &
    & im, jm, nm=>nn, xy_w, w_xy, &
    & w_StreamPotential2Vector,   &
    & w_Vector2VorDiv, w_VectorCosLat2VorDiv
  implicit none
  private
  !> 同時に処理する最大データ数(層の数): SNPACK 用ルーチンとの互換性のため.
  integer(8) :: km=16
  public km                            ! 層数
  public wa_base_Initial               ! 初期化サブルーチン
  public wa_base_Finalize              ! 終了処理
  public xya_wa                        ! 変換関数
  public wa_xya                        ! 変換関数
  public wa_StreamPotential2Vector     ! 流線ポテンシャルから速度場計算
  public wa_Vector2VorDiv              ! 速度場から渦度発散を計算
  public wa_VectorCosLat2VorDiv        ! 速度場から渦度発散を計算
  save km                              ! 最大データ数(層数)を記憶
contains
  !----------------------------------------------------------
  !> SNPACK 用 wa_base_initial の互換のためのダミーサブルーチン
  subroutine wa_base_initial(k_in)
    !> 最大データ数(層数)
    integer,intent(in) :: k_in

    km = k_in
    call MessageNotify('M','wa_base_initial: dummy',&
      'No need to set maximum level number and in wa_base_module')

  end subroutine wa_base_Initial

  !----------------------------------------------------------
  !> 球面調和関数スペクトル -> 格子点:
  !> スペクトルデータから格子点へと逆変換する(多層用)
  !>
  function xya_wa(wa_data,ipow,iflag)
    !> スペクトルデータ((nm+1)*(nm+1),:)
    real(DP), intent(in)   :: wa_data(:,:)
    !> 格子点データ(0:im-1,1:jm,:)
    real(DP)               :: xya_wa(0:im-1,1:jm,size(wa_data,2))
    !> 作用させる 1/cosφ の次数. 省略時は 0.
    integer, intent(in), optional  :: ipow
    !> 変換の種類. 省略時は 0
    !> *  0 : 通常の正変換
    !> * -1 : 経度微分を作用させた逆変換
    !> *  1 : 緯度微分 \f$cosφ・∂/∂φ \f$ を作用させた逆変換
    !> *  2 : \f$sinφ\f$を作用させた逆変換
    integer, intent(in), optional  :: iflag

    integer, parameter  :: ipow_default  = 0
    integer, parameter  :: iflag_default = 0
    integer ipval, ifval
    integer k

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

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

    do k=1,size(wa_data,2)
      xya_wa(:,:,k) = xy_w(wa_data(:,k),iflag=ifval,ipow=ipval)
    enddo

  end function xya_wa

  !----------------------------------------------------------
  !> 格子点 -> 球面調和関数スペクトル:
  !> 格子データからスペクトルデータへ(正)変換する(多層用).
  !>
  function wa_xya(xya_data,ipow,iflag) !
    !> 格子点データ(0:im-1,1:jm,:)
    real(DP), intent(in)   :: xya_data(0:,:,:)
    !> スペクトルデータ((nm+1)*(nm+1),:)
    real(DP)               :: wa_xya((nm+1)*(nm+1),size(xya_data,3))
    !> 変換時に同時に作用させる \f$1/cosφ\f$ の次数. 省略時は 0.
    integer, intent(in), optional  :: ipow
    !> 変換の種類. 省略時は 0
    !> *   0 : 通常の正変換
    !> *  -1 : 経度微分を作用させた正変換
    !> *   1 : 緯度微分 \f$1/cosφ・∂(f \cos^2φ)/∂φ\f$ を作用させた正変換
    !> *   2 : \f$sinφ\f$を作用させた正変換
    integer, intent(in), optional  :: iflag

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

    integer ipval
    integer ifval
    integer k

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

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

    do k=1,size(xya_data,3)
      wa_xya(:,k) = w_xy(xya_data(:,:,k),iflag=ifval,ipow=ipval)
    enddo

  end function wa_xya

  !----------------------------------------------------------
  !> 流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に
  !> (逆)変換する(多層用)
  !>
  !> スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
  !>\f{align*}{
  !>   u \cosφ & =       ∂χ/∂λ - \cosφ∂ψ/∂φ, \\
  !>   v \cosφ   = \cosφ∂χ/∂φ +       ∂ψ/∂λ
  !>\f}
  subroutine wa_StreamPotential2Vector(wa_Psi, wa_Chi, xya_U, xya_V)
    !> 流線関数((nm+1)*(nm+1),:)
    real(DP), intent(in)   :: wa_Psi(:,:)
    !> 速度ポテンシャル((nm+1)*(nm+1),:)
    real(DP), intent(in)   :: wa_Chi(:,:)
    !> 速度経度成分(0:im-1,1:jm,:)
    real(DP), intent(out)   :: xya_U(0:im-1,1:jm,size(wa_Psi,2))
    !> 速度緯度成分(0:im-1,1:jm,:)
    real(DP), intent(out)   :: xya_V(0:im-1,1:jm,size(wa_Psi,2))

    integer :: k

    do k=1,size(wa_Psi,2)
      call w_StreamPotential2Vector( &
        & wa_Psi(:,k), wa_Chi(:,k), xya_U(:,:,k), xya_V(:,:,k) )
    enddo

  end subroutine wa_StreamPotential2Vector

  !----------------------------------------------------------
  !> 速度場(格子データ)から渦度・発散(スペクトルデータ)に
  !> (正)変換する(多層用)
  !>
  !> スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
  !> \f{align*}{
  !>   ζ & = 1/\cosφ ∂v/∂λ - 1/\cosφ ∂(u \cosφ)/∂φ, \\
  !>    D   = 1/\cosφ ∂u/∂λ + 1/\cosφ ∂(v \cosφ)/∂φ
  !>\f}
  subroutine wa_Vector2VorDiv(xya_U, xya_V, wa_Vor, wa_Div)
    !> 速度経度成分(0:im-1,1:jm,:)
    real(DP), intent(in)   :: xya_U(0:,:,:)
    !> 速度緯度成分(0:im-1,1:jm,:)
    real(DP), intent(in)   :: xya_V(0:,:,:)
    !> 流線関数
    real(DP), intent(out)   :: wa_Vor((nm+1)*(nm+1),size(xya_U,3))
    !> 速度ポテンシャル
    real(DP), intent(out)   :: wa_Div((nm+1)*(nm+1),size(xya_U,3))

    integer :: k

    do k=1,size(xya_U,3)
      call w_Vector2VorDiv( &
        & xya_U(:,:,k), xya_V(:,:,k), wa_Vor(:,k), wa_Div(:,k))
    enddo
  end subroutine wa_Vector2VorDiv

  !----------------------------------------------------------
  !> 速度場(格子データ)から渦度・発散(スペクトルデータ)に
  !> (正)変換する(多層用)
  !>
  !> スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
  !> \f{align*}{
  !>   ζ & = 1/\cosφ∂v/∂λ - 1/\cosφ ∂(u \cosφ)/∂φ, \\
  !>    D & = 1/\cosφ∂u/∂λ + 1/\cosφ ∂(v \cosφ)/∂φ
  !> \f}
  subroutine wa_VectorCosLat2VorDiv(xya_UCosLat, xya_VCosLat, wa_Vor, wa_Div)
    !> 速度経度成分 * cos(lat) (0:im-1,1:jm,:)
    real(DP), intent(in)   :: xya_UCosLat(0:,:,:)
    !> 速度緯度成分 * cos(lat) (0:im-1,1:jm,:)
    real(DP), intent(in)   :: xya_VCosLat(0:,:,:)
    !> 流線関数
    real(DP), intent(out)   :: wa_Vor((nm+1)*(nm+1),size(xya_UCosLat,3))
    !> 速度ポテンシャル
    real(DP), intent(out)   :: wa_Div((nm+1)*(nm+1),size(xya_UCosLat,3))

    integer :: k

    do k=1,size(xya_UCosLat,3)
      call w_VectorCosLat2VorDiv &
        (xya_UCosLat(:,:,k), xya_VCosLat(:,:,k), wa_Vor(:,k), wa_Div(:,k))
    enddo
  end subroutine wa_VectorCosLat2VorDiv

  !----------------------------------------------------------
  !> 後方互換性のためのダミールーチン.
  !> 上位サブルーチン wa_Finalize を使用すること.
  !>
  subroutine wa_base_Finalize
    call MessageNotify('M','wa_base_finalize: dummy', &
      & 'wa_base_module is finalized')

  end subroutine wa_base_Finalize

end module wa_module_base
