!--
!----------------------------------------------------------------------
! Copyright (c) 2015 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!表題  w_plane_module
!
!   spml/w_module モジュールは 2 次元無限平面上での 2 次元流体運動を
!   球面調和函数を用いたスペクトル法によって数値計算するための
!   Fortran90 関数を提供する. 
!
!   内部で w_module を呼んでいる. スペクトルデータおよび格子点データの
!   格納方法や変換の詳しい計算法については w_module の下部ルーチンである
!   ISPACK/SNPACK,SPPACK のマニュアルを参照されたい.
!
!履歴  2015/12/26  竹広真一  新規作成
!
!++
module w_plane_module2
  !
  != w_plane_module
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id:$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/w_plane_module モジュールは 2 次元無限平面上での 2 次元流体運動を
  ! 球面調和函数を用いたスペクトル法によって数値計算するための Fortran90 
  ! 関数を提供する. 
  !
  ! 内部で w_module を呼んでいる. スペクトルデータおよび格子点データの
  ! 格納方法や変換の詳しい計算法については w_module の下部ルーチンである
  ! ISPACK/SNPACK,SPPACK のマニュアルを参照されたい.
  !
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !
  ! * 関数名の先頭 (w_, nm_, n_, xy_, x_, y_) は, 返す値の形を示している.
  !   w_  :: スペクトルデータ
  !   xy_ :: 2 次元格子点データ
  !
  ! * 関数名の間の文字列(Dx, Dy, Lapla, LaplaInv)は, その関数の作用を表している.
  !
  ! * 関数名の最後 (_w_w, _w, _xy) は, 入力変数の形スペクトルデータ
  !   および格子点データであることを示している.
  !   _w   :: スペクトルデータ
  !   _xy  :: 2 次元格子点データ
  !
  !=== 各データの種類の説明
  !
  ! * xy : 2 次元格子点データ.
  !   * 変数の種類と次元は real(8), dimension(0:im-1,1:jm). 
  !   * im, jm はそれぞれ経度, 緯度座標の格子点数であり, サブルーチン
  !     w_plane_Initial にてあらかじめ設定しておく.
  !
  ! * w : スペクトルデータ.
  !   * 変数の種類と次元は real(8), dimension((nm+1)*(nm+1)). 
  !   * nm は球面調和函数の最大全波数であり, サブルーチン w_Initial にて
  !     あらかじめ設定しておく. 
  !   * スペクトルデータの格納のされ方は関数 l_nm, nm_l によって
  !     調べることができる.
  !
  ! * w_ で始まる関数が返す値はスペクトルデータに同じ.
  !
  ! * xy_ で始まる関数が返す値は 2 次元格子点データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したものことである.
  !
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化 
  !
  ! w_plane_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
  ! 
  !==== 終了処理 
  !
  ! w_plane_Finalize :: モジュールの終了処理(割り付け配列の解放)をおこなう. 
  ! 
  !==== 座標変数
  !
  ! sp_S, sp_Phi   :: 格子点データの 2 次元極座標(格子点データ型 2 次元配列)
  ! sp_X, sp_Y     :: 格子点データの 2 次元 X, Y 座標(格子点データ型 2 次元配列)
  !
  !==== 基本変換
  !
  ! sp_w :: スペクトルデータから格子データへの変換
  ! w_sp :: 格子データからスペクトルデータへの変換
  !
  !==== 微分
  !
  ! w_Lapla_w       :: スペクトルデータにラプラシアンを作用させる
  ! w_LaplaInv_w    :: スペクトルデータにラプラシアンの逆変換を作用させる
  ! sp_Dx_w         :: スペクトルデータに d/dx を作用させる
  ! sp_Dy_w         :: スペクトルデータに d/dy を作用させる
  ! w_Dx_sp         :: 格子データに d/dx を作用させる
  ! w_Dy_sp         :: 格子データに d/dy を作用させる
  !
  !==== 補間
  !
  ! Interpolate_w :: スペクトルデータから任意の点での値を求める. 
  !
  !==== 積分・平均
  !
  ! IntSPhi_sp, AvrSPhi_sp :: 2 次元格子点データの全領域積分および平均
  !
  use dc_message, only : MessageNotify
  use w_module, only : p_Phi=>x_Lon, s_Lat=>y_Lat, &
                       ps_Non_w=>xy_w, w_Non_ps=>w_xy, &
                       w_Initial, w_Finalize, l_nm, nm_l, &  
                       w_Lapla_w, w_LaplaInv_w,w_DPhi_w => w_DLon_w, &
                       xy_GradLambda_w, xy_GradMu_w, w_DivMu_xy, &
                       InterpolateOnLonLat_w => Interpolate_w

  implicit none

  integer               :: im=64            ! 格子点の設定(東西)
  integer               :: jm=32            ! 格子点の設定(南北)
  integer               :: nm=21            ! 切断波数の設定
  real(8)               :: Radius           ! 球半径           

  real(8), allocatable  :: s_S(:)             ! 2 次元極座標動径
  real(8), allocatable  :: ps_S(:,:), ps_Phi(:,:)
  real(8), allocatable  :: s_Mu(:)            ! 球面緯度(sin φ)
  real(8), allocatable  :: ps_Mu(:,:)         ! 球面緯度(sin φ)
  real(8), allocatable  :: ps_X(:,:), ps_Y(:,:) ! 2 次元直交座標
  
  private

  public w_plane_Initial                      ! 初期化
  public w_plane_Finalize                     ! 終了処理

  public s_S, p_Phi                           ! 格子座標(im,jm)
  public ps_S, ps_Phi                         ! 格子座標(im,jm)
  public ps_X, ps_Y                           ! 格子座標(im,jm)
  public ps_w, w_ps                           ! 変換関数

  public l_nm, nm_l                           ! スペクトル成分位置情報

  public ps_Lapla_w, w_LaplaInv_ps            ! ラプラシアンと逆演算
  public w_DPhi_w                             ! 方位角微分
  public ps_DS_w, ps_DPhi_w                   ! 微分(極座標)
  public w_DPhi_ps                            ! 微分(極座標)
  public ps_DX_w, ps_DY_w                     ! 微分(直交座標)
  public w_Lapla_w, w_LaplaInv_w              ! 球面上でのラプラシアンと逆演算

!!$  public w_DX_ps, w_DY_ps                     ! 微分(直交座標)

  public Interpolate_w                        ! ps 座標での補間
  public InterpolateOnXY_w                    ! xy 座標での補間
  public xy_w                                 ! xy 座標での補間
  
  save im, jm, nm, Radius, s_S, ps_S, ps_Phi, s_Mu, ps_Mu, ps_X, ps_Y
  
contains

  !--------------- 初期化 -----------------
  subroutine w_plane_initial(n_in,i_in,j_in,r_in,np_in)
    !
    ! スペクトル変換の格子点数, 波数および OPENMP 使用時の
    ! 最大スレッド数を設定する.
    !
    ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を
    ! しなければならない. 
    !
    ! np_in に 1 より大きな値を指定すれば ISPACK の球面調和函数変換 
    ! OPENMP 並列計算ルーチンが用いられる. 並列計算を実行するには, 
    ! 実行時に環境変数 OMP_NUM_THREADS を np_in 以下の数字に設定する等の
    ! システムに応じた準備が必要となる. 
    !
    ! np_in に 1 より大きな値を指定しなければ並列計算ルーチンは呼ばれない.
    !
    integer,intent(in) :: i_in              !(in) 格子点数(方位角)
    integer,intent(in) :: j_in              !(in) 格子点数(動径)
    integer,intent(in) :: n_in              !(in) 切断波数の設定
    real(8),intent(in) :: r_in              !(in) 変換パラメター(球の半径)
    integer,intent(in), optional :: np_in   !(in) OPENMP での最大スレッド数

    im = i_in  ; jm = j_in  ; nm = n_in ; Radius = r_in

    if ( present (np_in) )then
       call w_initial(n_in,i_in,j_in,np_in)
    else
       call w_initial(n_in,i_in,j_in)
    endif

    allocate(s_S(jm))
    allocate(s_Mu(jm),ps_Mu(0:im-1,jm))
    allocate(ps_S(0:im-1,jm),ps_Phi(0:im-1,jm))
    allocate(ps_X(0:im-1,jm),ps_Y(0:im-1,jm))

    s_Mu = sin(s_Lat) ; ps_Mu = spread(s_Mu,1,im)

    s_S = 2*Radius*sqrt((1+s_Mu)/(1-s_Mu))

    ps_S   = spread(s_S,1,im)
    ps_Phi = spread(p_Phi,2,jm)
    ps_X = ps_S * cos(ps_Phi) ; ps_Y = ps_S * sin(ps_Phi)

    call MessageNotify('M','w_plane_initial',&
         'w_plane_module2 (2016/01/15) is initialized')

  end subroutine w_plane_initial

  !--------------- 終了処理 -----------------
  subroutine w_plane_Finalize
    !
    ! モジュールの終了処理(割り付け配列の解放)をおこなう. 
    !
    ! 解像度を変更する際にはこのサブルーチンを呼んで終了処理を
    ! おこなったのちに再度 w_Initial で初期設定しなければ
    ! ならない. 
    !
    call w_Finalize

    call MessageNotify('M','w_plane_Finalize',&
         'w_plane_module2 (2015/12/26) is finalized')

  end subroutine w_plane_Finalize

  !--------------- 基本変換 -----------------
  function w_ps(ps_Data)
    real(8), intent(IN) :: ps_Data(0:im-1,jm)
    real(8)             :: w_ps((nm+1)**2)

    w_ps = w_Non_ps(ps_Data/(1-ps_Mu))

  end function w_ps
    
  function ps_w(w_Data)
    real(8), intent(IN) :: w_Data((nm+1)**2)
    real(8)             :: ps_w(0:im-1,jm)

    ps_w = (1-ps_Mu)*ps_Non_w(w_Data)

  end function ps_w
    
  !--------------- 微分計算 -----------------
  function ps_Lapla_w(w_data)
    !
    ! 入力スペクトルデータにラプラシアン
    !
    !    ▽^2 = (1/s)(d/ds)[s(d/ds)] + 1/s^2 (d^2/d\lambda^2)
    !
    ! を作用する.
    !
    ! スペクトルデータのラプラシアンとは, 対応する格子点データに
    ! ラプラシアンを作用させたデータのスペクトル変換のことである. 
    !
    real(8)              :: ps_Lapla_w(0:im-1,jm)
    !(out) 入力スペクトルデータのラプラシアン

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

    ps_Lapla_w = (1-ps_Mu)**2/(4*Radius**2) &
         * ( ps_Non_w(w_Lapla_w(w_Non_ps(ps_w(w_data))))) 
!!$    ps_Lapla_w = (1-ps_Mu)**3/(4*Radius**2) * ps_Non_w(w_Lapla_w(w_data)) &
!!$         - (1-ps_Mu)**2/(4*Radius**2) &
!!$           * ps_Non_w(w_DivMu_xy((1-ps_Mu**2)*ps_Non_w(w_data)))

  end function ps_Lapla_w

  function w_LaplaInv_ps(ps_data)
    !
    ! 入力スペクトルデータに逆ラプラシアン
    !
    !    ▽^{-2}
    !      =[(1/s)(d/ds)[s(d/ds)] + 1/s^2 (d^2/d\lambda^2)]^{-1}
    !
    ! を作用する. 
    !
    ! スペクトルデータの逆ラプラシアンとは, 対応する格子点データに
    ! 逆ラプラシアンを作用させたデータのスペクトル変換のことである. 
    !
    real(8)              :: w_LaplaInv_ps((nm+1)*(nm+1))
    !(out) スペクトルデータの逆ラプラシアン

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

    w_LaplaInv_ps = w_ps(ps_Non_w(w_LaplaInv_w(w_Non_ps(ps_data*4*Radius**2/(1-ps_Mu)**2))))

  end function w_LaplaInv_ps

  !--------------- 微分計算 (S,Phi座標系用) -----------------

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

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

    ps_Ds_w = xy_GradMu_w(w_data)/ps_S * (1-ps_Mu) &
            - ps_Non_w(w_data) * (1-ps_Mu**2)/ps_S

  end function ps_DS_w

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

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

    ps_DPhi_w = xy_GradLambda_w(w_data) * (1-ps_Mu)

  end function ps_DPhi_w

  function w_DPhi_ps(ps_data)
    !
    ! 格子点データに方位角微分 d/d\lambda
    ! を作用させてスペクトルデータに変換して返す(1 層用).
    !
    real(8)              :: w_DPhi_ps((nm+1)*(nm+1))
    !(out) 格子点データを発散型経度微分したスペクトルデータ

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

    w_DPhi_ps = w_Non_ps(ps_data/(1-ps_Mu),ipow=0,iflag=-1)

  end function w_DPhi_ps

!!$    function w_sDS_ps(ps_data)
!!$      !
!!$      ! 格子点データに動径微分 s(d/ds) を作用させて
!!$      ! スペクトルデータに変換して返す.
!!$      !
!!$      real(8)              :: w_sDS_ps((nm+1)*(nm+1))
!!$      !(out) 格子点データを発散型緯度微分したスペクトルデータ
!!$
!!$      real(8), intent(in)  :: ps_data(0:im-1,1:jm)
!!$      !(in) 入力格子点データ
!!$
!!$      w_sDS_ps = w_ps(ps_data,ipow=2,iflag=1)
!!$
!!$    end function w_sDS_ps

  !--------------- 微分計算 (xy 座標系用) -----------------
  
  function ps_DX_w(w_data)
    !
    ! スペクトルデータに微分 d/dx
    ! を作用させて格子点データに変換して返す. 
    !
    real(8)              :: ps_DX_w(0:im-1,1:jm)
    !(out) スペクトルデータを勾配型緯度微分した格子点データ

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

    ps_DX_w = cos(ps_Phi)*ps_DS_w(w_data) &
         - sin(ps_Phi)/ps_S * ps_DPhi_w(w_data)

  end function ps_DX_w

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

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

    ps_DY_w = sin(ps_Phi)*ps_DS_w(w_data) &
         + cos(ps_Phi)/ps_S * ps_DPhi_w(w_data)

  end function ps_DY_w

  !--------------- 補間計算 -----------------
  function  Interpolate_w(w_data,pval,sval)
    real(8), intent(IN) :: w_data((nm+1)*(nm+1))  ! 
    real(8), intent(IN) :: pval                   ! 補間する位置(方位角)
    real(8), intent(IN) :: sval                   ! 補間する位置(動径)
    real(8)             :: Interpolate_w          ! 補間値

    real(8) :: lat, t

    t = (sval/(2*Radius))**2
    lat = asin((t-1)/(t+1))

    Interpolate_w = InterpolateOnLonLat_w(w_data,pval,lat)*(1-sin(lat))

  end function Interpolate_w

  function  InterpolateOnXY_w(w_data,xval,yval)
    real(8), intent(IN) :: w_data((nm+1)*(nm+1))  ! 
    real(8), intent(IN) :: xval                  ! 補間する位置(X)
    real(8), intent(IN) :: yval                   ! 補間する位置(Y)
    real(8)             :: InterpolateOnXY_w      ! 補間値

    real(8) :: lon, lat
    real(8) :: sval, t
    real(8), parameter  :: PI=3.141592653589793

    sval = sqrt(xval**2+yval**2)
    t = (sval/(2*Radius))**2
    lat = asin((t-1)/(t+1))

    if ( xval .eq. 0.0D0 .AND. yval .eq. 0.0D0 ) then
       lon = 0.0D0
    else if ( xval .eq. 0.0D0 .AND. yval .gt. 0.0D0 ) then
       lon = PI/2.0D0
    else if ( xval .eq. 0.0D0 .AND. yval .lt. 0.0D0 ) then
       lon = 3.0D0*PI/2.0D0
    else if ( xval .gt. 0.0D0 .AND. yval .ne. 0.0D0 ) then
       lon = atan(yval/xval) 
    else if ( xval .lt. 0.0D0 .AND. yval .ne. 0.0D0 ) then
       lon = atan(yval/xval) + PI
    endif

    if ( -1.0D0/2.0D0*PI .lt. lon .AND. lon .lt. 0.0 ) then
       lon = lon + PI
    endif

    InterpolateOnXY_w = InterpolateOnLonLat_w(w_data,lon,lat)*(1-sin(lat))

  end function InterpolateOnXY_w

  subroutine xy_w( w_data, imx, jmy, xmin, xmax, ymin, ymax, &
       xy_Data, x_X, y_Y)

    real(8), intent(IN) :: w_data((nm+1)*(nm+1)) ! 入力スペクトルデータ
    integer, intent(IN) :: imx                   ! X 格子点数
    integer, intent(IN) :: jmy                   ! Y 格子点数
    real(8), intent(IN) :: xmin, xmax            ! X 座標範囲
    real(8), intent(IN) :: ymin, ymax            ! Y 座標

    real(8), intent(OUT) :: xy_Data(imx,jmy)     ! XY データ
    real(8), intent(OUT) :: x_X(imx)             ! X 格子点座標値
    real(8), intent(OUT) :: y_Y(jmy)             ! Y 格子点座標値

    integer :: i,j

    do i=1,imx
       x_X(i) = xmin + (xmax-xmin)/(imx-1)*(i-1)
    end do

    do j=1,jmy
       y_Y(j) = ymin + (ymax-ymin)/(jmy-1)*(j-1)
    end do

    do j=1,jmy
       do i=1,imx
          xy_Data(i,j) = InterpolateOnXY_w(w_data,x_X(i),y_Y(j))
       enddo
    enddo

  end subroutine xy_w
    
end module w_plane_module2
