Class: Turbulence

モデルの物理過程を計算するために必要となる関数群を束ねたモジュール具体的には以下の項を計算するための関数を格納する.

 * 乱流拡散項
 * 乱流エネルギーの時間発展方程式に含まれる各項
 * 散逸加熱項

Methods

turbulence_KW1978_forcing turbulence_kw1978_init

Included Modules

dc_types dc_iounit dc_message gtool_historyauto mpi_wrapper timeset gridset basicset constants composition ChemCalc axesset xyz_deriv_c4_module xyz_deriv_module namelist_util setmargin DExnerDt

Public Instance methods

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 )

Subroutine :

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

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

xyz_KmAl(imin:imax,jmin:jmax,kmin:kmax) :

real(DP),intent(out)

乱流拡散係数

xyz_KhAl(imin:imax,jmin:jmax,kmin:kmax) :

real(DP),intent(out)

乱流拡散係数

[Source]

  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_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_Buoy = &
!        &  - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) &
!        &       * xyz_avr_xyr( xyr_dz_xyz( xyz_PTempBlAll ) ) &
!        &       / ( 2.0d0 * xyz_PTempBZ )
    
    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
    
    ! 時間積分
    !
    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_KhBl ) * xyr_dz_xyz( xyz_CDensBl ) )   

    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_KhBl ) * xyr_dz_xyz( xyz_PTempBlAll ) )  
    
    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_KhBl ) * xyr_dz_xyz( xyzf_QMixBlAll(:,:,:,f) ) )    
      
      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_KmBl ) * pyr_dx_xyr( xyr_VelZBl ) + pyr_avr_xyz( xyz_KmBl ) * pyr_dz_pyz( pyz_VelXBl ) ) - 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_KmBl ) * xqr_dy_xyr( xyr_VelZBl ) + xqr_avr_xyz( xyz_KmBl ) * xqr_dz_xqz( xqz_VelYBl ) ) - 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_KmBl * xyz_dz_xyr( xyr_VelZBl ) ) + 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_KmBl ** 2.0d0 ) ) / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )

    xyr_DVelZDt = xyr_DVelZDt0 + xyr_Turb

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

    !--------------------
    ! Exner function
    !
    if ( .not. 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

turbulence_kw1978_init( )

Subroutine :

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

This procedure input/output NAMELIST#turbulence_kw1978_nml .

[Source]

  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)

    !混合距離
    MinDelX = minval(x_dx)
    MinDelY = minval(y_dy)
    MinDelZ = minval(z_dz)

    MixLen = min( MinDelZ, min( MinDelX, MinDelY ) )

    ! 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

  end subroutine turbulence_kw1978_init