!----------------------------------------------------------------------
! Copyright (c) 2008-2013 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!表題  w_base_mpi_module
!
!  spml/w_base_mpi_module モジュールは球面上での 2 次元流体運動を
!  球面調和函数を用いたスペクトル法と MPI によって数値計算するための 
!  モジュール w_mpi_module の下部モジュールであり, スペクトル計算の
!  基本的な Fortran90 関数を提供する. 
!
!  内部で ISPACK の SPPACK と SNPACK の Fortran77 サブルーチンを呼んでいる. 
!  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!  ついては ISPACK/SNPACK,SPPACK のマニュアルを参照されたい.
!
!
!履歴  2008/05/26  竹広真一  w_base_module を MPI 化
!      2010/01/07  佐々木洋平  RDoc 用のドキュメント修正, 
!      2012/03/30  竹広真一  コメント修正
!      2013/02/12  竹広真一  w_StreamPotential2VectorMPI,  
!                            w_Vector2VorDivMPI 導入
!      2013/02/15  竹広真一  w_VectorCosLat2VorDivMPI 導入
!      2013/02/23  竹広真一  w_base_mpi_Finalize 導入
!
module w_base_mpi_module
  !
  ! w_base_mpi_module
  !
  !  spml/w_base_mpi_module モジュールは球面上での 2 次元流体運動を
  !  球面調和函数を用いたスペクトル法と MPI によって数値計算するための 
  !  モジュール w_mpi_module の下部モジュールであり, スペクトル法の
  !  基本的なな Fortran90 関数を提供する. 
  !
  !  内部で ISPACK の SPPACK と SNPACK の Fortran77 サブルーチンを呼んでいる. 
  !  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  !  ついては ISPACK/SNPACK,SPPACK のマニュアルを参照されたい.
  !
  use dc_message
  use w_base_module, only : im, jm, nm, x_Lon

  implicit none

  integer               :: it(6)            ! 変換用配列(分散格子点用)
  real(8), allocatable  :: t(:)             ! 変換用配列(分散格子点用)
  integer, allocatable  :: ip(:)            ! 変換用配列(分散格子点用)
  real(8), allocatable  :: p(:), r(:)       ! 変換用配列(分散格子点用)
  integer, allocatable  :: ia(:)            ! 変換用配列(分散格子点用)
  real(8), allocatable  :: a(:)             ! 変換用配列(分散格子点用)
  real(8), allocatable  :: y(:)             ! 変換用配列(分散格子点用)

  integer               :: jc               ! 分散配置用変数
  real(8), allocatable  :: yy(:,:)          ! 変換用配列
  
  real(8), allocatable  :: q(:)             ! 作業配列
  real(8), allocatable  :: ww(:), ws(:)     ! 作業配列
  real(8), allocatable  :: w(:)             ! 作業配列

  real(8), allocatable  :: v_Lat(:),v_Lat_Weight(:)      ! 緯度経度

  real(8), allocatable  :: xv_Lon(:,:), xv_Lat(:,:)

  real(8), allocatable  :: xv_work(:,:)     ! w_xv,xv_w 変換用配列

  integer               :: id=65, jd=33     ! xv_work の大きさ

  real(8), parameter    :: pi=3.1415926535897932385D0

  logical               :: w_base_initialize=.false.   ! 初期化フラッグ

  private
  private im, jm, nm                          ! Intel Fortran 対策

  public it, t, y, ip, p, r, ia, a            ! 変換用作業配列
  public id, jd                               ! 作業用配列の大きさ
  public jc                                   ! 分散配置情報

  public w_base_mpi_Initial                   ! 初期化サブルーチン
  public w_base_mpi_Finalize                  ! 終了サブルーチン

  public v_Lat, v_Lat_Weight                  ! 経度分散格子座標・重み
  public xv_Lon, xv_Lat                       ! 分散格子座標(im,jc)
  public xv_w, w_xv                           ! 変換関数
  public w_StreamPotential2VectorMPI          ! 流線ポテンシャルから速度場計算
  public w_Vector2VorDivMPI                   ! 速度場から渦度発散を計算
  public w_VectorCosLat2VorDivMPI             ! 速度場から渦度発散を計算

  save it, t, y, ip, p, r, ia, a              ! 変換用配列を記憶
  save jc                                     ! 分散格子点配列の大きさ
  save id, jd                                 ! 変換用配列の大きさ
  save w_base_initialize                      ! 初期化フラグ

  contains
  !--------------- 初期化 -----------------
    subroutine w_base_mpi_Initial
      !
      ! スペクトル変換の格子点数, 波数を設定する.
      !
      ! 実際の使用には上位サブルーチン w_mpi_Initial を用いること.
      !
      integer :: iw, i, j

      allocate(t(im*2))                       ! 変換用配列(分散配置)
      allocate(ip(((nm+1)/2+nm+1)*2))         ! 変換用配列(分散配置)
      allocate(p(((nm+1)/2+nm+1)*jm))         ! 変換用配列(分散配置)
      allocate(r(((nm+1)/2*2+3)*(nm/2+1)))    ! 変換用配列(分散配置)
      allocate(ia((nm+1)*(nm+1)*4))           ! 変換用配列(分散配置)
      allocate(a((nm+1)*(nm+1)*6))            ! 変換用配列(分散配置)
      allocate(y(jm*2))                       ! 変換用配列(分散配置)

      ! 注意 : 別ルーチンによって w_base_Initial が呼んであることを仮定
      call snmini(nm,im,jm,jc,it,t,y,ip,p,r,ia,a)

      if ( im/2*2 .eq. im ) then
         id = im+1 
      else
         id = im
      endif
      if ( jc/2*2 .eq. jc ) then
         jd = jc+1
      else
         jd = jc
      endif
      allocate(xv_work(id,jd))                ! 変換用配列

      allocate(q(((nm+1)/2+nm+1)*jm))         ! 作業配列
 
      iw=max((nm+4)*(nm+3),jd*3*(nm+1),jd*im)
      allocate(ws(iw),ww(iw), w((nm+1)*(nm+1)))    ! 作業用配列
      allocate(yy(jc/2,4))                         ! 変換用配列

      allocate(v_Lat(jc),v_Lat_Weight(jc))             ! 格子点座標格納配列

      allocate(xv_Lon(0:im-1,jc),xv_Lat(0:im-1,jc))   ! 格子点座標格納配列

      yy = reshape(y(1:2*jc),(/jc/2,4/))

      do j=1,jc/2
         v_Lat(jc/2+j)   =  asin(yy(j,1))        ! 緯度座標
         v_Lat(jc/2-j+1) = -asin(yy(j,1))        ! 緯度座標
         v_Lat_Weight(jc/2+j)   = 2*yy(j,2)      ! 緯度重み(Gauss grid)
         v_Lat_Weight(jc/2-j+1) = 2*yy(j,2)      ! 緯度重み(Gauss grid)
      enddo
  
      do j=1,jc
         xv_Lon(:,j) = x_Lon
      enddo

      do i=0,im-1
         xv_Lat(i,:) = v_Lat
      enddo

      w_base_initialize = .true.

      call MessageNotify('M','w_base_mpi_initial',&
                         'w_base_mpi_module (2013/02/23) is initialized')
    end subroutine w_base_mpi_Initial

  !--------------- 基本変換(格子点分散配置) -----------------

    function xv_w(w_data,ipow,iflag)
      !
      ! スペクトルデータから分散格子データへ変換する(1 層用).
      !
      real(8)               :: xv_w(0:im-1,jc)
      !(out) 格子点データ

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

      integer, intent(in), optional  :: ipow      
      !(in) 作用させる 1/cosφ の次数. 省略時は 0. 

      integer, intent(in), optional  :: iflag
      !(in) 変換の種類
      !     0 : 通常の正変換
      !     1 : 経度微分を作用させた正変換
      !    -1 : 緯度微分を作用させた正変換
      !     2 : sinφを作用させた正変換
      !     省略時は 0.
      !
      integer, parameter  :: ipow_default  = 0
      integer, parameter  :: iflag_default = 0

      integer ipval, ifval

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

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

      call snts2g(nm,im,id,jc,jd,1,w_data,xv_work,&
                  it,t,y,ip,p,r,ia,a,q,ws,ww,ipval,ifval)

      xv_w=xv_work(1:im,1:jc)

    end function xv_w

    function w_xv(xv_data,ipow,iflag)
      !
      ! 格子データからスペクトルデータへ(正)変換する(1 層用).
      !
      real(8)               :: w_xv((nm+1)*(nm+1))
      !(out) スペクトルデータ

      real(8), intent(in)   :: xv_data(0:im-1,jc)
      !(in) 格子点データ

      integer, intent(in), optional  :: ipow
      !(in) 変換時に同時に作用させる 1/cosφ の次数. 省略時は 0.

      integer, intent(in), optional  :: iflag
      ! 変換の種類
      !     0 : 通常の正変換
      !     1 : 経度微分を作用させた正変換
      !    -1 : 緯度微分を作用させた正変換
      !     2 : sinφを作用させた正変換
      !   省略時は 0.


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

      integer ipval, ifval

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

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

      xv_work(1:im,1:jc)=xv_data

      call sntgms(nm,im,id,jc,jd,1,xv_work,w_xv,&
                 it,t,y,ip,p,r,ia,a,q,ws,ww,ipval,ifval,w)

    end function w_xv

  !----------- 速度, 渦度・発散, 流線・ポテンシャル計算 -------------

    subroutine w_StreamPotential2VectorMPI(w_Psi, w_Chi, xv_U, xv_V)
      !
      ! 流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に
      ! (逆)変換する(1 層用, MPI)
      !
      ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. 
      !
      !   u cosφ =      ∂χ/∂λ - cosφ∂ψ/∂φ, 
      !   v cosφ = cosφ∂χ/∂φ +      ∂ψ/∂λ 
      !
      real(8), intent(in)   :: w_Psi((nm+1)*(nm+1))
      !(in) 流線関数
      real(8), intent(in)   :: w_Chi((nm+1)*(nm+1))
      !(in) 速度ポテンシャル

      real(8), intent(out)   :: xv_U(0:im-1,1:jc)
      !(out) 速度経度成分
      real(8), intent(out)   :: xv_V(0:im-1,1:jc)
      !(out) 速度緯度成分

      if ( .not. w_base_initialize ) then
         call MessageNotify('E','w_StreamPotential2VectorMPI',&
              'w_base_module not initialize yet.')
      endif
      !
      !   u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ の計算
      !
      xv_U = xv_w(w_Chi,ipow=1,iflag=-1) - xv_w(w_Psi,ipow=1,iflag=1) 
      !
      !   v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ の計算
      !
      xv_V = xv_w(w_Chi,ipow=1,iflag=1) + xv_w(w_Psi,ipow=1,iflag=-1) 

    end subroutine w_StreamPotential2VectorMPI

    subroutine w_Vector2VorDivMPI(xv_U, xv_V, w_Vor, w_Div)
      !
      ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に
      ! (正)変換する(1 層用, MPI)
      ! 
      ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. 
      !
      !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ 
      !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
      !
      real(8), intent(in)   :: xv_U(0:im-1,1:jm)
      !(in) 速度経度成分
      real(8), intent(in)   :: xv_V(0:im-1,1:jm)
      !(in) 速度緯度成分

      real(8), intent(out)   :: w_Vor((nm+1)*(nm+1))
      !(out) 流線関数
      real(8), intent(out)   :: w_Div((nm+1)*(nm+1))
      !(out) 速度ポテンシャル

      if ( .not. w_base_initialize ) then
         call MessageNotify('E','w_Vector2VorDivMPI',&
              'w_base_module not initialize yet.')
      endif

      !
      !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ 
      !
      w_Vor = w_xv(xv_V,ipow=1,iflag=-1) - w_xv(xv_U,ipow=1,iflag=1)
      !
      !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
      !
      w_Div = w_xv(xv_U,ipow=1,iflag=-1) + w_xv(xv_V,ipow=1,iflag=1)

    end subroutine w_Vector2VorDivMPI

    subroutine w_VectorCosLat2VorDivMPI(xv_UCosLat, xv_VCosLat, w_Vor, w_Div)
      !
      ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に
      ! (正)変換する(1 層用, MPI)
      ! 
      ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. 
      !
      !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ 
      !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
      !
      real(8), intent(in)   :: xv_UCosLat(0:im-1,1:jm)
      !(in) 速度経度成分 * cos(lat)
      real(8), intent(in)   :: xv_VCosLat(0:im-1,1:jm)
      !(in) 速度緯度成分 * cos(lat)

      real(8), intent(out)   :: w_Vor((nm+1)*(nm+1))
      !(out) 流線関数
      real(8), intent(out)   :: w_Div((nm+1)*(nm+1))
      !(out) 速度ポテンシャル

      if ( .not. w_base_initialize ) then
         call MessageNotify('E','w_VectorCosLat2VorDivMPI',&
              'w_base_module not initialize yet.')
      endif

      !
      !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ 
      !
      w_Vor =   w_xv(xv_VCosLat,ipow=2,iflag=-1) &
              - w_xv(xv_UCosLat,ipow=2,iflag=1)
      !
      !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
      !
      w_Div =   w_xv(xv_UCosLat,ipow=2,iflag=-1) &
              + w_xv(xv_VCosLat,ipow=2,iflag=1)

    end subroutine w_VectorCosLat2VorDivMPI

  !--------------- 終了処理 -----------------
    subroutine w_base_mpi_Finalize
      !
      ! モジュールの終了処理(割り付け配列の解放)をおこなう. 
      !
      ! 実際の使用には上位サブルーチン w_Finalize を用いること.
      !
      if ( .not. w_base_initialize ) then
         call MessageNotify('W','w_base_mpi_Finalize',&
              'w_base_mpi_module not initialized yet')
         return
      endif

      deallocate(t)            ! 変換用配列(分散配置)
      deallocate(ip)           ! 変換用配列(分散配置)
      deallocate(p)            ! 変換用配列(分散配置)
      deallocate(r)            ! 変換用配列(分散配置)
      deallocate(ia)           ! 変換用配列(分散配置)
      deallocate(a)            ! 変換用配列(分散配置)
      deallocate(y)            ! 変換用配列(分散配置)

      deallocate(xv_work)      ! 変換用配列
      deallocate(q)            ! 作業配列
      deallocate(ws,ww,w)      ! 作業用配列
      deallocate(yy)           ! 変換用配列

      deallocate(v_Lat,v_Lat_Weight)   ! 格子点座標格納配列
      deallocate(xv_Lon,xv_Lat)        ! 格子点座標格納配列

      w_base_initialize = .false.

      call MessageNotify('M','w_base_mpi_Finalize',&
           'w_base_mpi_module (2013/02/23) is finalized')

    end subroutine w_base_mpi_Finalize

  end module w_base_mpi_module
