!----------------------------------------------------------------------
! Copyright (c) 2012--2016 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!表題  w_base_mpi_module_sjpack
!
!  spml/w_base_mpi_module_sjpack モジュールは球面上での 2 次元流体運動を
!  球面調和函数を用いたスペクトル法と MPI によって数値計算するための 
!  モジュール w_mpi_module_sjpack の下部モジュールであり, スペクトル計算の
!  基本的な Fortran90 関数を提供する. 
!
!  内部で ISPACK の SJPACK-MPI の Fortran77 サブルーチンを呼んでいる. 
!  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!  ついては ISPACK/SJPACK-MPI のマニュアルを参照されたい.
!
!
!履歴  2012/03/31  竹広真一  w_base_mpi_module を sjpack 化
!      2012/04/02  竹広真一  格子点数等を private 変数に
!      2013/02/12  竹広真一  w_StreamPotential2VectorMPI,  
!                            w_Vector2VorDivMPI 導入
!      2013/02/15  竹広真一  w_VectorCosLat2VorDivMPI 導入
!      2013/02/23  竹広真一  w_base_mpi_Finalize 導入
!      2016/07/07  竹広真一  AVX 用 32 bit alignment 対応
!
module w_base_mpi_module_sjpack
  !
  ! w_base_mpi_module_sjpack
  !
  !  spml/w_base_mpi_module_sjpack モジュールは球面上での 2 次元流体運動を
  !  球面調和函数を用いたスペクトル法と MPI によって数値計算するための 
  !  モジュール w_mpi_module_sjpack の下部モジュールであり, スペクトル法の
  !  基本的なな Fortran90 関数を提供する. 
  !
  !  内部で ISPACK の SJPACK-MPI の Fortran77 サブルーチンを呼んでいる. 
  !  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  !  ついては ISPACK/SJPACK のマニュアルを参照されたい.
  !
  use dc_message
  use iso_c_binding
  use w_base_module_sjpack, only : im, jm, nn, nm, mm, openmp, np, x_Lon

  implicit none

  real(8), allocatable  :: p(:), r(:)       ! 変換用配列
  integer               :: it(4)            ! 変換用配列

  real(8), allocatable  :: c(:)             ! 作業配列

  integer               :: jc               ! 分散配置用変数

  real(8), allocatable  :: y(:,:)           ! 作業配列

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

  real(8), pointer      :: xv_pdata(:,:),q(:),t(:),w(:),wg(:)
  type(C_PTR)           :: pg, pq, pt, pw, pwg

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

  private
  private im, jm, nn, mm, nm                  ! 格子点数, 切断波数

  public it, t, p, r                          ! 変換用作業配列
  public openmp, np                           ! OPENMP 用変数
  public jc                                   ! 分散配置用変数

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

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

  save it, p, r                               ! 変換用配列を記憶
  save c                                      ! 変換用配列の大きさ
  save jc                                     ! 分散格子点配列の大きさ
  save xv_pdata, q, t, w, wg
  save w_base_initialize                      ! 初期化フラグ


  contains
  !--------------- 初期化 -----------------
    subroutine w_base_mpi_Initial
      !
      ! スペクトル変換の格子点数, 波数および OPENMP 使用時の
      ! 最大スレッド数を設定する.
      !
      ! 実際の使用には上位サブルーチン w_Initial を用いること.
      !
      integer :: i, j

      allocate(p(jm/2*(mm+4)))                ! 変換用配列
      allocate(r((mm+1)*(2*nm-mm-1)+1))       ! 変換用配列

      allocate(c((mm+1)*(mm+1)))              ! 変換用作業配列

      call mlallc(pt,im*6)
      call c_f_pointer(pt, t, [im*6])

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

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

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

      allocate(y(jc/2,mm+4))
      y = reshape(p(1:jc/2*(mm+4)),(/jc/2,mm+4/))

      do j=1,jc/2
         v_Lat(jc/2+j)   =  asin(y(j,1))        ! 緯度座標
         v_Lat(jc/2-j+1) = -asin(y(j,1))        ! 緯度座標
         v_Lat_Weight(jc/2+j)   = 2*y(j,2)      ! 緯度重み(Gauss grid)
         v_Lat_Weight(jc/2-j+1) = 2*y(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

      call mlallc(pg,jc*im)
      call mlallc(pq,jm/2*7*np)
      call mlallc(pw,(jc+1)*im)
      call mlallc(pwg,(jc+2)*im)
      call c_f_pointer(pg, xv_pdata, [IM,JC])
      call c_f_pointer(pq, q, [jm/2*7*np])
      call c_f_pointer(pw, w, [(jc+1)*im])
      call c_f_pointer(pwg, wg, [(jc+2)*im])      

      w_base_initialize = .true.

      call MessageNotify('M','w_base_mpi_initial',&
                         'w_base_mpi_module_sjpack (2016/07/07) is initialized')
    end subroutine w_base_mpi_Initial

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

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

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

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

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

      integer ipval, ifval

      real(8)             :: w_Rdata((2*nn+1-mm)*mm+nn+1)
      ! 作業用スペクトルデータ(SJTS2G 出力用)
      real(8)             :: w_Xdata((mm+1)*(mm+1))
      ! 作業用スペクトルデータ(SJCS2X 出力用)
      real(8)             :: w_Ydata((mm+4)*mm+2)
      ! 作業用スペクトルデータ(SJCS2Y 出力用)

      real(8)  :: ws(2*(nn+1)*np)            ! 変換用作業配列
      real(8)  :: ws2(2*(nm+1)*np)           ! 変換用作業配列

      logical :: first=.true.                    ! 初回判定スイッチ
      save first

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

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

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

      if ( openmp .and. first ) then
         call MessageNotify('M','xy_w', &
              'OpenMP routine SJTSOG/SJPACK-MPI is used for spherical harmonic transformation.')
      endif

      if ( ifval==0 ) then
         call sjcrup(mm,nn,w_data,w_Rdata)
         if ( openmp ) then
            call sjtsog(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         else
            call sjts2g(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         endif
         xv_w = xv_pdata
      else if( ifval==-1 ) then
         call sjcs2x(mm,w_data,w_Xdata)
         call sjcrup(mm,nn,w_Xdata,w_Rdata)
         if ( openmp ) then
            call sjtsog(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         else
            call sjts2g(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         endif
         xv_w = xv_pdata
      else if( ifval==1 ) then
         call sjcs2y(mm,w_data,w_Ydata,c)
         if ( openmp ) then
            call sjtsog(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
                        it,t,p,q,r,ws2,wg,w,ipval)
         else
            call sjts2g(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
                        it,t,p,q,r,ws2,wg,w,ipval)
         endif
         xv_w = xv_pdata
      else if( ifval==2 ) then
         call sjcrup(mm,nn,w_data,w_Rdata)
         if ( openmp ) then
            call sjtsog(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         else
            call sjts2g(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         endif
         xv_w = xv_pdata * sin(xv_Lat)
      else
         call MessageNotify('E','xv_w','invalid value of iflag')
      endif

      first = .false.

    end function xv_w

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

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

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

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


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

      integer ipval, ifval

      real(8)             :: w_Rdata((2*nn+1-mm)*mm+nn+1)
      ! 作業用スペクトルデータ(SJTS2G 出力用)
      real(8)             :: w_Xdata((mm+1)*(mm+1))
      ! 作業用スペクトルデータ(SJCS2X 出力用)
      real(8)             :: w_Ydata((mm+4)*nm+2)
      ! 作業用スペクトルデータ(SJCY2S 出力用)

      real(8)  :: ws(2*(nn+1)*np+(2*NN+1-MM)*MM+NN+1)  ! 変換用作業配列
      real(8)  :: ws2(2*(nm+1)*np+(2*NM+1-MM)*MM+NM+1) ! 変換用作業配列

      logical :: first=.true.                     ! 初回判定スイッチ
      save first

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

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

      if (present(iflag)) then
         ifval = iflag
      else
         ifval = iflag_default
      endif
      
      if ( openmp .and. first ) then
         call MessageNotify('M','w_xv', &
              'OpenMP routine SJPGOS/SJPACK-MPI is used for spherical harmonic transformation.')
      endif

      if ( ifval == 0 ) then
         xv_pdata = xv_data
         if ( openmp ) then
            call sjpgos(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         else
            call sjpg2s(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         endif
         call sjcrdn(mm,nn,w_Rdata,w_xv)
      else if ( ifval == -1 ) then
         xv_pdata = xv_data
         if ( openmp ) then
            call sjpgos(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         else
            call sjpg2s(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         endif
         call sjcrdn(mm,nn,w_Rdata,w_Xdata)
         call sjcs2x(mm,w_Xdata,w_xv)
      else if ( ifval == 1 ) then
         xv_pdata = xv_data
         if ( openmp ) then
            call sjpgos(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
                        it,t,p,q,r,ws2,wg,w,ipval)
         else
            call sjpg2s(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
                        it,t,p,q,r,ws2,wg,w,ipval)
         endif
         call sjcy2s(mm,w_Ydata,w_xv,c)
      else if ( ifval == 2 ) then
         xv_pdata = xv_data * sin(xv_Lat)
         if ( openmp ) then
            call sjpgos(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         else
            call sjpg2s(mm,nm,nn,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws,wg,w,ipval)
         endif
         call sjcrdn(mm,nn,w_Rdata,w_xv)
      end if

      first = .false.

    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((mm+1)*(mm+1))
      !(in) 流線関数
      real(8), intent(in)   :: w_Chi((mm+1)*(mm+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) 速度緯度成分

      real(8)             :: w_Rdata((mm+4)*mm+2)
      ! 作業用スペクトルデータ(SJTS2G 出力用)
      real(8)             :: w_Xdata((mm+1)*(mm+1))
      ! 作業用スペクトルデータ(SJCS2X 出力用)
      real(8)             :: w_Ydata((mm+4)*mm+2)
      ! 作業用スペクトルデータ(SJCS2Y 出力用)

      real(8)  :: ws2(2*(nm+1)*np)           ! 変換用作業配列

      logical :: first=.true.                     ! 初回判定スイッチ
      save first

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

      if ( openmp .and. first ) then
         call MessageNotify('M','w_StreamPotential2Vector', &
              'OpenMP routine SJTSOG/SNPACK-MPI is used for spherical harmonic transformation.')
      endif
      first = .false.

      !
      !   u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ の計算
      !
      call sjcs2x(mm,w_Chi,w_Xdata)
      call sjcs2y(mm,w_Psi,w_Ydata,c)
      call sjcrup(mm,nm,w_Xdata,w_Rdata)
      w_Rdata = w_Rdata - w_Ydata
      !
      !   u の計算
      !
      if ( openmp ) then
         call sjtsog(mm,nm,nm,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws2,wg,w,1)
      else
         call sjts2g(mm,nm,nm,im,jc,w_Rdata,xv_pdata,&
              it,t,p,q,r,ws2,wg,w,1)
      endif
      xv_U = xv_pdata
      !
      !   v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ の計算
      !
      call sjcs2y(mm,w_Chi,w_Ydata,c)
      call sjcs2x(mm,w_Psi,w_Xdata)
      call sjcrup(mm,nm,w_Xdata,w_Rdata)
      w_Rdata= w_Rdata + w_Ydata
      !
      !   v の計算
      !
      if ( openmp ) then
         call sjtsog(mm,nm,nm,im,jc,w_Rdata,xv_pdata,&
                        it,t,p,q,r,ws2,wg,w,1)
      else
         call sjts2g(mm,nm,nm,im,jc,w_Rdata,xv_pdata,&
              it,t,p,q,r,ws2,wg,w,1)
      endif
      xv_V = xv_pdata

    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:jc)
      !(in) 速度経度成分
      real(8), intent(in)   :: xv_V(0:im-1,1:jc)
      !(in) 速度緯度成分

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

      real(8)             :: w_Xdata((mm+1)*(mm+1))
      ! 作業用スペクトルデータ(SJCS2X 出力用)
      real(8)             :: w_Ydata((mm+4)*nm+2)
      ! 作業用スペクトルデータ(SJCY2S 出力用)

      real(8)  :: w_Data1((mm+1)*(mm+1))
      real(8)  :: w_Data2((mm+1)*(mm+1))

      real(8)  :: ws2(2*(nm+1)*np+(2*NM+1-MM)*MM+NM+1) ! 変換用作業配列

      logical :: first=.true.                     ! 初回判定スイッチ
      save first

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

      if ( openmp .and. first ) then
         call MessageNotify('M','w_Vector2VorDiv', &
              'OpenMP routine SJPGOS/SNPACK-MPI is used for spherical harmonic transformation.')
      endif
      first = .false.

      !
      ! 1/cosφ∂u/∂λ, 1/cosφ ∂(u cosφ)/∂φ の計算
      !
      xv_pdata = xv_U
      if ( openmp ) then
         call sjpgos(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
              it,t,p,q,r,ws2,wg,w,1)
      else
         call sjpg2s(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
              it,t,p,q,r,ws2,wg,w,1)
      endif
      call sjcrdn(mm,nm,w_Ydata,w_Xdata)
      call sjcs2x(mm,w_Xdata,w_Div)
      call sjcy2s(mm,w_Ydata,w_Data1,c)
      !
      ! 1/cosφ∂v/∂λ, 1/cosφ ∂(v cosφ)/∂φ の計算
      !
      xv_pdata = xv_V
      if ( openmp ) then
         call sjpgos(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
              it,t,p,q,r,ws2,wg,w,1)
      else
         call sjpg2s(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
              it,t,p,q,r,ws2,wg,w,1)
      endif
      call sjcrdn(mm,nm,w_Ydata,w_Xdata)
      call sjcs2x(mm,w_Xdata,w_Vor)
      call sjcy2s(mm,w_Ydata,w_Data2,c)
      !
      !  渦度・発散の計算
      !
      !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ 
      !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
      !
      w_Vor = w_Vor - w_Data1
      w_Div = w_Div + w_Data2

    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:jc)
      !(in) 速度経度成分 * cos(lat)
      real(8), intent(in)   :: xv_VCosLat(0:im-1,1:jc)
      !(in) 速度緯度成分 * cos(lat)

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

      real(8)             :: w_Xdata((mm+1)*(mm+1))
      ! 作業用スペクトルデータ(SJCS2X 出力用)
      real(8)             :: w_Ydata((mm+4)*nm+2)
      ! 作業用スペクトルデータ(SJCY2S 出力用)

      real(8)  :: w_Data1((mm+1)*(mm+1))
      real(8)  :: w_Data2((mm+1)*(mm+1))

      real(8)  :: ws2(2*(nm+1)*np+(2*NM+1-MM)*MM+NM+1) ! 変換用作業配列

      logical :: first=.true.                     ! 初回判定スイッチ
      save first

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

      if ( openmp .and. first ) then
         call MessageNotify('M','w_VectorCosLat2VorDiv', &
              'OpenMP routine SJPGOS/SNPACK-MPI is used for spherical harmonic transformation.')
      endif
      first = .false.

      !
      ! 1/cosφ∂u/∂λ, 1/cosφ ∂(u cosφ)/∂φ の計算
      !
      xv_pdata = xv_UCosLat
      if ( openmp ) then
         call sjpgos(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
              it,t,p,q,r,ws2,wg,w,2)
      else
         call sjpg2s(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
              it,t,p,q,r,ws2,wg,w,2)
      endif
      call sjcrdn(mm,nm,w_Ydata,w_Xdata)
      call sjcs2x(mm,w_Xdata,w_Div)
      call sjcy2s(mm,w_Ydata,w_Data1,c)
      !
      ! 1/cosφ∂v/∂λ, 1/cosφ ∂(v cosφ)/∂φ の計算
      !
      xv_pdata = xv_VCosLat
      if ( openmp ) then
         call sjpgos(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
              it,t,p,q,r,ws2,wg,w,2)
      else
         call sjpg2s(mm,nm,nm,im,jc,w_Ydata,xv_pdata,&
              it,t,p,q,r,ws2,wg,w,2)
      endif
      call sjcrdn(mm,nm,w_Ydata,w_Xdata)
      call sjcs2x(mm,w_Xdata,w_Vor)
      call sjcy2s(mm,w_Ydata,w_Data2,c)
      !
      !  渦度・発散の計算
      !
      !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ 
      !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
      !
      w_Vor = w_Vor - w_Data1
      w_Div = w_Div + w_Data2

    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_sjpack not initialized yet')
         return
      endif

      deallocate(p)              ! 変換用配列
      deallocate(r)              ! 変換用配列
      deallocate(c)              ! 変換用作業配列

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

      nullify(xv_pdata)              ! 変換用作業配列
      nullify(q)                     ! 変換用作業配列
      nullify(t)                     ! 変換用作業配列
      nullify(w)                     ! 変換用作業配列
      nullify(wg)                    ! 変換用作業配列
      call mlfree(pg)                ! 変換用作業配列
      call mlfree(pq)                ! 変換用作業配列
      call mlfree(pt)                ! 変換用作業配列
      call mlfree(pw)                ! 変換用作業配列
      call mlfree(pwg)               ! 変換用作業配列      

      w_base_initialize = .false.

      call MessageNotify('M','w_base_mpi_Finalize',&
           'w_base_mpi_module_sjpack (2016/07/07) for AVX is finalized')

    end subroutine w_base_mpi_Finalize

end module w_base_mpi_module_sjpack
