!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2006. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine Dynamic_exner
!
! Developer:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: dynamic_exner.f90,v 1.1 2011-04-26 04:08:10 sugiyama 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 Dynamic_Exner
  !
  !陰解法を用いた力学過程の各項の計算モジュール.
  !エクスナー関数を陰解法で解く際に必要となる, 係数行列の要素の決定, 
  !
  !本モジュールでは科学計算ライブラリとして, 以下をサポートする. 
  !  * 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 constants,only: CpDry,            &  !乾燥成分の比熱
    &                 Nu                   !音波の減衰係数
  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 dynamic_exner_implicit_init    !初期化ルーチン
  public calc_dynamic_exner_implicit    !エクスナー関数の計算
  
  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
  real(DP), save:: DTS = 1.0

  !値の保存
  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 dynamic_exner_implicit_init()
    !
    !エクスナー関数を陰解法で解く際に必要となる, 係数行列の要素を決め, 
    !LU 分解を行う. 
    !

    !暗黙の型宣言禁止
    implicit none

    ! 整合性のために名前の付け替え
    RegXMin = 1
    RegXMax = nx
    RegYMin = 1
    RegYMax = ny
    RegZMin = 1
    RegZMax = nz
    DTS = DelTimeShort

    !配列の割り付け
    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 dynamic_exner_implicit_init
  

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

    !暗黙の型宣言禁止
    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), intent(out)  :: xyz_ExnerAs &
      &                     (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)

    !変数の初期化
    xyz_ExnerAs = 0.0d0
    
    !速度の収束を計算
    xyz_DivVelNs =  xyz_dx_pyz( pyz_VelXNs )   &
      &             + xyz_dy_xqz( xqz_VelYNs ) &
      &             + xyz_dz_xyr( xyr_VelZNs )
  
    !行列計算のための係数
    E = xyr_dz_xyz( Nu * 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( Nu * 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_ExnerAs(ix,jy,kz) = D(ix + (RegXMax-RegXMin+1) * (jy - 1 ), kz)
        end do 
      end do
    end do	
    
  end subroutine calc_dynamic_Exner_implicit

!!!--------------------------------------------------------------------!!!
  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

  
end module Dynamic_Exner
