subroutine GWDM1987( xyz_U, xyz_V, xyz_Temp, xyz_Press, xyr_Press, xyz_Exner, xyz_Height, xy_SurfHeightStd, xyz_DUDt, xyz_DVDt )
!
!
!
! Tendency of gravity wave drag based on McFarlane (1987)
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav, GasRDry
! $ R $ [J kg-1 K-1].
! 乾燥大気の気体定数.
! Gas constant of air
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 宣言文 ; Declaration statements
!
real(DP), intent(in ) :: xyz_U (0:imax-1, 1:jmax, 1:kmax)
! Zonal wind
real(DP), intent(in ) :: xyz_V (0:imax-1, 1:jmax, 1:kmax)
! Meridional wind
real(DP), intent(in ) :: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! Temperature
real(DP), intent(in ) :: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
! Pressure
real(DP), intent(in ) :: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! Pressure
real(DP), intent(in ) :: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
! Exner function
real(DP), intent(in ) :: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
! Height
real(DP), intent(in ) :: xy_SurfHeightStd(0:imax-1, 1:jmax)
real(DP), intent(out) :: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
! 東西風変化率.
! Zonal wind tendency
real(DP), intent(out) :: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
! 南北風変化率.
! Meridional wind tendency
! 作業変数
! Work variables
!
real(DP) :: xy_OrogEffHeight(0:imax-1, 1:jmax)
real(DP) :: xyz_DelPress (0:imax-1, 1:jmax, 1:kmax)
integer :: xy_KIndexRef (0:imax-1, 1:jmax)
real(DP) :: xy_URef (0:imax-1, 1:jmax)
real(DP) :: xy_VRef (0:imax-1, 1:jmax)
real(DP) :: xyz_WindSpeed (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_AbsWindSpeedRef(0:imax-1, 1:jmax)
real(DP) :: xy_RhoRef (0:imax-1, 1:jmax)
real(DP) :: xy_NRef (0:imax-1, 1:jmax)
real(DP) :: xyz_Rho (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_PotTemp (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_N (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_Amp (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_ZeroOne (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_MomFluxRef (0:imax-1, 1:jmax)
real(DP) :: xyz_MomFlux (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: MomFluxSaturated
real(DP) :: xyz_DMomFluxDU(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DWindSpeedDt(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyr_MomFluxA (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyr_MomFluxXA(0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyr_MomFluxYA(0:imax-1, 1:jmax, 0:kmax)
real(DP) :: WindSpeedTentative
!!$ real(DP) :: MomFluxTentative
integer :: mmax
integer :: ms
integer :: me
real(DP) :: az_A(0:imax*jmax-1, 1:kmax)
real(DP) :: az_B(0:imax*jmax-1, 1:kmax)
real(DP) :: az_C(0:imax*jmax-1, 1:kmax)
real(DP) :: az_D(0:imax*jmax-1, 1:kmax)
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 :: kp
integer :: kn
!!$ real(DP) :: xyz_WindSpeedTentative(0:imax-1, 1:jmax, 1:kmax)
!!$ integer :: itr
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. gwd_m1987_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 additional variables
!
xy_OrogEffHeight = 2.0_DP * xy_SurfHeightStd
!
do k = 1, kmax
xyz_DelPress(:,:,k) = xyr_Press(:,:,k-1) - xyr_Press(:,:,k)
end do
!
! Determine reference level
!
if ( FlagDetermineRefLevByStd ) then
xy_KIndexRef = 2
do k = 1+1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_Height(i,j,k) < xy_OrogEffHeight(i,j) ) then
xy_KIndexRef(i,j) = k
end if
end do
end do
end do
else
xy_KIndexRef = 2
do k = 1+1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_Press(i,j,k) / xyr_Press(i,j,0) > SigmaRef ) then
xy_KIndexRef(i,j) = k
end if
end do
end do
end do
end if
!
! Set reference level wind velocity
!
do j = 1, jmax
do i = 0, imax-1
xy_URef = xyz_U(i,j,xy_KIndexRef(i,j))
xy_VRef = xyz_V(i,j,xy_KIndexRef(i,j))
end do
end do
xy_AbsWindSpeedRef = sqrt( xy_URef**2 + xy_VRef**2 )
! Calculation of additional variables
!
xyz_Rho = xyz_Press / ( GasRDry * xyz_Temp )
do j = 1, jmax
do i = 0, imax-1
xy_RhoRef = xyz_Rho(i,j,xy_KIndexRef(i,j))
end do
end do
!
xyz_PotTemp = xyz_Temp / xyz_Exner
!
do k = 1, kmax
kp = max( k - 1, 1 )
kn = min( k + 1, kmax )
xyz_N(:,:,k) = Grav / xyz_PotTemp(:,:,k) * ( xyz_PotTemp(:,:,kn) - xyz_PotTemp(:,:,kp) ) / ( xyz_Height (:,:,kn) - xyz_Height (:,:,kp) )
end do
!!$ xyz_N = max( xyz_N, 1.0e-6_DP )
xyz_N = max( xyz_N, 0.0_DP )
xyz_N = sqrt( xyz_N )
do j = 1, jmax
do i = 0, imax-1
xy_NRef = xyz_N(i,j,xy_KIndexRef(i,j))
end do
end do
!
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xy_AbsWindSpeedRef(i,j) /= 0.0_DP ) then
xyz_WindSpeed(i,j,k) = ( xyz_U(i,j,k) * xy_URef(i,j) + xyz_V(i,j,k) * xy_VRef(i,j) ) / xy_AbsWindSpeedRef(i,j)
else
xyz_WindSpeed(i,j,k) = 0.0_DP
end if
end do
end do
end do
! Negative wind speed is inrelevant in the current formulation.
xyz_WindSpeed = max( xyz_WindSpeed, 0.0_DP )
! Wave amplitude
!
! Momentum flux parallel to the reference level flow at reference level
!
xy_MomFluxRef = - MomFluxFactor * xy_OrogEffHeight**2 * xy_RhoRef * xy_NRef * xy_AbsWindSpeedRef
!
! Momentum flux parallel to the reference level flow
!
do j = 1, jmax
do i = 0, imax-1
! Region at and below the reference level
do k = 1, xy_KIndexRef(i,j)
xyz_MomFlux(i,j,k) = xy_MomFluxRef(i,j)
xyz_ZeroOne(i,j,k) = 0.0_DP
end do
! Region above the reference level and below the highest model level
do k = xy_KIndexRef(i,j)+1, kmax-1
! momentum flux is same as that in lower level tentatively
xyz_MomFlux(i,j,k) = xyz_MomFlux(i,j,k-1)
! calculate momentum flux in the case of saturation
if ( xyz_N(i,j,k) > 0.0_DP ) then
MomFluxSaturated = - MomFluxFactor * CrtlFNumSq * xyz_Rho(i,j,k) / xyz_N(i,j,k) * xyz_WindSpeed(i,j,k)**3
else
MomFluxSaturated = 0.0_DP
end if
! comparison of momentum flux
! check whether saturationed or not
! It should be noted that momentum flux here is negative, since
! a direction of the momentum flux is parallel to reference level
! flow.
if ( MomFluxSaturated > xyz_MomFlux(i,j,k) ) then
! saturation region
xyz_MomFlux(i,j,k) = MomFluxSaturated
xyz_ZeroOne(i,j,k) = 1.0_DP
else
! non-saturation region
xyz_ZeroOne(i,j,k) = 0.0_DP
end if
end do
! highest model level
do k = kmax, kmax
! momentum flux is same as that in lower level
xyz_MomFlux(i,j,k) = xyz_MomFlux(i,j,k-1)
xyz_ZeroOne(i,j,k) = 0.0_DP
end do
end do
end do
!
! Momentum flux derivative with reference level flow
!
!!$ xyz_DMomFluxDU = &
!!$ & - 3.0_DP * MomFluxFactor * CrtlFNumSq * xyz_Rho / xyz_N &
!!$ & * xyz_WindSpeed**2 &
!!$ & * xyz_ZeroOne
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_N(i,j,k) > 0.0_DP ) then
xyz_DMomFluxDU(i,j,k) = - 3.0_DP * MomFluxFactor * CrtlFNumSq * xyz_Rho(i,j,k) / xyz_N(i,j,k) * xyz_WindSpeed(i,j,k)**2 * xyz_ZeroOne(i,j,k)
else
xyz_DMomFluxDU(i,j,k) = 0.0_DP
end if
end do
end do
end do
!
! calculation of wind velocity tendency
!
do j = 1, jmax
do i = 0, imax-1
! No deceleration is assumed in the highest level
k = kmax
xyz_DWindSpeedDt(i,j,k) = 0.0_DP
do k = kmax-1, 1+1, -1
xyz_DWindSpeedDt(i,j,k) = Grav / xyz_DelPress(i,j,k) / 2.0_DP * ( xyz_MomFlux(i,j,k-1) - xyz_MomFlux(i,j,k+1) - xyz_DMomFluxDU(i,j,k+1) * xyz_DWindSpeedDt(i,j,k+1) * ( 2.0_DP * DelTime ) ) / ( 1.0_DP - Grav / xyz_DelPress(i,j,k) / 2.0_DP * xyz_DMomFluxDU(i,j,k) * ( 2.0_DP * DelTime ) )
! Wind speed tendency at k level and momentum flux at k-1 level
! are estimated again.
if ( k >= xy_KIndexRef(i,j) ) then
! Region above reference level
WindSpeedTentative = xyz_WindSpeed(i,j,k) + xyz_DWindSpeedDt(i,j,k) * ( 2.0_DP * DelTime )
if ( WindSpeedTentative < 0.0_DP ) then
xyz_DWindSpeedDt(i,j,k) = ( 0.0_DP - xyz_WindSpeed(i,j,k) ) / ( 2.0_DP * DelTime )
xyz_MomFlux(i,j,k-1) = ( 2.0_DP * xyz_DelPress(i,j,k) / Grav - xyz_DMomFluxDU(i,j,k) * ( 2.0_DP * DelTime ) ) * xyz_DWindSpeedDt(i,j,k) + xyz_MomFlux(i,j,k+1) + xyz_DMomFluxDU(i,j,k+1) * xyz_DWindSpeedDt(i,j,k+1) * ( 2.0_DP * DelTime )
end if
else
! below the reference level
xyz_DWindSpeedDt(i,j,k) = 0.0_DP
xyz_MomFlux(i,j,k-1) = ( 2.0_DP * xyz_DelPress(i,j,k) / Grav - xyz_DMomFluxDU(i,j,k) * ( 2.0_DP * DelTime ) ) * xyz_DWindSpeedDt(i,j,k) + xyz_MomFlux(i,j,k+1) + xyz_DMomFluxDU(i,j,k+1) * xyz_DWindSpeedDt(i,j,k+1) * ( 2.0_DP * DelTime )
end if
end do
! No deceleration is assumed in the lowest level
k = 1
xyz_DWindSpeedDt(i,j,k) = 0.0_DP
end do
end do
!!$ ! No deceleration is assumed in the lowest level
!!$ k = 1
!!$ xyz_DWindSpeedDt(:,:,k) = 0.0_DP
!!$ do k = 1+1, kmax-1
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ !
!!$ ! calculation of wind velocity tendency
!!$ !
!!$ xyz_DWindSpeedDt(i,j,k) = &
!!$ & Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$ & * ( xyz_MomFlux(i,j,k-1) &
!!$ & + xyz_DMomFluxDU(i,j,k-1) * xyz_DWindSpeedDt(i,j,k-1) &
!!$ & * ( 2.0_DP * DelTime ) &
!!$ & - xyz_MomFlux(i,j,k+1) ) &
!!$ & / ( 1.0_DP + Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$ & * xyz_DMomFluxDU(i,j,k) * ( 2.0_DP * DelTime ) )
!!$ xyz_DWindSpeedDt(i,j,k) = 0.0_DP
!!$
!!$ do itr = 1, 10
!!$
!!$ xyz_WindSpeedTentative(i,j,k) = &
!!$ & xyz_WindSpeed(i,j,k) + xyz_DWindSpeedDt(i,j,k) * ( 2.0_DP * DelTime )
!!$
!!$ if ( xyz_N(i,j,k) > 0.0_DP ) then
!!$ xyz_DMomFluxDU(i,j,k) = &
!!$ & - 3.0_DP * MomFluxFactor * CrtlFNumSq &
!!$ & * xyz_Rho(i,j,k) / xyz_N(i,j,k) &
!!$ & * xyz_WindSpeedTentative(i,j,k)**2 &
!!$ & * xyz_ZeroOne(i,j,k)
!!$ else
!!$ xyz_DMomFluxDU(i,j,k) = 0.0_DP
!!$ end if
!!$
!!$ xyz_DWindSpeedDt(i,j,k) = &
!!$ & Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$ & * ( xyz_MomFlux(i,j,k-1) &
!!$ & + xyz_DMomFluxDU(i,j,k-1) * xyz_DWindSpeedDt(i,j,k-1) &
!!$ & * ( 2.0_DP * DelTime ) &
!!$ & - xyz_MomFlux(i,j,k+1) ) &
!!$ & / ( 1.0_DP + Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$ & * xyz_DMomFluxDU(i,j,k) * ( 2.0_DP * DelTime ) )
!!$
!!$ end do
!!$ !
!!$ ! preparation for next level
!!$ !
!!$ ! estimated momentum flux at k and next time step
!!$ MomFluxTentative = &
!!$ & xyz_MomFlux(i,j,k) &
!!$ & + xyz_DMomFluxDU(i,j,k) * xyz_DWindSpeedDt(i,j,k) &
!!$ & * ( 2.0_DP * DelTime )
!!$ xyz_MomFlux(i,j,k+1) = MomFluxTentative
!!$ ! calculate momentum flux at k+1 in the case of saturation
!!$ MomFluxSaturated = &
!!$ & - MomFluxFactor * CrtlFNumSq &
!!$ & * xyz_Rho(i,j,k+1) / xyz_N(i,j,k+1) * xyz_WindSpeed(i,j,k+1)**3
!!$ ! comparison of momentum flux
!!$ ! check whether saturationed or not
!!$ if ( abs( MomFluxSaturated ) < abs( xyz_MomFlux(i,j,k+1) ) ) then
!!$ ! saturation region
!!$ ! set saturated momentum flux
!!$ xyz_MomFlux(i,j,k+1) = MomFluxSaturated
!!$ ! derivative of momentum flux with reference level flow
!!$ xyz_DMomFluxDU(i,j,k+1) = &
!!$ & - 3.0_DP * MomFluxFactor * CrtlFNumSq &
!!$ & * xyz_Rho(i,j,k+1) / xyz_N(i,j,k+1) &
!!$ & * xyz_WindSpeed(i,j,k+1)**2
!!$ else
!!$ ! non-saturation region
!!$ ! derivative of momentum flux with reference level flow
!!$ xyz_DMomFluxDU(i,j,k+1) = 0.0_DP
!!$ end if
!!$ end do
!!$ end do
!!$ end do
!!$ ! No deceleration is assumed in the highest level
!!$ k = kmax
!!$ xyz_DWindSpeedDt(:,:,k) = 0.0_DP
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xy_AbsWindSpeedRef(i,j) /= 0.0_DP ) then
xyz_DUDt(i,j,k) = xyz_DWindSpeedDt(i,j,k) * xy_URef(i,j) / xy_AbsWindSpeedRef(i,j)
xyz_DVDt(i,j,k) = xyz_DWindSpeedDt(i,j,k) * xy_VRef(i,j) / xy_AbsWindSpeedRef(i,j)
else
xyz_DUDt(i,j,k) = 0.0_DP
xyz_DVDt(i,j,k) = 0.0_DP
end if
end do
end do
end do
!!$ if ( FlagGWDDamp ) then
!!$
!!$ GWDDampCoef = 1.0_DP / DelTime * ( GWDDampPeriod - TimeN ) / GWDDampPeriod
!!$
!!$ if ( GWDDampCoef < 0.0_DP ) GWDDampCoef = 0.0_DP
!!$
!!$ xyz_DUDt = xyz_DUDt / ( 1.0_DP + 2.0_DP * DelTime * DivDampCoef )
!!$
!!$ end if
! estimated momentum flux at layer interface and at next time step
do k = 0+1, kmax-1
xyr_MomFluxA(:,:,k) = ( xyz_MomFlux(:,:,k) + xyz_MomFlux(:,:,k+1) + xyz_DMomFluxDU(:,:,k+1) * xyz_DWindSpeedDt(:,:,k+1) * ( 2.0_DP * DelTime ) ) / 2.0_DP
!!$ xyr_MomFluxA(:,:,k) = &
!!$ & ( xyz_MomFlux(:,:,k) &
!!$ & + xyz_DMomFluxDU(:,:,k) * xyz_DWindSpeedDt(:,:,k) &
!!$ & * ( 2.0_DP * DelTime ) &
!!$ & + xyz_MomFlux(:,:,k+1) ) &
!!$ & / 2.0_DP
end do
k = 0
xyr_MomFluxA(:,:,k) = xyr_MomFluxA(:,:,k+1)
k = kmax
xyr_MomFluxA(:,:,k) = xyr_MomFluxA(:,:,k-1)
do k = 0, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xy_AbsWindSpeedRef(i,j) /= 0.0_DP ) then
xyr_MomFluxXA(i,j,k) = xyr_MomFluxA(i,j,k) * xy_URef(i,j) / xy_AbsWindSpeedRef(i,j)
xyr_MomFluxYA(i,j,k) = xyr_MomFluxA(i,j,k) * xy_VRef(i,j) / xy_AbsWindSpeedRef(i,j)
else
xyr_MomFluxXA(i,j,k) = 0.0_DP
xyr_MomFluxYA(i,j,k) = 0.0_DP
end if
end do
end do
end do
!!$ ! Set up of simultaneously linear equations
!!$ do k = 1, kmax
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ az_A(i+imax*(j-1),k) = &
!!$ & Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$ & * xyz_DMomFluxDU(i,j,k-1) &
!!$ & * ( 2.0_DP * DelTime )
!!$ az_A(i+imax*(j-1),k) = - 1.0_DP
!!$ az_C(i+imax*(j-1),k) = &
!!$ & - Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$ & * xyz_DMomFluxDU(i,j,k+1) &
!!$ & * ( 2.0_DP * DelTime )
!!$ az_D(i+imax*(j-1),k) = &
!!$ & - Grav / xyz_DelPress(i,j,k) / 2.0_DP &
!!$ & * ( xyz_MomFlux(i,j,k-1) - xyz_MomFlux(i,j,k+1) )
!!$ end do
!!$ end do
!!$ end do
!!$
!!$ mmax = imax*jmax
!!$ ms = 1
!!$ me = imax*jmax
!!$ call tridiag( mmax, kmax, az_A, az_B, az_C, az_D, ms, me )
!!$
!!$ do k = 1, kmax
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ xyz_DUDt = az_D(i+imax*(j-1),k) * xy_URef(i,j) / xy_WindSpeedRef(i,j)
!!$ xyz_DVDt = az_D(i+imax*(j-1),k) * xy_VRef(i,j) / xy_WindSpeedRef(i,j)
!!$ end do
!!$ end do
!!$ end do
! ヒストリデータ出力
! History data output
!
call HistoryAutoPut( TimeN, 'GWMomFlux' , xyr_MomFluxA )
call HistoryAutoPut( TimeN, 'GWMomFluxX', xyr_MomFluxXA )
call HistoryAutoPut( TimeN, 'GWMomFluxY', xyr_MomFluxYA )
call HistoryAutoPut( TimeN, 'DUDtGWD' , xyz_DUDt )
call HistoryAutoPut( TimeN, 'DVDtGWD' , xyz_DVDt )
call HistoryAutoPut( TimeN, 'DWSDtGWD' , xyz_DWindSpeedDt )
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine GWDM1987