!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module Linlib
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: linlib.f90.g95,v 1.2.2.2 2005/11/14 01:55:12 kitamo Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!線形計算ライブラリ SSL II と LAPACK の g95 用ラッパー
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
                                                                 !=end

module linlib
  implicit none
  
  private
                                                                 !=begin
  !== Public Interface
  public linlib_init       !初期化ルーチン
  public linsolv           !実 3 項行列の連立 1 次方程式(倍精度)

  integer       :: LN = 10 !配列のサイズ. 値はダミー
  character(1)  :: LLib    !用いるライブラリのフラグ
  real(8), allocatable  ::  A(:)              ! 係数行列
  real(8), allocatable  ::  B(:)              ! 係数行列
  real(8), allocatable  ::  C(:)              ! 係数行列
  real(8), allocatable  ::  B2(:)             ! 係数行列(LU 分解後)
  real(8), allocatable  ::  ss_F1BasicZ(:,:)  ! 係数行列の計算に用いる配列
  real(8), allocatable  ::  sf_F2BasicZ(:,:)  ! 係数行列の計算に用いる配列
  real(8), allocatable  ::  IPIV(:)           ! 

  save LLib, LN, A, B, C, B2, ss_F1BasicZ, sf_F2BasicZ, IPIV
                                                                 !=end    
contains
                                                                 !=begin
  !== Procedure Interface 
  !
  !=== Initialize module and acquire NAMELIST
  !
  !利用する線形計算パッケージを NAMELIST から設定し, 
  !配列の大きさを決定する
  !
  subroutine linlib_init(cfgfile, N)

    !=== Dependency
    use dc_trace, only : BeginSub, EndSub
    use dc_message, only: MessageNotify
    use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, &
         &              RegXMax, RegXMin, RegZMin, RegZMax, DelZ
    use timeset,  only: DelTimeShort
    use basicset, only: ss_VelSoundBasicZ, ss_CpBasicZ, &
         &              ss_DensBasicZ, ss_PotTempBasicZ
    use arareset, only: beta
    use ssl2_linear, only: ssl2_ltx_init
    use lapack_linear, only: lapack_dgtsv_init, lapack_dgttrf
    use average, only : sf_avr_ss

    !=== Input
    integer, intent(in)      :: N       !配列サイズ
    character(*), intent(in) :: cfgfile !設定ファイル (NAMELIST)
    
    !=== NAMELIST
    NAMELIST /linlib/ LLib
                                                                 !=end
    open (10, FILE=cfgfile)
    read(10, NML=linlib)
    close(10)

    call BeginSub("linlib_init", &
&                 fmt="%c",       &
&                 c1="Initialize linear algebra library.")

    !--- 配列サイズを保管
    LN = N

    !--- 配列の割り付け
	allocate(  &
	& A(RegZMin+1:RegZMax),   &
	& B(RegZMin+2:RegZMax),   &  
	& C(RegZMin+1:RegZMax-1), &  
        & B2(RegZMin+3:RegZMax),  &
	& ss_F1BasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &  
	& sf_F2BasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
	& IPIV(RegZMin+1:RegZMax) )

    !--- 用いるライブラリの初期化
    if (LLib == 's') then 
       call MessageNotify("Message", "linlib.f90", &
&                         "Can't Call SSL II in the calculation with g95")

!       call ssl2_ltx_init(LN)

    elseif (LLib == 'l') then 
       call MessageNotify("Message", "linlib.f90", &
&                         "Call LAPACK in the calculation")

       call lapack_dgtsv_init(LN)

    else
      call MessageNotify("Error", "linlib.f90", &
&                        "Unknown LLib ", c1=LLib)
    end if

    !--- 係数行列の計算

    ss_F1BasicZ = ( ss_VelSoundBasicZ ** 2.0d0 ) &
       &        / (ss_CpBasicZ * ss_DensBasicZ *(ss_PotTempBasicZ ** 2.0d0))

    sf_F2BasicZ = sf_avr_ss( &
         & ss_CpBasicZ * ss_DensBasicZ * (ss_PotTempBasicZ ** 2.0d0) )

    A(RegZMin+2: RegZMax-1) = &
         & 1.0d0 + (beta ** 2.0d0) &
         &   * ss_F1BasicZ(RegXMax, RegZMin+2: RegZMax-1)  &  
         &   * ( sf_F2BasicZ(RegXMax, RegZMin+2: RegZMax-1)  &
         &       +  sf_F2BasicZ(RegXMax, RegZMin+1: RegZMax-2) ) &
         &   * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0)
  
    A(RegZMin+1) = 1.0d0 &
         & + (beta ** 2.0d0) &
         &   * ss_F1BasicZ(RegXMax, RegZMin+1)  &
         &   * sf_F2BasicZ(RegXMax, RegZMin+1) &
         &   * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0)

    A(RegZMax) = 1.0d0 &
         & + (beta ** 2.0d0) &
         &   * ss_F1BasicZ(RegXMax, RegZMax)  &
         &   * sf_F2BasicZ(RegXMax, RegZMax-1) &
         &   * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0)
  
    B(RegZMin+2: RegZMax) = &
         & - ( beta ** 2.0d0 ) &
         &   * ss_F1BasicZ(RegXMax, RegZMin+1: RegZMax-1) &
         &   * sf_F2BasicZ(RegXMax, RegZMin+1: RegZMax-1) &
         &   * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0 )
  
    C(RegZMin+1: RegZMax-1) = &
         & - ( beta ** 2.0d0 ) &
         &   * ss_F1BasicZ(RegXMax, RegZMin+2: RegZMax) &
         &   * sf_F2BasicZ(RegXMax, RegZMin+1: RegZMax-1) &
         &   * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0 )

    !--- 係数行列を LU 分解 (LAPACK)
    if (LLib == 'l') then
       call lapack_dgttrf(A, B, C, B2, IPIV)
    endif

    call EndSub("linlib_init")

  end subroutine linlib_init
  
                                                                 !=begin
  !== Procedure Interface 
  !
  !=== Initialize module
  !
  !利用する線形計算パッケージに配列を渡す
  !
  subroutine LinSolv(D)

    !=== Dependency
    use dc_trace, only : BeginSub, EndSub
    use dc_message, only: MessageNotify
!    use ssl2_linear, only: ssl2_ltx
    use lapack_linear, only: lapack_dgttrs

    !=== In/Out
    real(8), intent(inout) :: D(LN)
                                                                 !=end
    !=== Work
    real(8)                :: AA(LN)
    real(8)                :: BB(LN-1)
    real(8)                :: CC(LN-1) 


    call BeginSub("linsolv", &
&                 fmt="%c",       &
&                 c1="call linear algebra library.")

    !--- 連立 1 次方程式を解く
    if (LLib == 's') then 
!       call ssl2_ltx(A, B, C, D)
       call MessageNotify("Error", "linlib", &
&                        "unknown LLib for g95")

    elseif (LLib == 'l') then 
       !=== 属性が inout なので, 内部変数を使うことに.
       AA = A; BB = B; CC = C
       call lapack_dgttrs(AA, BB, CC, D, B2, IPIV)
       
    else
      call MessageNotify("Error", "linlib", &
&                        "unknown LLib")
    end if

    call EndSub("linsolv")

  end subroutine LinSolv
  
end module linlib
