成分の質量を補正します. xyzf_QMixRef が与えられた場合には,
全球積分値が xyzf_QMixRef のそれと 同じになるように補正します.
xyzf_QMixRef が与えられない場合には, 全球積分値が補正前のそれと
同じになるように補正します. xyzf_DQMixDt には xyz_QMix
の変化量が返ります.
subroutine MassFixer( xyr_Press, xyzf_QMix, xyr_PressRef, xyzf_QMixRef, xyzf_DQMixDt )
!
! 成分の質量を補正します.
! *xyzf_QMixRef* が与えられた場合には, 全球積分値が *xyzf_QMixRef* のそれと
! 同じになるように補正します.
! *xyzf_QMixRef* が与えられない場合には, 全球積分値が補正前のそれと
! 同じになるように補正します.
! *xyzf_DQMixDt* には *xyz_QMix* の変化量が返ります.
!
! Fix masses of constituents
! If *xyzf_QMixRef* is given, the mass is fixed to match its global integrated
! value is the same as that of *xyzf_QMixRef*.
! If *xyzf_QMixRef* is not given, the mass is fixed to match its global integrated
! value is the same as that of pre-fixed value.
! Variation of *xyzf_QMix* is returned to *xyz_DQMixDt*.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav ! $ g $ [m s-2].
! 重力加速度.
! Gravitational acceleration
! 積分と平均の操作
! Operation for integral and average
!
use intavr_operate, only: IntLonLat_xy
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimesetClockStart, TimesetClockStop
! 宣言文 ; 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(inout) :: xyzf_QMix (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q $ . 比湿. Specific humidity
real(DP), intent(in ), optional:: xyr_PressRef(0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
real(DP), intent(in ), optional:: xyzf_QMixRef(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ q \Delta p / g $ . 積分値を合わせる層内の成分の質量.
! Reference specific mass of constituent in a layer
real(DP), intent(out ), optional:: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \DP{q}{t} $ . 比湿補正率.
! Specific humidity correction
! 作業変数
! Work variables
!
real(DP):: xyz_QMixBefCor (0:imax-1, 1:jmax, 1:kmax)
! 修正前の比湿.
! Specific humidity before correction.
real(DP):: xyz_DelMass (0:imax-1, 1:jmax, 1:kmax)
! $ \Delta p / g $
!
real(DP):: xyz_DelMassRef (0:imax-1, 1:jmax, 1:kmax)
! $ \Delta p / g $ of reference
!
real(DP):: xyz_DelConsMass (0:imax-1, 1:jmax, 1:kmax)
! 各層内の成分の質量.
! Mass of constituents in a layer.
real(DP):: xyz_DelConsMassRef(0:imax-1, 1:jmax, 1:kmax)
! 積分値を合わせる各層内の成分の質量.
! Reference mass of constituents.
real(DP):: xy_ConsMass (0:imax-1, 1:jmax)
! 成分のカラム質量.
! Mass of constituents in a layer.
real(DP):: xy_ConsMassRef (0:imax-1, 1:jmax)
! 積分値を合わせる成分のカラム質量.
! Reference mass of constituents in a layer.
real(DP):: ConsMass
! 全球の各成分の質量
! Total mass of constituents
real(DP):: ConsMassRef
! 積分値を合わせる全球の各成分の質量
! Reference total mass of constituents.
!
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitudinal direction
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitudinal direction
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
!
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! 初期化
! Initialization
!
if ( .not. mass_fixer_inited ) call MassFixerInit
! Check arguments
!
if ( present( xyr_PressRef ) .or. present( xyzf_QMixRef ) ) then
if ( .not. ( present( xyr_PressRef ) .and. present( xyzf_QMixRef ) ) ) then
call MessageNotify( 'E', module_name, 'If xyr_PressRef or xyzf_QMixRef is given, both have to be given.' )
end if
end if
! $ \Delta p / g $ の計算
! Calculate $ \Delta p / g $
!
do k = 1, kmax
xyz_DelMass(:,:,k) = ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
end do
if ( present( xyr_PressRef ) ) then
do k = 1, kmax
xyz_DelMassRef(:,:,k) = ( xyr_PressRef(:,:,k-1) - xyr_PressRef(:,:,k) ) / Grav
end do
end if
do n = 1, ncmax
! Calculate mass of constituents
!
xyz_DelConsMass = xyzf_QMix(:,:,:,n) * xyz_DelMass
if ( present( xyzf_QMixRef ) ) then
xyz_DelConsMassRef = xyzf_QMixRef(:,:,:,n) * xyz_DelMassRef
else
xyz_DelConsMassRef = xyz_DelConsMass
end if
if ( present( xyzf_DQMixDt ) ) then
xyz_QMixBefCor = xyzf_QMix(:,:,:,n)
end if
! 負の質量を直下の層の質量で埋め合わせ.
! Negative mass is removed by filling it with the mass in a layer just below.
!
do k = kmax, 2, -1
do j = 1, jmax
do i = 0, imax-1
if ( xyz_DelConsMass(i,j,k) < 0.0_DP ) then
xyz_DelConsMass(i,j,k-1) = xyz_DelConsMass(i,j,k-1) + xyz_DelConsMass(i,j,k)
xyz_DelConsMass(i,j,k ) = 0.0_DP
end if
end do
end do
end do
k = 1
do j = 1, jmax
do i = 0, imax-1
if ( xyz_DelConsMass(i,j,k) < 0.0_DP ) then
xyz_DelConsMass(i,j,k) = 0.0_DP
end if
end do
end do
! 全球での補正
! Correction in globe
! 質量保存のために全体の質量を減少させる.
! Total mass is decreased to conserve mass.
!
xy_ConsMass = 0.0d0
xy_ConsMassRef = 0.0d0
do k = kmax, 1, -1
xy_ConsMass = xy_ConsMass + xyz_DelConsMass (:,:,k)
xy_ConsMassRef = xy_ConsMassRef + xyz_DelConsMassRef(:,:,k)
end do
ConsMass = IntLonLat_xy( xy_ConsMass )
ConsMassRef = IntLonLat_xy( xy_ConsMassRef )
if ( ConsMassRef < 0.0_DP ) then
call MessageNotify( 'M', module_name, 'ConsMassRef is negative. ConsMassRef is reset to zero, n = %d, ConsMassRef = %f.', i = (/ n /), d = (/ ConsMassRef /) )
ConsMassRef = 0.0_DP
!!$ call MessageNotify( 'E', module_name, 'ConsMassRef is negative, n = %d.', i = (/ n /) )
end if
if ( ConsMass /= 0.0_DP ) then
xyz_DelConsMass = ConsMassRef / ConsMass * xyz_DelConsMass
else
xyz_DelConsMass = 0.0_DP
end if
xyzf_QMix(:,:,:,n) = xyz_DelConsMass / xyz_DelMass
! 比湿変化の算出
! Calculate specific humidity variance
!
if ( present( xyzf_DQMixDt ) ) then
xyzf_DQMixDt(:,:,:,n) = xyzf_DQMixDt(:,:,:,n) + ( xyzf_QMix(:,:,:,n) - xyz_QMixBefCor ) / ( 2.0_DP * DelTime )
end if
end do
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyzf_QMix(i,j,k,1) < 0.0_DP ) then
write( 6, * ) 'NEGATIVE: ', i, j, k, xyzf_QMix(i,j,k,1)
end if
end do
end do
end do
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine MassFixer