!= Module DynamicsHEVI
!
! Authors::   杉山耕一朗(SUGIYAMA Ko-ichiro), ODAKA Masatsugu 
! Version::   $Id: dynamicshevi.f90,v 1.3 2011-06-17 19:06:17 sugiyama Exp $ 
! Tag Name::  $Name: arare5-20110623 $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
!モデルの力学過程を計算するために必要となる関数群を束ねたモジュール
!具体的には以下の項を計算するための関数を格納する.  
!  * 移流項
!  * 浮力項
!  * 気圧傾度力項
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * エクスナー関数の空間方向の離散化において, 2 次精度の離散化を陽に利
!    用しているため, 気圧傾度力項の計算プログラムにおいて
!    differentiate_center4 モジュールを指定することはできないので注意.
!
!== Future Plans
!

module DynamicsHEVI
  !
  !陽解法を用いた力学過程の各項の計算モジュール. 
  !具体的には以下の項を計算するための関数を格納する.  
  !  * 移流項
  !  * 浮力項
  !  * 気圧傾度力項
  !

  !モジュール読み込み
  use dc_types,   only : DP, STRING
  use dc_iounit,  only : FileOpen
  use dc_message, only : MessageNotify
  use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut 

  use mpi_wrapper,only: myrank
  use gridset, only: &
    &                 imin,            &! x 方向の配列の下限
    &                 imax,            &! x 方向の配列の上限
    &                 jmin,            &! y 方向の配列の下限
    &                 jmax,            &! y 方向の配列の上限
    &                 kmin,            &! z 方向の配列の下限
    &                 kmax,            &! z 方向の配列の上限
    &                 nx,              &! x 方向の物理領域の上限
    &                 ny,              &! x 方向の物理領域の上限
    &                 nz,              &! y 方向の物理領域の上限
    &                 ncmax                ! 物質数
  use constants,only: CpDry,           &! 乾燥成分の比熱
    &                 MolWtDry,        &! 乾燥成分の分子量
    &                 Grav              ! 重力加速度
  use composition, only: SpcWetSymbol,    &!
    &                 GasNum,          &! 
    &                 IdxG,            &!
    &                 MolWtWet          ! 湿潤成分の分子量
  use timeset, only:  DelTimeShort, DelTimeLong, TimeN
  use axesset, only:  x_dx, y_dy, z_dz, r_dz,                &! 格子間隔
    &                 xyz_avr_pyz, xyr_avr_pyr, xqz_avr_pqz, &
    &                 pyz_avr_xyz, pyr_avr_xyr, pqz_avr_xqz, &
    &                 xyz_avr_xqz, pyz_avr_pqz, xyr_avr_xqr, &
    &                 xqz_avr_xyz, pqz_avr_pyz, xqr_avr_xyr, &
    &                 xyz_avr_xyr, pyz_avr_pyr, xqz_avr_xqr, &
    &                 xyr_avr_xyz, pyr_avr_pyz, xqr_avr_xqz 
  use basicset, only: xyz_VelSoundBZ, &!基本場の音速 
    &                 xyz_DensBZ,     &!基本場の密度
    &                 xyz_PTempBZ,    &!基本場の温位
    &                 xyz_VPTempBZ,    &!基本場の温位
    &                 xyzf_QMixBZ,    &
    &                 xyr_QMixBZPerMolWt,  &
    &                 xyr_QMixBZ   
  use xyz_deriv_module
  use xyz_deriv_c4_module, only: pyz_c4dx_xyz, xqz_c4dy_xyz, xyr_c4dz_xyz, &
    &                            xyz_c4dx_pyz, pqz_c4dy_pyz, pyr_c4dz_pyz, &
    &                            pqz_c4dx_xqz, xyz_c4dy_xqz, xqr_c4dz_xqz, &
    &                            pyr_c4dx_xyr, xqr_c4dy_xyr, xyz_c4dz_xyr
  use setmargin,only: SetMargin_xyzf, SetMargin_xyz, &
    &                 SetMargin_pyz, SetMargin_xqz, SetMargin_xyr
  use fillnegative, only: FillNegativeQMix
  use namelist_util, only: namelist_filename

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private
  
  real(DP), save         :: beta  = 5.0d-1   !クランクニコルソン法なら 0.5
                                             !完全陰解法なら 1
  real(DP), allocatable, save :: xyz_F1BZ(:,:,:)  !係数行列の計算に用いる配列
  real(DP), allocatable, save :: xyr_F2BZ(:,:,:)  !係数行列の計算に用いる配列

  integer, save         :: N = 10            !係数行列/改行列の次数, 整合寸法
  integer, save         :: M = 10            !方程式の組数
  integer, save         :: NUD = 1           !係数行列の上三角部分の帯幅
  integer, save         :: NLD = 1           !係数行列の下三角部分の帯幅
  integer, save         :: NAL = 1           !LU 分解の結果 L の整合寸法
  integer, save         :: NA = 3            !NUD + NLD + 1

  real(DP), allocatable, save :: A(:)              !係数行列の対角成分
  real(DP), allocatable, save :: B(:)              !係数行列の上三角部分
  real(DP), allocatable, save :: C(:)              !係数行列の下三角部分
!  real(DP), allocatable, save :: AU2(:,:)          !LU 分解の結果 U (2 次元配列)
  real(DP), allocatable, save :: AL1(:)            !LU 分解の結果 L (1 次元配列)
!  real(DP), allocatable, save :: AL2(:,:)          !LU 分解の結果 L (2 次元配列)
  integer, allocatable, save  :: 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 :: Nu  = 0.0d0           !音波減衰項の減衰係数
  real(DP), save :: NuHh   = 0.0d0        !数値粘性の係数 (水平方向)
  real(DP), save :: NuVh   = 0.0d0        !数値粘性の係数 (鉛直方向)
  real(DP), save :: NuHm  = 0.0d0         !数値粘性の係数 (水平方向)
  real(DP), save :: NuVm  = 0.0d0         !数値粘性の係数 (鉛直方向)

  !public 
  public Dynamics_Init
  public Dynamics_Long_forcing
  public Dynamics_Short_forcing


contains

  subroutine Dynamics_Init

    !暗黙の型宣言禁止
    implicit none
    
    real(DP)  :: DelXMin, DelYMin, DelZMin
    real(DP)  :: AlphaSound = 1.0e-7  !音波減衰項の係数
    real(DP)  :: AlphaNDiff = 1.0d-4
    real(DP)  :: NDiffRatio = 1.0d0  !速度に対する粘性を上げる場合は数字を 1 以上にする. 
    integer   :: unit            !装置番号
    integer   :: l

    NAMELIST /Dynamics_nml/ AlphaSound, AlphaNDiff, NDiffRatio

    !ファイルオープン. 情報取得. 
    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=dynamics_nml)
    close(unit)
    
    !-------------------------------------------------------------------
    ! 音波減衰項の減衰率   Min(DelX, DelZ) ** 2.0 に比例
    !
    DelXMin = minval(x_dx)
    DelYMin = minval(y_dy)
    DelZMin = minval(z_dz)
    Nu = AlphaSound * ( Min(DelXMin, DelZMin) ** 2.0d0 ) / DelTimeShort
    
    ! CReSS マニュアルでは, 2 次精度中心差分の場合, Alpha < 1/8 くらいが適当と述べている. 
    ! deepconv では NuH として 500 以上の値が入っていたそうだ,
    NuHh = AlphaNDiff * ( SQRT(DelXMin*DelYMin) ** 2.0d0 ) / DelTimeLong
    NuVh = AlphaNDiff * ( DelZMin ** 2.0d0 ) / DelTimeLong
    NuHm = NuHh * NDiffRatio
    NuVm = NuVh * NDiffRatio

    if (myrank == 0) then 
      call MessageNotify( "M", "DynamicsHEVI_init", "Nu = %f", d=(/Nu/) )
      call MessageNotify( "M", "DynamicsHEVI_init", "NuHh = %f", d=(/NuHh/) )
      call MessageNotify( "M", "DynamicsHEVI_init", "NuVh = %f", d=(/NuVh/) )
      call MessageNotify( "M", "DynamicsHEVI_init", "NuHm = %f", d=(/NuHm/) )
      call MessageNotify( "M", "DynamicsHEVI_init", "NuVm = %f", d=(/NuVm/) )
    end if

    ! 陰解法の計算設定の初期化
    !
    call DynamicsVI_init()


    call HistoryAutoAddVariable(  &
      & varname='PTempAdv', &
      & dims=(/'x','y','z','t'/),     &
      & longname='Advection term of potential temperature',  &
      & units='K.s-1',    &
      & xtype='double')
    
    call HistoryAutoAddVariable(  &
      & varname='PTempNDiff',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Numerical diffusion term of potential temperature',&
      & units='K.s-1',    &
      & xtype='double')

    call HistoryAutoAddVariable(  &
      & varname='CDensAdv', &
      & dims=(/'x','y','z','t'/),     &
      & longname='Advection term of cloud density',  &
      & units='K.m-3.s-1',    &
      & xtype='double')
    
    call HistoryAutoAddVariable(  &
      & varname='CDensNDiff',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Numerical diffusion term of cloud density',&
      & units='K.m-3.s-1',    &
      & xtype='double')

    do l = 1, ncmax
      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_Adv', &
        & dims=(/'x','y','z','t'/),     &
        & longname='Advection term of '          &
        &           //trim(SpcWetSymbol(l))//' mixing ratio',  &
        & units='kg.kg-1.s-1',    &
        & xtype='double')
      
      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_NDiff', & 
        & dims=(/'x','y','z','t'/),     &
        & longname='Diffusion term of '          &
        &           //trim(SpcWetSymbol(l))//' mixing ratio',  &
        & units='kg.kg-1.s-1',    &
        & xtype='double')
    end do

    call HistoryAutoAddVariable(  &
      & varname='VelXAdv', &
      & dims=(/'x','y','z','t'/),     &
      & longname='Advection term of velocity (x)',  &
      & units='K.m-3.s-1',    &
      & xtype='double')
    
    call HistoryAutoAddVariable(  &
      & varname='VelXNDiff',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Numerical diffusion term of velocity (x)',&
      & units='K.m-3.s-1',    &
      & xtype='double')

    call HistoryAutoAddVariable(  &
      & varname='VelXPGrad', &
      & dims=(/'x','y','z','t'/),     &
      & longname='Pressure gradient term of velocity (x)',  &
      & units='K.m-3.s-1',    &
      & xtype='double')

    call HistoryAutoAddVariable(  &
      & varname='VelYAdv', &
      & dims=(/'x','y','z','t'/),     &
      & longname='Advection term of velocity (y)',  &
      & units='K.m-3.s-1',    &
      & xtype='double')
    
    call HistoryAutoAddVariable(  &
      & varname='VelYNDiff',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Numerical diffusion term of velocity (y)',&
      & units='K.m-3.s-1',    &
      & xtype='double')

    call HistoryAutoAddVariable(  &
      & varname='VelYPGrad', &
      & dims=(/'x','y','z','t'/),     &
      & longname='Pressure gradient term of velocity (y)',  &
      & units='K.m-3.s-1',    &
      & xtype='double')

    call HistoryAutoAddVariable(  &
      & varname='VelZAdv', &
      & dims=(/'x','y','z','t'/),     &
      & longname='Advection term of velocity (z)',  &
      & units='K.m-3.s-1',    &
      & xtype='double')
    
    call HistoryAutoAddVariable(  &
      & varname='VelZNDiff',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Numerical diffusion term of Velocity (z)',&
      & units='K.m-3.s-1',    &
      & xtype='double')

    call HistoryAutoAddVariable(  &
      & varname='VelZBuoyT',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Buoyancy (Temperature)',&
      & units='K.m-3.s-1',    &
      & xtype='double')

    call HistoryAutoAddVariable(  &
      & varname='VelZBuoyM',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Buoyancy (MolWt)',&
      & units='K.m-3.s-1',    &
      & xtype='double')

    call HistoryAutoAddVariable(  &
      & varname='VelZBuoyD',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Buoyancy (Drag)',&
      & units='K.m-3.s-1',    &
      & xtype='double')

    call HistoryAutoAddVariable(  &
      & varname='VelZPGrad', &
      & dims=(/'x','y','z','t'/),     &
      & longname='Pressure gradient term of velocity (z)',  &
      & units='K.m-3.s-1',    &
      & xtype='double')

  end subroutine Dynamics_Init


  subroutine Dynamics_Long_forcing(        &
    & pyz_VelXBl,  pyz_VelXNl,    & ! (in)
    & xqz_VelYBl,  xqz_VelYNl,    & ! (in)
    & xyr_VelZBl,  xyr_VelZNl,    & ! (in)
    & xyz_PTempBl, xyz_PTempNl,   & ! (in)
    & xyzf_QMixBl, xyzf_QMixNl,   & ! (in)
    & xyz_KmBl,    xyz_KmNl,      & ! (in)
    & xyz_CDensBl, xyz_CDensNl,   & ! (in)
    & pyz_DVelXDtNl,              & ! (inout)
    & xqz_DVelYDtNl,              & ! (inout)
    & xyr_DVelZDtNl,              & ! (inout)
    & xyz_DPTempDtNl,             & ! (inout)
    & xyzf_DQMixDtNl,             & ! (inout)
    & xyz_DKmDtNl,                & ! (inout)
    & xyz_DCDensDtNl,             & ! (inout)
    & xyz_PTempAl,                & ! (out)
    & xyzf_QMixAl                 & ! (out)
    & )

    implicit none

    real(DP), intent(in) :: pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(in) :: xyz_PTempBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyz_PTempNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyzf_QMixBl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP), intent(in) :: xyzf_QMixNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP), intent(in) :: xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyz_KmNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: pyz_DVelXDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xqz_DVelYDtNl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(inout) :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(inout) :: xyz_DPTempDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xyzf_DQMixDtNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP), intent(inout) :: xyz_DKmDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyz_PTempAl(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP), intent(out) :: xyzf_QMixAl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP), intent(in) :: xyz_CDensBl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyz_CDensNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout) :: xyz_DCDensDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: pyz_Orig(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: pyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: pyz_NDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_Orig(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_NDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_Orig(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_NDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_BuoyT(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_BuoyM(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_BuoyD(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_Orig(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_NDiff(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_PTempAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyzf_Orig(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyzf_Adv(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyzf_NDiff(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyzf_QMixPerMolWt(imin:imax,jmin:jmax,kmin:kmax, 1:GasNum)
    integer              :: f

    !------------------------------
    ! tendency of cloud density    
    ! 

    ! initialize 
    xyz_Orig = xyz_DCDensDtNl

    ! フラックス項の計算. 4 次精度中心差分を用いて計算
    !
    xyz_Adv = &
      & - xyz_c4dx_pyz(pyz_VelXNl * pyz_avr_xyz(xyz_CDensNl)) &
      & - xyz_c4dy_xqz(xqz_VelYNl * xqz_avr_xyz(xyz_CDensNl)) &
      & - xyz_c4dz_xyr(xyr_VelZNl * xyr_avr_xyz(xyz_CDensNl)) 

    ! 数値粘性項の計算
    xyz_NDiff = &
      &  + NuHh * (xyz_dx_pyz(pyz_dx_xyz(xyz_CDensBl))) &
      &  + NuHh * (xyz_dy_xqz(xqz_dy_xyz(xyz_CDensBl))) &
      &  + NuVh * (xyz_dz_xyr(xyr_dz_xyz(xyz_CDensBl))) 
    
    xyz_DCDensDtNl = xyz_Orig + xyz_Adv + xyz_NDiff

    call HistoryAutoPut(TimeN, 'CDensAdv',  xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'CDensNDiff', xyz_NDiff(1:nx,1:ny,1:nz))

    !----------------------------------
    ! 拡散係数
    !

    ! Initialize
    !
    xyz_Orig = xyz_DKmDtNl

    ! Advection term
    !
    xyz_Adv  =                                              &
      & - xyz_avr_pyz(pyz_VelXNl * pyz_c4dx_xyz(xyz_KmNl))  &
      & - xyz_avr_xqz(xqz_VelYNl * xqz_c4dy_xyz(xyz_KmNl))  &
      & - xyz_avr_xyr(xyr_VelZNl * xyr_c4dz_xyz(xyz_KmNl))    

    ! Numerical diffusion term 
    !
    xyz_NDiff =                                        &
      &    NuHm * (xyz_dx_pyz(pyz_dx_xyz( xyz_KmBl ))) &
      &  + NuHm * (xyz_dy_xqz(xqz_dy_xyz( xyz_KmBl ))) &
      &  + NuVm * (xyz_dz_xyr(xyr_dz_xyz( xyz_KmBl ))) 

    xyz_DKmDtNl = xyz_Orig + xyz_Adv + xyz_NDiff


    ! tendency of potential temperature
    ! 

    ! initialize 
    xyz_Orig = xyz_DPTempDtNl
    xyz_PTempAll = xyz_PTempNl + xyz_PTempBZ

    ! Advection term
    ! xyz_AdvScalar( xyz_PTempNl + xyz_PTempBZ, pyz_VelXNl, pyz_VelXNl, xyr_VelZNl) 
    !
    xyz_Adv = &
      & - xyz_avr_pyz(pyz_VelXNl * pyz_c4dx_xyz(xyz_PTempAll))  &
      & - xyz_avr_xqz(xqz_VelYNl * xqz_c4dy_xyz(xyz_PTempAll))  &
      & - xyz_avr_xyr(xyr_VelZNl * xyr_c4dz_xyz(xyz_PTempAll))  

    ! numerical diffusion term
    ! xyz_Num = xyz_NumDiffScalar( xyz_PTempBl)
    !
    xyz_NDiff = &
      &    NuHh * (xyz_dx_pyz(pyz_dx_xyz( xyz_PTempBl ))) &
      &  + NuHh * (xyz_dy_xqz(xqz_dy_xyz( xyz_PTempBl ))) &
      &  + NuVh * (xyz_dz_xyr(xyr_dz_xyz( xyz_PTempBl ))) 

    ! sum
    !
    xyz_DPTempDtNl = xyz_Orig + xyz_Adv + xyz_NDiff
    
    ! output
    !
    call HistoryAutoPut(TimeN, 'PTempAdv',   xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'PTempNDiff', xyz_NDiff(1:nx,1:ny,1:nz))

    xyz_PTempAl = xyz_PTempBl + 2.0d0 * DelTimeLong * xyz_DPTempDtNl

    ! Set Margin
    !
    call SetMargin_xyz(xyz_PTempAl)

    
    !------------------------------
    ! tendency of mixing ratio
    ! 

    ! initialize
    xyzf_Orig = xyzf_DQMixDtNl
    xyzf_QMixAll = xyzf_QMixNl + xyzf_QMixBZ

    do f = 1, ncmax
      ! Advection term
      !xyzf_Adv  = xyzf_AdvScalar(xyzf_QMixNl + xyzf_QMixBZ, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
      !
      xyzf_Adv(:,:,:,f) = &
        & - xyz_avr_pyz(pyz_VelXNl * pyz_c4dx_xyz(xyzf_QMixAll(:,:,:,f))) &
        & - xyz_avr_xqz(xqz_VelYNl * xqz_c4dy_xyz(xyzf_QMixAll(:,:,:,f))) &
        & - xyz_avr_xyr(xyr_VelZNl * xyr_c4dz_xyz(xyzf_QMixAll(:,:,:,f)))    

      ! numerical diffusion term
      ! xyzf_Diff = xyzf_NumDiffScalar(xyzf_QMixBl) 
      !
      xyzf_NDiff(:,:,:,f) = &
        &    NuHh * (xyz_dx_pyz(pyz_dx_xyz( xyzf_QMixBl(:,:,:,f) ))) &
        &  + NuHh * (xyz_dy_xqz(xqz_dy_xyz( xyzf_QMixBl(:,:,:,f) ))) &
        &  + NuVh * (xyz_dz_xyr(xyr_dz_xyz( xyzf_QMixBl(:,:,:,f) ))) 
    end do

    ! sum
    !
    xyzf_DQMixDtNl = xyzf_Orig + xyzf_Adv + xyzf_NDiff

    ! output
    !
    do f = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_Adv',   xyzf_Adv(1:nx,1:ny,1:nz,f))
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_NDiff', xyzf_NDiff(1:nx,1:ny,1:nz,f))
    end do

    ! time integration
    !
    xyzf_QMixAl = xyzf_QMixBl + 2.0d0 * DelTimeLong * xyzf_DQMixDtNl

    ! Set Margin
    ! 
    call SetMargin_xyzf(xyzf_QMixAl)

    ! 負の値を埋める
    !
    call FillNegativeQMix(xyzf_QMixAl)

    ! Set Margin
    ! 
    call SetMargin_xyzf(xyzf_QMixAl)

    !------------------------------
    ! tendency of VelX
    ! 

    ! initializa
    !
    pyz_Orig = pyz_DVelXDtNl

    ! Advection term
    !pyz_Adv  = pyz_AdvVelX(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) 
    !
    pyz_Adv  = &
      & - pyz_VelXNl * pyz_avr_xyz( xyz_c4dx_pyz( pyz_VelXNl ) )                &
      & - pyz_avr_pqz( pqz_avr_xqz( xqz_VelYNl ) * pqz_c4dy_pyz( pyz_VelXNl ) ) &
      & - pyz_avr_pyr( pyr_avr_xyr( xyr_VelZNl ) * pyr_c4dz_pyz( pyz_VelXNl ) )

    ! Numerical diffusion term 
    !pyz_Diff = pyz_NumDiffVelX(pyz_VelXBl)
    pyz_NDiff = &
      &   NuHm * ( pyz_dx_xyz( xyz_dx_pyz( pyz_VelXBl ) ) ) &
      & + NuHm * ( pyz_dy_pqz( pqz_dy_pyz( pyz_VelXBl ) ) ) &
      & + NuVm * ( pyz_dz_pyr( pyr_dz_pyz( pyz_VelXBl ) ) )

    ! sum
    !
    pyz_DVelXDtNl = pyz_Orig + pyz_Adv + pyz_NDiff

    call HistoryAutoPut(TimeN, 'VelXAdv',   pyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXNDiff', pyz_NDiff(1:nx,1:ny,1:nz))

    !------------------------------
    ! tendency of VelY
    !     

    ! ininitalize
    xqz_Orig = xqz_DVelYDtNl

    ! Advection term
    ! xqz_Adv  = xqz_AdvVelY(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
    !
    xqz_Adv  = &
      & - xqz_avr_pqz( pqz_avr_pyz( pyz_VelXNl ) * pqz_c4dx_xqz( xqz_VelYNl ) ) &
      & - xqz_VelYNl * xqz_avr_xyz( xyz_c4dy_xqz( xqz_VelYNl ) ) &
      & - xqz_avr_xqr( xqr_avr_xyr( xyr_VelZNl ) * xqr_c4dz_xqz( xqz_VelYNl ) )

    ! Numerical diffusion term
    ! xqz_Diff = xqz_NumDiffVelY(xqz_VelYBl)
    !
    xqz_NDiff = &
      &   NuHm * ( xqz_dx_pqz( pqz_dx_xqz( xqz_VelYBl ) ) ) &
      & + NuHm * ( xqz_dy_xyz( xyz_dy_xqz( xqz_VelYBl ) ) ) &
      & + NuVm * ( xqz_dz_xqr( xqr_dz_xqz( xqz_VelYBl ) ) )

    ! sum
    !
    xqz_DVelYDtNl = xqz_Orig + xqz_Adv + xqz_NDiff
    
    call HistoryAutoPut(TimeN, 'VelYAdv',   xqz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYNDiff', xqz_NDiff(1:nx,1:ny,1:nz))

    !------------------------------
    ! tendency of VelZ
    ! 

    ! Initialization
    !
    xyr_Orig = xyr_DVelZDtNl
    
    do f = 1, GasNum
      xyzf_QMixPerMolWt(:,:,:,f) = xyzf_QMixNl(:,:,:,IdxG(f)) / MolWtWet(IdxG(f))
    end do

    ! Advection term
    ! xyr_Adv  = xyr_AdvVelZ(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
    !
    xyr_Adv  = &
      & - xyr_avr_pyr( pyr_avr_pyz( pyz_VelXNl ) * pyr_c4dx_xyr( xyr_VelZNl ) ) &
      & - xyr_avr_xqr( xqr_avr_xqz( xqz_VelYNl ) * xqr_c4dy_xyr( xyr_VelZNl ) ) &
      & - xyr_VelZNl * xyr_avr_xyz( xyz_c4dz_xyr( xyr_VelZNl ) )

    ! Numerical diffusion term
    !xyr_Diff = xyr_NumDiffVelZ(xyr_VelZBl)
    !
    xyr_NDiff = &
      &   NuHm * ( xyr_dx_pyr( pyr_dx_xyr( xyr_VelZBl ) ) ) &
      & + NuHm * ( xyr_dy_xqr( xqr_dy_xyr( xyr_VelZBl ) ) ) & 
      & + NuVm * ( xyr_dz_xyz( xyz_dz_xyr( xyr_VelZBl ) ) )

    ! Buoyancy due to temperature disturbunce
    !xyr_BuoyT = xyr_Buoy(xyz_PTempNl)
    !
    xyr_BuoyT = Grav * xyr_avr_xyz( xyz_PTempNl / xyz_PTempBZ)

    ! Buoyancy due to molecular weight
    !
    ! 全量でなくて良いのかな??
    xyr_BuoyM =                                                   &
      & + Grav * xyr_avr_xyz( sum(xyzf_QMixPerMolWt, 4) )        &
      &    / ( 1.0d0 / MolWtDry + xyr_QMixBZPerMolWt )       &
      & - Grav * xyr_avr_xyz( sum(xyzf_QMixNl(:,:,:,1:GasNum), 4) )&
      &    / ( 1.0d0 + xyr_QmixBZ )

    ! Buoyancy due to loading
    !
    xyr_BuoyD =                                                           &
      & - Grav * xyr_avr_xyz( sum(xyzf_QMixNl(:,:,:,GasNum+1:ncmax), 4) ) &
      &    / ( 1.0d0 + xyr_QMixBZ )

    ! sum
    !
    xyr_DVelZDtNl = xyr_Orig + xyr_Adv + xyr_NDiff + xyr_BuoyT + xyr_BuoyM + xyr_BuoyD

    call HistoryAutoPut(TimeN, 'VelZAdv',   xyr_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelZNDiff', xyr_NDiff(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'VelZBuoyT', xyr_BuoyT(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'VelZBuoyM', xyr_BuoyM(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'VelZBuoyD', xyr_BuoyD(1:nx,1:ny,1:nz))    

  end subroutine Dynamics_Long_forcing
 

  subroutine Dynamics_Short_forcing(  &
        &  pyz_VelXNs,          & ! (in)
        &  xqz_VelYNs,          & ! (in)
        &  xyr_VelZNs,          & ! (in)
        &  xyz_ExnerNs,         & ! (in)
        &  pyz_DVelXDtNl,       & ! (in)
        &  xqz_DVelYDtNl,       & ! (in)
        &  xyr_DVelZDtNl,       & ! (in)
        &  xyz_DExnerDtNs,      & ! (in)
        &  pyz_VelXAs,          & ! (out)
        &  xqz_VelYAs,          & ! (out)
        &  xyr_VelZAs,          & ! (out)
        &  xyz_ExnerAs          & ! (out)
        & )

    real(DP), intent(in)  :: pyz_VelXNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xqz_VelYNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyr_VelZNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyz_ExnerNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: pyz_DVelXDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xqz_DVelYDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyz_DExnerDtNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: pyz_VelXAs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xqz_VelYAs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyr_VelZAs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(out) :: xyz_ExnerAs(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: pyz_DVelXDtNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xqz_DVelYDtNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyr_DVelZDtNs(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_DivVel(imin:imax,jmin:jmax,kmin:kmax)

    ! initialize: Divergence of velocity
    !
    xyz_DivVel = &
      &   xyz_dx_pyz( pyz_VelXNs ) &
      & + xyz_dy_xqz( xqz_VelYNs ) &
      & + xyz_dz_xyr( xyr_VelZNs )
    
    !--------------------------------------
    ! VelX
    !
    pyz_DVelXDtNs =  &
      &   - pyz_avr_xyz( CpDry * xyz_PTempBZ ) &
      &     * ( pyz_dx_xyz( xyz_ExnerNs ) - pyz_dx_xyz( Nu * xyz_DivVel ) )  
    
    ! Time integration
    !
    pyz_VelXAs    = pyz_VelXNs + DelTimeShort * (pyz_DVelXDtNl + pyz_DVelXDtNs)

    call HistoryAutoPut(TimeN, 'VelXPGrad', pyz_DVelXDtNs(1:nx,1:ny,1:nz))

    ! Set Margin
    !
    call SetMargin_pyz( pyz_VelXAs ) ! (inout)

    !--------------------------------------
    ! VelY
    !    
    xqz_DVelYDtNs =  &
      &   - xqz_avr_xyz( CpDry * xyz_PTempBZ ) &
      &     * ( xqz_dy_xyz( xyz_ExnerNs ) - xqz_dy_xyz( Nu * xyz_DivVel ) )
    
    ! Time integration
    !
    xqz_VelYAs = xqz_VelYNs + DelTimeShort * (xqz_DVelYDtNl + xqz_DVelYDtNs)

    call HistoryAutoPut(TimeN, 'VelYPGrad', xqz_DVelYDtNs(1:nx,1:ny,1:nz))

    ! Set Margin
    !
    call SetMargin_xqz( xqz_VelYAs ) ! (inout)
    
    !--------------------------------------
    ! Exner function
    !
    xyz_ExnerAs = xyz_Exner( &
      & pyz_VelXNs,          &
      & pyz_VelXAs,          &
      & xqz_VelYNs,          &
      & xqz_VelYAs,          &
      & xyr_VelZNs,          &
      & xyz_ExnerNs,         &
      & xyr_DVelZDtNl,       &
      & xyz_DExnerDtNs       &
      & )

    ! Set Margin
    !
    call SetMargin_xyz( xyz_ExnerAs ) ! (inout)

    !--------------------------------------
    ! VelZ
    !
    xyr_DVelZDtNs =                                         &
      & - xyr_avr_xyz(CpDry * xyz_PTempBZ )                   &
      &   * (                                               &
      &       beta * xyr_dz_xyz( xyz_ExnerAs )              &
      &       + (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerNs )  &
      &       - xyr_dz_xyz( Nu * xyz_DivVel )               &
      &     )                                                  
    
    ! Time integration
    !
    xyr_VelZAs = xyr_VelZNs + DelTimeShort * (xyr_DVelZDtNl + xyr_DVelZDtNs)

    call HistoryAutoPut(TimeN, 'VelZPGrad', xyr_DVelZDtNs(1:nx,1:ny,1:nz))

    ! Set Margin
    !
    call SetMargin_xyr( xyr_VelZAs ) ! (inout)

  end subroutine Dynamics_Short_forcing

  
!!!--------------------------------------------------------------------!!!
  subroutine DynamicsVI_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) )
    allocate( B(RegZMin+1:RegZMax) )
    allocate( C(RegZMin:RegZMax-1) )
    allocate( xyz_F1BZ(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyr_F2BZ(imin:imax,jmin:jmax,kmin:kmax) )

    !----------------------------------------------------------------
    ! 係数行列と共通して利用される配列の値を決める
    !----------------------------------------------------------------

    !係数行列の計算
    !  A, B, C を求める際, F1BZ と F2BZ は X 方向に一様なので. 
    !  RegXMax, RegYMax の値で代表させることとした. 
    xyz_F1BZ =                                                &
      &  ( xyz_VelSoundBZ ** 2.0d0 )                          &
      &   / (CpDry * xyz_DensBZ * (xyz_VPTempBZ ** 2.0d0))

    xyr_F2BZ =                                                 &
      &  xyr_avr_xyz(                                               &
      &    CpDry * xyz_DensBZ * ( xyz_VPTempBZ ** 2.0d0 ) &
      &   )
    
    A(RegZMin+1: RegZMax-1) =                                &
      & (beta ** 2.0d0)                                    &
      &    * xyz_F1BZ(RegXMax,RegYMax,RegZMin+1: RegZMax-1)      &  
      &    * (                                               &
      &          xyr_F2BZ(RegXMax,RegYMax,RegZMin+1: RegZMax-1)  &
      &            / r_dz(RegZMin+1: RegZMax-1)              &
      &        + xyr_F2BZ(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_F1BZ(RegXMax,RegYMax,RegZMin)       &
      &   * xyr_F2BZ(RegXMax,RegYMax,RegZMin)       &
      &     / r_dz(RegZMin)                             &
      &   * (DTS ** 2.0d0)                              &
      &   / z_dz(RegZMin)                               &
      & + 1.0d0                                         

    A(RegZMax) =                                        &
      & (beta ** 2.0d0)                               &
      &   * xyz_F1BZ(RegXMax,RegYMax,RegZMax)       &
      &   * xyr_F2BZ(RegXMax,RegYMax,RegZMax-1)     &
      &     / r_dz(RegZMax-1)                           &
      &   * (DTS ** 2.0d0)                              &
      &   / z_dz(RegZMax)                               &
      & + 1.0d0                                         
    
    B(RegZMin+1:RegZMax) =                             &
      & - (beta ** 2.0d0)                               &
      &   * xyz_F1BZ(RegXMax,RegYMax,RegZMin:RegZMax-1)    &
      &   * xyr_F2BZ(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_F1BZ(RegXMax,RegYMax,RegZMin+1:RegZMax)    &
      &   * xyr_F2BZ(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) )
    allocate( AL1(N), IP(N) )

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

   
  end subroutine DynamicsVI_init
  

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

    !暗黙の型宣言禁止
    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_DVelZDtNl &
      &                     (imin:imax,jmin:jmax,kmin:kmax) 
                                                           !Z 方向の外力項[t]
    real(DP), intent(in)   :: xyz_DExnerDtNs &
      &                     (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 ! 短い時間格子間隔
    integer   :: ix, jy, kz

    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( Nu * xyz_DivVelNs )                       &
      & - ( 1.0d0 - beta ) * xyr_dz_xyz( xyz_ExnerNs )        &
      & + xyr_DVelZDtNl / xyr_avr_xyz( CpDry * xyz_VPTempBZ ) 
    
    F = - beta * xyz_F1BZ * DTS                                     &
      & * xyz_dz_xyr(                                               &
      &    xyr_avr_xyz( xyz_DensBZ * xyz_VPTempBZ)                  &
      &    * (                                                      &
      &         xyr_VelZNs                                          &
      &       - xyr_avr_xyz(CpDry * xyz_VPTempBZ) * DTS             &
      &         * (                                                 &
      &               (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerNs )    &
      &             - xyr_dz_xyz( Nu * xyz_DivVelNs )               &
      &            )                                                &
      &       + xyr_DVelZDtNl * DTS                                 &
      &      )                                                      &
      &   )                                                         &
      & + xyz_DExnerDtNs * DTS

    D1 = xyz_ExnerNs                                                &
      & - (1.0d0 - beta)                                            &
      &   * xyz_F1BZ * DTS                                          &
      &   * xyz_dz_xyr(                                             &
      &       xyr_avr_xyz(xyz_DensBZ * xyz_VPTempBZ) * xyr_VelZNs   &
      &     )                                                       &
      & - (xyz_VelSoundBZ ** 2.0d0) * DTS                           &
      &   / (CpDry * xyz_VPTempBZ) * xyz_dx_pyz( pyz_VelXAs )       &
      & - (xyz_VelSoundBZ ** 2.0d0) * DTS                           & 
      &   / (CpDry * xyz_VPTempBZ) * xyz_dy_xqz( xqz_VelYAs )       &
      & + F

    D1(:,:,RegZMin) = D1(:,:,RegZMin)                           &
      & - beta * xyz_F1BZ(:,:,RegZMin) * (DTS ** 2.0d0)         &
      &   * xyr_F2BZ(:,:, RegZMin-1) * E(:,:,RegZMin-1)         &
      &   / z_dz(RegZMin)
    
    D1(:,:,RegZMax) = D1(:,:,RegZMax)                           &
      & + beta * xyz_F1BZ(:,:,RegZMax) * (DTS ** 2.0d0)         &
      &   * xyr_F2BZ(:,:,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

end module DynamicsHEVI

