!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2006. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine DynImpFunc
!
! Developer:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: dynimpfunc.f90.lapack,v 1.7 2008-06-26 10:56:17 odakker2 Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2003. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!陰解法を用いた力学過程の各項の計算モジュール.
!エクスナー関数を陰解法で解く際に必要となる, 係数行列の要素の決定,
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * 離散化する際, 上下境界条件として鉛直速度が零を与えている.
!  * 空間方向に 2 次精度の離散化を陽に利用しているため, differentiate_center4 モジュールを指定することはできない.
!
!== Future Plans
!

module DynImpFunc
  !
  !陰解法を用いた力学過程の各項の計算モジュール.
  !エクスナー関数を陰解法で解く際に必要となる, 係数行列の要素の決定,
  !
  !本モジュールでは科学計算ライブラリとして, 以下をサポートする.
  !  * MATRIX/MPP (HITACH 最適化 Fortran)
  !  * LAPACK
  !

  !モジュール読み込み
  use dc_types, only: DP
  use dc_iounit, only: FileOpen
  use dc_message, only: MessageNotify     !メッセージ出力

  use gridset, only: &
    &                imin,            &! x 方向の配列の下限
    &                imax,            &! x 方向の配列の上限
    &                jmin,            &! y 方向の配列の下限
    &                jmax,            &! y 方向の配列の上限
    &                kmin,            &! z 方向の配列の下限
    &                kmax,            &! z 方向の配列の上限
    &                nx,              &! x 方向の物理領域の下限
    &                ny,              &! x 方向の物理領域の上限
    &                nz                ! y 方向の物理領域の下限
  use axesset, only: z_dz, r_dz, xyz_dx, xyr_avr_xyz ! 格子間隔
  use timeset, only: DelTimeShort         !短い時間ステップ
  use damping, only: DampSound            !音波の減衰係数
  use constants,only: CpDry               !乾燥成分の比熱
  use basicset, only: xyz_VelSoundBasicZ, &!基本場の音速
    &                 xyz_DensBasicZ,     &!基本場の密度
    &                 xyz_PotTempBasicZ    !基本場の温位
  use xyz_deriv_module, only: xyr_dz_xyz, xyz_dz_xyr, xyz_dx_pyz, xyz_dy_xqz

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public に設定
  public xyz_Exner_init     !初期化ルーチン
  public xyz_Exner          !エクスナー関数の計算
  public xyr_GradPi         !陰解法を用いた圧力傾度

  real(DP)               :: beta  = 5.0d-1    !クランクニコルソン法なら 0.5
                                             !完全陰解法なら 1
  real(DP), allocatable  :: xyz_F1BasicZ(:,:,:)  !係数行列の計算に用いる配列
  real(DP), allocatable  :: xyr_F2BasicZ(:,:,:)  !係数行列の計算に用いる配列
  real(DP), allocatable  :: xyz_VPotTempBasicZ(:,:,:)
                                             !基本場の仮温位

  integer               :: N = 10            !係数行列/改行列の次数, 整合寸法
  integer               :: M = 10            !方程式の組数
  integer               :: NUD = 1           !係数行列の上三角部分の帯幅
  integer               :: NLD = 1           !係数行列の下三角部分の帯幅
  integer               :: NAL = 1           !LU 分解の結果 L の整合寸法
  integer               :: NA = 3            !NUD + NLD + 1
  integer               :: unit              !装置番号
  integer               :: ix, jy, kz        !ループ変数

  real(DP), allocatable :: A(:)              !係数行列の対角成分
  real(DP), allocatable :: B(:)              !係数行列の上三角部分
  real(DP), allocatable :: C(:)              !係数行列の下三角部分
  real(DP), allocatable :: AU2(:,:)          !LU 分解の結果 U (2 次元配列)
  real(DP), allocatable :: AL1(:)            !LU 分解の結果 L (1 次元配列)
  real(DP), allocatable :: AL2(:,:)          !LU 分解の結果 L (2 次元配列)
  integer, allocatable  :: IP(:)             !部分ピボット交換の情報を格納

  integer, save :: RegXMin = 1
  integer, save :: RegXMax = 1
  integer, save :: RegYMin = 1
  integer, save :: RegYMax = 1
  integer, save :: RegZMin = 1
  integer, save :: RegZMax = 1

  !値の保存
  save beta, xyz_F1BasicZ, xyr_F2BasicZ, xyz_VPotTempBasicZ
  save N, M, NUD, NLD, NAL, NA
  save A, B, C, AU2, AL1, AL2, IP

contains

!!!--------------------------------------------------------------------!!!
  subroutine xyz_Exner_init()
    !
    !エクスナー関数を陰解法で解く際に必要となる, 係数行列の要素を決め,
    !LU 分解を行う.
    !

    !暗黙の型宣言禁止
    implicit none

    real(DP)  :: DTS ! 短い時間格子

    DTS = DelTimeShort

    RegXMin = 1
    RegXMax = nx
    RegYMin = 1
    RegYMax = ny
    RegZMin = 1
    RegZMax = nz

    !配列の割り付け
    allocate(  &
      & A(RegZMin:RegZMax),                                  &
      & B(RegZMin+1:RegZMax),                                &
      & C(RegZMin:RegZMax-1),                                &
      & xyz_F1BasicZ(imin:imax,jmin:jmax,kmin:kmax), &
      & xyr_F2BasicZ(imin:imax,jmin:jmax,kmin:kmax),&
      & xyz_VPotTempBasicZ(imin:imax,jmin:jmax,kmin:kmax) )

    !----------------------------------------------------------------
    ! 係数行列と共通して利用される配列の値を決める
    !----------------------------------------------------------------
    !仮温位の定義
    xyz_VPotTempBasicZ = xyz_PotTempBasicZ

    !係数行列の計算
    !  A, B, C を求める際, F1BasicZ と F2BasicZ は X 方向に一様なので.
    !  RegXMax, RegYMax の値で代表させることとした.

    xyz_F1BasicZ =                                                &
      &  ( xyz_VelSoundBasicZ ** 2.0d0 )                          &
      &   / (CpDry * xyz_DensBasicZ * (xyz_VPotTempBasicZ ** 2.0d0))

    xyr_F2BasicZ =                                                 &
      &  xyr_avr_xyz(                                               &
      &    CpDry * xyz_DensBasicZ * ( xyz_VPotTempBasicZ ** 2.0d0 ) &
      &   )

    A(RegZMin+1: RegZMax-1) =                                &
      & (beta ** 2.0d0)                                    &
      &    * xyz_F1BasicZ(RegXMax,RegYMax,RegZMin+1: RegZMax-1)      &
      &    * (                                               &
      &          xyr_F2BasicZ(RegXMax,RegYMax,RegZMin+1: RegZMax-1)  &
      &            / r_dz(RegZMin+1: RegZMax-1)              &
      &        + xyr_F2BasicZ(RegXMax,RegYMax,RegZMin  : RegZMax-2)  &
      &            / r_dz(RegZMin: RegZMax-2)                &
      &       )                                              &
      &    * (DTS ** 2.0d0)                                  &
      &    / z_dz(RegZMin+1: RegZMax-1)                      &
      & + 1.0d0

    A(RegZMin) =                                        &
      & (beta ** 2.0d0)                               &
      &   * xyz_F1BasicZ(RegXMax,RegYMax,RegZMin)       &
      &   * xyr_F2BasicZ(RegXMax,RegYMax,RegZMin)       &
      &     / r_dz(RegZMin)                             &
      &   * (DTS ** 2.0d0)                              &
      &   / z_dz(RegZMin)                               &
      & + 1.0d0

    A(RegZMax) =                                        &
      & (beta ** 2.0d0)                               &
      &   * xyz_F1BasicZ(RegXMax,RegYMax,RegZMax)       &
      &   * xyr_F2BasicZ(RegXMax,RegYMax,RegZMax-1)     &
      &     / r_dz(RegZMax-1)                           &
      &   * (DTS ** 2.0d0)                              &
      &   / z_dz(RegZMax)                               &
      & + 1.0d0

    B(RegZMin+1:RegZMax) =                             &
      & - (beta ** 2.0d0)                               &
      &   * xyz_F1BasicZ(RegXMax,RegYMax,RegZMin:RegZMax-1)    &
      &   * xyr_F2BasicZ(RegXMax,RegYMax,RegZMin:RegZMax-1)    &
      &     / r_dz(RegZmin:RegZMax-1)                   &
      &   * (DTS ** 2.0d0)                     &
      &   / z_dz(RegZMin:RegZMax-1)

    C(RegZMin: RegZMax-1) =                           &
      & - ( beta ** 2.0d0 )                           &
      &   * xyz_F1BasicZ(RegXMax,RegYMax,RegZMin+1:RegZMax)    &
      &   * xyr_F2BasicZ(RegXMax,RegYMax,RegZMin  :RegZMax-1)  &
      &     / r_dz(RegZmin:RegZMax-1)                   &
      &   * (DTS ** 2.0d0) &
      &   / z_dz(RegZMin+1:RegZMax)

    !デバッグ出力ファイルオープン. 情報取得.
!    call FileOpen(unit, file="dynimpfunc_log", mode='w')
!    write(unit,*) "kz, A_k"
!    write(unit,*) "------------------------------------"
!    do kz = RegZmin, RegZmax
!      write(unit,*) kz, A(kz)
!    end do
!    write(unit,*) "------------------------------------"
!    write(unit,*) "kz, B_k"
!    write(unit,*) "------------------------------------"
!    do kz = RegZmin+1, RegZmax
!      write(unit,*) kz, B(kz)
!    end do
!    write(unit,*) "------------------------------------"
!    write(unit,*) "kz, C_k"
!    write(unit,*) "------------------------------------"
!    do kz = RegZmin, RegZmax-1
!      write(unit,*) kz, C(kz)
!    end do
!    write(unit,*) "------------------------------------"
!    close(unit)

    !----------------------------------------------------------------
    ! 係数行列を LU 分解
    !----------------------------------------------------------------
    !配列の大きさを保管
    N   = RegZMax - RegZMin +1 !係数行列/改行列の次数, 整合寸法
    M   = (RegXMax - RegXMin +1)*(RegYMax - RegYMin +1)
                               !方程式の組数
    NUD = 1                    !係数行列の上三角部分の帯幅
    NLD = 1                    !係数行列の下三角部分の帯幅
    NAL = NLD                  !LU 分解の結果 L の整合寸法
    NA  = NUD + NLD + 1

    !配列の割り当て
    allocate( AL1(N), AL2(NAL, N), AU2(NA, N), IP(N) )

    !LU 分解の実行
    !  LAPACK の利用
    call ResolvLU_Lapack( )

  end subroutine xyz_Exner_init

!!!--------------------------------------------------------------------!!!
  function xyz_Exner(xyr_FzNl, &
    &                pyz_VelXNs, pyz_VelXAs, &
    &                xqz_VelYNs, xqz_VelYAs, &
    &                xyr_VelZNs, xyz_ExnerNs)
    !
    !陰解法を用いたエクスナー関数の計算.
    !

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(DP), intent(in)   :: pyz_VelXNs &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                                           !速度 u [τ]
    real(DP), intent(in)   :: pyz_VelXAs &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                                           !速度 u [τ+Δτ]
    real(DP), intent(in)   :: xqz_VelYNs &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                                           !速度 v [τ]
    real(DP), intent(in)   :: xqz_VelYAs &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                                           !速度 v [τ+Δτ]
    real(DP), intent(in)   :: xyr_VelZNs &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                                           !速度 w [τ]
    real(DP), intent(in)   :: xyr_FzNl &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                                           !Z 方向の外力項[t]
    real(DP), intent(in)   :: xyz_ExnerNs &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                                           !無次元圧力
    real(DP)               :: xyz_Exner &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                                           !無次元圧力[τ+Δτ]

    !変数定義
    real(DP)  :: D1(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)  :: D2(RegXMin:RegXMax,RegYMin:RegYMax,RegZMin:RegZMax)
    real(DP)  :: D((RegXMax-RegXMin+1)*(RegYMax-RegYMin+1),(RegZMax-RegZMin+1))
    real(DP)  :: E(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)  :: F(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)  :: xyz_DivVelNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: DTS ! 短い時間格子間隔

    DTS = DelTimeShort

    !変数の初期化
    xyz_Exner = 0.0d0

    !速度の収束を計算
    xyz_DivVelNs =  xyz_dx_pyz( pyz_VelXNs )   &
      &             + xyz_dy_xqz( xqz_VelYNs ) &
      &             + xyz_dz_xyr( xyr_VelZNs )

    !行列計算のための係数
    E = xyr_dz_xyz( DampSound * xyz_DivVelNs )                      &
      & - ( 1.0d0 - beta ) * xyr_dz_xyz( xyz_ExnerNs )              &
      & + xyr_FzNl / xyr_avr_xyz( CpDry * xyz_VPotTempBasicZ )

    F = - beta * xyz_F1BasicZ * DTS &
      & * xyz_dz_xyr(                                               &
      &    xyr_avr_xyz( xyz_DensBasicZ * xyz_VPotTempBasicZ)        &
      &    * (                                                      &
      &         xyr_VelZNs                                          &
      &       - xyr_avr_xyz(CpDry * xyz_VPotTempBasicZ) * DTS       &
      &         * (                                                 &
      &               (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerNs )    &
      &             - xyr_dz_xyz( DampSound * xyz_DivVelNs )        &
      &            )                                                &
      &       + xyr_FzNl * DTS                                      &
      &      )                                                      &
      &   )

    D1 = xyz_ExnerNs                                           &
      & - (1.0d0 - beta)                                            &
      &   * xyz_F1BasicZ * DTS                                      &
      &   * xyz_dz_xyr(                                             &
      &       xyr_avr_xyz(xyz_DensBasicZ * xyz_VPotTempBasicZ) * xyr_VelZNs   &
      &     )                                                       &
      & - (xyz_VelSoundBasicZ ** 2.0d0) *DTS        &
      &   / (CpDry * xyz_VPotTempBasicZ)                            &
      &   * xyz_dx_pyz( pyz_VelXAs )                                &
      & - (xyz_VelSoundBasicZ ** 2.0d0) *DTS        &
      &   / (CpDry * xyz_VPotTempBasicZ)                            &
      &   * xyz_dy_xqz( xqz_VelYAs )                                &
      & + F

    D1(:,:,RegZMin) = D1(:,:,RegZMin)                       &
      & - beta * xyz_F1BasicZ(:,:,RegZMin) * (DTS ** 2.0d0)         &
      &   * xyr_F2BasicZ(:,:, RegZMin-1) * E(:,:,RegZMin-1)         &
      &   / z_dz(RegZMin)

    D1(:,:,RegZMax) = D1(:,:,RegZMax)                       &
      & + beta * xyz_F1BasicZ(:,:,RegZMax) * (DTS ** 2.0d0)         &
      &   * xyr_F2BasicZ(:,:,RegZMax) * E(:,:,RegZMax)              &
      &   / z_dz(RegZMax)

    D2 = D1(RegXMin:RegXMax,RegYMin:RegYMax,RegZMin:RegZMax)

    do kz = RegZMin, RegZMax
      do jy = RegYMin, RegYMax
        do ix = RegXMin, RegXMax
	  D(ix + (RegXMax-RegXMin+1) * (jy - 1), kz) =  D2(ix,jy,kz)
        end do
      end do
    end do

    !-----------------------------------------------------------
    !連立一次方程式の解を求める
    !------------------------------------------------------------

    !解の計算
    !  LAPACK 利用
    call LinSolv_Lapack( D )

    !戻り値を出力

    do kz = RegZMin, RegZMax
      do jy = RegYMin, RegYMax
        do ix = RegXMin, RegXMax
          xyz_Exner(ix,jy,kz) = D(ix + (RegXMax-RegXMin+1) * (jy - 1 ), kz)
        end do
      end do
    end do

  end function xyz_Exner

!!!--------------------------------------------------------------------!!!
  subroutine ResolvLU_Lapack(  )
    !
    !実 3 項行列の LU 分解(倍精度). LAPACK 利用
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer    :: INFO  !解のコンディションチェック

    !変数の初期化
    INFO = 0

    !解行列の計算. LAPACK を使用.
    call DGTTRF(N, C, A, B, AL1, IP, INFO)

!    !解のコンディションをチェック.
!    if (INFO /= 0) then
!      call MessageNotify("Error", "lapack_linear", "INFO is not 0")
!      stop
!    end if

  end subroutine ResolvLU_Lapack

!!!--------------------------------------------------------------------!!!
  subroutine LinSolv_Lapack( X )
    !
    !LU 分解された実 3 項行列の連立 1 次方程式(倍精度). LAPACK 利用
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(inout) :: X(M, N)     !定数/解行列
    real(DP)                :: TX(N, M)    !解行列を転置したもの
    integer                :: NRHS         !
    integer                :: INFO
    integer                :: LDB
    character(1),parameter :: TRANS = 'N'

    !変数の初期化
    NRHS = M
    INFO = 0
    LDB  = N
    TX = transpose( X )

    !解行列の計算. LAPACK を使用.
    call DGTTRS(TRANS, N, NRHS, C, A, B, AL1, IP, TX, LDB, INFO)

    !解の出力
    X = transpose( TX )

    !解のコンディションをチェック.
!    if (INFO /= 0) then
!      call MessageNotify("Error", "lapack_linear", "INFO is not 0")
!      stop
!    end if

  end subroutine LinSolv_Lapack

!!!--------------------------------------------------------------------!!!
  function xyr_GradPi(xyz_ExnerA, xyz_ExnerN, pyz_VelX, xqz_VelY, xyr_VelZ)
    !
    ! x 方向に半格子ずれた点での圧力傾度力項の計算.
    ! クランク・ニコルソン法を用いるために, 時刻 τ と τ+Δτでの
    ! エクスナー関数の値を引数として取る.
    ! 音波減衰項を含めた形式で定式化してあることに注意.
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(in)  :: xyz_ExnerA &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !エクスナー関数の擾乱
    real(DP), intent(in)  :: xyz_ExnerN &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !エクスナー関数の擾乱
    real(DP), intent(in)  :: pyz_VelX &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !水平速度
    real(DP), intent(in)  :: xqz_VelY &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !水平速度
    real(DP), intent(in)  :: xyr_VelZ &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !鉛直速度
    real(DP)              :: xyr_GradPi &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !圧力傾度力
    real(DP)              :: xyz_DivVel &
      &                     (imin:imax,jmin:jmax,kmin:kmax)
                                               !速度の収束

    !速度の収束
    xyz_DivVel =  xyz_dx_pyz( pyz_VelX )   &
      &           + xyz_dy_xqz( xqz_VelY ) &
      &           + xyz_dz_xyr( xyr_VelZ )

    !速度 w の圧力勾配
!    xyr_GradPi = 0.0d0
    xyr_GradPi =   &
      & xyr_avr_xyz(CpDry * xyz_PotTempBasicZ )   &
      &   * (                                                      &
      &       beta * xyr_dz_xyz( xyz_ExnerA )                         &
      &       + (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerN )             &
      &       - xyr_dz_xyz( DampSound * xyz_DivVel )                  &
      &     )

  end function xyr_GradPi

end module DynImpFunc