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, xqz_VelYBl,  xqz_VelYNl, xyr_VelZBl,  xyr_VelZNl, xyz_PTempBl, xyz_PTempNl, xyzf_QMixBl, xyzf_QMixNl, xyz_KmBl,    xyz_KmNl, xyz_CDensBl, xyz_CDensNl, pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, xyz_DPTempDtNl, xyzf_DQMixDtNl, xyz_DKmDtNl, 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(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, xqz_VelYNs, xyr_VelZNs, xyz_ExnerNs, pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, 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(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