!= 鉛直拡散フラックス (Mellor and Yamada, 1974, 1982) ! != Vertical diffusion flux (Mellor and Yamada, 1974, 1982) ! ! Authors:: Yasuhiro MORIKAWA, Yukiko YAMADA, Yoshiyuki O. Takahashi, Shin-ichi Takehiro ! Version:: $Id: vdiffusion_my.f90,v 1.5 2025/09/17 22:00:00 takepiro Exp $ ! Tag Name:: $Name: $ ! Copyright:: Copyright (C) GFD Dennou Club, 2008-2025. All rights reserved. ! License:: See COPYRIGHT[link:../../../COPYRIGHT] ! module vdiffusion_my ! != 鉛直拡散フラックス (Mellor and Yamada, 1974, 1982) ! != Vertical diffusion flux (Mellor and Yamada, 1974, 1982) ! ! Note that Japanese and English are described in parallel. ! ! 鉛直拡散フラックスを計算します. ! ! Vertical diffusion flux is calculated. ! !== Procedures List ! ! VDiffusion :: 鉛直拡散フラックスの計算 ! VDiffusionOutPut :: フラックスの出力 ! ------------ :: ------------ ! VDiffusion :: Calculate vertical diffusion fluxes ! VDiffusionOutPut :: Output fluxes ! !== NAMELIST ! ! NAMELIST#vdiffusion_my_nml ! ! モジュール引用 ; USE statements ! !== References ! ! * Mellor, G. L., and T. Yamada, 1974: ! A hierarchy of turbulence closure models for ! planetary boundary layers. ! J. Atmos. Sci., 31, 1791--1806. ! * Mellor, G. L., and T. Yamada, 1982: ! Development of a turbulence closure model for geophysical fluid problems, ! Rev. Geophys. Space Phys., 20, 851--875. ! * Galperin, B., L. H. Kantha, S. Hassid, and T. Rosati, 1988: ! A quasi-equilibrium turbulent energy model for geophysical flows, ! J. Atmos. Sci., 45, 55--62. ! 格子点設定 ! Grid points settings ! use gridset, only: imax, & ! 経度格子点数. ! Number of grid points in longitude & jmax, & ! 緯度格子点数. ! Number of grid points in latitude & kmax ! 鉛直層数. ! Number of vertical level ! 組成に関わる配列の設定 ! Settings of array for atmospheric composition ! use composition, only: ncmax, IndexH2OVap ! 種別型パラメタ ! Kind type parameter ! use dc_types, only: DP, & ! 倍精度実数型. Double precision. & STRING ! 文字列. Strings. ! メッセージ出力 ! Message output ! use dc_message, only: MessageNotify ! 宣言文 ; Declaration statements ! implicit none private ! 公開手続き ! Public procedure ! public :: VDiffusionMY25 public :: VDiffusionMY251DWrapper3D public :: VDiffusionMY25GBT94 public :: VDiffusion public :: VDiffusionExpTendency public :: VDiffusionOutPut public :: VDiffusionInit ! 公開変数 ! Public variables ! ! 非公開変数 ! Private variables ! logical, save :: vdiffusion_my_inited = .false. ! 初期設定フラグ. ! Initialization flag ! 鉛直拡散フラックス算出のための変数 ! Variables for calculation of vertical diffusion flux ! logical , save:: FlagConstDiffCoef ! Flag for use of constant diffusion coefficient real(DP), save:: ConstDiffCoefM ! Diffusion coefficient for momentum that is used in the ! case of FlagConstDiffCoef == .true. real(DP), save:: ConstDiffCoefH ! Diffusion coefficient for heat that is used in the ! case of FlagConstDiffCoef == .true. real(DP), save:: SquareVelMin ! 風二乗差最小値. ! Minimum value of square of velocity real(DP), save:: BulkRiNumMin ! バルク $ R_i $ 数最小値. ! Minimum value of bulk $ R_i $ ! 拡散係数算出のための変数 ! Veriables for calculation of diffusion coefficients ! real(DP), save:: MixLengthMax ! 最大混合距離. ! Maximum mixing length real(DP), save:: ShMin ! $ S_h $ 最小値. ! Minimum $ S_h $ real(DP), save:: SmMin ! $ S_m $ 最小値. ! Minimum $ S_m $ real(DP), save:: VelDiffCoefMin ! $ \Dvect{u} $ 拡散係数最小値. ! Minimum diffusion coefficient of $ \Dvect{u} $ real(DP), save:: TempDiffCoefMin ! $ T $ 拡散係数最小値. ! Minimum diffusion coefficient of $ T $ real(DP), save:: VelDiffCoefMax ! $ \Dvect{u} $ 拡散係数最大値. ! Maximum diffusion coefficient of $ \Dvect{u} $ real(DP), save:: TempDiffCoefMax ! $ T $ 拡散係数最大値. ! Maximum diffusion coefficient of $ T $ ! Mellor Yamada Level 2 定数 ! Constants for Mellor Yamada Level 2 ! real(DP), save:: MYConstA1 real(DP), save:: MYConstB1 real(DP), save:: MYConstA2 real(DP), save:: MYConstB2 real(DP), save:: MYConstC1 character(*), parameter:: module_name = 'vdiffusion_my' ! モジュールの名称. ! Module name character(*), parameter:: version = & & '$Name: $' // & & '$Id: vdiffusion_my.f90,v 1.5 2025/09/17 22:00:00 takepiro Exp $' ! モジュールのバージョン ! Module version contains !-------------------------------------------------------------------------------------- subroutine VDiffusionMY25( & & xyz_U, xyz_V, xyzf_QMix, & ! (in) & xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, & ! (in) & xy_SurfHeight, & ! (in) & xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, & ! (in) & xyz_TurKinEne, & ! (in) & xy_SurfMomFluxX, xy_SurfMomFluxY, & ! (in) & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (out) & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, & ! (in) & xyz_DTurKinEneDt & ! (out) & ) ! ! 鉛直拡散フラックスを計算します. ! ! Vertical diffusion flux is calculated by use of MY2.5 model. ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: & & FKarm, & ! $ k $ . ! カルマン定数. ! Karman constant & Grav, & ! $ g $ [m s-2]. ! 重力加速度. ! Gravitational acceleration & GasRDry, & ! $ R $ [J kg-1 K-1]. ! 乾燥大気の気体定数. ! Gas constant of air & CpDry ! $ C_p $ [J kg-1 K-1]. ! 乾燥大気の定圧比熱. ! Specific heat of air at constant pressure ! 時刻管理 ! Time control ! use timeset, only: & & DelTime, & ! $ \Delta t $ [s] & TimeN, & ! ステップ $ t $ の時刻. Time of step $ t $. & TimesetClockStart, TimesetClockStop ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! MPI 関連ルーチン ! MPI related routines ! use mpi_wrapper, only : MPIWrapperChkTrue ! 陰解法による時間積分のためのルーチン ! Routines for time integration with implicit scheme ! use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3 ! 宣言文 ; Declaration statements ! implicit none real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax) ! $ u $ . 東西風速. Eastward wind real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax) ! $ v $ . 南北風速. Northward wind real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) ! $ q $ . 質量混合比. Mass mixing ratio real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax) ! $ T $ . 温度. Temperature real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{T} $ . 温度 (半整数レベル). ! Temperature (half level) real(DP), intent(in):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax) ! $ T_v $ . 仮温度. Virtual temperature real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{T}_v $ . 仮温度 (半整数レベル). ! Virtual temperature (half level) real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{p} $ . 気圧 (半整数レベル). ! Air pressure (half level) real(DP), intent(in):: xy_SurfHeight (0:imax-1,1:jmax) ! $ z_s $ . 地表面高度. ! Surface height. real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax) ! 高度 (整数レベル). ! Height (full level) real(DP), intent(in):: xyr_Height (0:imax-1, 1:jmax, 0:kmax) ! 高度 (半整数レベル). ! Height (half level) real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax) ! Exner 関数 (整数レベル). ! Exner function (full level) real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax) ! Exner 関数 (半整数レベル). ! Exner function (half level) real(DP), intent(in):: xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax) ! ! Turbulent kinetic energy (m2 s-2) real(DP), intent(in):: xy_SurfMomFluxX (0:imax-1, 1:jmax) ! ! Eastward momentum flux at surface real(DP), intent(in):: xy_SurfMomFluxY (0:imax-1, 1:jmax) ! ! Northward momentum flux at surface real(DP), intent(out):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax) ! 東西方向運動量フラックス. ! Eastward momentum flux real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax) ! 南北方向運動量フラックス. ! Northward momentum flux real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax) ! 熱フラックス. ! Heat flux real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) ! 質量フラックス. ! Mass flux of compositions real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：温度. ! Diffusion coefficient: temperature real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：比湿. ! Diffusion coefficient: specific humidity real(DP), intent(out):: xyz_DTurKinEneDt (0:imax-1, 1:jmax, 1:kmax) ! ! Tendency of turbulent kinetic energy ! 作業変数 ! Work variables ! real(DP) :: xyz_MixLength(0:imax-1, 1:jmax, 1:kmax) ! 混合距離. ! Mixing length real(DP) :: xyz_DVelDzSq(0:imax-1, 1:jmax, 1:kmax) ! ! Vertical shear squared (s-2) real(DP) :: xyz_StatStab(0:imax-1, 1:jmax, 1:kmax) ! ! Static stability (s-2) real(DP) :: GhMin ! ! Minimum of G_h real(DP) :: GhMax ! ! Maximum of G_h real(DP) :: xyz_Gm(0:imax-1, 1:jmax, 1:kmax) ! ! G_m real(DP) :: xyz_Gh(0:imax-1, 1:jmax, 1:kmax) ! ! G_h !!$ real(DP) :: xyz_SmHat(0:imax-1, 1:jmax, 1:kmax) !!$ ! !!$ ! \hat{S}_M !!$ real(DP) :: xyz_ShHat(0:imax-1, 1:jmax, 1:kmax) !!$ ! !!$ ! \hat{S}_h !!$ real(DP) :: xyz_DSmHatDTKE(0:imax-1, 1:jmax, 1:kmax) !!$ ! !!$ ! derivative of \hat{S}_M !!$ real(DP) :: xyz_DShHatDTKE(0:imax-1, 1:jmax, 1:kmax) !!$ ! !!$ ! derivative of \hat{S}_h real(DP) :: xyz_Sm(0:imax-1, 1:jmax, 1:kmax) ! ! S_M real(DP) :: xyz_Sh(0:imax-1, 1:jmax, 1:kmax) ! ! S_h real(DP), parameter :: Stke = 0.2_DP ! ! S_{TKE} = 0.2 real(DP) :: xyz_VelDiffCoef (0:imax-1, 1:jmax, 1:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP) :: xyz_TempDiffCoef (0:imax-1, 1:jmax, 1:kmax) ! 拡散係数：温度. ! Transfer coefficient: temperature real(DP) :: xyr_TurKinEneDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! ! Diffusion coefficient: turbulent kinetic energy real(DP) :: xyz_TurKinEneDiffCoef (0:imax-1, 1:jmax, 1:kmax) ! ! Diffusion coefficient: turbulent kinetic energy real(DP) :: xyr_TurKinEneTransCoef(0:imax-1, 1:jmax, 0:kmax) ! ! Transfer coefficient: turbulent kinetic energy real(DP) :: xyr_TurKinEneFlux(0:imax-1, 1:jmax, 0:kmax) ! ! Turbulent energy flux real(DP) :: xyz_CShe1(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_CShe2(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_CBuo1(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_CBuo2(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_CDis1(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_CDis2(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_TurKinEneProShear(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_TurKinEneProBuoya(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_FricVelSq (0:imax-1, 1:jmax) real(DP) :: xy_TurKinEneAtLB(0:imax-1, 1:jmax) real(DP) :: xyza_TurKinEneMtx(0:imax-1, 1:jmax, 1:kmax, -1:1) ! ! Implicit matrix for turbulent kinetic energy real(DP) :: xyz_TurKinEneVec(0:imax-1, 1:jmax, 1:kmax) ! ! Implicit vector for turbulent kinetic energy real(DP) :: xyza_TurKinEneLUMtx (0:imax-1, 1:jmax, 1:kmax, -1:1) ! LU 行列. ! LU matrix real(DP) :: xyz_DelTurKinEneLUVec(0:imax-1, 1:jmax, 1:kmax) ! ! Tendency of turbulent kinetic energy real(DP) :: xyz_TurKinEneDiss(0:imax-1, 1:jmax, 1:kmax) ! ! Dissipation rate of turbulent kinetic energy (m2 s-3) real(DP) :: xyz_TurKinEneNonZero(0:imax-1, 1:jmax, 1:kmax) ! ! Turbulent kinetic energy with offset (m2 s-2) real(DP), parameter :: TurKinEneOffset = ( 1.0e-3_DP )**2 / 2.0_DP logical :: FlagReCalc ! ! Flag for recalculation logical :: a_FlagReCalcLocal (1) logical :: a_FlagReCalcGlobal(1) integer :: iloop integer :: nloop integer:: i ! 経度方向に回る DO ループ用作業変数 ! Work variables for DO loop in longitude integer:: j ! 緯度方向に回る DO ループ用作業変数 ! Work variables for DO loop in latitude integer:: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer:: n ! 組成方向に回る DO ループ用作業変数 ! Work variables for DO loop in dimension of constituents ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. vdiffusion_my_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 計算時間計測開始 ! Start measurement of computation time ! call TimesetClockStart( module_name ) ! Calculate turbulent kinetic energy with offset ! xyz_TurKinEneNonZero = xyz_TurKinEne + TurKinEneOffset ! ! Calculation of vertical shear squared do k = 1, kmax if ( k == 1 ) then xyz_DVelDzSq(:,:,k) = & & ( ( xyz_U(:,:,k+1) - xyz_U(:,:,k ) )**2 & & + ( xyz_V(:,:,k+1) - xyz_V(:,:,k ) )**2 ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k ) )**2 !!$ xyz_DVelDzSq(:,:,k) = & !!$ & ( ( xyz_U(:,:,k+1) - 0.0_DP )**2 & !!$ & + ( xyz_V(:,:,k+1) - 0.0_DP )**2 ) & !!$ & / ( xyz_Height(:,:,k+1) - xy_SurfHeight )**2 else if ( k == kmax ) then xyz_DVelDzSq(:,:,k) = & & ( ( xyz_U(:,:,k ) - xyz_U(:,:,k-1) )**2 & & + ( xyz_V(:,:,k ) - xyz_V(:,:,k-1) )**2 ) & & / ( xyz_Height(:,:,k ) - xyz_Height(:,:,k-1) )**2 else xyz_DVelDzSq(:,:,k) = & & ( ( xyz_U(:,:,k+1) - xyz_U(:,:,k-1) )**2 & & + ( xyz_V(:,:,k+1) - xyz_V(:,:,k-1) )**2 ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) )**2 end if end do ! Calculation of static stability do k = 1, kmax if ( k == 1 ) then xyz_StatStab(:,:,k) = & & Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) & & * ( xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1) & & - xyz_VirTemp(:,:,k ) / xyz_Exner(:,:,k ) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k ) ) else if ( k == kmax ) then xyz_StatStab(:,:,k) = & & Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) & & * ( xyz_VirTemp(:,:,k ) / xyz_Exner(:,:,k ) & & - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) & & / ( xyz_Height(:,:,k ) - xyz_Height(:,:,k-1) ) else xyz_StatStab(:,:,k) = & & Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) & & * ( xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1) & & - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) ) end if end do ! 混合距離の算出 ! Calculate mixing length ! do k = 1, kmax xyz_MixLength(:,:,k) = & & FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) & & / (1.0_DP + FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) / MixLengthMax ) end do ! Limit mixing length (Galperin et al., 1988) and avoid zero xyz_MixLength = & & min( xyz_MixLength, & & 0.53_DP & & * sqrt( 2.0_DP * xyz_TurKinEneNonZero / max( xyz_StatStab, 1.0e-10_DP ) ) ) & & + 1.0e-10_DP xyz_Gh = - xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) & & * xyz_StatStab ! Actually, xyz_Gm is not used below. xyz_Gm = xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) & & * xyz_DVelDzSq ! Limit Gh (Galperin et al., 1988) GhMin = - 0.53_DP**2 GhMax = 1.0_DP & & / ( MYConstA2 & & * ( 12.0_DP * MYConstA1 + MYConstB1 + 3.0_DP * MYConstB2 ) ) xyz_Gh = max( GhMin, min( xyz_Gh, GhMax ) ) xyz_Sh = & & MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) & & / ( 1.0_DP & & - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * xyz_Gh ) xyz_Sm = & & ( & & MYConstA1 * ( 1.0_DP - 3.0_DP * MYConstC1 & & - 6.0_DP * MYConstA1 / MYConstB1 ) & & + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) & & * xyz_Gh * xyz_Sh & & ) & & / ( 1.0_DP & & - 9.0_DP * MYConstA1 * MYConstA2 * xyz_Gh ) !!$ xyz_DShHatDTKE = & !!$ & - 2.0_DP * MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) & !!$ & / ( 2.0_DP * xyz_TurKinEneNonZero & !!$ & - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * xyz_GhPrime )**2 !!$ xyz_DSmHatDTKE = & !!$ & ( & !!$ & 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) & !!$ & * xyz_GhPrime * xyz_DShHatDTKE & !!$ & * ( 2.0_DP * xyz_TurKinEneNonZero & !!$ & - 9.0_DP * MYConstA1 * MYConstA2 * xyz_GhPrime ) & !!$ & - 2.0_DP & !!$ & * ( & !!$ & MYConstA1 * ( 1.0_DP - 3.0_DP * MYConstC1 & !!$ & - 6.0_DP * MYConstA1 / MYConstB1 ) & !!$ & + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) & !!$ & * xyz_GhPrime * xyz_ShHat & !!$ & ) & !!$ & ) & !!$ & / ( 2.0_DP * xyz_TurKinEneNonZero & !!$ & - 9.0_DP * MYConstA1 * MYConstA2 * xyz_GhPrime )**2 ! 拡散係数の計算 ! Calculation of diffusion coefficient ! xyz_VelDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sm xyz_TempDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sh ! do k = 0, kmax if ( ( k == 0 ) .or. ( k == kmax ) ) then xyr_VelDiffCoef (:,:,k) = 0.0_DP xyr_TempDiffCoef(:,:,k) = 0.0_DP else xyr_VelDiffCoef (:,:,k) = ( xyz_VelDiffCoef (:,:,k) + xyz_VelDiffCoef (:,:,k+1) ) / 2.0_DP xyr_TempDiffCoef(:,:,k) = ( xyz_TempDiffCoef(:,:,k) + xyz_TempDiffCoef(:,:,k+1) ) / 2.0_DP end if end do ! do k = 1, kmax-1 do j = 1, jmax do i = 0, imax-1 xyr_VelDiffCoef(i,j,k) = & & max( min( xyr_VelDiffCoef(i,j,k), VelDiffCoefMax ), & & VelDiffCoefMin & & ) xyr_TempDiffCoef(i,j,k) = & & max( min( xyr_TempDiffCoef(i,j,k), TempDiffCoefMax ), & & TempDiffCoefMin & & ) end do end do end do ! xyr_QMixDiffCoef = xyr_TempDiffCoef ! 輸送係数とフラックスの計算 ! Calculate transfer coefficient and flux ! call VDiffusionCalcFlux( & & xyz_U, xyz_V, xyzf_QMix, & ! (in) & xyz_Temp, xyr_VirTemp, xyr_Press, & ! (in) & xyz_Height, xyz_Exner, xyr_Exner, & ! (in) & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, & ! (in) & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux & ! (out) & ) ! Calculate tendency of turbulent kinetic energy ! Set diffusion coefficient for turbulent kinetic energy xyz_TurKinEneDiffCoef = & & xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * Stke ! do k = 0, kmax if ( k == 0 ) then xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,1) else if ( k == kmax ) then xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,kmax) else xyr_TurKinEneDiffCoef(:,:,k) = & & ( xyz_TurKinEneDiffCoef(:,:,k) + xyz_TurKinEneDiffCoef(:,:,k+1) ) & & / 2.0_DP end if end do ! Calculate turbulent kinetic energy at lower boundary ! xy_FricVelSq = sqrt( xy_SurfMomFluxX**2 + xy_SurfMomFluxY**2 ) & & / ( xyr_Press(:,:,0) / ( GasRDry * xyr_VirTemp(:,:,0) ) ) xy_TurKinEneAtLB = MYConstB1**(2.0_DP/3.0_DP) / 2.0_DP * xy_FricVelSq xy_TurKinEneAtLB = xy_TurKinEneAtLB + TurKinEneOffset ! Calculate transfer coefficient and flux of turbulent kinetic energy ! ! When transfer coefficient at lower boundary is calculated, ! diffusion coefficient at mid-point of 1st layer is used. ! In addition, transfer coefficient at upper boundary is assumed ! to be zero. k = 0 xyr_TurKinEneTransCoef(:,:,k) = & & xyz_TurKinEneDiffCoef(:,:,1) & & * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) & & / ( xyz_Height(:,:,1) - xy_SurfHeight ) do k = 1, kmax-1 xyr_TurKinEneTransCoef(:,:,k) = & & xyr_TurKinEneDiffCoef(:,:,k) & & * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) end do k = kmax xyr_TurKinEneTransCoef(:,:,k) = 0.0_DP ! do k = 1, kmax-1 xyr_TurKinEneFlux(:,:,k) = & & - xyr_TurKinEneTransCoef(:,:,k) & & * ( xyz_TurKinEneNonZero(:,:,k+1) - xyz_TurKinEneNonZero(:,:,k) ) end do k = 0 xyr_TurKinEneFlux(:,:,k) = & & - xyr_TurKinEneTransCoef(:,:,k) & & * ( xyz_TurKinEneNonZero(:,:,k+1) - xy_TurKinEneAtLB ) k = kmax xyr_TurKinEneFlux(:,:,k) = 0.0_DP !!$ xyz_CShe1 = 2.0d0**1.5 * xyz_MixLength * xyz_SmHat & !!$ & * xyz_DVelDzSq & !!$! & * xyz_TurKinEne**1.5 !!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5 !!$ xyz_CShe2 = 1.5_DP * 2.0d0**1.5 * xyz_MixLength * xyz_SmHat & !!$ & * xyz_DVelDzSq & !!$! & * xyz_TurKinEne**0.5 !!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5 !!$ xyz_CBuo1 = - 2.0d0**1.5 * xyz_MixLength * xyz_ShHat & !!$ & * xyz_StatStab & !!$! & * xyz_TurKinEne**1.5 !!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5 !!$ xyz_CBuo2 = - 1.5_DP * 2.0d0**1.5 * xyz_MixLength * xyz_ShHat & !!$ & * xyz_StatStab & !!$! & * xyz_TurKinEne**0.5 !!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5 !!$ xyz_CDis1 = 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) & !!$! & * xyz_TurKinEne**1.5 !!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5 !!$ xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) & !!$! & * xyz_TurKinEne**0.5 !!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5 !!$! xyz_CShe1 = 2.0d0**1.5 * xyz_MixLength * xyz_SmHat & !!$! & * xyz_DVelDzSq & !!$!! & * xyz_TurKinEne**1.5 !!$! & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5 !!$ xyz_CShe1 = 2.0d0**1.5 * xyz_MixLength * xyz_Sm & !!$ & * xyz_DVelDzSq & !!$! & * xyz_TurKinEne**1.5 !!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5 !!$ xyz_CShe2 = 0.0_DP !!$! xyz_CBuo1 = - 2.0d0**1.5 * xyz_MixLength * xyz_ShHat & !!$! & * xyz_StatStab & !!$!! & * xyz_TurKinEne**1.5 !!$! & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5 !!$ xyz_CBuo1 = - 2.0d0**1.5 * xyz_MixLength * xyz_Sh & !!$ & * xyz_StatStab & !!$! & * xyz_TurKinEne**1.5 !!$ & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5 !!$ xyz_CBuo2 = 0.0_DP !!$ xyz_CDis1 = 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) & !!$ & * xyz_TurKinEne**1.5 !!$! & * ( xyz_TurKinEne + TurKinEneNonZero )**1.5 !!$ xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) & !!$ & * xyz_TurKinEne**0.5 !!$! & * ( xyz_TurKinEne + TurKinEneNonZero )**0.5 !!$ xyz_CShe1 = 2.0d0**1.5 * xyz_MixLength * xyz_DVelDzSq & !!$ & * xyz_SmHat & !!$ & * xyz_TurKinEneNonZero**1.5 !!$ xyz_CShe2 = 2.0d0**1.5 * xyz_MixLength * xyz_DVelDzSq & !!$ & * ( xyz_DSmHatDTKE & !!$ & * xyz_TurKinEneNonZero**1.5 & !!$ & + 1.5_DP * xyz_SmHat & !!$ & * xyz_TurKinEneNonZero**0.5 & !!$ & ) !!$ xyz_CBuo1 = - 2.0d0**1.5 * xyz_MixLength * xyz_StatStab & !!$ & * xyz_ShHat & !!$ & * xyz_TurKinEneNonZero**1.5 !!$ xyz_CBuo2 = - 2.0d0**1.5 * xyz_MixLength * xyz_StatStab & !!$ & * ( xyz_DShHatDTKE & !!$ & * xyz_TurKinEneNonZero**1.5 & !!$ & + 1.5_DP * xyz_ShHat & !!$ & * xyz_TurKinEneNonZero**0.5 & !!$ & ) !!$ xyz_CDis1 = 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) & !!$ & * xyz_TurKinEneNonZero**1.5 !!$ xyz_CDis2 = 1.5_DP * 2.0_DP**1.5 / ( MYConstB1 * xyz_MixLength ) & !!$ & * xyz_TurKinEneNonZero**0.5 xyz_CShe1 = sqrt( 2.0_DP ) * xyz_MixLength * xyz_DVelDzSq & & * xyz_Sm & & * sqrt( xyz_TurKinEneNonZero ) !!$ xyz_CShe2 = 0.5_DP * sqrt( 2.0_DP ) * xyz_MixLength * xyz_DVelDzSq & !!$ & * xyz_Sm & !!$ & / sqrt( xyz_TurKinEneNonZero ) xyz_CShe2 = 0.0_DP xyz_CBuo1 = - sqrt( 2.0_DP ) * xyz_MixLength * xyz_StatStab & & * xyz_Sh & & * sqrt( xyz_TurKinEneNonZero ) !!$ xyz_CBuo2 = - 0.5_DP * sqrt( 2.0_DP ) * xyz_MixLength * xyz_StatStab & !!$ & * xyz_Sh & !!$ & / sqrt( xyz_TurKinEneNonZero ) xyz_CBuo2 = 0.0_DP xyz_CDis1 = 2.0_DP**1.5_DP / ( MYConstB1 * xyz_MixLength ) & & * xyz_TurKinEneNonZero**1.5_DP xyz_CDis2 = 1.5_DP * 2.0_DP**1.5_DP / ( MYConstB1 * xyz_MixLength ) & & * xyz_TurKinEneNonZero**0.5_DP nloop = kmax loop_solve : do iloop = 1, nloop ! ! Construct implicit matrix from transfer coefficient of vertical ! diffusion scheme (turbulent kinetic energy) ! k = 1 xyza_TurKinEneMtx(:,:,k,-1) = 0.0_DP xyza_TurKinEneMtx(:,:,k, 0) = & & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) & & + xyr_TurKinEneTransCoef(:,:,k-1) & & + xyr_TurKinEneTransCoef(:,:,k ) & & + ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav & & * ( xyz_CShe2(:,:,k) + xyz_CBuo2(:,:,k) - xyz_CDis2(:,:,k) ) xyza_TurKinEneMtx(:,:,k, 1) = & & - xyr_TurKinEneTransCoef(:,:,k ) ! do k = 2, kmax-1 xyza_TurKinEneMtx(:,:,k,-1) =& & - xyr_TurKinEneTransCoef(:,:,k-1) xyza_TurKinEneMtx(:,:,k, 0) = & & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) & & + xyr_TurKinEneTransCoef(:,:,k-1) & & + xyr_TurKinEneTransCoef(:,:,k ) & & + ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav & & * ( xyz_CShe2(:,:,k) + xyz_CBuo2(:,:,k) - xyz_CDis2(:,:,k) ) xyza_TurKinEneMtx(:,:,k, 1) = & & - xyr_TurKinEneTransCoef(:,:,k ) end do ! k = kmax xyza_TurKinEneMtx(:,:,k,-1) =& & - xyr_TurKinEneTransCoef(:,:,k-1) xyza_TurKinEneMtx(:,:,k, 0) = & & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) & & + xyr_TurKinEneTransCoef(:,:,k-1) & & + xyr_TurKinEneTransCoef(:,:,k ) & & + ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav & & * ( xyz_CShe2(:,:,k) + xyz_CBuo2(:,:,k) - xyz_CDis2(:,:,k) ) xyza_TurKinEneMtx(:,:,k, 1) = 0.0_DP do k = 1, kmax xyz_TurKinEneVec(:,:,k) = & & - ( xyr_TurKinEneFlux(:,:,k) - xyr_TurKinEneFlux(:,:,k-1) ) & & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav & & * ( xyz_CShe1(:,:,k) + xyz_CBuo1(:,:,k) - xyz_CDis1(:,:,k) ) end do ! ! Solve simultaneous linear equations by use of LU decomposition technique ! xyza_TurKinEneLUMtx = xyza_TurKinEneMtx ! call PhyImplLUDecomp3( & & xyza_TurKinEneLUMtx, & ! (inout) & imax * jmax, kmax & ! (in) & ) xyz_DelTurKinEneLUVec = xyz_TurKinEneVec ! call PhyImplLUSolve3( & & xyz_DelTurKinEneLUVec, & ! (inout) & xyza_TurKinEneLUMtx, & ! (in) & 1, imax * jmax , kmax & ! (in) & ) xyz_DTurKinEneDt = xyz_DelTurKinEneLUVec / ( 2.0_DP * DelTime ) ! Calculation of dissipation rate of turbulent kinetic energy ! ! Calculate production rate of turbulent kinetic energy ! by shear and buoyancy xyz_TurKinEneProShear = xyz_CShe1 + xyz_CShe2 * xyz_DTurKinEneDt * 2.0_DP * DelTime xyz_TurKinEneProBuoya = xyz_CBuo1 + xyz_CBuo2 * xyz_DTurKinEneDt * 2.0_DP * DelTime xyz_TurKinEneDiss = xyz_CDis1 + xyz_CDis2 * xyz_DTurKinEneDt * 2.0_DP * DelTime ! Check of turbulent kinetic energy dissipation rate ! If it is negative, tendency is recalculated without dissipation. ! FlagReCalc = .false. do k = 1, kmax do j = 1, jmax do i = 0, imax-1 if ( xyz_TurKinEneDiss(i,j,k) < 0.0_DP ) then xyz_CDis1(i,j,k) = 0.0_DP xyz_CDis2(i,j,k) = 0.0_DP FlagReCalc = .true. end if end do end do end do ! Check convergence a_FlagReCalcLocal = FlagReCalc call MPIWrapperChkTrue( & & 1, a_FlagReCalcLocal, & ! (in) & a_FlagReCalcGlobal & ! (out) & ) if ( .not. a_FlagReCalcGlobal(1) ) exit loop_solve end do loop_solve !!$ write( 6, * ) TimeN, iloop !!$ write( 6, * ) xyz_TurKinEne(0,jmax/2+1,1:4) !!$ write( 6, * ) xyz_TempDiffCoef(0,jmax/2+1,1:4) !!$ write( 6, * ) xyz_TurKinEneProShear(0,jmax/2+1,1:4) !!$ write( 6, * ) xyz_TurKinEneProBuoya(0,jmax/2+1,1:4) !!$ write( 6, * ) xyz_TurKinEneDiss(0,jmax/2+1,1:4) !!$ write( 6, * ) xyz_DTurKinEneDt(0,jmax/2+1,1:4) ! 拡散係数の出力 ! Output diffusion coefficients ! ! 拡散係数出力 ! Diffusion coeffficients output ! call HistoryAutoPut( TimeN, 'VelDiffCoef', xyr_VelDiffCoef ) call HistoryAutoPut( TimeN, 'TempDiffCoef', xyr_TempDiffCoef ) call HistoryAutoPut( TimeN, 'QVapDiffCoef', xyr_QMixDiffCoef ) call HistoryAutoPut( TimeN, 'TKEPShear', xyz_TurKinEneProShear ) call HistoryAutoPut( TimeN, 'TKEPBuoy' , xyz_TurKinEneProBuoya ) call HistoryAutoPut( TimeN, 'TKEDiss' , xyz_TurKinEneDiss ) call HistoryAutoPut( TimeN, 'MixLength' , xyz_MixLength ) ! 計算時間計測一時停止 ! Pause measurement of computation time ! call TimesetClockStop( module_name ) end subroutine VDiffusionMY25 !-------------------------------------------------------------------------------------- subroutine VDiffusionMY251DWrapper3D( & & xyz_U, xyz_V, xyzf_QMix, & ! (in) & xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, & ! (in) & xy_SurfHeight, & ! (in) & xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, & ! (in) & xyz_TurKinEne, & ! (in) & xy_SurfMomFluxX, xy_SurfMomFluxY, & ! (in) & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (out) & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, & ! (in) & xyz_DTurKinEneDt & ! (out) & ) ! ! 鉛直拡散フラックスを計算します. ! ! Vertical diffusion flux is calculated by use of MY2.5 model. ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: & & FKarm, & ! $ k $ . ! カルマン定数. ! Karman constant & Grav, & ! $ g $ [m s-2]. ! 重力加速度. ! Gravitational acceleration & GasRDry, & ! $ R $ [J kg-1 K-1]. ! 乾燥大気の気体定数. ! Gas constant of air & CpDry ! $ C_p $ [J kg-1 K-1]. ! 乾燥大気の定圧比熱. ! Specific heat of air at constant pressure ! 時刻管理 ! Time control ! use timeset, only: & & DelTime, & ! $ \Delta t $ [s] & TimeN, & ! ステップ $ t $ の時刻. Time of step $ t $. & TimesetClockStart, TimesetClockStop ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! 陰解法による時間積分のためのルーチン ! Routines for time integration with implicit scheme ! use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3 ! 宣言文 ; Declaration statements ! implicit none real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax) ! $ u $ . 東西風速. Eastward wind real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax) ! $ v $ . 南北風速. Northward wind real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) ! $ q $ . 質量混合比. Mass mixing ratio real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax) ! $ T $ . 温度. Temperature real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{T} $ . 温度 (半整数レベル). ! Temperature (half level) real(DP), intent(in):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax) ! $ T_v $ . 仮温度. Virtual temperature real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{T}_v $ . 仮温度 (半整数レベル). ! Virtual temperature (half level) real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{p} $ . 気圧 (半整数レベル). ! Air pressure (half level) real(DP), intent(in):: xy_SurfHeight (0:imax-1,1:jmax) ! $ z_s $ . 地表面高度. ! Surface height. real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax) ! 高度 (整数レベル). ! Height (full level) real(DP), intent(in):: xyr_Height (0:imax-1, 1:jmax, 0:kmax) ! 高度 (半整数レベル). ! Height (half level) real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax) ! Exner 関数 (整数レベル). ! Exner function (full level) real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax) ! Exner 関数 (半整数レベル). ! Exner function (half level) real(DP), intent(in):: xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax) ! ! Turbulent kinetic energy (m2 s-2) real(DP), intent(in):: xy_SurfMomFluxX (0:imax-1, 1:jmax) ! ! Eastward momentum flux at surface real(DP), intent(in):: xy_SurfMomFluxY (0:imax-1, 1:jmax) ! ! Northward momentum flux at surface real(DP), intent(out):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax) ! 東西方向運動量フラックス. ! Eastward momentum flux real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax) ! 南北方向運動量フラックス. ! Northward momentum flux real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax) ! 熱フラックス. ! Heat flux real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) ! 質量フラックス. ! Mass flux of compositions real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：温度. ! Diffusion coefficient: temperature real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：比湿. ! Diffusion coefficient: specific humidity real(DP), intent(out):: xyz_DTurKinEneDt (0:imax-1, 1:jmax, 1:kmax) ! ! Tendency of turbulent kinetic energy ! 作業変数 ! Work variables ! real(DP) :: z_U (1:kmax) real(DP) :: z_V (1:kmax) real(DP) :: zf_QMix(1:kmax, 1:ncmax) real(DP) :: z_Temp (1:kmax) real(DP) :: r_Temp (0:kmax) real(DP) :: z_VirTemp (1:kmax) real(DP) :: r_VirTemp (0:kmax) real(DP) :: r_Press (0:kmax) real(DP) :: SurfHeight real(DP) :: z_Height (1:kmax) real(DP) :: r_Height (0:kmax) real(DP) :: z_Exner (1:kmax) real(DP) :: r_Exner (0:kmax) real(DP) :: z_TurKinEne(1:kmax) real(DP) :: SurfMomFluxX real(DP) :: SurfMomFluxY real(DP) :: r_MomFluxX (0:kmax) real(DP) :: r_MomFluxY (0:kmax) real(DP) :: r_HeatFlux (0:kmax) real(DP) :: rf_QMixFlux(0:kmax, 1:ncmax) real(DP) :: r_VelDiffCoef (0:kmax) real(DP) :: r_TempDiffCoef(0:kmax) real(DP) :: r_QMixDiffCoef(0:kmax) real(DP) :: z_DTurKinEneDt (1:kmax) ! ! Tendency of turbulent kinetic energy integer:: i ! 経度方向に回る DO ループ用作業変数 ! Work variables for DO loop in longitude integer:: j ! 緯度方向に回る DO ループ用作業変数 ! Work variables for DO loop in latitude ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. vdiffusion_my_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 計算時間計測開始 ! Start measurement of computation time ! call TimesetClockStart( module_name ) do j = 1, jmax do i = 0, imax-1 z_U = xyz_U (i,j,:) z_V = xyz_V (i,j,:) zf_QMix = xyzf_QMix (i,j,:,:) z_Temp = xyz_Temp (i,j,:) r_Temp = xyr_Temp (i,j,:) z_VirTemp = xyz_VirTemp (i,j,:) r_VirTemp = xyr_VirTemp (i,j,:) r_Press = xyr_Press (i,j,:) SurfHeight = xy_SurfHeight (i,j) z_Height = xyz_Height (i,j,:) r_Height = xyr_Height (i,j,:) z_Exner = xyz_Exner (i,j,:) r_Exner = xyr_Exner (i,j,:) z_TurKinEne = xyz_TurKinEne (i,j,:) SurfMomFluxX = xy_SurfMomFluxX(i,j) SurfMomFluxY = xy_SurfMomFluxY(i,j) call VDiffusionMY251D( & & z_U, z_V, zf_QMix, & ! (in) & z_Temp, r_Temp, z_VirTemp, r_VirTemp, r_Press, & ! (in) & SurfHeight, & ! (in) & z_Height, r_Height, z_Exner, r_Exner, & ! (in) & z_TurKinEne, & ! (in) & SurfMomFluxX, SurfMomFluxY, & ! (in) & r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux, & ! (out) & r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef, & ! (out) & z_DTurKinEneDt & ! (out) & ) xyr_MomFluxX (i,j,:) = r_MomFluxX xyr_MomFluxY (i,j,:) = r_MomFluxY xyr_HeatFlux (i,j,:) = r_HeatFlux xyrf_QMixFlux (i,j,:,:) = rf_QMixFlux xyr_VelDiffCoef (i,j,:) = r_VelDiffCoef xyr_TempDiffCoef(i,j,:) = r_TempDiffCoef xyr_QMixDiffCoef(i,j,:) = r_QMixDiffCoef xyz_DTurKinEneDt (i,j,:) = z_DTurKinEneDt end do end do !!$ ! 拡散係数の出力 !!$ ! Output diffusion coefficients !!$ ! !!$ !!$ ! 拡散係数出力 !!$ ! Diffusion coeffficients output !!$ ! !!$ call HistoryAutoPut( TimeN, 'VelDiffCoef', xyr_VelDiffCoef ) !!$ call HistoryAutoPut( TimeN, 'TempDiffCoef', xyr_TempDiffCoef ) !!$ call HistoryAutoPut( TimeN, 'QVapDiffCoef', xyr_QMixDiffCoef ) !!$ !!$ call HistoryAutoPut( TimeN, 'TKEPShear', xyz_TurKinEneProShear ) !!$ call HistoryAutoPut( TimeN, 'TKEPBuoy' , xyz_TurKinEneProBuoya ) !!$ call HistoryAutoPut( TimeN, 'TKEDiss' , xyz_TurKinEneDiss ) !!$ !!$ call HistoryAutoPut( TimeN, 'MixLength' , xyz_MixLength ) ! 計算時間計測一時停止 ! Pause measurement of computation time ! call TimesetClockStop( module_name ) end subroutine VDiffusionMY251DWrapper3D !-------------------------------------------------------------------------------------- subroutine VDiffusionMY251D( & & z_U, z_V, zf_QMix, & ! (in) & z_Temp, r_Temp, z_VirTemp, r_VirTemp, r_Press, & ! (in) & SurfHeight, & ! (in) & z_Height, r_Height, z_Exner, r_Exner, & ! (in) & z_TurKinEne, & ! (in) & SurfMomFluxX, SurfMomFluxY, & ! (in) & r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux, & ! (out) & r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef, & ! (out) & z_DTurKinEneDt & ! (out) & ) ! ! 鉛直拡散フラックスを計算します. ! ! Vertical diffusion flux is calculated by use of MY2.5 model. ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: & & FKarm, & ! $ k $ . ! カルマン定数. ! Karman constant & Grav, & ! $ g $ [m s-2]. ! 重力加速度. ! Gravitational acceleration & GasRDry, & ! $ R $ [J kg-1 K-1]. ! 乾燥大気の気体定数. ! Gas constant of air & CpDry ! $ C_p $ [J kg-1 K-1]. ! 乾燥大気の定圧比熱. ! Specific heat of air at constant pressure ! 時刻管理 ! Time control ! use timeset, only: & & DelTime, & ! $ \Delta t $ [s] & TimeN, & ! ステップ $ t $ の時刻. Time of step $ t $. & TimesetClockStart, TimesetClockStop ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! 陰解法による時間積分のためのルーチン ! Routines for time integration with implicit scheme ! use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3 ! 宣言文 ; Declaration statements ! implicit none real(DP), intent(in):: z_U (1:kmax) ! $ u $ . 東西風速. Eastward wind real(DP), intent(in):: z_V (1:kmax) ! $ v $ . 南北風速. Northward wind real(DP), intent(in):: zf_QMix(1:kmax, 1:ncmax) ! $ q $ . 質量混合比. Mass mixing ratio real(DP), intent(in):: z_Temp (1:kmax) ! $ T $ . 温度. Temperature real(DP), intent(in):: r_Temp (0:kmax) ! $ \hat{T} $ . 温度 (半整数レベル). ! Temperature (half level) real(DP), intent(in):: z_VirTemp (1:kmax) ! $ T_v $ . 仮温度. Virtual temperature real(DP), intent(in):: r_VirTemp (0:kmax) ! $ \hat{T}_v $ . 仮温度 (半整数レベル). ! Virtual temperature (half level) real(DP), intent(in):: r_Press (0:kmax) ! $ \hat{p} $ . 気圧 (半整数レベル). ! Air pressure (half level) real(DP), intent(in):: SurfHeight ! $ z_s $ . 地表面高度. ! Surface height. real(DP), intent(in):: z_Height (1:kmax) ! 高度 (整数レベル). ! Height (full level) real(DP), intent(in):: r_Height (0:kmax) ! 高度 (半整数レベル). ! Height (half level) real(DP), intent(in):: z_Exner (1:kmax) ! Exner 関数 (整数レベル). ! Exner function (full level) real(DP), intent(in):: r_Exner (0:kmax) ! Exner 関数 (半整数レベル). ! Exner function (half level) real(DP), intent(in):: z_TurKinEne(1:kmax) ! ! Turbulent kinetic energy (m2 s-2) real(DP), intent(in):: SurfMomFluxX ! ! Eastward momentum flux at surface real(DP), intent(in):: SurfMomFluxY ! ! Northward momentum flux at surface real(DP), intent(out):: r_MomFluxX (0:kmax) ! 東西方向運動量フラックス. ! Eastward momentum flux real(DP), intent(out):: r_MomFluxY (0:kmax) ! 南北方向運動量フラックス. ! Northward momentum flux real(DP), intent(out):: r_HeatFlux (0:kmax) ! 熱フラックス. ! Heat flux real(DP), intent(out):: rf_QMixFlux(0:kmax, 1:ncmax) ! 質量フラックス. ! Mass flux of compositions real(DP), intent(out):: r_VelDiffCoef (0:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP), intent(out):: r_TempDiffCoef(0:kmax) ! 拡散係数：温度. ! Diffusion coefficient: temperature real(DP), intent(out):: r_QMixDiffCoef(0:kmax) ! 拡散係数：比湿. ! Diffusion coefficient: specific humidity real(DP), intent(out):: z_DTurKinEneDt (1:kmax) ! ! Tendency of turbulent kinetic energy ! 作業変数 ! Work variables ! real(DP) :: z_MixLength(1:kmax) ! 混合距離. ! Mixing length real(DP) :: z_DVelDzSq(1:kmax) ! ! Vertical shear squared (s-2) real(DP) :: z_StatStab(1:kmax) ! ! Static stability (s-2) real(DP) :: GhMin ! ! Minimum of G_h real(DP) :: GhMax ! ! Maximum of G_h real(DP) :: z_Gm(1:kmax) ! ! G_m real(DP) :: z_Gh(1:kmax) ! ! G_h real(DP) :: z_Sm(1:kmax) ! ! S_M real(DP) :: z_Sh(1:kmax) ! ! S_h real(DP), parameter :: Stke = 0.2_DP ! ! S_{TKE} = 0.2 real(DP) :: z_VelDiffCoef (1:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP) :: z_TempDiffCoef (1:kmax) ! 拡散係数：温度. ! Transfer coefficient: temperature real(DP) :: r_TurKinEneDiffCoef (0:kmax) ! ! Diffusion coefficient: turbulent kinetic energy real(DP) :: z_TurKinEneDiffCoef (1:kmax) ! ! Diffusion coefficient: turbulent kinetic energy real(DP) :: r_TurKinEneTransCoef(0:kmax) ! ! Transfer coefficient: turbulent kinetic energy real(DP) :: r_TurKinEneFlux(0:kmax) ! ! Turbulent energy flux real(DP) :: z_CShe1(1:kmax) real(DP) :: z_CShe2(1:kmax) real(DP) :: z_CBuo1(1:kmax) real(DP) :: z_CBuo2(1:kmax) real(DP) :: z_CDis1(1:kmax) real(DP) :: z_CDis2(1:kmax) real(DP) :: z_TurKinEneProShear(1:kmax) real(DP) :: z_TurKinEneProBuoya(1:kmax) real(DP) :: FricVelSq real(DP) :: TurKinEneAtLB real(DP) :: za_TurKinEneMtx(1:kmax, -1:1) ! ! Implicit matrix for turbulent kinetic energy real(DP) :: z_TurKinEneVec(1:kmax) ! ! Implicit vector for turbulent kinetic energy real(DP) :: aaza_TurKinEneLUMtx (1:1, 1:1, 1:kmax, -1:1) ! LU 行列. ! LU matrix real(DP) :: aaz_DelTurKinEneLUVec(1:1, 1:1, 1:kmax) ! ! Tendency of turbulent kinetic energy real(DP) :: z_TurKinEneDiss(1:kmax) ! ! Dissipation rate of turbulent kinetic energy (m2 s-3) real(DP) :: z_TurKinEneNonZero(1:kmax) ! ! Turbulent kinetic energy with offset (m2 s-2) real(DP), parameter :: TurKinEneOffset = ( 1.0e-3_DP )**2 / 2.0_DP logical :: FlagReCalc ! ! Flag for recalculation logical :: a_FlagReCalcLocal (1) logical :: a_FlagReCalcGlobal(1) integer :: iloop integer :: nloop integer:: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer:: n ! 組成方向に回る DO ループ用作業変数 ! Work variables for DO loop in dimension of constituents ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. vdiffusion_my_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if !!$ ! 計算時間計測開始 !!$ ! Start measurement of computation time !!$ ! !!$ call TimesetClockStart( module_name ) ! Calculate turbulent kinetic energy with offset ! z_TurKinEneNonZero = z_TurKinEne + TurKinEneOffset ! ! Calculation of vertical shear squared do k = 1, kmax if ( k == 1 ) then z_DVelDzSq(k) = & & ( ( z_U(k+1) - z_U(k ) )**2 & & + ( z_V(k+1) - z_V(k ) )**2 ) & & / ( z_Height(k+1) - z_Height(k ) )**2 else if ( k == kmax ) then z_DVelDzSq(k) = & & ( ( z_U(k ) - z_U(k-1) )**2 & & + ( z_V(k ) - z_V(k-1) )**2 ) & & / ( z_Height(k ) - z_Height(k-1) )**2 else z_DVelDzSq(k) = & & ( ( z_U(k+1) - z_U(k-1) )**2 & & + ( z_V(k+1) - z_V(k-1) )**2 ) & & / ( z_Height(k+1) - z_Height(k-1) )**2 end if end do ! Calculation of static stability do k = 1, kmax if ( k == 1 ) then z_StatStab(k) = & & Grav / ( z_VirTemp(k) / z_Exner(k) ) & & * ( z_VirTemp(k+1) / z_Exner(k+1) & & - z_VirTemp(k ) / z_Exner(k ) ) & & / ( z_Height(k+1) - z_Height(k ) ) else if ( k == kmax ) then z_StatStab(k) = & & Grav / ( z_VirTemp(k) / z_Exner(k) ) & & * ( z_VirTemp(k ) / z_Exner(k ) & & - z_VirTemp(k-1) / z_Exner(k-1) ) & & / ( z_Height(k ) - z_Height(k-1) ) else z_StatStab(k) = & & Grav / ( z_VirTemp(k) / z_Exner(k) ) & & * ( z_VirTemp(k+1) / z_Exner(k+1) & & - z_VirTemp(k-1) / z_Exner(k-1) ) & & / ( z_Height(k+1) - z_Height(k-1) ) end if end do ! 混合距離の算出 ! Calculate mixing length ! do k = 1, kmax z_MixLength(k) = & & FKarm * ( z_Height(k) - SurfHeight ) & & / (1.0_DP + FKarm * ( z_Height(k) - SurfHeight ) / MixLengthMax ) end do ! Limit mixing length (Galperin et al., 1988) and avoid zero z_MixLength = & & min( z_MixLength, & & 0.53_DP & & * sqrt( 2.0_DP * z_TurKinEneNonZero / max( z_StatStab, 1.0e-10_DP ) ) ) & & + 1.0e-10_DP z_Gh = - z_MixLength**2 / ( 2.0_DP * z_TurKinEneNonZero ) & & * z_StatStab ! Actually, xyz_Gm is not used below. z_Gm = z_MixLength**2 / ( 2.0_DP * z_TurKinEneNonZero ) & & * z_DVelDzSq ! Limit Gh (Galperin et al., 1988) GhMin = - 0.53_DP**2 GhMax = 1.0_DP & & / ( MYConstA2 & & * ( 12.0_DP * MYConstA1 + MYConstB1 + 3.0_DP * MYConstB2 ) ) z_Gh = max( GhMin, min( z_Gh, GhMax ) ) z_Sh = & & MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) & & / ( 1.0_DP & & - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * z_Gh ) z_Sm = & & ( & & MYConstA1 * ( 1.0_DP - 3.0_DP * MYConstC1 & & - 6.0_DP * MYConstA1 / MYConstB1 ) & & + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) & & * z_Gh * z_Sh & & ) & & / ( 1.0_DP & & - 9.0_DP * MYConstA1 * MYConstA2 * z_Gh ) ! 拡散係数の計算 ! Calculation of diffusion coefficient ! z_VelDiffCoef = z_MixLength * sqrt( 2.0_DP * z_TurKinEneNonZero ) * z_Sm z_TempDiffCoef = z_MixLength * sqrt( 2.0_DP * z_TurKinEneNonZero ) * z_Sh ! do k = 0, kmax if ( ( k == 0 ) .or. ( k == kmax ) ) then r_VelDiffCoef (k) = 0.0_DP r_TempDiffCoef(k) = 0.0_DP else r_VelDiffCoef (k) = ( z_VelDiffCoef (k) + z_VelDiffCoef (k+1) ) / 2.0_DP r_TempDiffCoef(k) = ( z_TempDiffCoef(k) + z_TempDiffCoef(k+1) ) / 2.0_DP end if end do ! do k = 1, kmax-1 r_VelDiffCoef(k) = & & max( min( r_VelDiffCoef(k), VelDiffCoefMax ), & & VelDiffCoefMin & & ) r_TempDiffCoef(k) = & & max( min( r_TempDiffCoef(k), TempDiffCoefMax ), & & TempDiffCoefMin & & ) end do ! r_QMixDiffCoef = r_TempDiffCoef ! 輸送係数とフラックスの計算 ! Calculate transfer coefficient and flux ! call VDiffusionCalcFlux1D( & & z_U, z_V, zf_QMix, & ! (in) & z_Temp, r_VirTemp, r_Press, & ! (in) & z_Height, z_Exner, r_Exner, & ! (in) & r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef, & ! (in) & r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux & ! (out) & ) ! Calculate tendency of turbulent kinetic energy ! Set diffusion coefficient for turbulent kinetic energy z_TurKinEneDiffCoef = & & z_MixLength * sqrt( 2.0_DP * z_TurKinEneNonZero ) * Stke ! do k = 0, kmax if ( k == 0 ) then r_TurKinEneDiffCoef(k) = z_TurKinEneDiffCoef(1) else if ( k == kmax ) then r_TurKinEneDiffCoef(k) = z_TurKinEneDiffCoef(kmax) else r_TurKinEneDiffCoef(k) = & & ( z_TurKinEneDiffCoef(k) + z_TurKinEneDiffCoef(k+1) ) & & / 2.0_DP end if end do ! Calculate turbulent kinetic energy at lower boundary ! FricVelSq = sqrt( SurfMomFluxX**2 + SurfMomFluxY**2 ) & & / ( r_Press(0) / ( GasRDry * r_VirTemp(0) ) ) TurKinEneAtLB = MYConstB1**(2.0_DP/3.0_DP) / 2.0_DP * FricVelSq TurKinEneAtLB = TurKinEneAtLB + TurKinEneOffset ! Calculate transfer coefficient and flux of turbulent kinetic energy ! ! When transfer coefficient at lower boundary is calculated, ! diffusion coefficient at mid-point of 1st layer is used. ! In addition, transfer coefficient at upper boundary is assumed ! to be zero. k = 0 r_TurKinEneTransCoef(k) = & & z_TurKinEneDiffCoef(1) & & * r_Press(k) / ( GasRDry * r_VirTemp(k) ) & & / ( z_Height(1) - SurfHeight ) do k = 1, kmax-1 r_TurKinEneTransCoef(k) = & & r_TurKinEneDiffCoef(k) & & * r_Press(k) / ( GasRDry * r_VirTemp(k) ) & & / ( z_Height(k+1) - z_Height(k) ) end do k = kmax r_TurKinEneTransCoef(k) = 0.0_DP ! do k = 1, kmax-1 r_TurKinEneFlux(k) = & & - r_TurKinEneTransCoef(k) & & * ( z_TurKinEneNonZero(k+1) - z_TurKinEneNonZero(k) ) end do k = 0 r_TurKinEneFlux(k) = & & - r_TurKinEneTransCoef(k) & & * ( z_TurKinEneNonZero(k+1) - TurKinEneAtLB ) k = kmax r_TurKinEneFlux(k) = 0.0_DP z_CShe1 = sqrt( 2.0_DP ) * z_MixLength * z_DVelDzSq & & * z_Sm & & * sqrt( z_TurKinEneNonZero ) z_CShe2 = 0.0_DP z_CBuo1 = - sqrt( 2.0_DP ) * z_MixLength * z_StatStab & & * z_Sh & & * sqrt( z_TurKinEneNonZero ) z_CBuo2 = 0.0_DP z_CDis1 = 2.0_DP**1.5_DP / ( MYConstB1 * z_MixLength ) & & * z_TurKinEneNonZero**1.5_DP z_CDis2 = 1.5_DP * 2.0_DP**1.5_DP / ( MYConstB1 * z_MixLength ) & & * z_TurKinEneNonZero**0.5_DP nloop = kmax loop_solve : do iloop = 1, nloop ! ! Construct implicit matrix from transfer coefficient of vertical ! diffusion scheme (turbulent kinetic energy) ! k = 1 za_TurKinEneMtx(k,-1) = 0.0_DP za_TurKinEneMtx(k, 0) = & & - ( r_Press(k) - r_Press(k-1) ) / Grav / ( 2.0_DP * DelTime ) & & + r_TurKinEneTransCoef(k-1) & & + r_TurKinEneTransCoef(k ) & & + ( r_Press(k) - r_Press(k-1) ) / Grav & & * ( z_CShe2(k) + z_CBuo2(k) - z_CDis2(k) ) za_TurKinEneMtx(k, 1) = & & - r_TurKinEneTransCoef(k ) ! do k = 2, kmax-1 za_TurKinEneMtx(k,-1) =& & - r_TurKinEneTransCoef(k-1) za_TurKinEneMtx(k, 0) = & & - ( r_Press(k) - r_Press(k-1) ) / Grav / ( 2.0_DP * DelTime ) & & + r_TurKinEneTransCoef(k-1) & & + r_TurKinEneTransCoef(k ) & & + ( r_Press(k) - r_Press(k-1) ) / Grav & & * ( z_CShe2(k) + z_CBuo2(k) - z_CDis2(k) ) za_TurKinEneMtx(k, 1) = & & - r_TurKinEneTransCoef(k ) end do ! k = kmax za_TurKinEneMtx(k,-1) =& & - r_TurKinEneTransCoef(k-1) za_TurKinEneMtx(k, 0) = & & - ( r_Press(k) - r_Press(k-1) ) / Grav / ( 2.0_DP * DelTime ) & & + r_TurKinEneTransCoef(k-1) & & + r_TurKinEneTransCoef(k ) & & + ( r_Press(k) - r_Press(k-1) ) / Grav & & * ( z_CShe2(k) + z_CBuo2(k) - z_CDis2(k) ) za_TurKinEneMtx(k, 1) = 0.0_DP do k = 1, kmax z_TurKinEneVec(k) = & & - ( r_TurKinEneFlux(k) - r_TurKinEneFlux(k-1) ) & & - ( r_Press(k) - r_Press(k-1) ) / Grav & & * ( z_CShe1(k) + z_CBuo1(k) - z_CDis1(k) ) end do ! ! Solve simultaneous linear equations by use of LU decomposition technique ! aaza_TurKinEneLUMtx(1,1,:,:) = za_TurKinEneMtx ! call PhyImplLUDecomp3( & & aaza_TurKinEneLUMtx, & ! (inout) & 1 * 1, kmax & ! (in) & ) aaz_DelTurKinEneLUVec(1,1,:) = z_TurKinEneVec ! call PhyImplLUSolve3( & & aaz_DelTurKinEneLUVec, & ! (inout) & aaza_TurKinEneLUMtx, & ! (in) & 1, 1 * 1 , kmax & ! (in) & ) z_DTurKinEneDt = aaz_DelTurKinEneLUVec(1,1,:) / ( 2.0_DP * DelTime ) ! Calculation of dissipation rate of turbulent kinetic energy ! ! Calculate production rate of turbulent kinetic energy ! by shear and buoyancy z_TurKinEneProShear = z_CShe1 + z_CShe2 * z_DTurKinEneDt * 2.0_DP * DelTime z_TurKinEneProBuoya = z_CBuo1 + z_CBuo2 * z_DTurKinEneDt * 2.0_DP * DelTime z_TurKinEneDiss = z_CDis1 + z_CDis2 * z_DTurKinEneDt * 2.0_DP * DelTime ! Check of turbulent kinetic energy dissipation rate ! If it is negative, tendency is recalculated without dissipation. ! FlagReCalc = .false. do k = 1, kmax if ( z_TurKinEneDiss(k) < 0.0_DP ) then z_CDis1(k) = 0.0_DP z_CDis2(k) = 0.0_DP FlagReCalc = .true. end if end do ! Check convergence a_FlagReCalcLocal = FlagReCalc !!$ call MPIWrapperChkTrue( & !!$ & 1, a_FlagReCalcLocal, & ! (in) !!$ & a_FlagReCalcGlobal & ! (out) !!$ & ) a_FlagReCalcGlobal = a_FlagReCalcLocal if ( .not. a_FlagReCalcGlobal(1) ) exit loop_solve end do loop_solve !!$ ! 計算時間計測一時停止 !!$ ! Pause measurement of computation time !!$ ! !!$ call TimesetClockStop( module_name ) end subroutine VDiffusionMY251D !-------------------------------------------------------------------------------------- subroutine VDiffusionMY25GBT94( & & xyz_U, xyz_V, xyzf_QMix, & ! (in) & xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, & ! (in) & xy_SurfHeight, & ! (in) & xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, & ! (in) & xyz_TurKinEne, & ! (in) & xy_SurfMomFluxX, xy_SurfMomFluxY, & ! (in) & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (out) & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, & ! (in) & xyz_DTurKinEneDt & ! (out) & ) ! ! 鉛直拡散フラックスを計算します. ! ! Vertical diffusion flux is calculated by use of MY2.5 model. ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: & & FKarm, & ! $ k $ . ! カルマン定数. ! Karman constant & Grav, & ! $ g $ [m s-2]. ! 重力加速度. ! Gravitational acceleration & GasRDry, & ! $ R $ [J kg-1 K-1]. ! 乾燥大気の気体定数. ! Gas constant of air & CpDry ! $ C_p $ [J kg-1 K-1]. ! 乾燥大気の定圧比熱. ! Specific heat of air at constant pressure ! 時刻管理 ! Time control ! use timeset, only: & & DelTime, & ! $ \Delta t $ [s] & TimeN, & ! ステップ $ t $ の時刻. Time of step $ t $. & TimesetClockStart, TimesetClockStop ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! MPI 関連ルーチン ! MPI related routines ! use mpi_wrapper, only : MPIWrapperChkTrue ! 陰解法による時間積分のためのルーチン ! Routines for time integration with implicit scheme ! use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3 ! 宣言文 ; Declaration statements ! implicit none real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax) ! $ u $ . 東西風速. Eastward wind real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax) ! $ v $ . 南北風速. Northward wind real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) ! $ q $ . 質量混合比. Mass mixing ratio real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax) ! $ T $ . 温度. Temperature real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{T} $ . 温度 (半整数レベル). ! Temperature (half level) real(DP), intent(in):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax) ! $ T_v $ . 仮温度. Virtual temperature real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{T}_v $ . 仮温度 (半整数レベル). ! Virtual temperature (half level) real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{p} $ . 気圧 (半整数レベル). ! Air pressure (half level) real(DP), intent(in):: xy_SurfHeight (0:imax-1,1:jmax) ! $ z_s $ . 地表面高度. ! Surface height. real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax) ! 高度 (整数レベル). ! Height (full level) real(DP), intent(in):: xyr_Height (0:imax-1, 1:jmax, 0:kmax) ! 高度 (半整数レベル). ! Height (half level) real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax) ! Exner 関数 (整数レベル). ! Exner function (full level) real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax) ! Exner 関数 (半整数レベル). ! Exner function (half level) real(DP), intent(in):: xyz_TurKinEne(0:imax-1, 1:jmax, 1:kmax) ! ! Turbulent kinetic energy (m2 s-2) real(DP), intent(in):: xy_SurfMomFluxX (0:imax-1, 1:jmax) ! ! Eastward momentum flux at surface real(DP), intent(in):: xy_SurfMomFluxY (0:imax-1, 1:jmax) ! ! Northward momentum flux at surface real(DP), intent(out):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax) ! 東西方向運動量フラックス. ! Eastward momentum flux real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax) ! 南北方向運動量フラックス. ! Northward momentum flux real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax) ! 熱フラックス. ! Heat flux real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) ! 質量フラックス. ! Mass flux of compositions real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：温度. ! Diffusion coefficient: temperature real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：比湿. ! Diffusion coefficient: specific humidity real(DP), intent(out):: xyz_DTurKinEneDt (0:imax-1, 1:jmax, 1:kmax) ! ! Tendency of turbulent kinetic energy ! 作業変数 ! Work variables ! real(DP) :: xyz_MixLength(0:imax-1, 1:jmax, 1:kmax) ! 混合距離. ! Mixing length real(DP) :: xyz_DVelDzSq(0:imax-1, 1:jmax, 1:kmax) ! ! Vertical shear squared (s-2) real(DP) :: xyz_StatStab(0:imax-1, 1:jmax, 1:kmax) ! ! Static stability (s-2) !!$ real(DP) :: GhMin !!$ ! !!$ ! Minimum of G_h !!$ real(DP) :: GhMax !!$ ! !!$ ! Maximum of G_h real(DP) :: xyz_Gm(0:imax-1, 1:jmax, 1:kmax) ! ! G_m real(DP) :: xyz_Gh(0:imax-1, 1:jmax, 1:kmax) ! ! G_h !!$ real(DP) :: xyz_SmHat(0:imax-1, 1:jmax, 1:kmax) !!$ ! !!$ ! \hat{S}_M !!$ real(DP) :: xyz_ShHat(0:imax-1, 1:jmax, 1:kmax) !!$ ! !!$ ! \hat{S}_h !!$ real(DP) :: xyz_DSmHatDTKE(0:imax-1, 1:jmax, 1:kmax) !!$ ! !!$ ! derivative of \hat{S}_M !!$ real(DP) :: xyz_DShHatDTKE(0:imax-1, 1:jmax, 1:kmax) !!$ ! !!$ ! derivative of \hat{S}_h real(DP) :: xyz_Sm(0:imax-1, 1:jmax, 1:kmax) ! ! S_M real(DP) :: xyz_Sh(0:imax-1, 1:jmax, 1:kmax) ! ! S_h real(DP), parameter :: Stke = 0.2_DP ! ! S_{TKE} = 0.2 real(DP) :: xyz_VelDiffCoef (0:imax-1, 1:jmax, 1:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP) :: xyz_TempDiffCoef (0:imax-1, 1:jmax, 1:kmax) ! 拡散係数：温度. ! Transfer coefficient: temperature real(DP) :: xyr_TurKinEneDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! ! Diffusion coefficient: turbulent kinetic energy real(DP) :: xyz_TurKinEneDiffCoef (0:imax-1, 1:jmax, 1:kmax) ! ! Diffusion coefficient: turbulent kinetic energy real(DP) :: xyr_TurKinEneTransCoef(0:imax-1, 1:jmax, 0:kmax) ! ! Transfer coefficient: turbulent kinetic energy real(DP) :: xyr_TurKinEneFlux(0:imax-1, 1:jmax, 0:kmax) ! ! Turbulent energy flux real(DP) :: xy_FricVelSq (0:imax-1, 1:jmax) real(DP) :: xy_TurKinEneAtLB(0:imax-1, 1:jmax) real(DP) :: xyza_TurKinEneMtx(0:imax-1, 1:jmax, 1:kmax, -1:1) ! ! Implicit matrix for turbulent kinetic energy real(DP) :: xyz_TurKinEneVec(0:imax-1, 1:jmax, 1:kmax) ! ! Implicit vector for turbulent kinetic energy real(DP) :: xyza_TurKinEneLUMtx (0:imax-1, 1:jmax, 1:kmax, -1:1) ! LU 行列. ! LU matrix real(DP) :: xyz_DelTurKinEneLUVec(0:imax-1, 1:jmax, 1:kmax) ! ! Tendency of turbulent kinetic energy real(DP) :: xyz_TurKinEneDiss(0:imax-1, 1:jmax, 1:kmax) ! ! Dissipation rate of turbulent kinetic energy (m2 s-3) real(DP) :: xyz_TurKinEneNonZero(0:imax-1, 1:jmax, 1:kmax) ! ! Turbulent kinetic energy with offset (m2 s-2) real(DP), parameter :: TurKinEneOffset = ( 1.0e-3_DP )**2 / 2.0_DP logical :: xyz_FlagUseRiNum (0:imax-1, 1:jmax, 1:kmax) logical :: xyz_FlagTKEAsymptToZero(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_RiNum(0:imax-1, 1:jmax, 1:kmax) real(DP) :: Beta1 real(DP) :: Beta2 real(DP) :: Beta3 real(DP) :: Beta4 real(DP) :: MixLength real(DP) :: xyz_A1 (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_A2 (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_R1 (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_R2 (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_SqrtA1(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_TKEInit (0:imax-1, 1:jmax, 1:kmax) ! ! Turbulent kinetic energy for current time step ! with offset (m2 s-2) real(DP) :: xyz_TKEx2 (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_SqrtTKEx2 (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_TKETentative(0:imax-1, 1:jmax, 1:kmax) real(DP) :: Alpha real(DP) :: BetaSq real(DP) :: StatStab real(DP) :: RiNum real(DP) :: DVelDzSq real(DP), parameter :: Epsilon = 1.0e-10_DP real(DP), parameter :: CrtlRiNum = 0.195_DP real(DP), parameter :: CrtlShear = 0.001_DP / 1000.0_DP !!$ logical :: FlagReCalc !!$ ! !!$ ! Flag for recalculation !!$ logical :: a_FlagReCalcLocal (1) !!$ logical :: a_FlagReCalcGlobal(1) !!$ integer :: iloop !!$ integer :: nloop integer:: i ! 経度方向に回る DO ループ用作業変数 ! Work variables for DO loop in longitude integer:: j ! 緯度方向に回る DO ループ用作業変数 ! Work variables for DO loop in latitude integer:: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction !!$ integer:: n ! 組成方向に回る DO ループ用作業変数 !!$ ! Work variables for DO loop in dimension of constituents ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. vdiffusion_my_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 計算時間計測開始 ! Start measurement of computation time ! call TimesetClockStart( module_name ) ! ! Calculation of vertical shear squared do k = 1, kmax if ( k == 1 ) then xyz_DVelDzSq(:,:,k) = & & ( ( xyz_U(:,:,k+1) - xyz_U(:,:,k ) )**2 & & + ( xyz_V(:,:,k+1) - xyz_V(:,:,k ) )**2 ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k ) )**2 !!$ xyz_DVelDzSq(:,:,k) = & !!$ & ( ( xyz_U(:,:,k+1) - 0.0_DP )**2 & !!$ & + ( xyz_V(:,:,k+1) - 0.0_DP )**2 ) & !!$ & / ( xyz_Height(:,:,k+1) - xy_SurfHeight )**2 else if ( k == kmax ) then xyz_DVelDzSq(:,:,k) = & & ( ( xyz_U(:,:,k ) - xyz_U(:,:,k-1) )**2 & & + ( xyz_V(:,:,k ) - xyz_V(:,:,k-1) )**2 ) & & / ( xyz_Height(:,:,k ) - xyz_Height(:,:,k-1) )**2 else xyz_DVelDzSq(:,:,k) = & & ( ( xyz_U(:,:,k+1) - xyz_U(:,:,k-1) )**2 & & + ( xyz_V(:,:,k+1) - xyz_V(:,:,k-1) )**2 ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) )**2 end if end do ! Calculation of static stability do k = 1, kmax if ( k == 1 ) then xyz_StatStab(:,:,k) = & & Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) & & * ( xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1) & & - xyz_VirTemp(:,:,k ) / xyz_Exner(:,:,k ) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k ) ) else if ( k == kmax ) then xyz_StatStab(:,:,k) = & & Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) & & * ( xyz_VirTemp(:,:,k ) / xyz_Exner(:,:,k ) & & - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) & & / ( xyz_Height(:,:,k ) - xyz_Height(:,:,k-1) ) else xyz_StatStab(:,:,k) = & & Grav / ( xyz_VirTemp(:,:,k) / xyz_Exner(:,:,k) ) & & * ( xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1) & & - xyz_VirTemp(:,:,k-1) / xyz_Exner(:,:,k-1) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k-1) ) end if end do ! 混合距離の算出 ! Calculate mixing length ! do k = 1, kmax xyz_MixLength(:,:,k) = & & FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) & & / (1.0_DP + FKarm * ( xyz_Height(:,:,k) - xy_SurfHeight ) / MixLengthMax ) end do !!$ ! Limit mixing length (Galperin et al., 1988) and avoid zero !!$ xyz_MixLength = & !!$ & min( xyz_MixLength, & !!$ & 0.53_DP & !!$ & * sqrt( 2.0_DP * xyz_TurKinEneNonZero / max( xyz_StatStab, 1.0d-10 ) ) ) & !!$ & + 1.0d-10 !******************************************************************** ! Gerrity et al. (1994) ! Set flag for using Richardson number do k = 1, kmax do j = 1, jmax do i = 0, imax-1 if ( xyz_DVelDzSq(i,j,k) < CrtlShear**2 ) then xyz_FlagUseRiNum(i,j,k) = .false. else xyz_FlagUseRiNum(i,j,k) = .true. end if end do end do end do ! Calculation of Richardson number do k = 1, kmax do j = 1, jmax do i = 0, imax-1 if ( xyz_FlagUseRiNum(i,j,k) ) then xyz_RiNum(i,j,k) = xyz_StatStab(i,j,k) / xyz_DVelDzSq(i,j,k) else xyz_RiNum(i,j,k) = - 1.0e100_DP end if end do end do end do ! Set flag for selecting an asymptotic equation do k = 1, kmax do j = 1, jmax do i = 0, imax-1 if ( xyz_FlagUseRiNum(i,j,k) ) then if ( CrtlRiNum <= xyz_RiNum(i,j,k) ) then ! Ric <= Ri xyz_FlagTKEAsymptToZero = .true. else ! Ri < Ric xyz_FlagTKEAsymptToZero = .false. end if else if ( xyz_StatStab(i,j,k) > 0.0_DP ) then xyz_FlagTKEAsymptToZero = .true. else xyz_FlagTKEAsymptToZero = .false. end if end if end do end do end do ! Calculation of roos for equations for steady state do k = 1, kmax do j = 1, jmax do i = 0, imax-1 if ( xyz_FlagUseRiNum(i,j,k) ) then ! Calculation with Richardson number MixLength = xyz_MixLength(i,j,k) DVelDzSq = xyz_DVelDzSq (i,j,k) RiNum = xyz_RiNum (i,j,k) ! Beta1 = & & - MixLength**2 * DVelDzSq & & * ( 6.53_DP - 49.0_DP * RiNum ) Beta2 = & & - MixLength**4 * DVelDzSq**2 & & * ( 51.2_DP - 262.7_DP * RiNum ) * RiNum Beta3 = & & MixLength**2 * DVelDzSq & & * ( 5.08_DP + 36.7_DP * RiNum ) Beta4 = & & MixLength**4 * DVelDzSq**2 & & * ( 88.8_DP + 187.4_DP * RiNum ) * RiNum ! xyz_A1(i,j,k) = - Beta1 / 2.0_DP & & + sqrt( Beta1**2 - 4.0_DP * Beta2 ) / 2.0_DP xyz_A2(i,j,k) = - Beta1 / 2.0_DP & & - sqrt( Beta1**2 - 4.0_DP * Beta2 ) / 2.0_DP xyz_R1(i,j,k) = - Beta3 / 2.0_DP & & + sqrt( Beta3**2 - 4.0_DP * Beta4 ) / 2.0_DP xyz_R2(i,j,k) = - Beta3 / 2.0_DP & & - sqrt( Beta3**2 - 4.0_DP * Beta4 ) / 2.0_DP else ! Calculation without Richardson number MixLength = xyz_MixLength(i,j,k) StatStab = xyz_StatStab (i,j,k) ! xyz_A1(i,j,k) = 24.49_DP * MixLength**2 * ( - StatStab ) & & + 18.36_DP * MixLength**2 * abs( - StatStab ) xyz_A2(i,j,k) = 24.49_DP * MixLength**2 * ( - StatStab ) & & - 18.36_DP * MixLength**2 * abs( - StatStab ) xyz_R1(i,j,k) = 18.35_DP * MixLength**2 * ( - StatStab ) & & + 12.22_DP * MixLength**2 * abs( - StatStab ) xyz_R2(i,j,k) = 18.35_DP * MixLength**2 * ( - StatStab ) & & - 12.22_DP * MixLength**2 * abs( - StatStab ) end if end do end do end do ! Set turbulent kinetic energy at current time step ! do k = 1, kmax do j = 1, jmax do i = 0, imax-1 if ( xyz_FlagUseRiNum(i,j,k) ) then if ( xyz_RiNum(i,j,k) < 0.0_DP ) then ! Ri < 0 (or CT > 0) xyz_TKEInit(i,j,k) = & & max( xyz_TurKinEne(i,j,k), xyz_R1(i,j,k) / 2.0_DP ) else if ( xyz_RiNum(i,j,k) < CrtlRiNum ) then ! 0 <= Ri < Ric xyz_TKEInit(i,j,k) = xyz_TurKinEne(i,j,k) else ! Ric <= Ri xyz_TKEInit(i,j,k) = xyz_TurKinEne(i,j,k) end if else if ( xyz_StatStab(i,j,k) < 0.0_DP ) then ! CT > 0 xyz_TKEInit(i,j,k) = & & max( xyz_TurKinEne(i,j,k), xyz_R1(i,j,k) / 2.0_DP ) else ! CT <= 0 xyz_TKEInit(i,j,k) = xyz_TurKinEne(i,j,k) end if end if xyz_TKEInit(i,j,k) = xyz_TKEInit(i,j,k) + TurKinEneOffset end do end do end do ! xyz_TKEx2 = 2.0_DP * xyz_TKEInit xyz_SqrtTKEx2 = sqrt( xyz_TKEx2 ) do k = 1, kmax do j = 1, jmax do i = 0, imax-1 if ( xyz_FlagTKEAsymptToZero(i,j,k) ) then xyz_SqrtA1(i,j,k) = 1.0e100_DP else xyz_SqrtA1(i,j,k) = sqrt( xyz_A1(i,j,k) ) end if end do end do end do ! do k = 1, kmax do j = 1, jmax do i = 0, imax-1 if ( xyz_FlagTKEAsymptToZero(i,j,k) ) then ! Use equations asymptotic to zero TKE, Eq. (10) ! Eq. (11) BetaSq = & & ( xyz_TKEx2(i,j,k) - xyz_A1(i,j,k)**2 ) & & * ( xyz_TKEx2(i,j,k) - xyz_A2(i,j,k)**2 ) & & / ( MYConstB1 * xyz_MixLength(i,j,k) & & * ( xyz_TKEx2(i,j,k) - xyz_R1(i,j,k) ) & & * ( xyz_TKEx2(i,j,k) - xyz_R2(i,j,k) ) ) ! Eq. (10) xyz_TKETentative(i,j,k) = & & xyz_TKEInit(i,j,k) / 2.0_DP & & / ( 1.0_DP & & + xyz_SqrtTKEx2(i,j,k) * BetaSq & & * ( 2.0_DP * DelTime ) )**2 else ! Use equations asymptotic to certain TKE, Eq. (8) ! Eq. (7b) BetaSq = & & xyz_TKEx2(i,j,k) & & * ( xyz_TKEx2(i,j,k) - xyz_A2(i,j,k)**2 ) & & / ( MYConstB1 * xyz_MixLength(i,j,k) & & * ( xyz_TKEx2(i,j,k) - xyz_R1(i,j,k) ) & & * ( xyz_TKEx2(i,j,k) - xyz_R2(i,j,k) ) ) ! Eq. (9) Alpha = & & ( xyz_SqrtTKEx2(i,j,k) - xyz_SqrtA1(i,j,k) ) & & / ( xyz_SqrtTKEx2(i,j,k) + xyz_SqrtA1(i,j,k) ) & & * exp( - 2.0_DP * xyz_SqrtA1(i,j,k) * BetaSq & & * ( 2.0_DP * DelTime ) ) Alpha = max( min( Alpha, 1.0_DP - Epsilon ), -1.0_DP ) ! Eq. (8) xyz_TKETentative(i,j,k) = & & xyz_A1(i,j,k) / 2.0_DP & & * ( ( 1.0_DP + Alpha ) / ( 1.0_DP - Alpha ) )**2 end if end do end do end do !******************************************************************** ! Set turbulent kinetic energy for diffusion calculation ! !!$ xyz_TurKinEneNonZero = xyz_TurKinEne + TurKinEneOffset xyz_TurKinEneNonZero = xyz_TKETentative + TurKinEneOffset xyz_Gh = - xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) & & * xyz_StatStab ! Actually, xyz_Gm is not used below. xyz_Gm = xyz_MixLength**2 / ( 2.0_DP * xyz_TurKinEneNonZero ) & & * xyz_DVelDzSq ! Limit Gh (Galperin et al., 1988) !!$ GhMin = - 0.53d0**2 !!$ GhMax = 1.0_DP & !!$ & / ( MYConstA2 & !!$ & * ( 12.0_DP * MYConstA1 + MYConstB1 + 3.0_DP * MYConstB2 ) ) !!$ xyz_Gh = max( GhMin, min( xyz_Gh, GhMax ) ) xyz_Sh = & & MYConstA2 * ( 1.0_DP - 6.0_DP * MYConstA1 / MYConstB1 ) & & / ( 1.0_DP & & - 3.0_DP * MYConstA2 * ( 6.0_DP * MYConstA1 + MYConstB2 ) * xyz_Gh ) xyz_Sm = & & ( & & MYConstA1 * ( 1.0_DP - 3.0_DP * MYConstC1 & & - 6.0_DP * MYConstA1 / MYConstB1 ) & & + 9.0_DP * MYConstA1 * ( 2.0_DP * MYConstA1 + MYConstA2 ) & & * xyz_Gh * xyz_Sh & & ) & & / ( 1.0_DP & & - 9.0_DP * MYConstA1 * MYConstA2 * xyz_Gh ) ! 拡散係数の計算 ! Calculation of diffusion coefficient ! xyz_VelDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sm xyz_TempDiffCoef = xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * xyz_Sh ! do k = 0, kmax if ( ( k == 0 ) .or. ( k == kmax ) ) then xyr_VelDiffCoef (:,:,k) = 0.0_DP xyr_TempDiffCoef(:,:,k) = 0.0_DP else xyr_VelDiffCoef (:,:,k) = ( xyz_VelDiffCoef (:,:,k) + xyz_VelDiffCoef (:,:,k+1) ) / 2.0_DP xyr_TempDiffCoef(:,:,k) = ( xyz_TempDiffCoef(:,:,k) + xyz_TempDiffCoef(:,:,k+1) ) / 2.0_DP end if end do ! do k = 1, kmax-1 do j = 1, jmax do i = 0, imax-1 xyr_VelDiffCoef(i,j,k) = & & max( min( xyr_VelDiffCoef(i,j,k), VelDiffCoefMax ), & & VelDiffCoefMin & & ) xyr_TempDiffCoef(i,j,k) = & & max( min( xyr_TempDiffCoef(i,j,k), TempDiffCoefMax ), & & TempDiffCoefMin & & ) end do end do end do ! xyr_QMixDiffCoef = xyr_TempDiffCoef ! 輸送係数とフラックスの計算 ! Calculate transfer coefficient and flux ! call VDiffusionCalcFlux( & & xyz_U, xyz_V, xyzf_QMix, & ! (in) & xyz_Temp, xyr_VirTemp, xyr_Press, & ! (in) & xyz_Height, xyz_Exner, xyr_Exner, & ! (in) & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, & ! (in) & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux & ! (out) & ) ! Calculate tendency of turbulent kinetic energy ! Set diffusion coefficient for turbulent kinetic energy xyz_TurKinEneDiffCoef = & & xyz_MixLength * sqrt( 2.0_DP * xyz_TurKinEneNonZero ) * Stke ! do k = 0, kmax if ( k == 0 ) then xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,1) else if ( k == kmax ) then xyr_TurKinEneDiffCoef(:,:,k) = xyz_TurKinEneDiffCoef(:,:,kmax) else xyr_TurKinEneDiffCoef(:,:,k) = & & ( xyz_TurKinEneDiffCoef(:,:,k) + xyz_TurKinEneDiffCoef(:,:,k+1) ) & & / 2.0_DP end if end do ! Calculate turbulent kinetic energy at lower boundary ! xy_FricVelSq = sqrt( xy_SurfMomFluxX**2 + xy_SurfMomFluxY**2 ) & & / ( xyr_Press(:,:,0) / ( GasRDry * xyr_VirTemp(:,:,0) ) ) xy_TurKinEneAtLB = MYConstB1**(2.0_DP/3.0_DP) / 2.0_DP * xy_FricVelSq xy_TurKinEneAtLB = xy_TurKinEneAtLB + TurKinEneOffset ! Calculate transfer coefficient and flux of turbulent kinetic energy ! ! When transfer coefficient at lower boundary is calculated, ! diffusion coefficient at mid-point of 1st layer is used. ! In addition, transfer coefficient at upper boundary is assumed ! to be zero. k = 0 xyr_TurKinEneTransCoef(:,:,k) = & & xyz_TurKinEneDiffCoef(:,:,1) & & * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) & & / ( xyz_Height(:,:,1) - xy_SurfHeight ) do k = 1, kmax-1 xyr_TurKinEneTransCoef(:,:,k) = & & xyr_TurKinEneDiffCoef(:,:,k) & & * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) end do k = kmax xyr_TurKinEneTransCoef(:,:,k) = 0.0_DP ! do k = 1, kmax-1 xyr_TurKinEneFlux(:,:,k) = & & - xyr_TurKinEneTransCoef(:,:,k) & & * ( xyz_TurKinEneNonZero(:,:,k+1) - xyz_TurKinEneNonZero(:,:,k) ) end do k = 0 xyr_TurKinEneFlux(:,:,k) = & & - xyr_TurKinEneTransCoef(:,:,k) & & * ( xyz_TurKinEneNonZero(:,:,k+1) - xy_TurKinEneAtLB ) k = kmax xyr_TurKinEneFlux(:,:,k) = 0.0_DP ! ! Construct implicit matrix from transfer coefficient of vertical ! diffusion scheme (turbulent kinetic energy) ! k = 1 xyza_TurKinEneMtx(:,:,k,-1) = 0.0_DP xyza_TurKinEneMtx(:,:,k, 0) = & & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) & & + xyr_TurKinEneTransCoef(:,:,k-1) & & + xyr_TurKinEneTransCoef(:,:,k ) xyza_TurKinEneMtx(:,:,k, 1) = & & - xyr_TurKinEneTransCoef(:,:,k ) ! do k = 2, kmax-1 xyza_TurKinEneMtx(:,:,k,-1) =& & - xyr_TurKinEneTransCoef(:,:,k-1) xyza_TurKinEneMtx(:,:,k, 0) = & & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) & & + xyr_TurKinEneTransCoef(:,:,k-1) & & + xyr_TurKinEneTransCoef(:,:,k ) xyza_TurKinEneMtx(:,:,k, 1) = & & - xyr_TurKinEneTransCoef(:,:,k ) end do ! k = kmax xyza_TurKinEneMtx(:,:,k,-1) =& & - xyr_TurKinEneTransCoef(:,:,k-1) xyza_TurKinEneMtx(:,:,k, 0) = & & - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) & & + xyr_TurKinEneTransCoef(:,:,k-1) & & + xyr_TurKinEneTransCoef(:,:,k ) xyza_TurKinEneMtx(:,:,k, 1) = 0.0_DP do k = 1, kmax xyz_TurKinEneVec(:,:,k) = & & - ( xyr_TurKinEneFlux(:,:,k) - xyr_TurKinEneFlux(:,:,k-1) ) end do ! ! Solve simultaneous linear equations by use of LU decomposition technique ! xyza_TurKinEneLUMtx = xyza_TurKinEneMtx ! call PhyImplLUDecomp3( & & xyza_TurKinEneLUMtx, & ! (inout) & imax * jmax, kmax & ! (in) & ) xyz_DelTurKinEneLUVec = xyz_TurKinEneVec ! call PhyImplLUSolve3( & & xyz_DelTurKinEneLUVec, & ! (inout) & xyza_TurKinEneLUMtx, & ! (in) & 1, imax * jmax , kmax & ! (in) & ) !!$ xyz_DTurKinEneDt = xyz_DelTurKinEneLUVec / ( 2.0_DP * DelTime ) xyz_TKETentative = xyz_TKETentative + xyz_DelTurKinEneLUVec xyz_DTurKinEneDt = & & ( xyz_TKETentative - xyz_TurKinEne ) / ( 2.0_DP * DelTime ) !!$ write( 6, * ) TimeN, iloop !!$ write( 6, * ) xyz_TurKinEne(0,jmax/2+1,1:4) !!$ write( 6, * ) xyz_TempDiffCoef(0,jmax/2+1,1:4) !!$ write( 6, * ) xyz_DTurKinEneDt(0,jmax/2+1,1:4) ! 拡散係数の出力 ! Output diffusion coefficients ! ! 拡散係数出力 ! Diffusion coeffficients output ! call HistoryAutoPut( TimeN, 'VelDiffCoef', xyr_VelDiffCoef ) call HistoryAutoPut( TimeN, 'TempDiffCoef', xyr_TempDiffCoef ) call HistoryAutoPut( TimeN, 'QVapDiffCoef', xyr_QMixDiffCoef ) call HistoryAutoPut( TimeN, 'MixLength' , xyz_MixLength ) ! 計算時間計測一時停止 ! Pause measurement of computation time ! call TimesetClockStop( module_name ) end subroutine VDiffusionMY25GBT94 !-------------------------------------------------------------------------------------- subroutine VDiffusion( & & xyz_U, xyz_V, xyzf_QMix, & ! (in) & xyz_Temp, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, & ! (in) & xy_SurfHeight, & ! (in) & xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, & ! (in) & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (out) & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef & ! (out) & ) ! ! 鉛直拡散フラックスを計算します. ! ! Vertical diffusion flux is calculated by use of MY2 model. ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: & & FKarm, & ! $ k $ . ! カルマン定数. ! Karman constant & Grav ! $ g $ [m s-2]. ! 重力加速度. ! Gravitational acceleration ! 時刻管理 ! Time control ! use timeset, only: & & TimeN, & ! ステップ $ t $ の時刻. Time of step $ t $. & TimesetClockStart, TimesetClockStop ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! 宣言文 ; Declaration statements ! implicit none real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax) ! $ u $ . 東西風速. Eastward wind real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax) ! $ v $ . 南北風速. Northward wind real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) ! $ q $ . 質量混合比. Mass mixing ratio real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax) ! $ T $ . 温度. Temperature real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{T} $ . 温度 (半整数レベル). ! Temperature (half level) real(DP), intent(in):: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax) ! $ T_v $ . 仮温度. Virtual temperature real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{T}_v $ . 仮温度 (半整数レベル). ! Virtual temperature (half level) real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{p} $ . 気圧 (半整数レベル). ! Air pressure (half level) real(DP), intent(in):: xy_SurfHeight (0:imax-1,1:jmax) ! $ z_s $ . 地表面高度. ! Surface height. real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax) ! 高度 (整数レベル). ! Height (full level) real(DP), intent(in):: xyr_Height (0:imax-1, 1:jmax, 0:kmax) ! 高度 (半整数レベル). ! Height (half level) real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax) ! Exner 関数 (整数レベル). ! Exner function (full level) real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax) ! Exner 関数 (半整数レベル). ! Exner function (half level) real(DP), intent(out):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax) ! 東西方向運動量フラックス. ! Eastward momentum flux real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax) ! 南北方向運動量フラックス. ! Northward momentum flux real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax) ! 熱フラックス. ! Heat flux real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) ! 質量フラックス. ! Mass flux of compositions real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP), intent(out):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：温度. ! Transfer coefficient: temperature real(DP), intent(out):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：比湿. ! Diffusion coefficient: specific humidity ! 作業変数 ! Work variables ! real(DP) :: xyr_DVelDz (0:imax-1, 1:jmax, 0:kmax) ! $ \DD{|\Dvect{v}|}{z} $ real(DP) :: xyr_BulkRiNum (0:imax-1, 1:jmax, 0:kmax) ! バルク $ R_i $ 数. ! Bulk $ R_i $ integer:: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer:: n ! 組成方向に回る DO ループ用作業変数 ! Work variables for DO loop in dimension of constituents ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. vdiffusion_my_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 計算時間計測開始 ! Start measurement of computation time ! call TimesetClockStart( module_name ) ! 拡散係数の計算 ! Calculation of diffusion coefficient ! if ( FlagConstDiffCoef ) then xyr_VelDiffCoef (:,:,0 ) = 0.0_DP xyr_VelDiffCoef (:,:,1:kmax-1) = ConstDiffCoefM xyr_VelDiffCoef (:,:,kmax ) = 0.0_DP xyr_TempDiffCoef(:,:,0 ) = 0.0_DP xyr_TempDiffCoef(:,:,1:kmax-1) = ConstDiffCoefH xyr_TempDiffCoef(:,:,kmax ) = 0.0_DP xyr_QMixDiffCoef(:,:,0 ) = 0.0_DP xyr_QMixDiffCoef(:,:,1:kmax-1) = ConstDiffCoefH xyr_QMixDiffCoef(:,:,kmax ) = 0.0_DP else ! バルク $ R_i $ 数算出 ! Calculate bulk $ R_i $ ! xyr_DVelDz(:,:,0) = 0.0_DP xyr_DVelDz(:,:,kmax) = 0.0_DP xyr_BulkRiNum(:,:,0) = 0.0_DP xyr_BulkRiNum(:,:,kmax) = 0.0_DP do k = 1, kmax-1 xyr_DVelDz(:,:,k) = & & sqrt( max( SquareVelMin , & & ( xyz_U(:,:,k+1) - xyz_U(:,:,k) )**2 & & + ( xyz_V(:,:,k+1) - xyz_V(:,:,k) )**2 & & ) & & ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) xyr_BulkRiNum(:,:,k) = & & Grav / ( xyr_VirTemp(:,:,k) / xyr_Exner(:,:,k) ) & & * ( xyz_VirTemp(:,:,k+1) / xyz_Exner(:,:,k+1) & & - xyz_VirTemp(:,:,k ) / xyz_Exner(:,:,k ) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) & & / xyr_DVelDz(:,:,k)**2 xyr_BulkRiNum(:,:,k) = max( xyr_BulkRiNum(:,:,k) , BulkRiNumMin ) end do ! 拡散係数の計算 ! Calculate diffusion coefficients ! call VDiffCoefficient( & & xy_SurfHeight, & ! (in) & xyr_Height, xyr_DVelDz, xyr_BulkRiNum, & ! (in) & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef & ! (out) & ) end if ! 浅い積雲対流 ! Shallow cumulus convection ! ! (AGCM5 から導入予定) ! 拡散係数の出力 ! Output diffusion coefficients ! ! (上記の「浅い積雲対流」導入後に作成) ! 拡散係数出力 ! Diffusion coeffficients output ! call HistoryAutoPut( TimeN, 'VelDiffCoef', xyr_VelDiffCoef ) call HistoryAutoPut( TimeN, 'TempDiffCoef', xyr_TempDiffCoef ) call HistoryAutoPut( TimeN, 'QVapDiffCoef', xyr_QMixDiffCoef ) ! 輸送係数とフラックスの計算 ! Calculate transfer coefficient and flux ! call VDiffusionCalcFlux( & & xyz_U, xyz_V, xyzf_QMix, & ! (in) & xyz_Temp, xyr_VirTemp, xyr_Press, & ! (in) & xyz_Height, xyz_Exner, xyr_Exner, & ! (in) & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, & ! (in) & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux & ! (out) & ) ! 計算時間計測一時停止 ! Pause measurement of computation time ! call TimesetClockStop( module_name ) end subroutine VDiffusion !-------------------------------------------------------------------------------------- subroutine VDiffusionCalcFlux( & & xyz_U, xyz_V, xyzf_QMix, & ! (in) & xyz_Temp, xyr_VirTemp, xyr_Press, & ! (in) & xyz_Height, xyz_Exner, xyr_Exner, & ! (in) & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, & ! (in) & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux & ! (out) & ) ! ! 鉛直拡散フラックスを計算します. ! ! Vertical diffusion flux is calculated. ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: & & FKarm, & ! $ k $ . ! カルマン定数. ! Karman constant & GasRDry, & ! $ R $ [J kg-1 K-1]. ! 乾燥大気の気体定数. ! Gas constant of air & CpDry ! $ C_p $ [J kg-1 K-1]. ! 乾燥大気の定圧比熱. ! Specific heat of air at constant pressure ! 時刻管理 ! Time control ! use timeset, only: & & TimeN, & ! ステップ $ t $ の時刻. Time of step $ t $. & TimesetClockStart, TimesetClockStop ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! 宣言文 ; Declaration statements ! implicit none real(DP), intent(in):: xyz_U (0:imax-1, 1:jmax, 1:kmax) ! $ u $ . 東西風速. Eastward wind real(DP), intent(in):: xyz_V (0:imax-1, 1:jmax, 1:kmax) ! $ v $ . 南北風速. Northward wind real(DP), intent(in):: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) ! $ q $ . 質量混合比. Mass mixing ratio real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax) ! $ T $ . 温度. Temperature real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{T}_v $ . 仮温度 (半整数レベル). ! Virtual temperature (half level) real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{p} $ . 気圧 (半整数レベル). ! Air pressure (half level) real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax) ! 高度 (整数レベル). ! Height (full level) real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax) ! Exner 関数 (整数レベル). ! Exner function (full level) real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax) ! Exner 関数 (半整数レベル). ! Exner function (half level) real(DP), intent(in):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP), intent(in):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：温度. ! Transfer coefficient: temperature real(DP), intent(in):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：比湿. ! Diffusion coefficient: specific humidity real(DP), intent(out):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax) ! 東西方向運動量フラックス. ! Eastward momentum flux real(DP), intent(out):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax) ! 南北方向運動量フラックス. ! Northward momentum flux real(DP), intent(out):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax) ! 熱フラックス. ! Heat flux real(DP), intent(out):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) ! 質量フラックス. ! Mass flux of compositions ! 作業変数 ! Work variables ! real(DP) :: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax) ! 輸送係数：運動量. ! Transfer coefficient: velocity real(DP) :: xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) ! 輸送係数：温度. ! Transfer coefficient: temperature real(DP) :: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) ! 輸送係数：質量. ! Transfer coefficient: mass (composition) integer:: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer:: n ! 組成方向に回る DO ループ用作業変数 ! Work variables for DO loop in dimension of constituents ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. vdiffusion_my_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 輸送係数の計算 ! Calculate transfer coefficient ! xyr_VelTransCoef (:,:,0) = 0.0_DP xyr_VelTransCoef (:,:,kmax) = 0.0_DP xyr_TempTransCoef(:,:,0) = 0.0_DP xyr_TempTransCoef(:,:,kmax) = 0.0_DP xyr_QMixTransCoef(:,:,0) = 0.0_DP xyr_QMixTransCoef(:,:,kmax) = 0.0_DP do k = 1, kmax-1 xyr_VelTransCoef(:,:,k) = & & xyr_VelDiffCoef(:,:,k) & & * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) xyr_TempTransCoef(:,:,k) = & & xyr_TempDiffCoef(:,:,k) & & * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) xyr_QMixTransCoef(:,:,k) = & & xyr_QMixDiffCoef(:,:,k) & & * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) end do ! フラックスの計算 ! Calculate fluxes ! xyr_MomFluxX(:,:,0) = 0.0_DP xyr_MomFluxX(:,:,kmax) = 0.0_DP xyr_MomFluxY(:,:,0) = 0.0_DP xyr_MomFluxY(:,:,kmax) = 0.0_DP xyr_HeatFlux(:,:,0) = 0.0_DP xyr_HeatFlux(:,:,kmax) = 0.0_DP do n = 1, ncmax xyrf_QMixFlux(:,:,0,n) = 0.0_DP xyrf_QMixFlux(:,:,kmax,n) = 0.0_DP end do do k = 1, kmax-1 xyr_MomFluxX(:,:,k) = & & - xyr_VelTransCoef(:,:,k) * ( xyz_U(:,:,k+1) - xyz_U(:,:,k) ) xyr_MomFluxY(:,:,k) = & & - xyr_VelTransCoef(:,:,k) * ( xyz_V(:,:,k+1) - xyz_V(:,:,k) ) xyr_HeatFlux(:,:,k) = & & - CpDry * xyr_TempTransCoef(:,:,k) * xyr_Exner(:,:,k) & & * ( xyz_Temp(:,:,k+1) / xyz_Exner(:,:,k+1) & & - xyz_Temp(:,:,k) / xyz_Exner(:,:,k) ) end do do n = 1, ncmax do k = 1, kmax-1 xyrf_QMixFlux(:,:,k,n) = & & - xyr_QMixTransCoef(:,:,k) & & * ( xyzf_QMix(:,:,k+1,n) - xyzf_QMix(:,:,k,n) ) end do end do end subroutine VDiffusionCalcFlux !-------------------------------------------------------------------------------------- subroutine VDiffusionCalcFlux1D( & & z_U, z_V, zf_QMix, & ! (in) & z_Temp, r_VirTemp, r_Press, & ! (in) & z_Height, z_Exner, r_Exner, & ! (in) & r_VelDiffCoef, r_TempDiffCoef, r_QMixDiffCoef, & ! (in) & r_MomFluxX, r_MomFluxY, r_HeatFlux, rf_QMixFlux & ! (out) & ) ! ! 鉛直拡散フラックスを計算します. ! ! Vertical diffusion flux is calculated. ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: & & FKarm, & ! $ k $ . ! カルマン定数. ! Karman constant & GasRDry, & ! $ R $ [J kg-1 K-1]. ! 乾燥大気の気体定数. ! Gas constant of air & CpDry ! $ C_p $ [J kg-1 K-1]. ! 乾燥大気の定圧比熱. ! Specific heat of air at constant pressure ! 時刻管理 ! Time control ! use timeset, only: & & TimeN, & ! ステップ $ t $ の時刻. Time of step $ t $. & TimesetClockStart, TimesetClockStop ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! 宣言文 ; Declaration statements ! implicit none real(DP), intent(in):: z_U (1:kmax) ! $ u $ . 東西風速. Eastward wind real(DP), intent(in):: z_V (1:kmax) ! $ v $ . 南北風速. Northward wind real(DP), intent(in):: zf_QMix(1:kmax, 1:ncmax) ! $ q $ . 質量混合比. Mass mixing ratio real(DP), intent(in):: z_Temp (1:kmax) ! $ T $ . 温度. Temperature real(DP), intent(in):: r_VirTemp (0:kmax) ! $ \hat{T}_v $ . 仮温度 (半整数レベル). ! Virtual temperature (half level) real(DP), intent(in):: r_Press (0:kmax) ! $ \hat{p} $ . 気圧 (半整数レベル). ! Air pressure (half level) real(DP), intent(in):: z_Height (1:kmax) ! 高度 (整数レベル). ! Height (full level) real(DP), intent(in):: z_Exner (1:kmax) ! Exner 関数 (整数レベル). ! Exner function (full level) real(DP), intent(in):: r_Exner (0:kmax) ! Exner 関数 (半整数レベル). ! Exner function (half level) real(DP), intent(in):: r_VelDiffCoef (0:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP), intent(in):: r_TempDiffCoef (0:kmax) ! 拡散係数：温度. ! Transfer coefficient: temperature real(DP), intent(in):: r_QMixDiffCoef (0:kmax) ! 拡散係数：比湿. ! Diffusion coefficient: specific humidity real(DP), intent(out):: r_MomFluxX (0:kmax) ! 東西方向運動量フラックス. ! Eastward momentum flux real(DP), intent(out):: r_MomFluxY (0:kmax) ! 南北方向運動量フラックス. ! Northward momentum flux real(DP), intent(out):: r_HeatFlux (0:kmax) ! 熱フラックス. ! Heat flux real(DP), intent(out):: rf_QMixFlux(0:kmax, 1:ncmax) ! 質量フラックス. ! Mass flux of compositions ! 作業変数 ! Work variables ! real(DP) :: r_VelTransCoef (0:kmax) ! 輸送係数：運動量. ! Transfer coefficient: velocity real(DP) :: r_TempTransCoef(0:kmax) ! 輸送係数：温度. ! Transfer coefficient: temperature real(DP) :: r_QMixTransCoef(0:kmax) ! 輸送係数：質量. ! Transfer coefficient: mass (composition) integer:: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer:: n ! 組成方向に回る DO ループ用作業変数 ! Work variables for DO loop in dimension of constituents ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. vdiffusion_my_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 輸送係数の計算 ! Calculate transfer coefficient ! r_VelTransCoef (0) = 0.0_DP r_VelTransCoef (kmax) = 0.0_DP r_TempTransCoef(0) = 0.0_DP r_TempTransCoef(kmax) = 0.0_DP r_QMixTransCoef(0) = 0.0_DP r_QMixTransCoef(kmax) = 0.0_DP do k = 1, kmax-1 r_VelTransCoef(k) = & & r_VelDiffCoef(k) & & * r_Press(k) / ( GasRDry * r_VirTemp(k) ) & & / ( z_Height(k+1) - z_Height(k) ) r_TempTransCoef(k) = & & r_TempDiffCoef(k) & & * r_Press(k) / ( GasRDry * r_VirTemp(k) ) & & / ( z_Height(k+1) - z_Height(k) ) r_QMixTransCoef(k) = & & r_QMixDiffCoef(k) & & * r_Press(k) / ( GasRDry * r_VirTemp(k) ) & & / ( z_Height(k+1) - z_Height(k) ) end do ! フラックスの計算 ! Calculate fluxes ! r_MomFluxX(0) = 0.0_DP r_MomFluxX(kmax) = 0.0_DP r_MomFluxY(0) = 0.0_DP r_MomFluxY(kmax) = 0.0_DP r_HeatFlux(0) = 0.0_DP r_HeatFlux(kmax) = 0.0_DP do n = 1, ncmax rf_QMixFlux(0,n) = 0.0_DP rf_QMixFlux(kmax,n) = 0.0_DP end do do k = 1, kmax-1 r_MomFluxX(k) = & & - r_VelTransCoef(k) * ( z_U(k+1) - z_U(k) ) r_MomFluxY(k) = & & - r_VelTransCoef(k) * ( z_V(k+1) - z_V(k) ) r_HeatFlux(k) = & & - CpDry * r_TempTransCoef(k) * r_Exner(k) & & * ( z_Temp(k+1) / z_Exner(k+1) & & - z_Temp(k) / z_Exner(k) ) end do do n = 1, ncmax do k = 1, kmax-1 rf_QMixFlux(k,n) = & & - r_QMixTransCoef(k) & & * ( zf_QMix(k+1,n) - zf_QMix(k,n) ) end do end do end subroutine VDiffusionCalcFlux1D !-------------------------------------------------------------------------------------- subroutine VDiffCoefficient( & & xy_SurfHeight, & ! (in) & xyr_Height, xyr_DVelDz, xyr_BulkRiNum, & ! (in) & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef & ! (out) & ) ! ! 鉛直拡散フラックスを計算します. ! ! Vertical diffusion flux is calculated. ! ! モジュール引用 ; USE statements ! ! 時刻管理 ! Time control ! use timeset, only: & & TimeN, & ! ステップ $ t $ の時刻. Time of step $ t $. & TimesetClockStart, TimesetClockStop ! 物理定数設定 ! Physical constants settings ! use constants, only: & & FKarm ! $ k $ . ! カルマン定数. ! Karman constant ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! 宣言文 ; Declaration statements ! implicit none real(DP), intent(in):: xy_SurfHeight (0:imax-1,1:jmax) ! $ z_s $ . 地表面高度. ! Surface height. real(DP), intent(in):: xyr_Height (0:imax-1, 1:jmax, 0:kmax) ! 高度 (半整数レベル). ! Height (half level) real(DP), intent(in):: xyr_DVelDz (0:imax-1, 1:jmax, 0:kmax) ! $ \DD{|\Dvect{v}|}{z} $ real(DP), intent(in):: xyr_BulkRiNum (0:imax-1, 1:jmax, 0:kmax) ! バルク $ R_i $ 数. ! Bulk $ R_i $ real(DP), intent(out):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP), intent(out):: xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：温度. ! Transfer coefficient: temperature real(DP), intent(out):: xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：質量 ! Diffusion coefficient: mass of constituents ! 作業変数 ! Work variables ! real(DP):: xyr_FluxRiNum (0:imax-1, 1:jmax, 0:kmax) ! フラックス $ R_i $ 数. ! Flux $ R_i $ number real(DP):: xyr_Sh (0:imax-1, 1:jmax, 0:kmax) ! $ S_h $ (温度, 比湿). ! $ S_h $ (temperature, specific humidity) real(DP):: xyr_Sm (0:imax-1, 1:jmax, 0:kmax) ! $ S_m $ (運動量). ! $ S_m $ (momentum) real(DP):: xyr_MixLength (0:imax-1, 1:jmax, 0:kmax) ! 混合距離. ! Mixing length real(DP):: Alpha1, Alpha2 real(DP):: Beta1, Beta2, Beta3, Beta4 real(DP):: Gamma1, Gamma2 real(DP):: CrtlFluxRiNum real(DP):: xyr_TurKinEne(0:imax-1, 1:jmax, 0:kmax) ! ! Turbulent kinetic energy integer:: i ! 経度方向に回る DO ループ用作業変数 ! Work variables for DO loop in longitude integer:: j ! 緯度方向に回る DO ループ用作業変数 ! Work variables for DO loop in latitude integer:: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. vdiffusion_my_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 定数計算 ! Calculate constants ! Gamma1 = ( 1.0_DP / 3.0_DP ) - ( 2.0_DP * MYConstA1 / MYConstB1 ) Gamma2 = ( MYConstB2 / MYConstB1 ) & & + ( 6.0_DP * MYConstA1 / MYConstB1 ) Alpha1 = 3.0_DP * MYConstA2 * Gamma1 Alpha2 = 3.0_DP * MYConstA2 * ( Gamma1 + Gamma2 ) Beta1 = MYConstA1 * MYConstB1 * ( Gamma1 - MYConstC1 ) Beta2 = MYConstA1 & & * ( MYConstB1 * ( Gamma1 - MYConstC1 ) & & + 6.0_DP * MYConstA1 & & + 3.0_DP * MYConstA2 ) Beta3 = MYConstA2 * MYConstB1 * Gamma1 Beta4 = MYConstA2 * ( MYConstB1 * ( Gamma1 + Gamma2 ) & & - 3.0_DP * MYConstA1 ) CrtlFluxRiNum = Gamma1 / ( Gamma1 + Gamma2 ) ! フラックス $ R_i $ 数の算出 ! Calculate flux $ R_i $ number ! xyr_FluxRiNum = & & ( Beta1 & & + Beta4 * xyr_BulkRiNum & & - sqrt( ( Beta1 + Beta4 * xyr_BulkRiNum )**2 & & - 4.0_DP * Beta2 * Beta3 * xyr_BulkRiNum & & ) & & ) / ( 2.0_DP * Beta2 ) ! $ \tilde{S_h} $ と $ \tilde{S_m} $ の算出 ! Calculate $ \tilde{S_h} $ and $ \tilde{S_m} $ ! xyr_Sh(:,:,kmax) = 0.0_DP xyr_Sm(:,:,kmax) = 0.0_DP do k = 0, kmax-1 do i = 0, imax-1 do j = 1, jmax if ( xyr_FluxRiNum(i,j,k) < CrtlFluxRiNum ) then xyr_Sh(i,j,k) = & & ( Alpha1 - Alpha2 * xyr_FluxRiNum(i,j,k) ) & & / ( 1.0_DP - 1.0_DP * xyr_FluxRiNum(i,j,k) ) xyr_Sm(i,j,k) = & & ( Beta1 - Beta2 * xyr_FluxRiNum(i,j,k) ) & & / ( Beta3 - Beta4 * xyr_FluxRiNum(i,j,k) ) & & * xyr_Sh(i,j,k) xyr_Sh(i,j,k) = max( xyr_Sh(i,j,k), ShMin ) xyr_Sm(i,j,k) = max( xyr_Sm(i,j,k), SmMin ) else xyr_Sh(i,j,k) = ShMin xyr_Sm(i,j,k) = SmMin end if end do end do end do ! 混合距離の算出 ! Calculate mixing length ! do k = 0, kmax xyr_MixLength(:,:,k) = & & FKarm * ( xyr_Height(:,:,k) - xy_SurfHeight(:,:) ) & & / (1.0_DP + FKarm * ( xyr_Height(:,:,k) - xy_SurfHeight(:,:) ) / MixLengthMax ) end do ! 拡散係数の算出 ! Calculate diffusion constants ! xyr_VelDiffCoef = & & xyr_MixLength**2 & & * xyr_DVelDz & & * sqrt ( MYConstB1 * ( 1.0_DP - xyr_FluxRiNum ) * xyr_Sm ) & & * xyr_Sm xyr_TempDiffCoef = & & xyr_MixLength ** 2 & & * xyr_DVelDz & & * sqrt ( MYConstB1 * ( 1.0_DP - xyr_FluxRiNum ) * xyr_Sm ) & & * xyr_Sh do k = 0, kmax-1 do i = 0, imax-1 do j = 1, jmax xyr_VelDiffCoef(i,j,k) = & & max( min( xyr_VelDiffCoef(i,j,k), VelDiffCoefMax ), & & VelDiffCoefMin & & ) xyr_TempDiffCoef(i,j,k) = & & max( min( xyr_TempDiffCoef(i,j,k), TempDiffCoefMax ), & & TempDiffCoefMin & & ) end do end do end do xyr_QMixDiffCoef = xyr_TempDiffCoef xyr_VelDiffCoef (:,:,0) = 0.0_DP xyr_VelDiffCoef (:,:,kmax) = 0.0_DP xyr_TempDiffCoef(:,:,0) = 0.0_DP xyr_TempDiffCoef(:,:,kmax) = 0.0_DP xyr_QMixDiffCoef(:,:,0) = 0.0_DP xyr_QMixDiffCoef(:,:,kmax) = 0.0_DP ! Calculation of turbulent kinetic energy ! Turbulent kinetic energy is diagnosed. xyr_TurKinEne = MYConstB1 * xyr_MixLength**2 & & * ( 1.0_DP - xyr_FluxRiNum ) * xyr_Sm * xyr_DVelDz**2 & & / 2.0_DP xyr_TurKinEne(:,:,0 ) = 0.0_DP xyr_TurKinEne(:,:,kmax) = 0.0_DP ! ヒストリデータ出力 ! History data output ! call HistoryAutoPut( TimeN, 'TurKinEne', xyr_TurKinEne ) end subroutine VDiffCoefficient !-------------------------------------------------------------------------------------- subroutine VDiffusionExpTendency( & & xyr_Press, & ! (in ) & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (in ) optional & xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt & ! (out) optional & ) ! ! 時間変化率の計算を行います. ! ! Calculate tendencies. ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: & & Grav, & ! $ g $ [m s-2]. ! 重力加速度. ! Gravitational acceleration & CpDry ! $ C_p $ [J kg-1 K-1]. ! 乾燥大気の定圧比熱. ! Specific heat of air at constant pressure ! 宣言文 ; Declaration statements ! implicit none real(DP), intent(in ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{p} $ . 気圧 (半整数レベル). ! Air pressure (half level) real(DP), intent(in ), optional :: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax) ! 東西方向運動量フラックス. ! Eastward momentum flux real(DP), intent(in ), optional :: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax) ! 南北方向運動量フラックス. ! Northward momentum flux real(DP), intent(in ), optional :: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax) ! 熱フラックス. ! Heat flux real(DP), intent(in ), optional :: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) ! 比湿フラックス. ! Specific humidity flux real(DP), intent(out), optional :: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax) ! $ \DP{u}{t} $ . 東西風速変化. ! Eastward wind tendency real(DP), intent(out), optional :: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax) ! $ \DP{v}{t} $ . 南北風速変化. ! Northward wind tendency real(DP), intent(out), optional :: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax) ! $ \DP{T}{t} $ . 温度変化. ! Temperature tendency real(DP), intent(out), optional :: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) ! $ \DP{q}{t} $ . 質量混合比変化. ! Mass mixing ratio tendency ! 作業変数 ! Work variables ! integer:: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer:: n ! 組成方向に回る DO ループ用作業変数 ! Work variables for DO loop in dimension of constituents ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. vdiffusion_my_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! Check arguments ! if ( present( xyz_DUDt ) ) then if ( .not. present( xyr_MomFluxX ) ) then call MessageNotify( 'E', module_name, 'xyr_MomFluxX has to be present.' ) end if end if if ( present( xyz_DVDt ) ) then if ( .not. present( xyr_MomFluxY ) ) then call MessageNotify( 'E', module_name, 'xyr_MomFluxY has to be present.' ) end if end if if ( present( xyz_DTempDt ) ) then if ( .not. present( xyr_HeatFlux ) ) then call MessageNotify( 'E', module_name, 'xyr_HeatFlux has to be present.' ) end if end if if ( present( xyzf_DQMixDt ) ) then if ( .not. present( xyrf_QMixFlux ) ) then call MessageNotify( 'E', module_name, 'xyrf_QMixFlux has to be present.' ) end if end if if ( present( xyz_DUDt ) ) then do k = 1, kmax xyz_DUDt(:,:,k) = & & + Grav * ( xyr_MomFluxX(:,:,k) - xyr_MomFluxX(:,:,k-1) ) & & / ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) end do end if if ( present( xyz_DVDt ) ) then do k = 1, kmax xyz_DVDt(:,:,k) = & & + Grav * ( xyr_MomFluxY(:,:,k) - xyr_MomFluxY(:,:,k-1) ) & & / ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) end do end if if ( present( xyz_DTempDt ) ) then do k = 1, kmax xyz_DTempDt(:,:,k) = & & + Grav / CpDry & & * ( xyr_HeatFlux(:,:,k) - xyr_HeatFlux(:,:,k-1) ) & & / ( xyr_Press (:,:,k) - xyr_Press (:,:,k-1) ) end do end if if ( present( xyzf_DQMixDt ) ) then do n = 1, ncmax do k = 1, kmax xyzf_DQMixDt(:,:,k,n) = & & + Grav * ( xyrf_QMixFlux(:,:,k,n) - xyrf_QMixFlux(:,:,k-1,n) ) & & / ( xyr_Press (:,:,k) - xyr_Press (:,:,k-1) ) end do end do end if end subroutine VDiffusionExpTendency !-------------------------------------------------------------------------------------- subroutine VDiffusionOutPut( & & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (in) & xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt, & ! (in) & xyr_Press, xyz_Exner, xyr_Exner, & ! (in) & xyr_VirTemp, xyz_Height, & ! (in) & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef & ! (in) & ) ! ! フラックス (xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux). ! について, その他の引数を用いて補正し, 出力を行う. ! ! Fluxes (xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux) are ! corrected by using other arguments, and the corrected values are output. ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: & & CpDry, & ! $ C_p $ [J kg-1 K-1]. ! 乾燥大気の定圧比熱. ! Specific heat of air at constant pressure & LatentHeat, & ! $ L $ [J kg-1] . ! 凝結の潜熱. ! Latent heat of condensation & GasRDry ! $ R $ [J kg-1 K-1]. ! 乾燥大気の気体定数. ! Gas constant of air ! 時刻管理 ! Time control ! use timeset, only: & & DelTime, & ! $ \Delta t $ [s] & TimeN, & ! ステップ $ t $ の時刻. Time of step $ t $. & TimesetClockStart, TimesetClockStop ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! 宣言文 ; Declaration statements ! implicit none real(DP), intent(in):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax) ! 東西方向運動量フラックス. ! Eastward wind flux real(DP), intent(in):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax) ! 南北方向運動量フラックス. ! Northward momentum flux real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax) ! 熱フラックス. ! Heat flux real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) ! 質量フラックス. ! Mass flux of constituents real(DP), intent(in):: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax) ! $ \DP{u}{t} $ . 東西風速変化. ! Eastward wind tendency real(DP), intent(in):: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax) ! $ \DP{v}{t} $ . 南北風速変化. ! Northward wind tendency real(DP), intent(in):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax) ! $ \DP{T}{t} $ . 温度変化. ! Temperature tendency real(DP), intent(in):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) ! $ \DP{q}{t} $ . 混合比変化. ! Mass mixing ratio tendency real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{p} $ . 気圧 (半整数レベル). ! Air pressure (half level) real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax) ! Exner 関数 (整数レベル). ! Exner function (full level) real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax) ! Exner 関数 (半整数レベル). ! Exner function (half level) real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax) ! $ \hat{T}_v $ . 仮温度 (半整数レベル). ! Virtual temperature (half level) real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax) ! 高度 (整数レベル). ! Height (full level) real(DP), intent(in):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：運動量. ! Diffusion coefficient: velocity real(DP), intent(in):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：温度. ! Transfer coefficient: temperature real(DP), intent(in):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax) ! 拡散係数：比湿. ! Diffusion coefficient: specific humidity ! 出力のための作業変数 ! Work variables for output ! real(DP):: xyr_MomFluxXCor (0:imax-1, 1:jmax, 0:kmax) ! 東西方向運動量フラックス. ! Eastward momentum flux real(DP):: xyr_MomFluxYCor (0:imax-1, 1:jmax, 0:kmax) ! 南北方向運動量フラックス. ! Northward momentum flux real(DP):: xyr_HeatFluxCor (0:imax-1, 1:jmax, 0:kmax) ! 熱フラックス. ! Heat flux real(DP):: xyrf_QMixFluxCor(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) ! 質量フラックス. ! Mass flux of constituents ! 作業変数 ! Work variables ! real(DP) :: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax) ! 輸送係数：運動量. ! Transfer coefficient: velocity real(DP) :: xyr_TempTransCoef (0:imax-1, 1:jmax, 0:kmax) ! 輸送係数：温度. ! Transfer coefficient: temperature real(DP) :: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) ! 輸送係数：質量. ! Transfer coefficient: mass of constituents integer:: i ! 経度方向に回る DO ループ用作業変数 ! Work variables for DO loop in longitude integer:: j ! 緯度方向に回る DO ループ用作業変数 ! Work variables for DO loop in latitude integer:: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer:: n ! 組成方向に回る DO ループ用作業変数 ! Work variables for DO loop in dimension of constituents real(DP):: LCp ! $ L / C_p $ [K]. ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. vdiffusion_my_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 計算時間計測開始 ! Start measurement of computation time ! call TimesetClockStart( module_name ) ! 輸送係数の計算 ! Calculate transfer coefficient ! xyr_VelTransCoef (:,:,0) = 0.0_DP xyr_VelTransCoef (:,:,kmax) = 0.0_DP xyr_TempTransCoef(:,:,0) = 0.0_DP xyr_TempTransCoef(:,:,kmax) = 0.0_DP xyr_QMixTransCoef(:,:,0) = 0.0_DP xyr_QMixTransCoef(:,:,kmax) = 0.0_DP do k = 1, kmax-1 xyr_VelTransCoef(:,:,k) = & & xyr_VelDiffCoef(:,:,k) & & * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) xyr_TempTransCoef(:,:,k) = & & xyr_TempDiffCoef(:,:,k) & & * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) xyr_QMixTransCoef(:,:,k) = & & xyr_QMixDiffCoef(:,:,k) & & * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) & & / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) ) end do ! 風速, 温度, 比湿フラックス補正 ! Correct fluxes of wind, temperature, specific humidity ! LCp = LatentHeat / CpDry do k = 1, kmax-1 do j = 1, jmax do i = 0, imax-1 xyr_MomFluxXCor( i,j,k ) = xyr_MomFluxX( i,j,k ) & & + ( xyz_DUDt( i,j,k ) - xyz_DUDt( i,j,k+1 ) ) & & * xyr_VelTransCoef( i,j,k ) * DelTime xyr_MomFluxYCor( i,j,k ) = xyr_MomFluxY( i,j,k ) & & + ( xyz_DVDt( i,j,k ) - xyz_DVDt( i,j,k+1 ) ) & & * xyr_VelTransCoef( i,j,k ) * DelTime xyr_HeatFluxCor( i,j,k ) = xyr_HeatFlux( i,j,k ) & & + ( xyz_DTempDt( i,j,k ) / xyz_Exner( i,j,k ) & & - xyz_DTempDt( i,j,k+1 ) / xyz_Exner( i,j,k+1 ) & & ) * CpDry * xyr_TempTransCoef( i,j,k ) & & * xyr_Exner( i,j,k ) * DelTime end do end do end do do n = 1, ncmax do k = 1, kmax-1 do j = 1, jmax do i = 0, imax-1 xyrf_QMixFluxCor( i,j,k,n ) = xyrf_QMixFlux( i,j,k,n ) & & + ( xyzf_DQMixDt( i,j,k,n ) - xyzf_DQMixDt( i,j,k+1,n ) ) & & * CpDry * xyr_QMixTransCoef( i,j,k ) * LCp * DelTime end do end do end do end do xyr_MomFluxXCor (:,:,0) = 0.0_DP xyr_MomFluxXCor (:,:,kmax) = 0.0_DP xyr_MomFluxYCor (:,:,0) = 0.0_DP xyr_MomFluxYCor (:,:,kmax) = 0.0_DP xyr_HeatFluxCor (:,:,0) = 0.0_DP xyr_HeatFluxCor (:,:,kmax) = 0.0_DP do n = 1, ncmax xyrf_QMixFluxCor(:,:,0,n) = 0.0_DP xyrf_QMixFluxCor(:,:,kmax,n) = 0.0_DP end do ! MEMO ! Output values of surface fluxes in MomFluxX, MomFluxY, HeatFlux, and QVapFlux ! are not correct. (YOT, 2009/08/14) ! Please refer to output variables, 'TauX', 'TauY', 'Sens', and 'Evap' for those ! values. (YOT, 2011/05/28) ! ヒストリデータ出力 ! History data output ! call HistoryAutoPut( TimeN, 'MomFluxX', xyr_MomFluxXCor ) call HistoryAutoPut( TimeN, 'MomFluxY', xyr_MomFluxYCor ) call HistoryAutoPut( TimeN, 'HeatFlux', xyr_HeatFluxCor ) call HistoryAutoPut( TimeN, 'QVapFlux', xyrf_QMixFluxCor ) call HistoryAutoPut( TimeN, 'DUDtVDiff' , xyz_DUDt ) call HistoryAutoPut( TimeN, 'DVDtVDiff' , xyz_DVDt ) call HistoryAutoPut( TimeN, 'DTempDtVDiff', xyz_DTempDt ) call HistoryAutoPut( TimeN, 'DQVapDtVDiff', xyzf_DQMixDt(:,:,:,IndexH2OVap) ) ! 計算時間計測一時停止 ! Pause measurement of computation time ! call TimesetClockStop( module_name ) end subroutine VDiffusionOutPut !-------------------------------------------------------------------------------------- subroutine VDiffusionInit ! ! vdiffusion_my モジュールの初期化を行います. ! NAMELIST#vdiffusion_my_nml の読み込みはこの手続きで行われます. ! ! "vdiffusion_my" module is initialized. ! "NAMELIST#vdiffusion_my_nml" is loaded in this procedure. ! ! モジュール引用 ; USE statements ! ! NAMELIST ファイル入力に関するユーティリティ ! Utilities for NAMELIST file input ! use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid ! ファイル入出力補助 ! File I/O support ! use dc_iounit, only: FileOpen ! 種別型パラメタ ! Kind type parameter ! use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output ! 文字列操作 ! Character handling ! use dc_string, only: StoA ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoAddVariable ! 座標データ設定 ! Axes data settings ! use axesset, only: & & AxnameX, & & AxnameY, & & AxnameZ, & & AxnameR, & & AxnameT ! メッセージ制御 ! Message control ! use mpi_messagecntl, only : DoesOutputMPIMessage ! 宣言文 ; Declaration statements ! implicit none integer:: unit_nml ! NAMELIST ファイルオープン用装置番号. ! Unit number for NAMELIST file open integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT. ! IOSTAT of NAMELIST read ! NAMELIST 変数群 ! NAMELIST group name ! namelist /vdiffusion_my_nml/ & & FlagConstDiffCoef, & & ConstDiffCoefM, ConstDiffCoefH, & ! & SquareVelMin, BulkRiNumMin, & ! & MixLengthMax, ShMin, SmMin, & & VelDiffCoefMin, TempDiffCoefMin, & & VelDiffCoefMax, TempDiffCoefMax, & ! & MYConstA1, MYConstB1, MYConstA2, MYConstB2, MYConstC1 ! ! デフォルト値については初期化手続 "vdiffusion_my#VDiffInit" ! のソースコードを参照のこと. ! ! Refer to source codes in the initialization procedure ! "vdiffusion_my#VDiffInit" for the default values. ! ! 実行文 ; Executable statement ! if ( vdiffusion_my_inited ) return ! デフォルト値の設定 ! Default values settings ! FlagConstDiffCoef = .false. ConstDiffCoefM = 0.0_DP ConstDiffCoefH = 0.0_DP SquareVelMin = 0.1_DP BulkRiNumMin = - 100.0_DP MixLengthMax = 300.0_DP ShMin = 0.0_DP SmMin = 0.0_DP VelDiffCoefMin = 0.1_DP TempDiffCoefMin = 0.1_DP VelDiffCoefMax = 10000.0_DP TempDiffCoefMax = 10000.0_DP ! Parameters proposed by Mellor and Yamada (1982). ! MYConstA1 = 0.92_DP MYConstB1 = 16.6_DP MYConstA2 = 0.74_DP MYConstB2 = 10.1_DP MYConstC1 = 0.08_DP ! NAMELIST の読み込み ! NAMELIST is input ! if ( trim(namelist_filename) /= '' ) then call FileOpen( unit_nml, & ! (out) & namelist_filename, mode = 'r' ) ! (in) rewind( unit_nml ) read( unit_nml, & ! (in) & nml = vdiffusion_my_nml, & ! (out) & iostat = iostat_nml ) ! (out) close( unit_nml ) call NmlutilMsg( iostat_nml, module_name ) ! (in) if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = vdiffusion_my_nml ) end if ! ヒストリデータ出力のためのへの変数登録 ! Register of variables for history data output ! call HistoryAutoAddVariable( 'VelDiffCoef', & & (/ AxNameX, AxNameY, AxNameR, AxNameT /), & & 'diffusion coef. momentum', 'm2 s-1' ) call HistoryAutoAddVariable( 'TempDiffCoef', & & (/ AxNameX, AxNameY, AxNameR, AxNameT /), & & 'diffusion coef. heat ', 'm2 s-1' ) call HistoryAutoAddVariable( 'QVapDiffCoef', & & (/ AxNameX, AxNameY, AxNameR, AxNameT /), & & 'diffusion coef. moisture', 'm2 s-1' ) call HistoryAutoAddVariable( 'MomFluxX', & & (/ AxNameX, AxNameY, AxNameR, AxNameT /), & & 'eastward momentum flux', 'N m-2' ) call HistoryAutoAddVariable( 'MomFluxY', & & (/ AxNameX, AxNameY, AxNameR, AxNameT /), & & 'northward momentum flux', 'N m-2' ) call HistoryAutoAddVariable( 'HeatFlux', & & (/ AxNameX, AxNameY, AxNameR, AxNameT /), & & 'heat flux', 'W m-2' ) call HistoryAutoAddVariable( 'QVapFlux', & & (/ AxNameX, AxNameY, AxNameR, AxNameT /), & & 'moisture flux', 'W m-2' ) call HistoryAutoAddVariable( 'DUDtVDiff', & & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), & & 'tendency of zonal wind by vertical diffusion', 'm s-2' ) call HistoryAutoAddVariable( 'DVDtVDiff', & & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), & & 'tendency of meridional wind by vertical diffusion', 'm s-2' ) call HistoryAutoAddVariable( 'DTempDtVDiff', & & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), & & 'tendency of temperature by vertical diffusion', 'K s-1' ) call HistoryAutoAddVariable( 'DQVapDtVDiff', & & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), & & 'tendency of specific humidity by vertical diffusion', 's-1' ) call HistoryAutoAddVariable( 'TurKinEne', & & (/ AxNameX, AxNameY, AxNameR, AxNameT /), & & 'turbulent kinetic energy', 'm2 s-2' ) call HistoryAutoAddVariable( 'TKEPShear', & & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), & & 'turbulent kinetic energy production rate by shear', 'm2 s-3' ) call HistoryAutoAddVariable( 'TKEPBuoy', & & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), & & 'turbulent kinetic energy production rate by buoyancy', 'm2 s-3' ) call HistoryAutoAddVariable( 'TKEDiss', & & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), & & 'turbulent kinetic energy dissipation rate', 'm2 s-3' ) call HistoryAutoAddVariable( 'MixLength', & & (/ AxNameX, AxNameY, AxNameZ, AxNameT /), & & 'mixing length', 'm' ) ! 印字 ; Print ! call MessageNotify( 'M', module_name, '----- Initialization Messages -----' ) call MessageNotify( 'M', module_name, 'For vertical diffusion flux:' ) call MessageNotify( 'M', module_name, ' FlagConstDiffCoef = %b', l = (/ FlagConstDiffCoef /) ) call MessageNotify( 'M', module_name, ' ConstDiffCoefM = %f', d = (/ ConstDiffCoefM /) ) call MessageNotify( 'M', module_name, ' ConstDiffCoefH = %f', d = (/ ConstDiffCoefH /) ) call MessageNotify( 'M', module_name, ' SquareVelMin = %f', d = (/ SquareVelMin /) ) call MessageNotify( 'M', module_name, ' BulkRiNumMin = %f', d = (/ BulkRiNumMin /) ) call MessageNotify( 'M', module_name, 'For diffusion coefficients:' ) call MessageNotify( 'M', module_name, ' MixLengthMax = %f', d = (/ MixLengthMax /) ) call MessageNotify( 'M', module_name, ' ShMin = %f', d = (/ ShMin /) ) call MessageNotify( 'M', module_name, ' SmMin = %f', d = (/ SmMin /) ) call MessageNotify( 'M', module_name, ' VelDiffCoefMin = %f', d = (/ VelDiffCoefMin /) ) call MessageNotify( 'M', module_name, ' TempDiffCoefMin = %f', d = (/ TempDiffCoefMin /) ) call MessageNotify( 'M', module_name, ' VelDiffCoefMax = %f', d = (/ VelDiffCoefMax /) ) call MessageNotify( 'M', module_name, ' TempDiffCoefMax = %f', d = (/ TempDiffCoefMax /) ) call MessageNotify( 'M', module_name, ' MYConstA1 = %f', d = (/ MYConstA1 /) ) call MessageNotify( 'M', module_name, ' MYConstB1 = %f', d = (/ MYConstB1 /) ) call MessageNotify( 'M', module_name, ' MYConstA2 = %f', d = (/ MYConstA2 /) ) call MessageNotify( 'M', module_name, ' MYConstB2 = %f', d = (/ MYConstB2 /) ) call MessageNotify( 'M', module_name, ' MYConstC1 = %f', d = (/ MYConstC1 /) ) call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) ) vdiffusion_my_inited = .true. end subroutine VDiffusionInit !-------------------------------------------------------------------------------------- end module vdiffusion_my