!= Module Turbulence_tendency
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA, Masatsugu
! Version::   $Id: turbulence_constKm.f90,v 1.22 2013-01-31 05:08:33 odakker Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
! 乱流拡散係数を一定にしたとき
!
!== Error Handling
!
!== Bugs
!
!== Note
!
!  * 定式化は CReSS のマニュアルを参照した
!
!== Future Plans
!
!

module Turbulence_utils
  !
  !  * 乱流拡散項
  !

  !モジュール読み込み 
  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 timeset, only:  DelTimeLong, TimeN
  use gridset, only:  imin,           &! x 方向の配列の下限
    &                 imax,           &! x 方向の配列の上限
    &                 jmin,           &! y 方向の配列の下限
    &                 jmax,           &! y 方向の配列の上限
    &                 kmin,           &! z 方向の配列の下限
    &                 kmax,           &! z 方向の配列の上限
    &                 nx,ny,nz,ncmax,FlagCalc3D  
  use basicset, only: xyz_PTempBZ,    &!基本場の温位
    &                 xyzf_QMixBZ,    &!基本場の混合比
    &                 xyz_ExnerBZ,     &!基本場のエクスナー関数
    &                 xyz_DensBZ       !基本場の密度
  use constants,only: CpDry
  use composition, only: SpcWetSymbol
  use axesset, only: &  
    &                 x_dx,           &! x 方向の格子点間隔
    &                 y_dy,           &! y 方向の格子点間隔
    &                 z_dz,           &! z 方向の格子点間隔
    &                 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, &
    &                 pqz_avr_xyz, pyr_avr_xyz, xqr_avr_xyz, &
    &                 xyz_avr_pqz, xyz_avr_pyr, xyz_avr_xqr
  use xyz_deriv_module, only: xyz_dx_pyz, xyz_dy_xqz, xyz_dz_xyr, &
    &                         pyz_dx_xyz, xqz_dy_xyz, xyr_dz_xyz, &
    &                         pqz_dx_xqz, pqz_dy_pyz, pyz_dy_pqz, &
    &                         pyr_dx_xyr, pyz_dz_pyr, pyr_dz_pyz, &
    &                         pyr_dz_pyz, pyr_dx_xyr, xyr_dx_pyr, &
    &                         xqr_dz_xqz, xqr_dy_xyr, xyr_dy_xqr, &
    &                         pqz_dy_pyz, pqz_dx_xqz, xqz_dx_pqz, &
    &                         xqr_dy_xyr, xqr_dz_xqz, xqz_dz_xqr  
  use namelist_util, only: namelist_filename
  use DExnerDt,  only: xyz_DExnerDt_xyz

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private
  
  !関数を public に設定
  public Turbulence_constKm_Init
  public Turbulence_constKm_Forcing

  !変数定義
  real(DP), save :: Cm     = 2.0d-1    !乱流エネルギー診断式の係数 
  real(DP), save :: MixLen = 0.0d0     !平均混合距離
  real(DP), save :: ConstKm = 0.0d0    !定数値
  real(DP), save :: ConstKh = 0.0d0    !定数値
  logical,  save :: FlagDispHeat = .false. 
  logical,  save :: FlagDExnerDtTurb =.true.

  character(*), parameter:: module_name = 'turbulence_constKm'

contains

!!!------------------------------------------------------------------------!!!

  subroutine turbulence_utils_tendency(        &
    & xyz_DispHeat,                            &
    & 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,               &
    )

    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

    
    !--------------------------------
    ! 雲密度の 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

    if ( FlagDispHeat ) then 
      xyz_DispHeat = (xyz_KmBl ** 3.0d0) &
        & / (xyz_ExnerBZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0))
    else 
      xyz_DispHeat = 0.0d0
    end if

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

  end subroutine Turbulence_tendency
  
end module Turbulence_utils
