subroutine turbulence_KW1978_forcing( pyz_VelXBl, xqz_VelYBl,   xyr_VelZBl, xyz_PTempBl, xyz_ExnerBl, xyzf_QMixBl, xyz_KmBl,    xyz_KhBl,    xyz_CDensBl, pyz_DVelXDt, xqz_DVelYDt,  xyr_DVelZDt, xyz_DPTempDt,xyz_DExnerDt, xyzf_DQMixDt, xyz_DKmDt,   xyz_DCDensDt, xyz_KmAl, xyz_KhAl )

    implicit none

    real(DP),intent(in) :: pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !水平速度
    real(DP),intent(in) :: xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !水平速度
    real(DP),intent(in) :: xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !鉛直速度
    real(DP),intent(in) :: xyz_PTempBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !温位
    real(DP),intent(in) :: xyz_ExnerBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !無次元圧力
    real(DP),intent(in) :: xyzf_QMixBl(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                                    !凝縮成分の混合比
    real(DP),intent(in) :: xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP),intent(in) :: xyz_KhBl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP),intent(in) :: xyz_CDensBl(imin:imax,jmin:jmax,kmin:kmax)

    real(DP),intent(inout):: pyz_DVelXDt(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP),intent(inout):: xqz_DVelYDt(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP),intent(inout):: xyr_DVelZDt(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP),intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax)

    real(DP),intent(inout):: xyz_DExnerDt(imin:imax,jmin:jmax,kmin:kmax)

    real(DP),intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax)

    real(DP),intent(inout):: xyz_DKmDt(imin:imax,jmin:jmax,kmin:kmax)

    real(DP),intent(inout) :: xyz_DCDensDt(imin:imax,jmin:jmax,kmin:kmax)
    
    real(DP),intent(out):: xyz_KmAl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP),intent(out):: xyz_KhAl(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流拡散係数
    real(DP)            :: xyz_Buoy(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数の
    real(DP)            :: xyz_BuoyT(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数の
    real(DP)            :: xyz_BuoyM(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数の
    real(DP)            :: xyz_Shear(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数の
    real(DP)            :: xyz_Diff(imin:imax,jmin:jmax,kmin:kmax)
                                                    !渦粘性係数の
    real(DP)            :: xyz_Disp(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流エネルギーの消散
    real(DP)            :: xyz_DispPI(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流エネルギーの消散
    real(DP)            :: xyz_DispHeat(imin:imax,jmin:jmax,kmin:kmax)
                                                    !乱流エネルギーの消散
    real(DP)            :: xyz_Turb(imin:imax,jmin:jmax,kmin:kmax)
                                                    !
    real(DP)            :: xyzf_Turb(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                                    !スカラー量の水平乱流拡散
    real(DP)            :: pyz_Turb(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xqz_Turb(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xyr_Turb(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: pyz_DVelXDt0(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP)            :: xqz_DVelYDt0(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP)            :: xyr_DVelZDt0(imin:imax,jmin:jmax,kmin:kmax)
                                                    !スカラー量の水平乱流拡散
    real(DP)            :: xyz_DPTempDt0(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xyz_DExnerDt0(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xyzf_DQMixDt0(imin:imax,jmin:jmax,kmin:kmax, ncmax)

    real(DP)            :: xyz_DKmDt0(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xyz_DCDensDt0(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)            :: xyz_PTempBlAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: xyzf_QMixBlAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    integer             :: f


    !----------------------------------
    ! 初期化
    !
    xyz_PTempBlAll = xyz_PTempBl + xyz_PTempBZ
    xyzf_QMixBlAll  = xyzf_QMixBl + xyzf_QMixBZ 
    pyz_DVelXDt0 = pyz_DVelXDt
    xqz_DVelYDt0 = xqz_DVelYDt
    xyr_DVelZDt0 = xyr_DVelZDt
    xyz_DPTempDt0 = xyz_DPTempDt
    xyz_DExnerDt0 = xyz_DExnerDt
    xyzf_DQMixDt0 = xyzf_DQMixDt
    xyz_DKmDt0    = xyz_DKmDt
    xyz_DCDensDt0 = xyz_DCDensDt
    
    !----------------------------------
    ! 拡散係数の時間発展 (エネルギー方程式を Km の式に変形したもの)
    !
    if ( .not. FlagConstTurbCoef ) then
    ! Buoyancy term
    !
    xyz_Buoy = xyz_BuoyMoistKm(xyz_PTempBl, xyz_ExnerBl, xyzf_QMixBl) 
    xyz_BuoyT = - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) * xyz_avr_xyr( xyr_dz_xyz( xyz_PTempBlAll ) ) / ( 2.0d0 * xyz_PTempBZ )
    xyz_BuoyM = xyz_Buoy - xyz_BuoyT


    xyz_Shear = ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) * ( ( xyz_dx_pyz( pyz_VelXBl ) ) ** 2.0d0 + ( xyz_dy_xqz( xqz_VelYBl ) ) ** 2.0d0 + ( xyz_dz_xyr( xyr_VelZBl ) ) ** 2.0d0 + 5.0d-1 * ( ( xyz_avr_pyr( pyr_dz_pyz( pyz_VelXBl ) ) + xyz_avr_pyr( pyr_dx_xyr( xyr_VelZBl ) ) ) ** 2.0d0 + ( xyz_avr_xqr( xqr_dy_xyr( xyr_VelZBl ) ) + xyz_avr_xqr( xqr_dz_xqz( xqz_VelYBl ) ) ) ** 2.0d0 + ( xyz_avr_pqz( pqz_dx_xqz( xqz_VelYBl ) ) + xyz_avr_pqz( pqz_dy_pyz( pyz_VelXBl ) ) ) ** 2.0d0 ) ) - xyz_KmBl * (  xyz_dx_pyz( pyz_VelXBl ) + xyz_dy_xqz( xqz_VelYBl ) + xyz_dz_xyr( xyr_VelZBl ) ) / 3.0d0

    xyz_Diff = 5.0d-1 * ( xyz_dx_pyz(pyz_dx_xyz(xyz_KmBl ** 2.0d0)) + xyz_dy_xqz(xqz_dy_xyz(xyz_KmBl ** 2.0d0)) + xyz_dz_xyr(xyr_dz_xyz(xyz_KmBl ** 2.0d0)) ) + ( (xyz_avr_pyz(pyz_dx_xyz(xyz_KmBl))) ** 2.0d0 + (xyz_avr_xqz(xqz_dy_xyz(xyz_KmBl))) ** 2.0d0 + (xyz_avr_xyr(xyr_dz_xyz(xyz_KmBl))) ** 2.0d0 )
    
    ! t - \Delta t で評価
    !
    xyz_Disp = - (xyz_KmBl ** 2.0d0) * 5.0d-1 / (MixLen ** 2.0d0)

    ! tendency
    !
    xyz_DKmDt = xyz_DKmDt0 + xyz_Buoy + xyz_Shear + xyz_Diff + xyz_Disp

    call HistoryAutoPut(TimeN, 'KmBuoyT', xyz_BuoyT(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'KmBuoyM', xyz_BuoyM(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'KmShear', xyz_Shear(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'KmDiff', xyz_Diff(1:nx,1:ny,1:nz))    
    call HistoryAutoPut(TimeN, 'KmDisp', xyz_Disp(1:nx,1:ny,1:nz))    
   
    ! 時間積分
    !
    xyz_KmAl = xyz_KmBl + (2.0d0 * DelTimeLong) * xyz_DKmDt

    ! 上限値の設定
    !
    xyz_KmAl = max( 0.0d0, min( xyz_KmAl, KmMax ) )

    ! Kh の設定
    !
    xyz_KhAl = 3.0d0 * xyz_KmAl

    else 
      xyz_KmAl = ConstKm
      xyz_KhAl = ConstKh
    end if
    
    !--------------------------------
    ! 雲密度の tendency
    !
    xyz_Turb = xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyz_CDensBl ) ) + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyz_CDensBl ) ) + xyz_dz_xyr( xyr_avr_xyz( xyz_DensBZ * xyz_KhBl ) * xyr_dz_xyz( xyz_CDensBl ) ) / xyz_DensBZ

    xyz_DCDensDt = xyz_DCDensDt0 + xyz_Turb

    ! output
    !
    call HistoryAutoPut(TimeN, 'CDensTurb', xyz_Turb(1:nx,1:ny,1:nz))    

    !--------------------------------
    ! 温位の tendency
    !
    xyz_Turb = xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyz_PTempBlAll ) ) + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyz_PTempBlAll ) ) + xyz_dz_xyr( xyr_avr_xyz( xyz_DensBZ * xyz_KhBl ) * xyr_dz_xyz( xyz_PTempBlAll ) ) / xyz_DensBZ
    
    xyz_DispHeat = (xyz_KmBl ** 3.0d0) / (xyz_ExnerBZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0))

    xyz_DPTempDt = xyz_DPTempDt0 + xyz_Turb + xyz_DispHeat

    call HistoryAutoPut(TimeN, 'PTempDisp', xyz_DispHeat(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'PTempTurb', xyz_Turb(1:nx, 1:ny, 1:nz))

    !--------------------------------
    ! 混合比の tendency
    !
    do f = 1, ncmax    
      xyzf_Turb(:,:,:,f) = xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyzf_QMixBlAll(:,:,:,f) ) ) + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyzf_QMixBlAll(:,:,:,f) ) ) + xyz_dz_xyr( xyr_avr_xyz( xyz_DensBZ * xyz_KhBl ) * xyr_dz_xyz( xyzf_QMixBlAll(:,:,:,f) ) ) / xyz_DensBZ
      
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_Turb', xyzf_Turb(1:nx,1:ny,1:nz,f))
    end do

    xyzf_DQMixDt = xyzf_DQMixDt0 + xyzf_Turb
    
    !--------------------------------
    ! 各速度成分の tendency
    !
    pyz_Turb = 2.0d0 * pyz_dx_xyz( xyz_KmBl * xyz_dx_pyz( pyz_VelXBl ) ) + pyz_dy_pqz( pqz_avr_xyz( xyz_KmBl ) * pqz_dx_xqz( xqz_VelYBl ) + pqz_avr_xyz( xyz_KmBl ) * pqz_dy_pyz( pyz_VelXBl ) ) + pyz_dz_pyr( pyr_avr_xyz( xyz_DensBZ * xyz_KmBl ) * pyr_dx_xyr( xyr_VelZBl ) + pyr_avr_xyz( xyz_DensBZ * xyz_KmBl ) * pyr_dz_pyz( pyz_VelXBl ) ) / xyz_DensBZ - 2.0d0 * pyz_dx_xyz( ( xyz_KmBl ** 2.0d0 ) ) / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )

    pyz_DVelXDt = pyz_DVelXDt0 + pyz_Turb

    call HistoryAutoPut(TimeN, 'VelXTurb', pyz_Turb(1:nx, 1:ny, 1:nz))

    xqz_Turb = 2.0d0 * xqz_dy_xyz( xyz_KmBl * xyz_dy_xqz( xqz_VelYBl ) ) + xqz_dx_pqz( pqz_avr_xyz( xyz_KmBl ) * pqz_dy_pyz( pyz_VelXBl ) + pqz_avr_xyz( xyz_KmBl ) * pqz_dx_xqz( xqz_VelYBl ) ) + xqz_dz_xqr( xqr_avr_xyz( xyz_DensBZ * xyz_KmBl ) * xqr_dy_xyr( xyr_VelZBl ) + xqr_avr_xyz( xyz_DensBZ * xyz_KmBl ) * xqr_dz_xqz( xqz_VelYBl ) ) / xyz_DensBZ - 2.0d0 * xqz_dy_xyz( ( xyz_KmBl ** 2.0d0 ) ) / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )
    
    xqz_DVelYDt = xqz_DVelYDt0 + xqz_Turb

    call HistoryAutoPut(TimeN, 'VelYTurb', xqz_Turb(1:nx, 1:ny, 1:nz))

    xyr_Turb = + 2.0d0 * xyr_dz_xyz( xyz_DensBZ * xyz_KmBl * xyz_dz_xyr( xyr_VelZBl ) ) / xyr_avr_xyz( xyz_DensBZ ) + xyr_dx_pyr( pyr_avr_xyz( xyz_KmBl ) * pyr_dz_pyz( pyz_VelXBl ) + pyr_avr_xyz( xyz_KmBl ) * pyr_dx_xyr( xyr_VelZBl ) ) + xyr_dy_xqr( xqr_avr_xyz( xyz_KmBl ) * xqr_dz_xqz( xqz_VelYBl ) + xqr_avr_xyz( xyz_KmBl ) * xqr_dy_xyr( xyr_VelZBl ) ) - 2.0d0 * xyr_dz_xyz( xyz_DensBZ * ( xyz_KmBl ** 2.0d0 ) ) / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) ) / xyr_avr_xyz( xyz_DensBZ ) 

    xyr_DVelZDt = xyr_DVelZDt0 + xyr_Turb

    call HistoryAutoPut(TimeN, 'VelZTurb', xyr_Turb(1:nx, 1:ny, 1:nz))

    !--------------------
    ! Exner function
    !
    if ( FlagDExnerDtTurb ) then
      xyz_DispPI = xyz_DExnerDt_xyz( xyz_DispHeat )
    else 
      xyz_DispPi = 0.0d0
    end if
    xyz_DExnerDt = xyz_DExnerDt0 + xyz_DispPI

    call HistoryAutoPut(TimeN, 'ExnerDisp', xyz_DispPI(1:nx, 1:ny, 1:nz))

    ! Set Margin
    !
    call SetMargin_xyz(xyz_KmAl)
    call SetMargin_xyz(xyz_KhAl)

  end subroutine Turbulence_KW1978_forcing

  subroutine turbulence_kw1978_init
    !
    ! Turbulence モジュールの初期化ルーチン
    ! 

    !暗黙の型宣言禁止
    implicit none

    real(DP):: MinDelX
    real(DP):: MinDelY
    real(DP):: MinDelZ
    integer :: l
    integer :: unit

    ! NAMELIST から情報を取得
    NAMELIST /turbulence_kw1978_nml/ Cm, KmMax, FlagConstTurbCoef, ConstKm, ConstKh, FlagDExnerDtTurb

    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=turbulence_kw1978_nml)
    close(unit)

    ! 混合距離
    ! 2 次元計算の場合には DelY に依存しないようにするために if 文を利用.
    ! 
    MinDelX = minval(x_dx)
    MinDelY = minval(y_dy)
    MinDelZ = minval(z_dz)

    if (jmin == jmax) then 
!       MixLen = min( MinDelZ, MinDelX)
       MixLen = sqrt(MinDelX * MinDelZ) 
    else
!       MixLen = min( MinDelZ, min( MinDelX, MinDelY ) )
       MixLen = ( MinDelZ * MinDelX * MinDelY ) ** (1.0 / 3.0)
    end if

    
    ! KmMax が設定されていない場合. 
    ! 安定性解析では, dt / l**2 < 0.5 を満たす必要がある. ここでは 0.1 にしておいた. 
    !
    if (KmMax == 0.0d0) then 
       KmMax = 0.1 * (MixLen ** 2.0d0) / (DelTimeLong * 2.0d0) !LeapFrog
    end if

    ! Output
    !
    if (myrank == 0) then 
      call MessageNotify( "M", module_name, "Cm = %f", d=(/Cm/))
      call MessageNotify( "M", module_name, "KmMax = %f", d=(/KmMax/))
      call MessageNotify( "M", module_name, "MixLen= %f", d=(/MixLen/))
      call MessageNotify( "M", module_name, "FlagConstTurbCoef= %b", l=(/ FlagConstTurbCoef /))
      call MessageNotify( "M", module_name, "ConstKm= %f", d=(/ConstKm/))
      call MessageNotify( "M", module_name, "ConstKh= %f", d=(/ConstKh/))
      call MessageNotify( "M", module_name, "FlagDExnerDtTurb= %b", l=(/ FlagDExnerDtTurb /))
    end if

    call HistoryAutoAddVariable( varname='VelXTurb', dims=(/'x','y','z','t'/), longname='Turbulence term of velocity (x)', units='K.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='VelYTurb', dims=(/'x','y','z','t'/), longname='Turbulence term of velocity (y)', units='K.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='VelZTurb', dims=(/'x','y','z','t'/), longname='Turbulence term of velocity (z)', units='K.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='PTempTurb', dims=(/'x','y','z','t'/), longname='Turbulence term of potential temperature', units='K.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='PTempDisp', dims=(/'x','y','z','t'/), longname='Dissipation term of potential temperature', units='K.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='ExnerDisp', dims=(/'x','y','z','t'/), longname='Dissipation term of exner function', units='K.s-1', xtype='float')

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

    do l = 1, ncmax
      call HistoryAutoAddVariable( varname=trim(SpcWetSymbol(l))//'_Turb', dims=(/'x','y','z','t'/), longname='Turbulence term of ' //trim(SpcWetSymbol(l))//' mixing ratio', units='kg.kg-1.s-1', xtype='float')
    end do

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

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

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

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

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


  end subroutine turbulence_kw1978_init