! -*- mode: f90; coding: utf-8 -*-
!-------------------------------------------------------------------------
! Copyright (c) 2019 SPMODEL Development Group. All rights reserved.
!
! 履歴  2019/03/16  竹広真一
!       2019/10/04  佐々木洋平
!       2019/10/09  佐々木洋平
!       2020/11/20  竹広真一
!-------------------------------------------------------------------------
!> @author Shin-ichi Takehiro, Youhei SASAKI
!> @copyright Copyright (C) SPMODEL Development Group, 2019.
!>            License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @brief @ref w_module の下部モジュール: 微分
!> @warning
!> 本モジュールを直接呼ぶ事は想定されていないことに注意されたい．
!> ユーザは @ref w_module を呼ぶこと．
!>
module w_module_deriv_mint
  use dc_message, only : MessageNotify
  use dc_types, only: DP
  use w_module_base_mint, only :   &
    & im, jm, mm, nn, mi,          &
    & w_base_Initial, xy_w, w_xy
  implicit none
  private
  public w_deriv_Initial
  public w_deriv_Finalize
  public w_Lapla_w
  public w_LaplaInv_w
  public w_DLon_w
  public xy_GradLon_w
  public xy_GradLat_w
  public w_DivLon_xy
  public w_DivLat_xy
  public w_Div_xy_xy
  public w_Jacobian_w_w
  public xy_GradLambda_w
  public xy_GradMu_w
  public w_DivLambda_xy
  public w_DivMu_xy
  public rn

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

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

  !> 初期化フラグ
  logical               :: w_deriv_initialize=.false.

  save D, rn

contains

  !----------------------------------------------------------
  !> スペクトル微分計算に必要となる作業領域を設定する.
  !>
  !> 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで
  !> 初期設定をしなければならない.
  !>
  !> @warning
  !> このサブルーチンを単独で用いるのでなく,
  !> 上位サブルーチン
  !> [w_Initial](@ref w_module#w_initial) を使用すること.
  !>
  subroutine w_deriv_initial
    allocate(D((MM/MI+1)*(2*(NN+1)-MM/MI*MI)*2_8))
    allocate(rn((MM/MI+1)*(2*(NN+1)-MM/MI*MI),2_8))

    call sxinmd(mm,nn,D,mi)

    rn = reshape(D,(/(MM/MI+1)*(2*(NN+1)-MM/MI*MI),2_8/))

    w_deriv_initialize=.true.

    call MessageNotify('M','w_deriv_initial',&
      'w_module_deriv_mint (2020/11/18) is initialized')

  end subroutine w_deriv_initial

  !----------------------------------------------------------
  !> 入力スペクトルデータにラプラシアン
  !> \f{align*}{
  !>    ▽^2 = 1/\cos^2φ・∂^2/∂λ^2 + 1/\cosφ・∂/∂φ(\cosφ∂/∂φ)
  !> \f}
  !> を作用する(1 層用).
  !>
  function w_Lapla_w(w_data)
    !> 入力スペクトルデータ
    real(DP), intent(in)  :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> 入力スペクトルデータのラプラシアン
    real(DP)              :: w_Lapla_w((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

    call sxclam(mm,nn,w_data,w_Lapla_w,D,1_8,mi)

  end function w_Lapla_w

  !----------------------------------------------------------
  !> 入力スペクトルデータに逆ラプラシアン
  !> \f{align*}{
  !>    ▽^{-2}
  !>      =[1/\cos^2φ・∂^2/∂λ^2 + 1/\cosφ・∂/∂φ(\cosφ∂/∂φ)]^{-1}
  !> \f}
  !> を作用する(1 層用).
  !>
  function w_LaplaInv_w(w_data)
    !>  入力スペクトルデータ
    real(DP), intent(in)  :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> スペクトルデータの逆ラプラシアン
    real(DP)              :: w_LaplaInv_w((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

    call sxclam(mm,nn,w_data,w_LaplaInv_w,D,2_8,mi)

  end function w_LaplaInv_w

  !----------------------------------------------------------
  !> スペクトルデータに経度微分 \f$∂/∂λ\f$ を作用させる(1 層用).
  !>
  function w_DLon_w(w_data)
    !> 入力スペクトルデータ
    real(DP), intent(in)  :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> スペクトルデータの経度微分
    real(DP)              :: w_DLon_w((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

    call sxcsmx(mm,nn,w_data,w_DLon_w,mi)

  end function w_DLon_w

  !-------------------------------------------------------------
  !> スペクトルデータに勾配型経度微分 \f$1/\cosφ・∂/∂λ\f$ を
  !> 作用させた格子点データを返す(1 層用).
  !>
  function xy_GradLon_w(w_data)
    !> 入力スペクトルデータ
    real(DP), intent(in)  :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> スペクトルデータを勾配型経度微分した格子点データ
    real(DP)              :: xy_GradLon_w(0:im-1,1:jm)

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

  end function xy_GradLon_w

  !---------------------------------------------------------------
  !> スペクトルデータに勾配型緯度微分 \f$∂/∂φ\f$ を作用させて
  !> 格子点データに変換して返す(1 層用).
  !>
  function xy_GradLat_w(w_data)
    !> 入力スペクトルデータ
    real(DP), intent(in)  :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> スペクトルデータを勾配型緯度微分した格子点データ
    real(DP)              :: xy_GradLat_w(0:im-1,1:jm)

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

  end function xy_GradLat_w

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

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

  end function w_DivLon_xy

  !-------------------------------------------------------------------
  !> 格子点データに発散型緯度微分 \f$1/\cosφ・∂(f \cosφ)/∂φ\f$
  !> を作用させてスペクトルデータに変換して返す(1 層用).
  !>
  function w_DivLat_xy(xy_data)
    !> 入力格子点データ
    real(DP), intent(in)  :: xy_data(0:im-1,1:jm)
    !> 格子点データを発散型緯度微分したスペクトルデータ
    real(DP)              :: w_DivLat_xy((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

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

  end function w_DivLat_xy

  !----------------------------------------------------------
  !> 2 つの入力格子点データをベクトル成分とする発散を計算し,
  !> スペクトルデータとして返す(1 層用).
  !>
  function w_Div_xy_xy(xy_u,xy_v)
    !> ベクトル経度成分の格子点データ
    real(DP), intent(in)  :: xy_u(0:im-1,1:jm)
    !> ベクトル緯度成分の格子点データ
    real(DP), intent(in)  :: xy_v(0:im-1,1:jm)
    !> 2 つの入力格子点データをベクトル成分とする発散のスペクトルデータ
    real(DP)              :: w_Div_xy_xy((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

    w_Div_xy_xy = w_Divlon_xy(xy_u) + w_Divlat_xy(xy_v)

  end function w_Div_xy_xy

  !-------------------------------------------------------------------------
  !> 2 つのスペクトルデータにヤコビアン
  !> \f{align*}{
  !>   J(f,g) & = ∂f/∂λ・∂g/∂μ - ∂g/∂λ・∂f/∂μ \\
  !>          & = ∂f/∂λ・1/\cos φ・∂g/∂φ
  !>               - ∂g/∂λ・1/\cos φ・∂f/∂φ
  !> \f}
  !> を作用させる(1 層用).
  !>
  function w_Jacobian_w_w(w_a,w_b)
    !>  1つ目の入力スペクトルデータ
    real(DP), intent(in) :: w_a((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !>  2つ目の入力スペクトルデータ
    real(DP), intent(in) :: w_b((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> 2 つのスペクトルデータのヤコビアン
    real(DP)             :: w_Jacobian_w_w((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

    w_Jacobian_w_w =                                          &
      & w_xy(  xy_w(w_DLon_w(w_a))*xy_w(w_b,ipow=2,iflag=1)   &
      &      - xy_w(w_DLon_w(w_b))*xy_w(w_a,ipow=2,iflag=1) )

  end function w_Jacobian_w_w

  !----------------------------------------------------------
  !> スペクトルデータに勾配型経度微分 ∂/∂λ を作用する(1 層用).
  !>
  function xy_GradLambda_w(w_data)
    !>  入力スペクトルデータ
    real(DP), intent(in)  :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> スペクトルデータを勾配型経度微分した格子点データ
    real(DP)              :: xy_GradLambda_w(0:im-1,1:jm)

    xy_GradLambda_w = xy_w(w_data,ipow=0,iflag=-1)

  end function xy_GradLambda_w

  !----------------------------------------------------------
  !> スペクトルデータに勾配型緯度微分
  !> \f$(1-μ^2)∂/∂μ ,\, (μ=\sinφ)\f$
  !> を作用させて格子点データに変換して返す(1 層用).
  !>
  function xy_GradMu_w(w_data)
    !>  入力スペクトルデータ
    real(DP), intent(in)  :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> スペクトルデータを勾配型緯度微分した格子点データ
    real(DP)              :: xy_GradMu_w(0:im-1,1:jm)

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

  end function xy_GradMu_w

  !----------------------------------------------------------
  !> 格子点データに発散型経度微分
  !> \f$1/(1-μ^2)・∂/∂λ, \, (μ=\sinφ)\f$
  !> を作用させてスペクトルデータに変換して返す(1 層用).
  !>
  function w_DivLambda_xy(xy_data)
    !>  入力格子点データ
    real(DP), intent(in)  :: xy_data(0:im-1,1:jm)
    !> 格子点データを発散型経度微分したスペクトルデータ
    real(DP)              :: w_DivLambda_xy((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

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

  end function w_DivLambda_xy

  !----------------------------------------------------------
  !> 格子点データに発散型緯度微分
  !> \f$∂/∂μ, \, (μ=\sinφ)\f$
  !> を作用させてスペクトルデータに変換して返す(1 層用).
  !>
  function w_DivMu_xy(xy_data)
    !>  入力格子点データ
    real(DP), intent(in)  :: xy_data(0:im-1,1:jm)
    !> 格子点データを発散型緯度微分したスペクトルデータ
    real(DP)              :: w_DivMu_xy((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

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

  end function w_DivMu_xy

  !----------------------------------------------------------
  !> モジュールの終了処理(割り付け配列の解放)をおこなう.
  !>
  !> @warning
  !> このサブルーチンを単独で用いるのでなく,
  !> 上位サブルーチン
  !> [w_Finalize](@ref w_module#w_finalize) を使用すること.
  !>
  subroutine w_deriv_finalize

    if ( .not. w_deriv_initialize ) then
      call MessageNotify('W','w_deriv_Finalize',&
        'w_module_deriv_mint not initialized yet')
      return
    endif

    deallocate(D)           ! ラプラシアン演算用配列
    deallocate(rn)          ! ラプラシアン演算用配列

    w_deriv_initialize = .false.

    call MessageNotify('M','w_deriv_finalize',&
      'w_module_deriv_mint (2020/11/18) is finalized')

  end subroutine w_deriv_finalize

end module w_module_deriv_mint
