!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
! physics_implicit.f90 
!
! History
!   2005/09/21 Yamada Yukiko     create
!

module physics_implicit_mod

  implicit none

  private
  public :: physics_implicit_init
  public :: physics_implicit_integrate

contains

  subroutine physics_implicit_init( &
         & xyr_VelLonFlux       , & ! (out) 速度経度成分フラックス
         & xyr_VelLatFlux       , & ! (out) 速度緯度成分フラックス
         & xyr_TempFlux         , & ! (out) 温度フラックス
         & xyr_QvapFlux         , & ! (out) 比湿フラックス
         & xyzo_VelMatrix       , & ! (out) 速度陰解行列
         & xyzo_TempMatrix      , & ! (out) 温度陰解行列
         & xyzo_QvapMatrix      , & ! (out) 比湿陰解行列
         & xyr_Press            , & ! (in) 気圧 (半整数)
         & DelTimePhy           , & ! (in) ２Δt
         & xy_SurfHeatCapacity  , & ! (in) 地表熱容量
         & xy_SurfCondition       ) ! (in) 地表状態

    !==== Dependency
    use type_mod,    only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod, only: im, jm, km
    use constants_mod, only: Cp, Grav
    use dc_trace,    only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    !==== Output
    !
    real(DBKIND), intent(out) :: &
         & xyr_VelLonFlux(im*jm,km+1)        , & ! (out) 速度経度成分フラックス
         & xyr_VelLatFlux(im*jm,km+1)        , & ! (out) 速度緯度成分フラックス
         & xyr_TempFlux(im*jm,km+1)          , & ! (out) 温度フラックス
         & xyr_QvapFlux(im*jm,km+1)          , & ! (out) 比湿フラックス
         & xyzo_VelMatrix(im*jm,km,-1:1)     , & ! (out) 速度陰解行列
         & xyzo_TempMatrix(im*jm,0:km,-1:1)  , & ! (out) 温度陰解行列
         & xyzo_QvapMatrix(im*jm,km,-1:1)        ! (out) 比湿陰解行列

    !==== Input
    !
    real(DBKIND), intent(in) :: &
         & xyr_Press(im*jm,km+1)             , & ! (in) 気圧 (半整数)
         & DelTimePhy                        , & ! (in) Δt
         & xy_SurfHeatCapacity(im*jm)            ! (in) 地表熱容量

    integer(INTKIND), intent(in) :: &
         & xy_SurfCondition(im*jm)               ! (in) 地表状態

    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "physics_implicit_init"

    ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    integer(INTKIND)    :: ij, k

    continue

    !----------------------------------------------------------------
    !   開始処理
    !----------------------------------------------------------------
    call BeginSub(subname)

    !----------------------------------------------------------------
    !   陰解行列初期化
    !----------------------------------------------------------------

    ! ---- 1. 質量, 熱容量の項 ----
    xyzo_VelMatrix  = 0.0 
    xyzo_TempMatrix = 0.0 
    xyzo_QvapMatrix = 0.0 

    do k = 1, km 
       xyzo_VelMatrix(:,k,0)  = ( xyr_Press(:,k) - xyr_Press(:,k+1) ) &
            &                   / Grav / DelTimePhy  
       xyzo_TempMatrix(:,k,0) = xyzo_VelMatrix(:,k,0) * Cp
       xyzo_QvapMatrix(:,k,0) = xyzo_VelMatrix(:,k,0) * Cp
    end do

    do ij = 1, im*jm
       if ( xy_SurfCondition(ij) .GE. 1 ) then 
          xyzo_TempMatrix(ij,0,0) = xy_SurfHeatCapacity(ij) / DelTimePhy 
       else
          xyzo_TempMatrix(ij,0,0) = 1.
       end if
    end do

    ! ---- 2. フラックスをリセット ----
    xyr_VelLonFlux = 0.0 
    xyr_VelLatFlux = 0.0 
    xyr_TempFlux   = 0.0 
    xyr_QvapFlux   = 0.0 

    !----------------------------------------------------------------
    !   終了処理
    !----------------------------------------------------------------
    call EndSub(subname)

  end subroutine physics_implicit_init


  subroutine physics_implicit_integrate( &
       & xyz_DVerdiffVelLonDt    , & ! (out) 経度成分 鉛直拡散加速度
       & xyz_DVerdiffVelLatDt    , & ! (out) 緯度成分 鉛直拡散加速度
       & xyz_DVerdiffTempDt      , & ! (out) 鉛直拡散加熱率
       & xyz_DVerdiffSurfTempDt  , & ! (out) 地表面 鉛直拡散加熱率
       & xyz_DVerdiffQvapDt      , & ! (out) 鉛直拡散加湿率
       & xyr_VelLonFlux          , & ! (in) 速度経度成分フラックス
       & xyr_VelLatFlux          , & ! (in) 速度緯度成分フラックス
       & xyr_TempFlux            , & ! (in) 温度フラックス
       & xyr_SurfRadSFlux        , & ! (in) 日射フラックス
       & xyr_SurfRadLFlux        , & ! (in) 長波フラックス
       & xy_GroundTempFlux       , & ! (in) 地中熱フラックス
       & xyr_QvapFlux            , & ! (in) 比湿フラックス
       & xyzo_VelMatrix          , & ! (in) 速度陰解行列
       & xyzo_TempMatrix         , & ! (in) 温度陰解行列
       & xyzo_QvapMatrix         , & ! (in) 比湿陰解行列
       & xy_SurfVelMatrix        , & ! (in) 速度陰解行列: 地表
       & xyoo_SurfTempMatrix     , & ! (in) 温度陰解行列: 地表
       & xyoo_SurfQvapMatrix     , & ! (in) 比湿陰解行列: 地表
       & xyo_SurfRadLMatrix      , & ! (in) Ｔ陰解行列：地表
       & DelTimePhy              , & ! (in) ２Δt
       & xy_SurfCondition          ) ! (in) 地表状態

    !==== Dependency
    use type_mod,    only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod, only: im, jm, km
    use constants_mod, only: EL, Cp 
    use dc_trace,    only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    !==== Output
    !
    real(DBKIND), intent(out) :: &
         & xyz_DVerdiffVelLonDt(im*jm,km)    , & ! (out) 経度成分 鉛直拡散加速度
         & xyz_DVerdiffVelLatDt(im*jm,km)    , & ! (out) 緯度成分 鉛直拡散加速度
         & xyz_DVerdiffTempDt(im*jm,km)      , & ! (out) 鉛直拡散加熱率
         & xyz_DVerdiffSurfTempDt(im*jm)     , & ! (out) 地表面 鉛直拡散加熱率
         & xyz_DVerdiffQvapDt(im*jm,km)          ! (out) 鉛直拡散加湿率
    
    !==== Input
    !
    real(DBKIND), intent(in) :: &
         & xyr_VelLonFlux(im*jm,km+1)          , & ! (in) 速度経度成分フラックス
         & xyr_VelLatFlux(im*jm,km+1)          , & ! (in) 速度緯度成分フラックス
         & xyr_TempFlux(im*jm,km+1)            , & ! (in) 温度フラックス
         & xyr_SurfRadSFlux(im*jm)             , & ! (in) 日射フラックス
         & xyr_SurfRadLFlux(im*jm)             , & ! (in) 長波フラックス
         & xy_GroundTempFlux(im*jm)            , & ! (in) 地中熱フラックス
         & xyr_QvapFlux(im*jm,km+1)            , & ! (in) 比湿フラックス
         & xyzo_VelMatrix(im*jm,km,-1:1)       , & ! (in) 速度陰解行列
         & xyzo_TempMatrix(im*jm,0:km,-1:1)    , & ! (in) 温度陰解行列
         & xyzo_QvapMatrix(im*jm,km,-1:1)      , & ! (in) 比湿陰解行列
         & xy_SurfVelMatrix(im*jm)             , & ! (in) 速度陰解行列: 地表
         & xyoo_SurfTempMatrix(im*jm,0:1,-1:1) , & ! (in) 温度陰解行列: 地表
         & xyoo_SurfQvapMatrix(im*jm,0:1,-1:1) , & ! (in) 比湿陰解行列: 地表
         & xyo_SurfRadLMatrix(im*jm,-1:1)      , & ! (in) Ｔ陰解行列：地表
         & DelTimePhy                              ! (in) ２Δt

    integer(INTKIND)    :: &
         & xy_SurfCondition(im*jm)               ! (in) 地表状態


    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "physics_implicit_integrate"

    ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    integer(INTKIND)    :: ij, k, l

    real(DBKIND) :: & 
         & xyz_DelTempQvap(im*jm,-km:km)            , & ! Ｔｑ時間変化
         & xyzo_TempQvapLUMatrix(im*jm,-km:km,-1:1) , & ! ＬＵ行列
         & xyzo_VelLUMatrix(im*jm,km,-1:1)              ! ＬＵ行列

    continue

    !----------------------------------------------------------------
    !   開始処理
    !----------------------------------------------------------------
    call BeginSub(subname)

    !----------------------------------------------------------------
    !   陰解行列計算
    !----------------------------------------------------------------

    ! ---- 1. 速度 (Vlon, Vlat) の解 ----

    xyzo_VelLUMatrix  = xyzo_VelMatrix
    xyzo_VelLUMatrix(:,1,0)  = xyzo_VelLUMatrix(:,1,0) + xy_SurfVelMatrix(:)

    call lu_decomposition_tridiagonal( &
         & xyzo_VelLUMatrix, im*jm, km )

    do k = 1, km
       xyz_DVerdiffVelLonDt(:,k) = xyr_VelLonFlux(:,k) - xyr_VelLonFlux(:,k+1)
       xyz_DVerdiffVelLatDt(:,k) = xyr_VelLatFlux(:,k) - xyr_VelLatFlux(:,k+1)
    end do

    call lu_solve_tridiagonal( &
         & xyz_DVerdiffVelLonDt , & 
         & xyzo_VelLUMatrix     , & 
         & 1, im*jm, km )

    call lu_solve_tridiagonal( &
         & xyz_DVerdiffVelLatDt , & 
         & xyzo_VelLUMatrix     , & 
         & 1, im*jm, km )

    ! ---- 2. 温度と比湿の解 ----

    do l = -1, 1

       do k = 1, km
          xyzo_TempQvapLUMatrix(:,k,l)   = xyzo_TempMatrix(:,k,l)
          xyzo_TempQvapLUMatrix(:,-k,-l) = xyzo_QvapMatrix(:,k,l)
       end do
       
       xyzo_TempQvapLUMatrix(:,1,l)   = xyzo_TempMatrix(:,1,l) &
            &                          + xyoo_SurfTempMatrix(:,1,l)
       xyzo_TempQvapLUMatrix(:,-1,-l) = xyzo_QvapMatrix(:,1,l) & 
            &                          + xyoo_SurfQvapMatrix(:,1,l)

    end do

    xyzo_TempQvapLUMatrix(:,0,0) = xyzo_TempMatrix(:,0,0) &
         & + xyoo_SurfTempMatrix(:,0,0) + xyoo_SurfQvapMatrix(:,0,0) &
         & + xyo_SurfRadLMatrix(:,0)

    xyzo_TempQvapLUMatrix(:,0,1) = &
         & + xyoo_SurfTempMatrix(:,0,1) + xyo_SurfRadLMatrix(:,1)

    xyzo_TempQvapLUMatrix(:,0,-1) =  xyoo_SurfQvapMatrix(:,0,1) 

    call lu_decomposition_tridiagonal( &
         & xyzo_TempQvapLUMatrix, im*jm, 2*km+1 )

    do k = 1, km
       xyz_DelTempQvap(:,k)  = xyr_TempFlux(:,k) - xyr_TempFlux(:,k+1)
       xyz_DelTempQvap(:,-k) = xyr_QvapFlux(:,k) - xyr_QvapFlux(:,k+1)
    end do

    xyz_DelTempQvap(:,0) = - xyr_SurfRadSFlux  - xyr_SurfRadLFlux  &
         &                 - xyr_TempFlux(:,1) - xyr_QvapFlux(:,1) &
         &                 + xy_GroundTempFlux


    call lu_solve_tridiagonal( &
         & xyz_DelTempQvap       , & 
         & xyzo_TempQvapLUMatrix , & 
         & 1, im*jm, 2*km+1 )

    ! ---- 2. 時間変化率 ----
    do k = 1, km
       xyz_DVerdiffVelLonDt(:,k) = xyz_DVerdiffVelLonDt(:,k) / DeltimePhy 
       xyz_DVerdiffVelLatDt(:,k) = xyz_DVerdiffVelLatDt(:,k) / DeltimePhy 
       xyz_DVerdiffTempDt(:,k)   = xyz_DelTempQvap(:,k) / DeltimePhy 
       xyz_DVerdiffQvapDt(:,k)   = xyz_DelTempQvap(:,-k) / DeltimePhy / EL * Cp
    end do
    
    do ij = 1, im*jm
       if ( xy_SurfCondition(ij) .GE. 1 ) then 
          xyz_DVerdiffSurfTempDt(ij) = xyz_DelTempQvap(ij,0) / DeltimePhy 
       else
          xyz_DVerdiffSurfTempDt(ij) = 0.
       end if
    end do

    !----------------------------------------------------------------
    !   終了処理
    !----------------------------------------------------------------
    call EndSub(subname)

  end subroutine physics_implicit_integrate


  subroutine lu_decomposition_tridiagonal( &
       & jno_LUMatrix, JDimMax, NDimMax )

    !==== Dependency
    use type_mod,    only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use dc_trace,    only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    !==== Parameter
    !
    integer(INTKIND), intent(in)    :: &
         JDimMax, NDimMax 

    !==== In/Out
    !
    real(DBKIND), intent(inout) :: &
         jno_LUMatrix(JDimMax, NDimMax, -1:1) ! 入力／ＬＵ行列

    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "lu_decomposition_tridiagonal"
    integer(INTKIND)    :: j, n ! do ループ用作業変数

    continue

    !----------------------------------------------------------------
    !   開始処理
    !----------------------------------------------------------------
    call BeginSub(subname)

    !----------------------------------------------------------------
    !   行列のＬＵ分解 [ ３重対角行列 ]
    !----------------------------------------------------------------

    do j = 1, JDimMax
       jno_LUMatrix(j,1,1) = jno_LUMatrix(j,1,1) / jno_LUMatrix(j,1,0)
    end do

    do n = 2, NDimMax-1
       do j = 1, JDimMax

       jno_LUMatrix(j,n,0) = jno_LUMatrix(j,n,0) &
            &               - jno_LUMatrix(j,n,-1) * jno_LUMatrix(j,n-1,1) 

       jno_LUMatrix(j,n,1) = jno_LUMatrix(j,n,1) /jno_LUMatrix(j,n,0) 

       end do
    end do

    do j = 1, JDimMax
       jno_LUMatrix(j,NDimMax,0) = jno_LUMatrix(j,NDimMax,0) &
            &    - jno_LUMatrix(j,NDimMax,-1) * jno_LUMatrix(j,NDimMax-1,1) 
    end do

    !----------------------------------------------------------------
    !   終了処理
    !----------------------------------------------------------------
    call EndSub(subname)

  end subroutine lu_decomposition_tridiagonal


  subroutine lu_solve_tridiagonal( &
         & ijn_Vector       , & 
         & jno_LUMatrix     , & 
         & IDimMax, JDimMax, NDimMax )

    !==== Dependency
    use type_mod,    only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use dc_trace,    only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    !==== Parameter
    !
    integer(INTKIND), intent(in)    :: &
         IDimMax, JDimMax, NDimMax 

    !==== In/Out
    !
    real(DBKIND), intent(inout) :: &
         & ijn_Vector(IDimMax, JDimMax, NDimMax) ! 右辺ベクトル／解

    !==== Input
    !
    real(DBKIND), intent(in) :: &
         jno_LUMatrix(JDimMax, NDimMax, -1:1) ! ＬＵ行列

    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "lu_solve_tridiagonal"
    integer(INTKIND)    :: i, j, n    ! do ループ用作業変数

    continue

    !----------------------------------------------------------------
    !   開始処理
    !----------------------------------------------------------------
    call BeginSub(subname)

    !----------------------------------------------------------------
    !   ＬＵ分解による解の計算 [ ３重対角行列 ]
    !----------------------------------------------------------------

    ! ---- 1. 前進代入 ----    

    do i = 1, IDimMax
       do j = 1, JDimMax
          ijn_Vector(i,j,1) = ijn_Vector(i,j,1) / jno_LUMatrix(j,1,0)
       end do
    end do

    do n = 2, NDimMax
       do i = 1, IDimMax
          do j = 1, JDimMax
             ijn_Vector(i,j,n) = ( ijn_Vector(i,j,n) &
                  & - ijn_Vector(i,j,n-1) * jno_LUMatrix(j,n,-1) ) &
                  &               / jno_LUMatrix(j,n,0)
          end do
       end do       
    end do

    ! ---- 2. 後退代入 ----    

    do n = NDimMax-1, 1, -1
       do i = 1, IDimMax
          do j = 1, JDimMax
             ijn_Vector(i,j,n) = ijn_Vector(i,j,n) &
                  & - ijn_Vector(i,j,n+1) * jno_LUMatrix(j,n,1)
          end do
       end do       
    end do

    !----------------------------------------------------------------
    !   終了処理
    !----------------------------------------------------------------
    call EndSub(subname)

  end subroutine lu_solve_tridiagonal


end module physics_implicit_mod
















