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)  :: AlphaNDiff2 = 0.0d0  !2次の数値拡散の係数. CReSS マニュアルより
    real(DP)  :: NDiffRatio = 1.0d0   !速度に対する粘性を上げる場合は数字を 1 以上にする. 
    integer   :: unit            !装置番号
    integer   :: l

    NAMELIST /Dynamics_nml/ AlphaSound, AlphaNDiff, AlphaNDiff2, NDiffRatio, beta, FactorBuoyTemp, FactorBuoyMolWt, FactorBuoyLoading

    !ファイルオープン. 情報取得. 
    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=dynamics_nml)
    close(unit)
   
    !-------------------------------------------------------------------
    ! 音波減衰項の減衰率   Min(DelX, DelZ) ** 2.0 に比例
    ! 2 次元計算の場合には DelY に依存しないようにするために if 文を利用.
    !
    DelXMin = minval(x_dx)
    DelYMin = minval(y_dy)
    DelZMin = minval(z_dz)
    if (jmin == jmax) then 
       Alpha = AlphaSound * ( Min(DelXMin, DelZMin) ** 2.0d0 ) / DelTimeShort
    else
       Alpha = AlphaSound * ( Min(DelXMin, DelYMin, DelZMin) ** 2.0d0 ) / DelTimeShort
    end if

    if (jmin == jmax) then     
       NuHh4 = AlphaNDiff * ( DelXMin ** 4.0d0 ) / (2.0d0 * DelTimeLong)
    else
       NuHh4 = AlphaNDiff * ( SQRT(DelXMin*DelYMin) ** 4.0d0 ) / (2.0d0 * DelTimeLong)
    end if
    NuVh4 = AlphaNDiff * ( DelZMin ** 4.0d0 ) / (2.0d0 * DelTimeLong)
    NuHm4 = NuHh4 * NDiffRatio
    NuVm4 = NuVh4 * NDiffRatio

    if (jmin == jmax) then     
       NuHh2 = AlphaNDiff2 * ( DelXMin ** 2.0d0 ) / (2.0d0 * DelTimeLong)
    else
       NuHh2 = AlphaNDiff2 * ( SQRT(DelXMin*DelYMin) ** 2.0d0 ) / (2.0d0 * DelTimeLong)
    end if
    NuVh2 = AlphaNDiff2 * ( DelZMin ** 2.0d0 ) / (2.0d0 * DelTimeLong)
    NuHm2 = NuHh2 * NDiffRatio
    NuVm2 = NuVh2 * NDiffRatio

    if (myrank == 0) then 
      if ( AlphaNDiff > 0.0d0 .AND. AlphaNDiff2 > 0.0d0 ) then 
        write(*,*) "CHECK! AlphaNDiff == 0.0d0 .OR. AlphaNDiff2 == 0.0d0"
        stop 
      end if
      call MessageNotify( "M", module_name, "Alpha = %f", d=(/Alpha/) )
      call MessageNotify( "M", module_name, "NuHh = %f", d=(/NuHh4/) )
      call MessageNotify( "M", module_name, "NuVh = %f", d=(/NuVh4/) )
      call MessageNotify( "M", module_name, "NuHm = %f", d=(/NuHm4/) )
      call MessageNotify( "M", module_name, "NuVm = %f", d=(/NuVm4/) )
      call MessageNotify( "M", module_name, "NuHh = %f", d=(/NuHh2/) )
      call MessageNotify( "M", module_name, "NuVh = %f", d=(/NuVh2/) )
      call MessageNotify( "M", module_name, "NuHm = %f", d=(/NuHm2/) )
      call MessageNotify( "M", module_name, "NuVm = %f", d=(/NuVm2/) )
    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='PTempNDiff2', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of potential temperature (2 order)', 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='ExnerNDiff2', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of exner function (2 order)', 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')

    call HistoryAutoAddVariable( varname='CDensNDiff2', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of cloud density (2 order)', 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')

      call HistoryAutoAddVariable( varname=trim(SpcWetSymbol(l))//'_NDiff2', dims=(/'x','y','z','t'/), longname='Diffusion term of ' //trim(SpcWetSymbol(l))//' mixing ratio (2 order)', 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='VelXNDiff2', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of velocity (x) (2 order)', 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='VelYNDiff2', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of velocity (y) (2 order)', 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='VelZNDiff2', dims=(/'x','y','z','t'/), longname='Numerical diffusion term of Velocity (z) (2 order)', 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='KmAdv', dims=(/'x','y','z','t'/), longname='Advection of Km', units='s-1', xtype='float')

    call HistoryAutoAddVariable( varname='KmNDiff', dims=(/'x','y','z','t'/), longname='Diffusion term of Km', units='s-1', xtype='float')

    call HistoryAutoAddVariable( varname='KmNDiff2', dims=(/'x','y','z','t'/), longname='Diffusion term of Km', units='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_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_NDiff2(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_NDiff4 = - NuHm4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyz_KmBl ))))) - NuHm4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_KmBl ))))) - NuVm4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyz_KmBl ))))) 

    xyz_NDiff2 = + NuHm2 * (xyz_dx_pyz(pyz_dx_xyz( xyz_KmBl ))) + NuHm2 * (xyz_dy_xqz(xqz_dy_xyz( xyz_KmBl ))) + NuVm2 * (xyz_dz_xyr(xyr_dz_xyz( xyz_KmBl ))) 
    
    xyz_DKmDtNl = xyz_Orig + xyz_Adv + xyz_NDiff4 + xyz_NDiff2

    call HistoryAutoPut(TimeN, 'KmAdv',    xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'KmNDiff',  xyz_NDiff4(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'KmNDiff2', xyz_NDiff2(1:nx,1:ny,1:nz))
    
    ! 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_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: pyz_NDiff2(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_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xqz_NDiff2(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_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyr_NDiff2(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_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_NDiff2(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_NDiff4(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax)
    real(DP)             :: xyzf_NDiff2(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_NDiff4 = - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_CDensBl))))) - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_CDensBl))))) - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_CDensBl))))) 

    xyz_NDiff2 = + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz(xyz_CDensBl))) + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz(xyz_CDensBl))) + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz(xyz_CDensBl))) 
    
    xyz_DCDensDtNl = xyz_Orig + xyz_Adv + xyz_NDiff4 + xyz_NDiff2

    call HistoryAutoPut(TimeN, 'CDensAdv',  xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'CDensNDiff', xyz_NDiff4(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'CDensNDiff2', xyz_NDiff2(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_NDiff4 = - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyz_PTempBl ))))) - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_PTempBl ))))) - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyz_PTempBl ))))) 

    xyz_NDiff2 = + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz( xyz_PTempBl ))) + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz( xyz_PTempBl ))) + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz( xyz_PTempBl ))) 

    ! sum
    !
    xyz_DPTempDtNl = xyz_Orig + xyz_Adv + xyz_NDiff4 + xyz_NDiff2
    
    ! output
    !
    call HistoryAutoPut(TimeN, 'PTempAdv',   xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'PTempNDiff', xyz_NDiff4(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'PTempNDiff2', xyz_NDiff2(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_NDiff4(:,:,:,f) = - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyzf_QMixBl(:,:,:,f) ))))) - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyzf_QMixBl(:,:,:,f) ))))) - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyzf_QMixBl(:,:,:,f) )))))

      xyzf_NDiff2(:,:,:,f) = + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz( xyzf_QMixBl(:,:,:,f) ))) + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz( xyzf_QMixBl(:,:,:,f) ))) + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz( xyzf_QMixBl(:,:,:,f) )))
    end do

    ! sum
    !
    xyzf_DQMixDtNl = xyzf_Orig + xyzf_Adv + xyzf_NDiff4 + xyzf_NDiff2

    ! 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_NDiff4(1:nx,1:ny,1:nz,f))
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_NDiff2', xyzf_NDiff2(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_NDiff4 = - NuHm4 * ( pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz( pyz_VelXBl ))))) - NuHm4 * ( pyz_dy_pqz(pqz_dy_pyz(pyz_dy_pqz(pqz_dy_pyz( pyz_VelXBl ))))) - NuVm4 * ( pyz_dz_pyr(pyr_dz_pyz(pyz_dz_pyr(pyr_dz_pyz( pyz_VelXBl )))))

    pyz_NDiff2 = + NuHm2 * ( pyz_dx_xyz(xyz_dx_pyz( pyz_VelXBl ))) + NuHm2 * ( pyz_dy_pqz(pqz_dy_pyz( pyz_VelXBl ))) + NuVm2 * ( pyz_dz_pyr(pyr_dz_pyz( pyz_VelXBl )))

    ! sum
    !
    pyz_DVelXDtNl = pyz_Orig + pyz_Adv + pyz_NDiff4 + pyz_NDiff2

    call HistoryAutoPut(TimeN, 'VelXAdv',   pyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXNDiff', pyz_NDiff4(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXNDiff2', pyz_NDiff2(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_NDiff4 = - NuHm4 * ( xqz_dx_pqz(pqz_dx_xqz(xqz_dx_pqz(pqz_dx_xqz( xqz_VelYBl ))))) - NuHm4 * ( xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz( xqz_VelYBl ))))) - NuVm4 * ( xqz_dz_xqr(xqr_dz_xqz(xqz_dz_xqr(xqr_dz_xqz( xqz_VelYBl )))))

    xqz_NDiff2 = + NuHm2 * ( xqz_dx_pqz(pqz_dx_xqz( xqz_VelYBl ))) + NuHm2 * ( xqz_dy_xyz(xyz_dy_xqz( xqz_VelYBl ))) + NuVm2 * ( xqz_dz_xqr(xqr_dz_xqz( xqz_VelYBl )))

    ! sum
    !
    xqz_DVelYDtNl = xqz_Orig + xqz_Adv + xqz_NDiff4 + xqz_NDiff2
    
    call HistoryAutoPut(TimeN, 'VelYAdv',   xqz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYNDiff', xqz_NDiff4(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYNDiff2', xqz_NDiff2(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_NDiff4 = - NuHm4 * ( xyr_dx_pyr(pyr_dx_xyr(xyr_dx_pyr(pyr_dx_xyr( xyr_VelZBl ))))) - NuHm4 * ( xyr_dy_xqr(xqr_dy_xyr(xyr_dy_xqr(xqr_dy_xyr( xyr_VelZBl ))))) - NuVm4 * ( xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr( xyr_VelZBl )))))

    xyr_NDiff2 = + NuHm2 * ( xyr_dx_pyr(pyr_dx_xyr( xyr_VelZBl ))) + NuHm2 * ( xyr_dy_xqr(xqr_dy_xyr( xyr_VelZBl ))) + NuVm2 * ( 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_NDiff4 + xyr_NDiff2 + xyr_BuoyT * FactorBuoyTemp + xyr_BuoyM * FactorBuoyMolWt + xyr_BuoyD * FactorBuoyLoading

    call HistoryAutoPut(TimeN, 'VelZAdv',   xyr_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelZNDiff', xyr_NDiff4(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'VelZNDiff2', xyr_NDiff2(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_NDiff4(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_NDiff2(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   =   Alpha * 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   =   Alpha * 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_NDiff4 = - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_ExnerNs))))) - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_ExnerNs))))) - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_ExnerNs))))) 

    xyz_NDiff2 = + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz(xyz_ExnerNs))) + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz(xyz_ExnerNs))) + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz(xyz_ExnerNs))) 

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

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

    call HistoryAutoPut(TimeN, 'ExnerAdv',   xyz_Adv(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerNDiff', xyz_NDiff4(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerNDiff2', xyz_NDiff2(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 =  Alpha * 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 

    ! 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)

  end subroutine Dynamics_Short_forcing