subroutine Dynamics_Init

    !暗黙の型宣言禁止
    implicit none
    
    real(DP)  :: DelXMin, DelYMin, DelZMin
    real(DP)  :: AlphaSound = 5.0d-2  !音波減衰項の係数 (気象庁数値予報課報告・別冊49 より)
    real(DP)  :: AlphaNDiff = 1.0d-3  !4次の数値拡散の係数. CReSS マニュアルより
    real(DP)  :: NDiffRatio = 1.0d0  !速度に対する粘性を上げる場合は数字を 1 以上にする. 
    integer   :: unit            !装置番号
    integer   :: l

    NAMELIST /Dynamics_nml/ AlphaSound, AlphaNDiff, NDiffRatio, beta

    !ファイルオープン. 情報取得. 
    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)
    AlphaH = AlphaSound * ( SQRT(DelXMin*DelYMin) ** 2.0d0 ) / DelTimeShort
    AlphaV = AlphaSound * ( DelZMin ** 2.0d0 ) / DelTimeShort
    
    NuHh = AlphaNDiff * ( SQRT(DelXMin*DelYMin) ** 4.0d0 ) / (2.0d0 * DelTimeLong)
    NuVh = AlphaNDiff * ( DelZMin ** 4.0d0 ) / (2.0d0 * DelTimeLong)
    NuHm = NuHh * NDiffRatio
    NuVm = NuVh * NDiffRatio

    if (myrank == 0) then 
      call MessageNotify( "M", module_name, "AlphaH = %f", d=(/AlphaH/) )
      call MessageNotify( "M", module_name, "AlphaV = %f", d=(/AlphaV/) )
      call MessageNotify( "M", module_name, "NuHh = %f", d=(/NuHh/) )
      call MessageNotify( "M", module_name, "NuVh = %f", d=(/NuVh/) )
      call MessageNotify( "M", module_name, "NuHm = %f", d=(/NuHm/) )
      call MessageNotify( "M", module_name, "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='float')
    
    call HistoryAutoAddVariable( varname='PTempNDiff', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of potential temperature', units='K.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='ExnerAdv', dims=(/'x','y','z','t'/), longname='Advection term of exner function', units='s-1', xtype='float')
    
    call HistoryAutoAddVariable( varname='ExnerNDiff', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of exner function', units='s-1', xtype='float')

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

    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='float')
      
      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='float')
    end do

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

    call HistoryAutoAddVariable( varname='VelXPGrad', dims=(/'x','y','z','t'/), longname='Pressure gradient term of velocity (x)', units='m.s-2', xtype='float')

    call HistoryAutoAddVariable( varname='VelXSWF', dims=(/'x','y','z','t'/), longname='Filter for acoustic mode (x)', units='m.s-2', xtype='float')

!    call HistoryAutoAddVariable(  &
!      & varname='VelXTndNs', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='Velocity Tendency (x)',  &
!      & units='m.s-2',    &
!      & xtype='float')

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

    call HistoryAutoAddVariable( varname='VelYPGrad', dims=(/'x','y','z','t'/), longname='Pressure gradient term of velocity (y)', units='m.s-2', xtype='float')

    call HistoryAutoAddVariable( varname='VelYSWF', dims=(/'x','y','z','t'/), longname='Filter for acoustic mode (y)', units='m.s-2', xtype='float')

!    call HistoryAutoAddVariable(  &
!      & varname='VelYTndNs', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='Velocity Tendency (y)',  &
!      & units='m.s-2',    &
!      & xtype='float')

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

    call HistoryAutoAddVariable( varname='VelZBuoyT', dims=(/'x','y','z','t'/), longname='Buoyancy (Temperature)', units='m.s-2', xtype='float')

    call HistoryAutoAddVariable( varname='VelZBuoyM', dims=(/'x','y','z','t'/), longname='Buoyancy (MolWt)', units='m.s-2', xtype='float')

    call HistoryAutoAddVariable( varname='VelZBuoyD', dims=(/'x','y','z','t'/), longname='Buoyancy (Drag)', units='m.s-2', xtype='float')

    call HistoryAutoAddVariable( varname='VelZPGrad', dims=(/'x','y','z','t'/), longname='Pressure gradient term of velocity (z)', units='m.s-2', xtype='float')

    call HistoryAutoAddVariable( varname='VelZSWF', dims=(/'x','y','z','t'/), longname='Filter for acoustic mode (z)', units='m.s-2', xtype='float')

!    call HistoryAutoAddVariable(  &
!      & varname='VelZTndNs', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='Velocity Tendency (z)',  &
!      & units='m.s-2',    &
!      & xtype='float')

!    call HistoryAutoAddVariable(  &
!      & varname='VelDiv', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='Velocity Divergence',  &
!      & units='s-1',    &
!      & xtype='float')

!    call HistoryAutoAddVariable(  &
!      & varname='VelLapla', &
!      & dims=(/'x','y','z','t'/),     &
!      & longname='Velocity Lapla',  &
!      & units='m-1.s-1',    &
!      & xtype='float')

  end subroutine Dynamics_Init

  subroutine Dynamics_Km_forcing( pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, xyz_KmBl,    xyz_KmNl, xyz_DKmDtNl )

    implicit none

    real(DP), intent(in) :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax) 
    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) :: xyz_DKmDtNl(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)

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

    ! 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_dx_pyz(pyz_dx_xyz( xyz_KmBl ))))) - NuHm * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_KmBl ))))) - NuVm * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyz_KmBl ))))) 
    
    xyz_DKmDtNl = xyz_Orig + xyz_Adv + xyz_NDiff
    
    ! Set Margin
    !
!    call SetMargin_xyz( xyz_DKmDtNl )

  end subroutine Dynamics_Km_forcing

  subroutine Dynamics_Long_forcing( pyz_VelXBl,  pyz_VelXNl, xqz_VelYBl,  xqz_VelYNl, xyr_VelZBl,  xyr_VelZNl, xyz_PTempBl, xyz_PTempNl, xyzf_QMixBl, xyzf_QMixNl, xyz_CDensBl, xyz_CDensNl, pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, xyz_DPTempDtNl, xyzf_DQMixDtNl, xyz_DCDensDtNl, xyz_PTempAl, xyzf_QMixAl )

    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(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(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_dx_pyz(pyz_dx_xyz(xyz_CDensBl))))) - NuHh * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_CDensBl))))) - NuVh * (xyz_dz_xyr(xyr_dz_xyz(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))


    ! 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_dx_pyz(pyz_dx_xyz( xyz_PTempBl ))))) - NuHh * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_PTempBl ))))) - NuVh * (xyz_dz_xyr(xyr_dz_xyz(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(xyz_dx_pyz(pyz_dx_xyz( xyzf_QMixBl(:,:,:,f) ))))) - NuHh * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyzf_QMixBl(:,:,:,f) ))))) - NuVh * (xyz_dz_xyr(xyr_dz_xyz(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_dx_xyz(xyz_dx_pyz( pyz_VelXBl ))))) - NuHm * ( pyz_dy_pqz(pqz_dy_pyz(pyz_dy_pqz(pqz_dy_pyz( pyz_VelXBl ))))) - NuVm * ( pyz_dz_pyr(pyr_dz_pyz(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_dx_pqz(pqz_dx_xqz( xqz_VelYBl ))))) - NuHm * ( xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz( xqz_VelYBl ))))) - NuVm * ( xqz_dz_xqr(xqr_dz_xqz(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_dx_pyr(pyr_dx_xyr( xyr_VelZBl ))))) - NuHm * ( xyr_dy_xqr(xqr_dy_xyr(xyr_dy_xqr(xqr_dy_xyr( xyr_VelZBl ))))) - NuVm * ( xyr_dz_xyz(xyz_dz_xyr(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))    

    ! Set Margin
    !
!    call SetMargin_pyz( pyz_DVelXDtNl )
!    call SetMargin_xqz( xqz_DVelYDtNl )
!    call SetMargin_xyr( xyr_DVelZDtNl )
!    call SetMargin_xyz( xyz_DPTempDtNl )
!    call SetMargin_xyz( xyz_DCDensDtNl )
!    call SetMargin_xyzf(xyzf_DQMixDtNl )
 
  end subroutine Dynamics_Long_forcing

  subroutine Dynamics_Short_forcing( pyz_VelXNs, xqz_VelYNs, xyr_VelZNs, xyz_ExnerNs, pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs, pyz_VelXAs, xqz_VelYAs, xyr_VelZAs, xyz_ExnerAs )

    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(inout)  :: xyz_DExnerDtNl(imin:imax,jmin:jmax,kmin:kmax) !test
    real(DP), intent(inout)  :: xyz_DExnerDtNs(imin:imax,jmin:jmax,kmin:kmax) !test
    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_VelDivNs(imin:imax,jmin:jmax,kmin:kmax)
!    real(DP) :: xyz_VelLaplaNs(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: pyz_VelXPGrad(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xqz_VelYPGrad(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyr_VelZPGrad(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: pyz_VelXSWF(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xqz_VelYSWF(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyr_VelZSWF(imin:imax,jmin:jmax,kmin:kmax)

    real(DP) :: xyz_ExnerAll(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)

    ! initialize: Divergence of velocity
    !
    xyz_VelDivNs = xyz_dx_pyz( pyz_VelXNs ) + xyz_dy_xqz( xqz_VelYNs ) + xyz_dz_xyr( xyr_VelZNs )
!    call HistoryAutoPut(TimeN, 'VelDiv', xyz_VelDivNs(1:nx,1:ny,1:nz))

!    xyz_VelLaplaNs = pyz_dx_xyz( xyz_VelDivNs ) +  xqz_dy_xyz( xyz_VelDivNs ) +  xyr_dz_xyz( xyz_VelDivNs )
!    call HistoryAutoPut(TimeN, 'VelLapla', xyz_VelLaplaNs(1:nx,1:ny,1:nz))

    !--------------------------------------
    ! VelX
    !
    pyz_VelXSWF   =   AlphaH * pyz_dx_xyz( xyz_VelDivNs ) 
    pyz_VelXPGrad = - pyz_avr_xyz( CpDry * xyz_VPTempBZ ) * pyz_dx_xyz( xyz_ExnerNs ) + pyz_VelXSWF
    pyz_DVelXDtNs =   pyz_VelXPGrad 

    ! Time integration
    !
    pyz_VelXAs    = pyz_VelXNs + DelTimeShort * (pyz_DVelXDtNl + pyz_DVelXDtNs)

    call HistoryAutoPut(TimeN, 'VelXPGrad', pyz_VelXPGrad(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXSWF',   pyz_VelXSWF(1:nx,1:ny,1:nz))
!    call HistoryAutoPut(TimeN, 'VelXTndNs', pyz_DVelXDtNs(1:nx,1:ny,1:nz))

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

    !--------------------------------------
    ! VelY
    !    
    xqz_VelYSWF   =   AlphaH * xqz_dy_xyz( xyz_VelDivNs ) 
    xqz_VelYPGrad = - xqz_avr_xyz( CpDry * xyz_VPTempBZ ) * xqz_dy_xyz( xyz_ExnerNs ) + xqz_VelYSWF
    xqz_DVelYDtNs =   xqz_VelYPGrad
    
    ! Time integration
    !
    xqz_VelYAs = xqz_VelYNs + DelTimeShort * (xqz_DVelYDtNl + xqz_DVelYDtNs)

    call HistoryAutoPut(TimeN, 'VelYPGrad', xqz_VelYPGrad(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYSWF',   xqz_VelYSWF(1:nx,1:ny,1:nz))
!    call HistoryAutoPut(TimeN, 'VelYTndNs', xqz_DVelYDtNs(1:nx,1:ny,1:nz))

    ! Set Margin
    !
    call SetMargin_xqz( xqz_VelYAs ) ! (inout)
    
    !--------------------------------------
    ! Exner function
    !

    ! フラックス項の計算. 4 次精度中心差分を用いて計算
    !
    xyz_Adv = - xyz_avr_pyz(pyz_VelXNs * pyz_c4dx_xyz(xyz_ExnerNs)) - xyz_avr_xqz(xqz_VelYNs * xqz_c4dy_xyz(xyz_ExnerNs)) - xyz_avr_xyr(xyr_VelZNs * xyr_c4dz_xyz(xyz_ExnerNs))  !&  
!      & + CpDry / CvDry * GasRDry * xyz_ExnerNs * xyz_VelDivNs

    ! 数値粘性項の計算
    xyz_NDiff = - NuHh * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_ExnerNs))))) - NuHh * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_ExnerNs))))) - NuVh * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_ExnerNs))))) 

    ! 短い時間ステップでのtendency 
    xyz_DExnerDtNs = xyz_DExnerDtNs + xyz_Adv + xyz_NDiff

    !!!!!!!!!!!!!!!!!!!!!!!
    !
    ! test 
    !
    !!!!!!!!!!!!!!!!!!!!!!!
    xyz_DExnerDtNs = 0.0d0

    call HistoryAutoPut(TimeN, 'ExnerAdv',   xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerNDiff', xyz_NDiff(1:nx,1:ny,1:nz))

    xyz_ExnerAs = xyz_Exner( pyz_VelXAs, xqz_VelYAs, xyr_VelZNs, xyz_VelDivNs, xyz_ExnerNs, xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs )

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

!    write(*,*) "+++ Exner +++", xyz_ExnerAs(imin:imax,1,kmin:kmax)

    !--------------------------------------
    ! VelZ
    !
    xyr_VelZSWF =  AlphaV * xyr_dz_xyz( xyz_VelDivNs ) 
    xyr_VelZPGrad = - xyr_avr_xyz(CpDry * xyz_VPTempBZ ) * ( beta * xyr_dz_xyz( xyz_ExnerAs ) + (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerNs ) ) + xyr_VelZSWF
    xyr_DVelZDtNs = xyr_VelZPGrad 

    !OK
    ! write(*,*) "+++ DVelZDtNl +++", xyr_DVelZDtNl(1:nx,1,1:nz)

    !OK
    !write(*,*) "+++ PGrad +++", xyr_VelZPGrad(1:nx,1,1:nz)
    
    ! Time integration
    !
    xyr_VelZAs = xyr_VelZNs + DelTimeShort * (xyr_DVelZDtNl + xyr_DVelZDtNs)

    call HistoryAutoPut(TimeN, 'VelZPGrad', xyr_VelZPGrad(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelZSWF',   xyr_VelZSWF(1:nx,1:ny,1:nz))
!    call HistoryAutoPut(TimeN, 'VelZTndNs', xyr_DVelZDtNs(1:nx,1:ny,1:nz))

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

    !OK
    !write(*,*) "+++ VelZ +++", xyr_VelZAs(1:nx,1,1:nz)
    !stop

  end subroutine Dynamics_Short_forcing