subroutine MoistConvAdjust( xyz_Temp, xyz_QVap, xyz_Press, xyr_Press, xyz_DQH2OLiqDt )
!
! 羚炊修絲丈�茯睡��鴻����若��������, 羝�墾���羚帥��茯睡�.
!
! Adjust temperature and specific humidity by moist convective adjustment
!
! �≪�吾�ャ�若����� ; USE statements
!
! ����絎��域┃絎�
! Physical constants settings
!
use constants, only: Grav, GasRDry, CpDry, LatentHeat
! $ L $ [J kg-1] .
! ��腟������.
! Latent heat of condensation
! ���紫���
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! ���鴻�������若�水�阪��
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 蕋遵��罸�羚帥�����
! Evaluate saturation specific humidity
!
use saturate, only: xyz_CalcQVapSat, xy_CalcDQVapSatDTemp
! 絎ｈ��� ; Declaration statements
!
implicit none
real(DP), intent(inout):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
! $ T $ . 羝�墾. Temperature
real(DP), intent(inout):: xyz_QVap (0:imax-1, 1:jmax, 1:kmax)
! $ q $ . 罸�羚�. Specific humidity
!!$ real(DP), intent(inout):: xy_Rain (0:imax-1, 1:jmax)
!!$ ! ��羂顔��.
!!$ ! Precipitation
real(DP), intent(in):: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
! $ p $ . 羂��� (�贋�違������).
! Air pressure (full level)
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 羂��� (���贋�違������).
! Air pressure (half level)
real(DP), intent(out):: xyz_DQH2OLiqDt(0:imax-1, 1:jmax, 1:kmax)
! 篏�罐�����
! Work variables
!
real(DP):: xy_RainCumulus (0:imax-1, 1:jmax)
! ��羂顔��.
! Precipitation
real(DP):: xyz_DTempDtCumulus (0:imax-1, 1:jmax, 1:kmax)
! 羝�墾紊�����.
! Temperature tendency
real(DP):: xyz_DQVapDtCumulus (0:imax-1, 1:jmax, 1:kmax)
! 罸�羚水�����.
! Specific humidity tendency
real(DP):: xyz_QVapB (0:imax-1, 1:jmax, 1:kmax)
! 茯睡������羚�.
! Specific humidity before adjustment
real(DP):: xyz_TempB (0:imax-1, 1:jmax, 1:kmax)
! 茯睡�����県綺�.
! Temperature before adjustment
logical:: xy_Adjust (0:imax-1, 1:jmax)
! 篁���茯睡�������������?.
! Whether it was adjusted this time or not?
logical:: xy_AdjustB (0:imax-1, 1:jmax)
! ����茯睡�������������?.
! Whether it was adjusted last time or not?
real(DP):: xyz_DelPress(0:imax-1, 1:jmax, 1:kmax)
! $ \Delta p $
!
real(DP):: xyz_DelMass (0:imax-1, 1:jmax, 1:kmax)
! $ \Delta m $
!
real(DP):: xyz_QVapSat (0:imax-1, 1:jmax, 1:kmax)
! 蕋遵��罸�羚�.
! Saturation specific humidity.
real(DP):: xyr_ConvAdjustFactor(0:imax-1, 1:jmax, 0:kmax)
! $ \frac{1}{2} \frac{ R }{Cp}
! \frac{ p_{k} - p_{k+1} }{ p_{k+1/2} } $
real(DP):: TempEquivToExcEne
! Temperature equivalent to the excess moist static energy
! (Moist static energy difference devided by specific heat)
real(DP):: DelQ
real(DP):: DelTempUppLev
! k+1 ������韻��������茯睡�������羝�墾�������.
! Temperature variation by adjustment at k+1 level
real(DP):: DelTempLowLev
! k ������韻��������茯睡�������羝�墾�������.
! Temperature variation by adjustment at k level
real(DP):: DQVapSatDTempUppLev
! $ \DD{q^{*}} (k+1)}{T} $
real(DP):: DQVapSatDTempLowLev
! $ \DD{q^{*}} (k)}{T} $
real(DP):: GamUppLev
! $ \gamma_{k+1} = \frac{L}{C_p} \DP{q^{*}}{T}_{k+1} $
real(DP):: GamLowLev
! $ \gamma_{k} = \frac{L}{C_p} \DP{q^{*}}{T}_{k} $
logical:: Adjust
! 篁������������������綺��с��茯睡�������������?.
! Whether it was adjusted even once in global
! this time or not?
real(DP):: TempLowLevBefAdj ! Variables for check routine
real(DP):: TempUppLevBefAdj
real(DP):: QVapLowLevBefAdj
real(DP):: QVapUppLevBefAdj
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:: itr ! �ゃ�����若�激�с�恰�劫�������� DO ���若�����罐�����
! Work variables for DO loop in iteration direction
real(DP):: xy_DQVapSatDTempUppLev(0:imax-1, 1:jmax)
real(DP):: xy_DQVapSatDTempLowLev(0:imax-1, 1:jmax)
real(DP):: ExchangeMass
!
! Mass transport
real(DP):: ExchangeMassDenom
!
! Variable for mass transport calculation
real(DP):: ExchangeMassLowLim
!
! Lower limit of mass transport calculation
real(DP), parameter :: ExchangeMassLowLimTempDiff = 1.0d-5
!
! Lower limit of temperature difference
! between two layers for mass transport
! calculation
real(DP):: xyz_DelQVapCond(0:imax-1, 1:jmax, 1:kmax)
real(DP):: DelQVapCondLowLev
real(DP):: DelQVapCondUppLev
real(DP):: DelQVapCond2Levs
real(DP):: MassCor
real(DP):: xyz_RainCumulus (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_NegDDelLWDt (0:imax-1, 1:jmax)
real(DP) :: xyz_DDelLWDtCCPLV(0:imax-1, 1:jmax, 1:kmax)
! 絎�茵��� ; Executable statement
!
! ������腆肴�
! Initialization check
!
if ( .not. moist_conv_adjust_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
!!$ if ( .not. FlagUse ) return
! 荐�膊�����荐�羝���紮�
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! 茯睡��� "QVap", "Temp" ���絖�
! Store "QVap", "Temp" before adjustment
!
xyz_QVapB = xyz_QVap
xyz_TempB = xyz_Temp
! 蕋遵��罸�羚粋�膊�
! Calculate saturation specific humidity
!
xyz_QVapSat = xyz_CalcQVapSat( xyz_Temp, xyz_Press )
! Calculate some values used for moist convective adjustment
!
do k = 1, kmax
xyz_DelPress(:,:,k) = xyr_Press(:,:,k-1) - xyr_Press(:,:,k)
end do
xyz_DelMass = xyz_DelPress / Grav
! \frac{1}{2} \frac{ R }{Cp} \frac{ p_{k} - p_{k+1} }{ p_{k+1/2} }
!
! The value at k = 0 is not used.
k = 0
xyr_ConvAdjustFactor(:,:,k) = 0.0_DP
!
do k = 1, kmax-1
xyr_ConvAdjustFactor(:,:,k) = GasRDry / CpDry * ( xyz_Press(:,:,k) - xyz_Press(:,:,k+1) ) / xyr_Press(:,:,k) / 2.0_DP
end do
! The value at k = kmax is not used.
k = kmax
xyr_ConvAdjustFactor(:,:,k) = 0.0_DP
!
! Initialization
!
xyz_DelQVapCond = 0.0_DP
! 茯睡�
! Adjustment
!
xy_AdjustB = .true.
! 膵違��菴���
! Iteration
!
do itr = 1, ItrtMax
xy_Adjust = .false.
do k = 1, kmax-1
xy_DQVapSatDTempUppLev = xy_CalcDQVapSatDTemp( xyz_Temp(:,:,k+1), xyz_QVapSat(:,:,k+1) )
xy_DQVapSatDTempLowLev = xy_CalcDQVapSatDTemp( xyz_Temp(:,:,k ), xyz_QVapSat(:,:,k ) )
do j = 1, jmax
do i = 0, imax-1
if ( xy_AdjustB(i,j) ) then
! Temperature equivalent to the excess moist static energy
! (Moist static energy difference devided by specific heat)
!
TempEquivToExcEne = xyz_Temp(i,j,k) - xyz_Temp(i,j,k+1) + LatentHeat / CpDry * ( xyz_QVapSat(i,j,k) - xyz_QVapSat(i,j,k+1) ) - xyr_ConvAdjustFactor(i,j,k) * ( xyz_Temp(i,j,k) + xyz_Temp(i,j,k+1) )
! Check vertical gradient of moist static energy
!
if ( TempEquivToExcEne > AdjustCriterion(itr) ) then
! Check relative humidity
!
if ( ( xyz_QVap(i,j,k+1) / xyz_QVapSat(i,j,k+1) >= CrtlRH ) .and. ( xyz_QVap(i,j,k ) / xyz_QVapSat(i,j,k ) >= CrtlRH ) ) then
DelQ = xyz_DelPress(i,j,k ) * ( xyz_QVap(i,j,k ) - xyz_QVapSat(i,j,k ) ) + xyz_DelPress(i,j,k+1) * ( xyz_QVap(i,j,k+1) - xyz_QVapSat(i,j,k+1) )
DQVapSatDTempUppLev = xy_DQVapSatDTempUppLev(i,j)
DQVapSatDTempLowLev = xy_DQVapSatDTempLowLev(i,j)
GamUppLev = LatentHeat / CpDry * DQVapSatDTempUppLev
GamLowLev = LatentHeat / CpDry * DQVapSatDTempLowLev
DelTempUppLev = ( xyz_DelPress(i,j,k) * ( 1.0_DP + GamLowLev ) * TempEquivToExcEne + ( 1.0_DP + GamLowLev - xyr_ConvAdjustFactor(i,j,k) ) * LatentHeat / CpDry * DelQ ) / ( xyr_ConvAdjustFactor(i,j,k) * ( xyz_DelPress(i,j,k ) * ( 1.0_DP + GamLowLev ) - xyz_DelPress(i,j,k+1) * ( 1.0_DP + GamUppLev ) ) + ( 1.0_DP + GamLowLev ) * ( 1.0_DP + GamUppLev ) * ( xyz_DelPress(i,j,k) + xyz_DelPress(i,j,k+1) ) )
DelTempLowLev = ( LatentHeat / CpDry * DelQ - xyz_DelPress(i,j,k+1) * ( 1.0_DP + GamUppLev ) * DelTempUppLev ) / ( ( 1.0_DP + GamLowLev ) * xyz_DelPress(i,j,k) )
!=========
TempLowLevBefAdj = xyz_Temp(i,j,k )
TempUppLevBefAdj = xyz_Temp(i,j,k+1)
QVapLowLevBefAdj = xyz_QVap(i,j,k )
QVapUppLevBefAdj = xyz_QVap(i,j,k+1)
!=========
! 羝�墾���膀�
! Adjust temperature
!
xyz_Temp(i,j,k ) = xyz_Temp(i,j,k ) + DelTempLowLev
xyz_Temp(i,j,k+1) = xyz_Temp(i,j,k+1) + DelTempUppLev
! 罸�羚帥���膀�
! Adjust specific humidity
!
xyz_QVap(i,j,k ) = xyz_QVapSat(i,j,k ) + DQVapSatDTempLowLev * DelTempLowLev
xyz_QVap(i,j,k+1) = xyz_QVapSat(i,j,k+1) + DQVapSatDTempUppLev * DelTempUppLev
!
! Mass exchange
! Denominator
ExchangeMassDenom = CpDry * ( TempLowLevBefAdj - TempUppLevBefAdj ) - GasRDry * ( TempLowLevBefAdj + TempUppLevBefAdj ) / 2.0_DP / xyr_Press(i,j,k) * ( xyz_Press(i,j,k) - xyz_Press(i,j,k+1) ) + LatentHeat * ( QVapLowLevBefAdj - QVapUppLevBefAdj )
ExchangeMassLowLim = CpDry * ExchangeMassLowLimTempDiff
! If a static energy difference between two layers is smaller
! than a specified lower limit, momentum and mixing ratio are
! not mixed to ensure numerical stability.
! If the lower limit is zero, some calculations are unstable.
! (yot, 2013/10/02)
if ( ExchangeMassDenom > ExchangeMassLowLim ) then
ExchangeMass = - ( CpDry * DelTempLowLev + LatentHeat * ( xyz_QVap(i,j,k) - QVapLowLevBefAdj ) ) / ExchangeMassDenom * xyz_DelMass(i,j,k)
else
ExchangeMass = 0.0_DP
end if
! Limitation of amount of mass exchange not to
! reverse a gradient
ExchangeMass = min( ExchangeMass, xyz_DelMass(i,j,k) * xyz_DelMass(i,j,k+1) / ( xyz_DelMass(i,j,k) + xyz_DelMass(i,j,k+1) ) )
DelQVapCondLowLev = ( QVapUppLevBefAdj - QVapLowLevBefAdj ) * ExchangeMass / xyz_DelMass(i,j,k ) - ( xyz_QVap(i,j,k ) - QVapLowLevBefAdj )
DelQVapCondUppLev = - ( QVapUppLevBefAdj - QVapLowLevBefAdj ) * ExchangeMass / xyz_DelMass(i,j,k+1) - ( xyz_QVap(i,j,k+1) - QVapUppLevBefAdj )
! Check
DelQVapCond2Levs = DelQVapCondLowLev * xyz_DelMass(i,j,k ) + DelQVapCondUppLev * xyz_DelMass(i,j,k+1)
if ( DelQVapCond2Levs < 0.0_DP ) then
!!$ call MessageNotify( 'M', module_name, &
!!$ & 'Condensation amount is negative, %f.', &
!!$ & d = (/ DelQVapCond2Levs /) )
else
if ( DelQVapCondLowLev < 0.0_DP ) then
MassCor = - DelQVapCondLowLev * xyz_DelMass(i,j,k )
DelQVapCondLowLev = 0.0_DP
DelQVapCondUppLev = ( DelQVapCondUppLev * xyz_DelMass(i,j,k+1) - MassCor ) / xyz_DelMass(i,j,k+1)
end if
if ( DelQVapCondUppLev < 0.0_DP ) then
MassCor = - DelQVapCondUppLev * xyz_DelMass(i,j,k+1)
DelQVapCondLowLev = ( DelQVapCondLowLev * xyz_DelMass(i,j,k ) - MassCor ) / xyz_DelMass(i,j,k )
DelQVapCondUppLev = 0.0_DP
end if
end if
xyz_DelQVapCond(i,j,k ) = xyz_DelQVapCond(i,j,k ) + DelQVapCondLowLev
xyz_DelQVapCond(i,j,k+1) = xyz_DelQVapCond(i,j,k+1) + DelQVapCondUppLev
!=========
! check routine
!---------
!!$ write( 6, * ) '====='
!!$ write( 6, * ) xyz_Temp(i,j,k), xyz_Temp(i,j,k+1), xyz_QVap(i,j,k), xyz_QVap(i,j,k+1)
!!$ write( 6, * ) DelTempLowLev, DelTempUppLev
!!$ write( 6, * ) 'Energy difference before and after adjustment and each energy'
!!$ write( 6, * ) &
!!$ & ( CpDry * TempLowLevBefAdj + LatentHeat * QVapLowLevBefAdj ) &
!!$ & * xyz_DelPress(i,j,k ) / Grav &
!!$ & + ( CpDry * TempUppLevBefAdj + LatentHeat * QVapUppLevBefAdj ) &
!!$ & * xyz_DelPress(i,j,k+1) / Grav &
!!$ & - ( CpDry * xyz_Temp(i,j,k ) + LatentHeat * xyz_QVap(i,j,k ) ) &
!!$ & * xyz_DelPress(i,j,k ) / Grav &
!!$ & - ( CpDry * xyz_Temp(i,j,k+1) + LatentHeat * xyz_QVap(i,j,k+1) ) &
!!$ & * xyz_DelPress(i,j,k+1) / Grav, &
!!$ & ( CpDry * TempLowLevBefAdj + LatentHeat * QVapLowLevBefAdj ) &
!!$ & * xyz_DelPress(i,j,k ) / Grav, &
!!$ & ( CpDry * TempUppLevBefAdj + LatentHeat * QVapUppLevBefAdj ) &
!!$ & * xyz_DelPress(i,j,k+1) / Grav, &
!!$ & ( CpDry * xyz_Temp(i,j,k ) + LatentHeat * xyz_QVap(i,j,k ) ) &
!!$ & * xyz_DelPress(i,j,k ) / Grav, &
!!$ & ( CpDry * xyz_Temp(i,j,k+1) + LatentHeat * xyz_QVap(i,j,k+1) ) &
!!$ & * xyz_DelPress(i,j,k+1) / Grav
!!$ write( 6, * ) 'Difference of moist static energy after adjustment'
!!$ write( 6, * ) &
!!$ & ( CpDry * xyz_Temp(i,j,k ) + LatentHeat * xyz_QVap(i,j,k ) ) &
!!$ & - ( CpDry * xyz_Temp(i,j,k+1) + LatentHeat * xyz_QVap(i,j,k+1) ) &
!!$ & - CpDry * xyr_ConvAdjustFactor(i,j,k) &
!!$ & * ( xyz_Temp(i,j,k) + xyz_Temp(i,j,k+1) ), &
!!$ & ( CpDry * xyz_Temp(i,j,k ) + LatentHeat * xyz_QVap(i,j,k ) ), &
!!$ & ( CpDry * xyz_Temp(i,j,k+1) + LatentHeat * xyz_QVap(i,j,k+1) ), &
!!$ & - CpDry * xyr_ConvAdjustFactor(i,j,k) &
!!$ & * ( xyz_Temp(i,j,k) + xyz_Temp(i,j,k+1) )
!!$ write( 6, * ) 'Difference of water vapor amount before and after adjustment'
!!$ write( 6, * ) &
!!$ & - LatentHeat * ( xyz_QVap(i,j,k ) - QVapLowLevBefAdj ) &
!!$ & * xyz_DelPress(i,j,k ) / Grav, &
!!$ & - LatentHeat * ( xyz_QVap(i,j,k+1) - QVapUppLevBefAdj ) &
!!$ & * xyz_DelPress(i,j,k+1) / Grav
!=========
xyz_QVapSat(i,j,k ) = xyz_QVap(i,j,k )
xyz_QVapSat(i,j,k+1) = xyz_QVap(i,j,k+1)
! 茯睡�����������?
! Whether it was adjusted or not?
!
xy_Adjust(i,j) = .true.
end if
end if
end if
end do
end do
end do
Adjust = .false.
do i = 0, imax-1
do j = 1, jmax
xy_AdjustB(i,j) = xy_Adjust(i,j)
Adjust = Adjust .or. xy_Adjust(i,j)
end do
end do
if ( .not. Adjust ) exit
end do
call MoistConvAdjustChkCons( xyz_TempB, xyz_Temp, xyz_QVapB, xyz_QVap, xyz_DelQVapCond, xyz_DelMass )
! 罸�羚水�����, 羝�墾紊�����, ��羂顔������
! Calculate specific humidity tendency, temperature tendency, precipitation
!
xyz_DQVapDtCumulus = ( xyz_QVap - xyz_QVapB ) / ( 2.0_DP * DelTime )
xyz_DTempDtCumulus = ( xyz_Temp - xyz_TempB ) / ( 2.0_DP * DelTime )
! old
!!$ xyz_DQH2OLiqDt = ( xyz_QVapB - xyz_QVap ) / ( 2.0_DP * DelTime )
! new (2014/12/04)
xyz_DQH2OLiqDt = xyz_DelQVapCond / ( 2.0_DP * DelTime )
! avoid negative cloud amount
xyz_DQH2OLiqDt = max( xyz_DQH2OLiqDt, 0.0_DP )
!!$ xyz_RainCumulus = xyz_DQH2OLiqDt * xyz_DelPress / Grav
!!$ xy_RainCumulus = 0.0d0
!!$ do k = kmax, 1, -1
!!$ xy_RainCumulus = xy_RainCumulus + xyz_RainCumulus(:,:,k)
!!$ end do
!!$ j = jmax/2+1
!!$ do i = 0, imax-1
!!$ if ( xy_RainCumulus(i,j) /= 0.0d0 ) then
!!$ write( 6, * ) xy_RainCumulus(i,j)
!!$ end if
!!$ end do
!!$ write( 6, * ) '---'
!!$ xy_Rain = xy_Rain + xy_RainCumulus
! calculation for output (tentative)
xyz_RainCumulus = xyz_DQH2OLiqDt * xyz_DelPress / Grav
xy_RainCumulus = 0.0_DP
do k = kmax, 1, -1
xy_RainCumulus = xy_RainCumulus + xyz_RainCumulus(:,:,k)
end do
! ���鴻�������若�水�阪��
! History data output
!
call HistoryAutoPut( TimeN, 'RainCumulus', xy_RainCumulus * LatentHeat )
call HistoryAutoPut( TimeN, 'DTempDtCumulus', xyz_DTempDtCumulus )
call HistoryAutoPut( TimeN, 'DQVapDtCumulus', xyz_DQVapDtCumulus )
!!$ if ( present( xyz_DQH2OLiqDt ) ) then
!!$
!!$ xyz_DDelLWDtCCPLV = &
!!$ & + ( xyz_QVapB - xyz_QVap ) &
!!$ & * xyz_DelPress / Grav / ( 2.0d0 * DelTime )
!!$
!!$ ! Negative cloud production rate is filled with values in lower layers.
!!$ !
!!$ xy_NegDDelLWDt = 0.0d0
!!$ do k = kmax, 1, -1
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ xyz_DDelLWDtCCPLV(i,j,k) = xyz_DDelLWDtCCPLV(i,j,k) + xy_NegDDelLWDt(i,j)
!!$ if ( xyz_DDelLWDtCCPLV(i,j,k) < 0.0d0 ) then
!!$ xy_NegDDelLWDt(i,j) = xyz_DDelLWDtCCPLV(i,j,k)
!!$ xyz_DDelLWDtCCPLV(i,j,k) = 0.0d0
!!$ end if
!!$ end do
!!$ end do
!!$ end do
!!$
!!$ xyz_DQH2OLiqDt = xyz_DDelLWDtCCPLV / ( xyz_DelPress / Grav )
!!$
!!$ end if
! 荐�膊�����荐�羝������罩�
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine MoistConvAdjust