Class | relaxed_arakawa_schubert |
In: |
cumulus/relaxed_arakawa_schubert.f90
|
Note that Japanese and English are described in parallel.
Change temperature and specific humidity by using the Relaxed Arakawa-Schubert scheme
Lord, S. J., W. C. Chao, and A. Arakawa, Interaction of a cumulus cloud ensemble with the large-scale environment. Part IV: The discrete model, J. Atmos. Sci., 39, 104-113, 1992. Moorthi, S., and M. J. Suarez, Relaxed Arakawa-Schubert: A parameterization of moist convection for general circulation models, Mon. Wea. Rev., 120, 978-1002, 1992.
RelaxedArakawaSchubert : | 温度と比湿の調節 |
———————- : | ———— |
RelaxedArakawaSchubert : | Change temperature and specific humidity |
Subroutine : | |||||
xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||||
xyz_Press(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||||
xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in )
| ||||
xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||||
xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in )
| ||||
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_QH2OVap(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_DQVapDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_DQH2OLiqDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out ) | ||||
xyz_MoistConvDetTend(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out ), optional | ||||
xyz_MoistConvSubsidMassFlux(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out ), optional |
relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化.
Change temperature and specific humidity by relaxed Arakawa-Schubert scheme
subroutine RAS1DWrapperTesting( xy_SurfTemp, xyz_Press, xyr_Press, xyz_Exner, xyr_Exner, xyz_Temp, xyz_QH2OVap, xyz_DTempDt, xyz_DQVapDt, xyz_DQH2OLiqDt, xyz_MoistConvDetTend, xyz_MoistConvSubsidMassFlux ) ! ! relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化. ! ! Change temperature and specific humidity by relaxed Arakawa-Schubert scheme ! ! モジュール引用 ; 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, xyz_CalcDQVapSatDTemp ! Arakawa-Schubert scheme by Lord et al. (1982) ! Arakawa-Schubert scheme by Lord et al. (1982) ! use arakawa_schubert_L1982, only : ArakawaSchubertL1982CalcCWFCrtl ! 宣言文 ; Declaration statements ! real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax) ! Pressure 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 ) :: xyr_Exner (0:imax-1, 1:jmax, 0:kmax) ! Exner function real(DP), intent(inout) :: xyz_Temp (0:imax-1, 1:jmax, 1:kmax) ! Temperature real(DP), intent(inout) :: xyz_QH2OVap (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(inout) :: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax) ! 温度変化率. ! Temperature tendency real(DP), intent(inout) :: xyz_DQVapDt (0:imax-1, 1:jmax, 1:kmax) ! 比湿変化率. ! Specific humidity tendency real(DP), intent(out ) :: xyz_DQH2OLiqDt(0:imax-1, 1:jmax, 1:kmax) real(DP), intent(out ), optional :: xyz_MoistConvDetTend (0:imax-1, 1:jmax, 1:kmax) real(DP), intent(out ), optional :: xyz_MoistConvSubsidMassFlux(0:imax-1, 1:jmax, 1:kmax) ! 作業変数 ! Work variables ! real(DP) :: SurfTemp ! Pressure real(DP) :: z_Press (1:kmax) ! Pressure real(DP) :: r_Press (0:kmax) ! Pressure real(DP) :: z_Exner (1:kmax) ! Exner function real(DP) :: r_Exner (0:kmax) ! Exner function real(DP) :: z_Temp (1:kmax) ! Temperature real(DP) :: z_QH2OVap (1:kmax) ! $ q $ . 比湿. Specific humidity !!$ real(DP), intent(inout) :: xy_Rain (0:imax-1, 1:jmax) !!$ ! 降水量. !!$ ! Precipitation real(DP) :: z_DTempDt (1:kmax) ! 温度変化率. ! Temperature tendency real(DP) :: z_DQVapDt (1:kmax) ! 比湿変化率. ! Specific humidity tendency real(DP) :: z_DQH2OLiqDt(1:kmax) real(DP) :: z_MoistConvDetTend (1:kmax) real(DP) :: z_MoistConvSubsidMassFlux(1:kmax) 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_DelPress(0:imax-1, 1:jmax, 1:kmax) ! $ \Delta p $ ! real(DP) :: xyz_PotTemp (0:imax-1, 1:jmax, 1:kmax) ! Potential temperature ! real(DP) :: xyz_QH2OVapSat (0:imax-1, 1:jmax, 1:kmax) ! 飽和比湿. ! Saturation specific humidity. ! Dry and moist static energy in environment (Env) and cloud (Cld) ! real(DP) :: xyz_EnvDryStaticEne (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyr_EnvDryStaticEne (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyz_EnvMoistStaticEne (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyr_EnvMoistStaticEne (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyz_EnvMoistStaticEneSat(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyr_EnvMoistStaticEneSat(0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyr_CldMoistStaticEne (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xy_Kernel (0:imax-1, 1:jmax) ! Tendency of cloud work function by cumulus convection, kernel real(DP) :: xy_CWF (0:imax-1, 1:jmax) ! Cloud work function real(DP) :: xyz_CWF (0:imax-1, 1:jmax, 1:kmax) ! Cloud work function ! (variable for output) real(DP) :: xy_DCWFDtLS (0:imax-1, 1:jmax) ! Tendency of cloud work function by large scale motion real(DP) :: xyz_DCWFDtLS (0:imax-1, 1:jmax, 1:kmax) ! Tendency of cloud work function by large scale motion ! (variable for output) real(DP) :: xy_CldMassFluxBottom (0:imax-1, 1:jmax) ! Cloud mass flux at cloud bottom real(DP) :: xyz_Beta (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_BetaCldTop (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_Gamma (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_GammaDSE (0:imax-1, 1:jmax, 1:kmax) ! Tendency of dry static energy per unit mass flux real(DP) :: xyz_GammaMSE (0:imax-1, 1:jmax, 1:kmax) ! Tendency of moist static energy per unit mass flux real(DP) :: xyz_Mu (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_Eps (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_PressCldBase (0:imax-1, 1:jmax) ! Pressure of cloud base real(DP) :: xyz_CWFCrtl (0:imax-1, 1:jmax, 1:kmax) ! "Critical value" of cloud work function real(DP) :: xyz_RainFactor (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_EntParam (0:imax-1, 1:jmax) ! Entrainment factor real(DP) :: xyz_EntParam (0:imax-1, 1:jmax, 1:kmax) ! Entrainment factor (variable for output) real(DP) :: xy_EntParamLL (0:imax-1, 1:jmax) ! Entrainment factor for a cloud with top at one layer ! higher level real(DP) :: xy_EntParamUL (0:imax-1, 1:jmax) ! Entrainment factor for a cloud with top at one layer ! lower level ! Difference of normalized mass flux between layer interface real(DP) :: xyz_DelNormMassFlux (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_DelNormMassFluxCldTop(0:imax-1, 1:jmax) ! Normalized mass flux at layer interface and cloud top real(DP) :: xyr_NormMassFlux (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xy_NormMassFluxCldTop (0:imax-1, 1:jmax) ! Liquid water at cloud top real(DP) :: xy_CldQH2OLiqCldTop (0:imax-1, 1:jmax) ! Mass flux distribution function real(DP) :: xyz_MassFluxDistFunc (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_DelH2OMass (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_H2OMassB (0:imax-1, 1:jmax) real(DP) :: xy_H2OMassA (0:imax-1, 1:jmax) real(DP) :: xyz_RainCumulus (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_DDelLWDtCCPLV(0:imax-1, 1:jmax, 1:kmax) logical :: xy_FlagCrossSatEquivPotTemp(0:imax-1, 1:jmax) ! ! Flag showing whether a parcel in cloud has moist static ! energy larger than environment's real(DP) :: xyr_QH2OVapSat (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyr_TempAdiabAscent (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xy_SurfPotTemp (0:imax-1, 1:jmax) !!$ real(DP) :: xyz_TempAdiabAscent (0:imax-1, 1:jmax, 1:kmax) ! Variables for looking for top of mixed layer ! logical :: xy_FlagMixLayTopFound (0:imax-1, 1:jmax) integer :: xy_IndexMixLayTop (0:imax-1, 1:jmax) ! Variables for modification of cloud mass flux ! real(DP) :: xyz_QH2OVapTentative (0:imax-1, 1:jmax, 1:kmax) real(DP) :: CldMassFluxCorFactor real(DP) :: xy_CldMassFluxCorFactor(0:imax-1, 1:jmax) real(DP) :: xyz_QH2OVapB(0:imax-1, 1:jmax, 1:kmax) ! 調節前の比湿. ! Specific humidity before adjustment ! Flags for modification of ! logical :: xy_FlagKernelNegative (0:imax-1, 1:jmax) logical :: xy_FlagNegH2OLiqCldTop(0:imax-1, 1:jmax) ! Variables for subsidence mass flux between updrafts ! real(DP) :: DelNormMassFluxHalfLayer real(DP) :: NormMassFlux ! Variables for debug ! !!$ real(DP) :: xyz_DelVal(0:imax-1, 1:jmax, 1:kmax) !!$ real(DP) :: xy_SumValB(0:imax-1, 1:jmax) !!$ real(DP) :: xy_SumValA(0:imax-1, 1:jmax) !!$ real(DP) :: Ratio real(DP) :: xy_SumTmp(0:imax-1, 1:jmax) character(STRING) :: VarName real(DP) :: xyrz_CldTemp (0:imax-1, 1:jmax, 0:kmax, 1:kmax) real(DP) :: xyrz_CldQH2OVap(0:imax-1, 1:jmax, 0:kmax, 1:kmax) real(DP) :: xyrz_CldQH2OLiq(0:imax-1, 1:jmax, 0:kmax, 1:kmax) real(DP) :: rz_CldTemp (0:kmax, 1:kmax) real(DP) :: rz_CldQH2OVap(0:kmax, 1:kmax) real(DP) :: rz_CldQH2OLiq(0:kmax, 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 :: l integer :: m integer :: n ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. relaxed_arakawa_schubert_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 SurfTemp = xy_SurfTemp(i,j) do k = 1, kmax z_Press (k) = xyz_Press (i,j,k) z_Exner (k) = xyz_Exner (i,j,k) z_Temp (k) = xyz_Temp (i,j,k) z_QH2OVap(k) = xyz_QH2OVap(i,j,k) z_DTempDt(k) = xyz_DTempDt(i,j,k) z_DQVapDt(k) = xyz_DQVapDt(i,j,k) end do do k = 0, kmax r_Press (k) = xyr_Press (i,j,k) r_Exner (k) = xyr_Exner (i,j,k) end do call RAS1DTesting( SurfTemp, z_Press, r_Press, z_Exner, r_Exner, z_Temp, z_QH2OVap, z_DTempDt, z_DQVapDt, z_DQH2OLiqDt, z_MoistConvDetTend, z_MoistConvSubsidMassFlux, rz_CldTemp, rz_CldQH2OVap, rz_CldQH2OLiq ) do k = 1, kmax xyz_Temp (i,j,k) = z_Temp (k) xyz_QH2OVap (i,j,k) = z_QH2OVap (k) xyz_DTempDt (i,j,k) = z_DTempDt (k) xyz_DQVapDt (i,j,k) = z_DQVapDt (k) xyz_DQH2OLiqDt(i,j,k) = z_DQH2OLiqDt(k) end do if ( present( xyz_MoistConvDetTend ) ) then do k = 1, kmax xyz_MoistConvDetTend(i,j,k) = z_MoistConvDetTend(k) end do end if if ( present( xyz_MoistConvSubsidMassFlux ) ) then do k = 1, kmax xyz_MoistConvSubsidMassFlux(i,j,k) = z_MoistConvSubsidMassFlux(k) end do end if do l = 1, kmax do k = 0, kmax xyrz_CldTemp (i,j,k,l) = rz_CldTemp (k,l) xyrz_CldQH2OVap(i,j,k,l) = rz_CldQH2OVap(k,l) xyrz_CldQH2OLiq(i,j,k,l) = rz_CldQH2OLiq(k,l) end do end do end do end do ! calculation for output do k = 1, kmax xyz_DelPress(:,:,k) = xyr_Press(:,:,k-1) - xyr_Press(:,:,k) end do 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 ! ヒストリデータ出力 ! History data output ! call HistoryAutoPut( TimeN, 'RainCumulus' , xy_RainCumulus * LatentHeat ) call HistoryAutoPut( TimeN, 'DTempDtCumulus' , xyz_DTempDtCumulus ) call HistoryAutoPut( TimeN, 'DQVapDtCumulus' , xyz_DQVapDtCumulus ) call HistoryAutoPut( TimeN, 'RASMassFluxDistFunc', xyz_MassFluxDistFunc ) call HistoryAutoPut( TimeN, 'RASEntParam' , xyz_EntParam ) call HistoryAutoPut( TimeN, 'RASCWF' , xyz_CWF ) call HistoryAutoPut( TimeN, 'RASCWFCrtl' , xyz_CWFCrtl ) call HistoryAutoPut( TimeN, 'RASDCWFDtLS' , xyz_DCWFDtLS ) !!$ call HistoryAutoPut( TimeN, 'RASMixLayTopIndex' , real( xy_IndexMixLayTop ) ) do l = 1, kmax do k = 0, kmax do j = 1, jmax do i = 0, imax-1 if ( xyrz_CldTemp (i,j,k,l) == 1.0d100 ) xyrz_CldTemp (i,j,k,l) = 0.0_DP if ( xyrz_CldQH2OVap(i,j,k,l) == 1.0d100 ) xyrz_CldQH2OVap(i,j,k,l) = 0.0_DP if ( xyrz_CldQH2OLiq(i,j,k,l) == 1.0d100 ) xyrz_CldQH2OLiq(i,j,k,l) = 0.0_DP end do end do end do end do do k = 1, kmax write( VarName, '(a,i3.3)' ) 'RASCldTemp', k call HistoryAutoPut( TimeN, VarName, xyrz_CldTemp (:,:,:,k) ) write( VarName, '(a,i3.3)' ) 'RASCldQH2OVap', k call HistoryAutoPut( TimeN, VarName, xyrz_CldQH2OVap(:,:,:,k) ) write( VarName, '(a,i3.3)' ) 'RASCldQH2OLiq', k call HistoryAutoPut( TimeN, VarName, xyrz_CldQH2OLiq(:,:,:,k) ) end do !!$ if ( present( xyz_DQH2OLiqDt ) ) then !!$ !!$ ! unit is kg m-2 s-1 !!$ xyz_DDelLWDtCCPLV = xyz_RainCumulus !!$ !!$ ! 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 !!$ !!$ ! unit is s-1 !!$ xyz_DQH2OLiqDt = xyz_DDelLWDtCCPLV / ( xyz_DelPress / Grav ) !!$ !!$ end if ! 計算時間計測一時停止 ! Pause measurement of computation time ! call TimesetClockStop( module_name ) end subroutine RAS1DWrapperTesting
Subroutine : | |||||
xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||||
xyz_Press(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||||
xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in )
| ||||
xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||||
xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in )
| ||||
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_QH2OVap(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_DQVapDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_DQH2OLiqDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out ) | ||||
xyz_MoistConvDetTend(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out ), optional | ||||
xyz_MoistConvSubsidMassFlux(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out ), optional |
relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化.
Change temperature and specific humidity by relaxed Arakawa-Schubert scheme
subroutine RelaxedArakawaSchubert( xy_SurfTemp, xyz_Press, xyr_Press, xyz_Exner, xyr_Exner, xyz_Temp, xyz_QH2OVap, xyz_DTempDt, xyz_DQVapDt, xyz_DQH2OLiqDt, xyz_MoistConvDetTend, xyz_MoistConvSubsidMassFlux ) ! ! relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化. ! ! Change temperature and specific humidity by relaxed Arakawa-Schubert scheme ! ! モジュール引用 ; 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, xyz_CalcDQVapSatDTemp ! Arakawa-Schubert scheme by Lord et al. (1982) ! Arakawa-Schubert scheme by Lord et al. (1982) ! use arakawa_schubert_L1982, only : ArakawaSchubertL1982CalcCWFCrtl ! 宣言文 ; Declaration statements ! real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax) ! Pressure 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 ) :: xyr_Exner (0:imax-1, 1:jmax, 0:kmax) ! Exner function real(DP), intent(inout) :: xyz_Temp (0:imax-1, 1:jmax, 1:kmax) ! Temperature real(DP), intent(inout) :: xyz_QH2OVap (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(inout) :: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax) ! 温度変化率. ! Temperature tendency real(DP), intent(inout) :: xyz_DQVapDt (0:imax-1, 1:jmax, 1:kmax) ! 比湿変化率. ! Specific humidity tendency real(DP), intent(out ) :: xyz_DQH2OLiqDt(0:imax-1, 1:jmax, 1:kmax) real(DP), intent(out ), optional :: xyz_MoistConvDetTend (0:imax-1, 1:jmax, 1:kmax) real(DP), intent(out ), optional :: xyz_MoistConvSubsidMassFlux(0:imax-1, 1:jmax, 1:kmax) ! 作業変数 ! Work variables ! real(DP) :: xyz_Height (0:imax-1, 1:jmax, 1:kmax) ! ! Height real(DP) :: xyr_Height (0:imax-1, 1:jmax, 0:kmax) ! ! Height 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_DelPress(0:imax-1, 1:jmax, 1:kmax) ! $ \Delta p $ ! real(DP) :: xyz_PotTemp (0:imax-1, 1:jmax, 1:kmax) ! Potential temperature ! real(DP) :: xyz_QH2OVapSat (0:imax-1, 1:jmax, 1:kmax) ! 飽和比湿. ! Saturation specific humidity. ! Dry and moist static energy in environment (Env) and cloud (Cld) ! real(DP) :: xyz_EnvDryStaticEne (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyr_EnvDryStaticEne (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyz_EnvMoistStaticEne (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyr_EnvMoistStaticEne (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyz_EnvMoistStaticEneSat(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyr_EnvMoistStaticEneSat(0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyr_CldMoistStaticEne (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xy_Kernel (0:imax-1, 1:jmax) ! Tendency of cloud work function by cumulus convection, kernel real(DP) :: xy_CWF (0:imax-1, 1:jmax) ! Cloud work function real(DP) :: xyz_CWF (0:imax-1, 1:jmax, 1:kmax) ! Cloud work function ! (variable for output) real(DP) :: xy_DCWFDtLS (0:imax-1, 1:jmax) ! Tendency of cloud work function by large scale motion real(DP) :: xyz_DCWFDtLS (0:imax-1, 1:jmax, 1:kmax) ! Tendency of cloud work function by large scale motion ! (variable for output) real(DP) :: xy_CldMassFluxBottom (0:imax-1, 1:jmax) ! Cloud mass flux at cloud bottom real(DP) :: xyz_Beta (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_BetaCldTop (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_Gamma (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_GammaDSE (0:imax-1, 1:jmax, 1:kmax) ! Tendency of dry static energy per unit mass flux real(DP) :: xyz_GammaMSE (0:imax-1, 1:jmax, 1:kmax) ! Tendency of moist static energy per unit mass flux real(DP) :: xyz_Mu (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_Eps (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_PressCldBase (0:imax-1, 1:jmax) ! Pressure of cloud base real(DP) :: xyz_CWFCrtl (0:imax-1, 1:jmax, 1:kmax) ! "Critical value" of cloud work function real(DP) :: xyz_RainFactor (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_EntParam (0:imax-1, 1:jmax) ! Entrainment factor real(DP) :: xyz_EntParam (0:imax-1, 1:jmax, 1:kmax) ! Entrainment factor (variable for output) real(DP) :: xy_EntParamLL (0:imax-1, 1:jmax) ! Entrainment factor for a cloud with top at one layer ! higher level real(DP) :: xy_EntParamUL (0:imax-1, 1:jmax) ! Entrainment factor for a cloud with top at one layer ! lower level ! Difference of normalized mass flux between layer interface real(DP) :: xyz_DelNormMassFlux (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_DelNormMassFluxCldTop(0:imax-1, 1:jmax) ! Normalized mass flux at layer interface and cloud top real(DP) :: xyr_NormMassFlux (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xy_NormMassFluxCldTop (0:imax-1, 1:jmax) ! Liquid water at cloud top real(DP) :: xy_CldQH2OLiqCldTop (0:imax-1, 1:jmax) ! Mass flux distribution function real(DP) :: xyz_MassFluxDistFunc (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_DelH2OMass (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_H2OMassB (0:imax-1, 1:jmax) real(DP) :: xy_H2OMassA (0:imax-1, 1:jmax) 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) logical :: xy_FlagCrossSatEquivPotTemp(0:imax-1, 1:jmax) ! ! Flag showing whether a parcel in cloud has moist static ! energy larger than environment's real(DP) :: xyr_QH2OVapSat (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyr_TempAdiabAscent (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xy_SurfPotTemp (0:imax-1, 1:jmax) !!$ real(DP) :: xyz_TempAdiabAscent (0:imax-1, 1:jmax, 1:kmax) ! Variables for looking for top of mixed layer ! logical :: xy_FlagMixLayTopFound (0:imax-1, 1:jmax) integer :: xy_IndexMixLayTop (0:imax-1, 1:jmax) ! Variables for modification of cloud mass flux ! real(DP) :: xyz_QH2OVapTentative (0:imax-1, 1:jmax, 1:kmax) real(DP) :: CldMassFluxCorFactor real(DP) :: xy_CldMassFluxCorFactor(0:imax-1, 1:jmax) real(DP) :: xyz_TempB (0:imax-1, 1:jmax, 1:kmax) ! 調節前の温度. ! Temperature before adjustment real(DP) :: xyz_QH2OVapB(0:imax-1, 1:jmax, 1:kmax) ! 調節前の比湿. ! Specific humidity before adjustment ! Flags for modification of ! logical :: xy_FlagKernelNegative (0:imax-1, 1:jmax) logical :: xy_FlagNegH2OLiqCldTop(0:imax-1, 1:jmax) ! Variables for subsidence mass flux between updrafts ! real(DP) :: DelNormMassFluxHalfLayer real(DP) :: NormMassFlux ! Variables for debug ! !!$ real(DP) :: xyz_DelVal(0:imax-1, 1:jmax, 1:kmax) !!$ real(DP) :: xy_SumValB(0:imax-1, 1:jmax) !!$ real(DP) :: xy_SumValA(0:imax-1, 1:jmax) !!$ real(DP) :: Ratio real(DP) :: xy_SumTmp(0:imax-1, 1:jmax) 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 :: l integer :: m integer :: n ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. relaxed_arakawa_schubert_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 計算時間計測開始 ! Start measurement of computation time ! call TimesetClockStart( module_name ) ! 調節前 "Temp", "QH2OVap" の保存 ! Store "Temp", "QH2OVap" before adjustment ! xyz_TempB = xyz_Temp xyz_QH2OVapB = xyz_QH2OVap ! Preparation of variables ! ! ! Auxiliary variables ! Pressure difference between upper and lower interface of the layer do k = 1, kmax xyz_DelPress(:,:,k) = xyr_Press(:,:,k-1) - xyr_Press(:,:,k) end do ! beta do k = 1, kmax xyz_Beta(:,:,k) = CpDry / Grav * ( xyr_Exner(:,:,k-1) - xyr_Exner(:,:,k) ) end do do k = 1, kmax xyz_BetaCldTop(:,:,k) = CpDry / Grav * ( xyr_Exner(:,:,k-1) - xyz_Exner(:,:,k) ) end do ! ! Search for top of mixed layer (lifting condensation level) based on ! a description in p.684 of Arakawa and Shubert (1974). ! call RelaxedArakawaSchubertHeight( xyz_Temp, xyz_Exner, xyz_Beta, xyz_BetaCldTop, xyz_Height, xyr_Height ) ! !==================================== ! !!$ xyz_TempAdiabAscent(:,:,1) = xyz_Temp(:,:,1) !!$ do k = 2, kmax !!$ xyz_TempAdiabAscent(:,:,k) = & !!$ & xyz_Temp(:,:,1) - Grav / CpDry * ( xyz_Height(:,:,k) - xyz_Height(:,:,1) ) !!$ end do !!$ xyz_TempAdiabAscent = max( xyz_TempAdiabAscent, 1.0_DP ) !!$ xyz_QH2OVapSat = xyz_CalcQVapSat( xyz_TempAdiabAscent, xyz_Press ) !!$ xy_IndexMixLayTop = 1 !!$ xy_FlagMixLayTopFound = .false. !!$ do k = 2, kmax !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( ( xyz_QH2OVap(i,j,1) >= xyz_QH2OVapSat(i,j,k) ) .and. & !!$ & ( .not. xy_FlagMixLayTopFound(i,j) ) ) then !!$ xy_IndexMixLayTop (i,j) = k - 1 !!$ xy_FlagMixLayTopFound(i,j) = .true. !!$ end if !!$ end do !!$ end do !!$ end do ! !------------------------------------ ! !!$ xyr_TempAdiabAscent(:,:,0) = xy_SurfTemp !!$ do k = 1, kmax !!$ xyr_TempAdiabAscent(:,:,k) = & !!$ & xy_SurfTemp - Grav / CpDry * ( xyr_Height(:,:,k) - 0.0_DP ) !!$ end do !!$ xyr_TempAdiabAscent = max( xyr_TempAdiabAscent, 1.0_DP ) !!$ xyr_TempAdiabAscent(:,:,0) = xy_SurfTemp xy_SurfPotTemp = xy_SurfTemp / xyr_Exner(:,:,0) do k = 1, kmax xyr_TempAdiabAscent(:,:,k) = xy_SurfPotTemp * xyr_Exner(:,:,k) end do ! xyr_QH2OVapSat(:,:,0 ) = 1.0d100 xyr_QH2OVapSat(:,:,1:kmax-1) = xyz_CalcQVapSat( xyr_TempAdiabAscent(:,:,1:kmax-1), xyr_Press(:,:,1:kmax-1) ) xyr_QH2OVapSat(:,:,kmax ) = xyr_QH2OVapSat(:,:,kmax-1) ! xy_IndexMixLayTop = 1 xy_FlagMixLayTopFound = .false. do k = 2, kmax do j = 1, jmax do i = 0, imax-1 if ( ( xyz_QH2OVap(i,j,1) >= xyr_QH2OVapSat(i,j,k) ) .and. ( .not. xy_FlagMixLayTopFound(i,j) ) ) then xy_IndexMixLayTop (i,j) = k - 1 xy_FlagMixLayTopFound(i,j) = .true. end if end do end do end do ! !==================================== ! do j = 1, jmax do i = 0, imax-1 if ( .not. xy_FlagMixLayTopFound(i,j) ) then xy_IndexMixLayTop(i,j) = kmax - 1 end if end do end do ! ! Critical cloud work function ! if ( FlagZeroCrtlCWF ) then xyz_CWFCrtl = 0.0_DP else do j = 1, jmax do i = 0, imax-1 xy_PressCldBase(i,j) = xyr_Press(i,j,xy_IndexMixLayTop(i,j)) end do end do call ArakawaSchubertL1982CalcCWFCrtl( xy_PressCldBase, xyz_Press, xyz_CWFCrtl ) end if ! ! Rain conversion factor ! if ( RainConversionFactor < 0.0_DP ) then do k = 1, kmax do j = 1, jmax do i = 0, imax-1 if ( xyz_Press(i,j,k) < 500.0d2 ) then xyz_RainFactor(i,j,k) = 1.0_DP else if ( xyz_Press(i,j,k) < 800.0d2 ) then xyz_RainFactor(i,j,k) = 0.8_DP + ( 800.0d2 - xyz_Press(i,j,k) ) / 1500.0d2 else xyz_RainFactor(i,j,k) = 0.8_DP end if end do end do end do else xyz_RainFactor = RainConversionFactor end if xyz_RainCumulus (:,:,1) = 0.0_DP xyz_EntParam (:,:,1) = 0.0_DP xyz_CWF (:,:,1) = 0.0_DP xyz_DCWFDtLS (:,:,1) = 0.0_DP xyz_MassFluxDistFunc(:,:,1) = 0.0_DP if ( present( xyz_MoistConvDetTend ) ) then xyz_MoistConvDetTend(:,:,1) = 0.0_DP end if if ( present( xyz_MoistConvSubsidMassFlux ) ) then ! Subsidence mass flux between the updrafts ! Initialization ! xyz_MoistConvSubsidMassFlux = 0.0_DP end if loop_cloud_top : do l = 2, kmax call RelaxedArakawaSchubertHeight( xyz_Temp, xyz_Exner, xyz_Beta, xyz_BetaCldTop, xyz_Height, xyr_Height ) ! Potential temperature ! xyz_PotTemp = xyz_Temp / xyz_Exner ! Saturation mixing ratio ! xyz_QH2OVapSat = xyz_CalcQVapSat( xyz_Temp, xyz_Press ) ! Calculation of dry and moist static energies ! xyz_EnvDryStaticEne = CpDry * xyz_Temp + Grav * xyz_Height xyz_EnvMoistStaticEne = xyz_EnvDryStaticEne + LatentHeat * xyz_QH2OVap ! k = 0 xyr_EnvDryStaticEne (:,:,k) = 1.0d100 xyr_EnvMoistStaticEne(:,:,k) = 1.0d100 do k = 1, kmax-1 xyr_EnvDryStaticEne (:,:,k) = ( xyz_EnvDryStaticEne (:,:,k) + xyz_EnvDryStaticEne (:,:,k+1) ) / 2.0_DP xyr_EnvMoistStaticEne(:,:,k) = ( xyz_EnvMoistStaticEne(:,:,k) + xyz_EnvMoistStaticEne(:,:,k+1) ) / 2.0_DP end do k = kmax xyr_EnvDryStaticEne (:,:,k) = xyz_EnvDryStaticEne (:,:,k) xyr_EnvMoistStaticEne(:,:,k) = xyz_EnvMoistStaticEne(:,:,k) ! Calculation of saturated moist static energy ! xyz_EnvMoistStaticEneSat = xyz_EnvDryStaticEne + LatentHeat * xyz_QH2OVapSat ! k = 0 xyr_EnvMoistStaticEneSat(:,:,k) = 1.0d100 do k = 1, kmax-1 xyr_EnvMoistStaticEneSat(:,:,k) = ( xyz_EnvMoistStaticEneSat(:,:,k) + xyz_EnvMoistStaticEneSat(:,:,k+1) ) / 2.0_DP end do k = kmax xyr_EnvMoistStaticEneSat(:,:,k) = xyz_EnvMoistStaticEneSat(:,:,k) ! Auxiliary variables ! xyz_Gamma = LatentHeat / CpDry * xyz_CalcDQVapSatDTemp( xyz_Temp, xyz_QH2OVapSat ) ! k = 1 xyz_Mu (:,:,k) = 1.0d100 xyz_Eps(:,:,k) = 1.0d100 do k = 2, kmax xyz_Mu (:,:,k) = ( xyz_Exner(:,:,k ) - xyr_Exner(:,:,k) ) / ( xyz_Exner(:,:,k) * ( 1.0_DP + xyz_Gamma(:,:,k) ) ) xyz_Eps(:,:,k) = ( xyr_Exner(:,:,k-1) - xyz_Exner(:,:,k) ) / ( xyz_Exner(:,:,k) * ( 1.0_DP + xyz_Gamma(:,:,k) ) ) end do ! Entrainment parameter ! call RelaxedArakawaSchubertEntParam( l, xyz_Beta, xyz_BetaCldTop, xyz_PotTemp, xyz_EnvMoistStaticEne, xyz_EnvMoistStaticEneSat, xy_IndexMixLayTop, xy_EntParam ) if ( l >= 3 ) then call RelaxedArakawaSchubertEntParam( l-1, xyz_Beta, xyz_BetaCldTop, xyz_PotTemp, xyz_EnvMoistStaticEne, xyz_EnvMoistStaticEneSat, xy_IndexMixLayTop, xy_EntParamLL ) else xy_EntParamLL = 1.0d100 end if if ( l <= kmax-1 ) then call RelaxedArakawaSchubertEntParam( l+1, xyz_Beta, xyz_BetaCldTop, xyz_PotTemp, xyz_EnvMoistStaticEne, xyz_EnvMoistStaticEneSat, xy_IndexMixLayTop, xy_EntParamUL ) else xy_EntParamUL = 1.0d100 end if ! for output xyz_EntParam(:,:,l) = xy_EntParam ! Difference of normalized mass flux ! ! difference of normalized mass flux between layer bottom and top ! xyz_DelNormMassFlux(:,:,1) = 1.0d100 do k = 2, l-1 xyz_DelNormMassFlux(:,:,k) = - xy_EntParam * xyz_Beta(:,:,k) * xyz_PotTemp(:,:,k) end do do k = l, kmax xyz_DelNormMassFlux(:,:,k) = 1.0d100 end do ! ! difference of normalized mass flux between layer bottom and mid-point ! xy_DelNormMassFluxCldTop = - xy_EntParam * xyz_BetaCldTop(:,:,l) * xyz_PotTemp(:,:,l) ! Normalized mass flux ! ! normalized mass flux at layer interface ! xyr_NormMassFlux(:,:,0) = 0.0_DP do k = 1, l-1 do j = 1, jmax do i = 0, imax-1 if ( k < xy_IndexMixLayTop(i,j) ) then xyr_NormMassFlux(i,j,k) = 0.0_DP else if ( k == xy_IndexMixLayTop(i,j) ) then xyr_NormMassFlux(i,j,k) = 1.0_DP else xyr_NormMassFlux(i,j,k) = xyr_NormMassFlux(i,j,k-1) - xyz_DelNormMassFlux(i,j,k) end if end do end do end do do k = l, kmax xyr_NormMassFlux(:,:,k) = 0.0_DP end do ! ! normalized mass flux at cloud top (at layer mid-point) ! xy_NormMassFluxCldTop = xyr_NormMassFlux(:,:,l-1) - xy_DelNormMassFluxCldTop ! Liquid water content at top of clouds ! If l is less than xy_IndexMixLayTop(i,j), i.e. the cloud top is below top of ! mixed layer, xy_SumTmp is zero, then, xy_CldQH2OLiqCldTop is also zero. ! do j = 1, jmax do i = 0, imax-1 if ( l > xy_IndexMixLayTop(i,j) ) then xy_SumTmp(i,j) = xyz_QH2OVap(i,j,xy_IndexMixLayTop(i,j)) do k = xy_IndexMixLayTop(i,j)+1, l-1 xy_SumTmp(i,j) = xy_SumTmp(i,j) - xyz_DelNormMassFlux(i,j,k) * xyz_QH2OVap(i,j,k) end do xy_SumTmp(i,j) = xy_SumTmp(i,j) - xy_DelNormMassFluxCldTop(i,j) * xyz_QH2OVap(i,j,l) else xy_SumTmp(i,j) = 0.0_DP end if end do end do xy_CldQH2OLiqCldTop = xy_SumTmp / ( xy_NormMassFluxCldTop + 1.0d-100 ) - xyz_QH2OVapSat(:,:,l) ! Check whether kernel is positive or negative. ! do j = 1, jmax do i = 0, imax-1 if ( xy_CldQH2OLiqCldTop(i,j) < 0.0_DP ) then xy_FlagNegH2OLiqCldTop(i,j) = .true. else xy_FlagNegH2OLiqCldTop(i,j) = .false. end if end do end do ! avoid negative value xy_CldQH2OLiqCldTop = max( xy_CldQH2OLiqCldTop, 0.0_DP ) ! Moist static energy in clouds ! xyr_CldMoistStaticEne(:,:,0) = 1.0d100 do k = 1, l-1 do j = 1, jmax do i = 0, imax-1 if ( k < xy_IndexMixLayTop(i,j) ) then xyr_CldMoistStaticEne(i,j,k) = 1.0d100 else if ( k == xy_IndexMixLayTop(i,j) ) then xyr_CldMoistStaticEne(i,j,k) = xyz_EnvMoistStaticEne(i,j,xy_IndexMixLayTop(i,j)) else xyr_CldMoistStaticEne(i,j,k) = ( xyr_NormMassFlux(i,j,k-1) * xyr_CldMoistStaticEne(i,j,k-1) - xyz_DelNormMassFlux(i,j,k) * xyz_EnvMoistStaticEne(i,j,k) ) / xyr_NormMassFlux(i,j,k) end if end do end do end do do k = l, kmax xyr_CldMoistStaticEne(:,:,k) = 1.0d100 end do !############################################### ! Check whether a parcel in cloud has moist static energy larger than environment's ! !!$ xy_FlagCrossSatEquivPotTemp = .false. !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ do k = xy_IndexMixLayTop(i,j), l-1 !!$ if ( xyr_EnvMoistStaticEneSat(i,j,k) < xyr_CldMoistStaticEne(i,j,k) ) then !!$ xy_FlagCrossSatEquivPotTemp(i,j) = .true. !!$ end if !!$ end do !!$ end do !!$ end do !############################################### ! Cloud work function ! xy_CWF = 0.0_DP do k = 2, l-1 do j = 1, jmax do i = 0, imax-1 if ( k > xy_IndexMixLayTop(i,j) ) then xy_CWF(i,j) = xy_CWF(i,j) + xyz_Mu (i,j,k) * xyr_NormMassFlux(i,j,k ) * ( xyr_CldMoistStaticEne(i,j,k ) - xyz_EnvMoistStaticEneSat(i,j,k) ) + xyz_Eps(i,j,k) * xyr_NormMassFlux(i,j,k-1) * ( xyr_CldMoistStaticEne(i,j,k-1) - xyz_EnvMoistStaticEneSat(i,j,k) ) end if end do end do end do k = l do j = 1, jmax do i = 0, imax-1 if ( k > xy_IndexMixLayTop(i,j) ) then xy_CWF(i,j) = xy_CWF(i,j) + xyz_Eps(i,j,k) * xyr_NormMassFlux(i,j,k-1) * ( xyr_CldMoistStaticEne(i,j,k-1) - xyz_EnvMoistStaticEneSat(i,j,k) ) end if end do end do ! for save xyz_CWF(:,:,l) = xy_CWF ! Time derivative of cloud work function by large scale motion ! xy_DCWFDtLS = ( xy_CWF - xyz_CWFCrtl(:,:,l) ) / ( 2.0_DP * DelTime ) ! for save xyz_DCWFDtLS(:,:,l) = xy_DCWFDtLS ! Tendency of dry static energy per unit mass flux ! do k = 1, l xyz_GammaDSE(:,:,k) = - Grav / xyz_DelPress(:,:,k) * ( xyr_NormMassFlux(:,:,k-1) * ( xyr_EnvDryStaticEne(:,:,k-1) - xyz_EnvDryStaticEne(:,:,k) ) + xyr_NormMassFlux(:,:,k ) * ( xyz_EnvDryStaticEne(:,:,k ) - xyr_EnvDryStaticEne(:,:,k) ) ) end do k = l xyz_GammaDSE(:,:,k) = xyz_GammaDSE(:,:,k) - Grav / xyz_DelPress(:,:,k) * LatentHeat * xy_CldQH2OLiqCldTop * xy_NormMassFluxCldTop * ( 1.0_DP - xyz_RainFactor(:,:,k) ) do k = l+1, kmax xyz_GammaDSE(:,:,k) = 0.0_DP end do ! Tendency of moist static energy per unit mass flux ! do k = 1, l xyz_GammaMSE(:,:,k) = - Grav / xyz_DelPress(:,:,k) * ( xyr_NormMassFlux(:,:,k-1) * ( xyr_EnvMoistStaticEne(:,:,k-1) - xyz_EnvMoistStaticEne(:,:,k) ) + xyr_NormMassFlux(:,:,k ) * ( xyz_EnvMoistStaticEne(:,:,k ) - xyr_EnvMoistStaticEne(:,:,k) ) ) end do k = l xyz_GammaMSE(:,:,k) = xyz_GammaMSE(:,:,k) + Grav / xyz_DelPress(:,:,k) * xy_NormMassFluxCldTop * ( xyz_EnvMoistStaticEneSat(:,:,k) - xyz_EnvMoistStaticEne(:,:,k) ) do k = l+1, kmax xyz_GammaMSE(:,:,k) = 0.0_DP end do ! Kernel, time derivative of cloud work function by cumulus convection per unit ! mass flux ! do j = 1, jmax do i = 0, imax-1 xy_Kernel(i,j) = xyz_Eps(i,j,xy_IndexMixLayTop(i,j)+1) * xyz_GammaMSE(i,j,xy_IndexMixLayTop(i,j)) - xyz_Eps(i,j,l) * xyr_NormMassFlux(i,j,l-1) * ( 1.0_DP + xyz_Gamma(i,j,l) ) * xyz_GammaDSE(i,j,l) do n = xy_IndexMixLayTop(i,j)+1, l-1 xy_SumTmp(i,j) = 0.0_DP do m = xy_IndexMixLayTop(i,j)+1, n xy_SumTmp(i,j) = xy_SumTmp(i,j) + xyz_DelNormMassFlux(i,j,m) * xyz_GammaMSE(i,j,m) end do xy_Kernel(i,j) = xy_Kernel(i,j) + ( xyz_Eps(i,j,n+1) + xyz_Mu(i,j,n) ) * ( xyz_GammaMSE(i,j,xy_IndexMixLayTop(i,j)) - xy_SumTmp(i,j) ) - ( xyz_Eps(i,j,n) * xyr_NormMassFlux(i,j,n-1) + xyz_Mu (i,j,n) * xyr_NormMassFlux(i,j,n ) ) * ( 1.0_DP + xyz_Gamma(i,j,n) ) * xyz_GammaDSE(i,j,n) end do end do end do ! Check whether kernel is positive or negative. ! do j = 1, jmax do i = 0, imax-1 if ( xy_Kernel(i,j) < 0.0_DP ) then xy_FlagKernelNegative(i,j) = .true. else xy_FlagKernelNegative(i,j) = .false. end if end do end do ! Load et al. (1982), p.108 xy_Kernel = min( xy_Kernel, -5.0d-3 ) ! Cloud mass flux at cloud bottom ! xy_CldMassFluxBottom = - xy_DCWFDtLS / xy_Kernel ! ! mass flux has to be zero or positive xy_CldMassFluxBottom = max( xy_CldMassFluxBottom, 0.0_DP ) ! mass flux is zero if entrainment parameter is zero or negative do j = 1, jmax do i = 0, imax-1 if ( xy_EntParam(i,j) <= 0.0_DP ) then xy_CldMassFluxBottom(i,j) = 0.0_DP end if end do end do !!$ ! mass flux is zero if it is below lifting condensation level !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( .not. xy_FlagCrossSatEquivPotTemp(i,j) ) then !!$ xy_CloudMassFluxBottom(i,j) = 0.0_DP !!$ end if !!$ end do !!$ end do ! mass flux is zero if the LNB is unstable for updrafts ! (i.e., if the parcel is positively buoyant just above the LNB). ! See Lord et al. (1982), p.112, for more details. ! Strictly speaking, the process below is different from that ! proposed by Lord et al. (1982). Lord et al. (1982) compare ! entrainment parameters at 3 levels. But, entrainment ! parameters at 2 levels are compared below, because comparison ! of values between 2 levels seems to be sufficient. !!$ if ( ( 3 <= l ) .and. ( l <= kmax-1 ) ) then !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( ( xy_EntParamLL(i,j) < xy_EntParam (i,j) ) .and. & !!$ & ( xy_EntParam (i,j) < xy_EntParamUL(i,j) ) ) then !!$ if ( ( xy_EntParamLL(i,j) > 0.0_DP ) .and. & !!$ & ( xy_EntParam (i,j) > 0.0_DP ) .and. & !!$ & ( xy_EntParamUL(i,j) > 0.0_DP ) ) then !!$ xy_CloudMassFluxBottom(i,j) = 0.0_DP !!$ end if !!$ end if !!$ end do !!$ end do !!$ end if do j = 1, jmax do i = 0, imax-1 !!$ if ( xy_IndexMixLayTop(i,j) == l ) then !!$ if ( ( xy_EntParam (i,j) > 0.0_DP ) .and. & !!$ & ( xy_EntParamUL(i,j) > 0.0_DP ) ) then !!$ if ( xy_EntParam (i,j) < xy_EntParamUL(i,j) ) then !!$ xy_CloudMassFluxBottom(i,j) = 0.0_DP !!$ end if !!$ end if !!$ else if ( ( xy_IndexMixLayTop(i,j) < l ) .and. ( l <= kmax-1 ) ) then !!$ if ( ( xy_EntParamLL(i,j) > 0.0_DP ) .and. & !!$ & ( xy_EntParam (i,j) > 0.0_DP ) .and. & !!$ & ( xy_EntParamUL(i,j) > 0.0_DP ) ) then !!$ if ( ( xy_EntParamLL(i,j) < xy_EntParam (i,j) ) .and. & !!$ & ( xy_EntParam (i,j) < xy_EntParamUL(i,j) ) ) then if ( ( xy_IndexMixLayTop(i,j) <= l ) .and. ( l <= kmax-1 ) ) then if ( ( xy_EntParam (i,j) > 0.0_DP ) .and. ( xy_EntParamUL(i,j) > 0.0_DP ) ) then if ( xy_EntParam (i,j) < xy_EntParamUL(i,j) ) then xy_CldMassFluxBottom(i,j) = 0.0_DP end if end if end if end do end do ! ! mass flux is zero unless kernel is negative ! do j = 1, jmax do i = 0, imax-1 if ( .not. xy_FlagKernelNegative(i,j) ) then xy_CldMassFluxBottom(i,j) = 0.0_DP end if end do end do ! ! mass flux is zero if liquid water at a cloud top is negative ! do j = 1, jmax do i = 0, imax-1 if ( xy_FlagNegH2OLiqCldTop(i,j) ) then xy_CldMassFluxBottom(i,j) = 0.0_DP end if end do end do ! ! multiply factor ! xy_CldMassFluxBottom = xy_CldMassFluxBottom * min( 2.0_DP * DelTime / AdjTimeConst, 1.0_DP ) ! ! for output xyz_MassFluxDistFunc(:,:,l) = xy_CldMassFluxBottom ! Check values of cloud mass flux ! If water vapor amount transported by convection is larger than that in a ! column, cloud mass flux is reduced. ! ! tendency of specific humidity is calculated tentatively do k = 1, kmax xyz_DQVapDtCumulus(:,:,k) = + xy_CldMassFluxBottom * ( xyz_GammaMSE(:,:,k) - xyz_GammaDSE(:,:,k) ) / LatentHeat end do ! total H2O mass in a vertical column after RAS xyz_QH2OVapTentative = xyz_QH2OVap + xyz_DQVapDtCumulus * 2.0_DP * DelTime xy_CldMassFluxCorFactor = 1.0_DP do k = 1, kmax do j = 1, jmax do i = 0, imax-1 if ( xyz_QH2OVapTentative(i,j,k) < 0.0_DP ) then CldMassFluxCorFactor = xyz_QH2OVap(i,j,k) / ( xyz_QH2OVap(i,j,k) - xyz_QH2OVapTentative(i,j,k) ) else CldMassFluxCorFactor = 1.0_DP end if if ( CldMassFluxCorFactor < xy_CldMassFluxCorFactor(i,j) ) then xy_CldMassFluxCorFactor(i,j) = CldMassFluxCorFactor end if end do end do end do ! modify cloud mass flux xy_CldMassFluxBottom = xy_CldMassFluxCorFactor * xy_CldMassFluxBottom !!$ do k = 1, kmax !!$ xyz_DQVapDtCumulus(:,:,k) = & !!$ & + xy_CloudMassFluxBottom * ( xyz_GammaMSE(:,:,k) - xyz_GammaDSE(:,:,k) ) & !!$ & / LatentHeat !!$ end do !!$ ! total H2O mass in a vertical column before RAS !!$ xyz_DelH2OMass = xyz_QH2OVap * xyz_DelPress / Grav !!$ xy_H2OMassB = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_H2OMassB = xy_H2OMassB + xyz_DelH2OMass(:,:,k) !!$ end do !!$ ! total H2O mass in a vertical column after RAS !!$ xyz_QH2OVapTentative = xyz_QH2OVap + xyz_DQVapDtCumulus * 2.0_DP * DelTime !!$ xyz_DelH2OMass = xyz_QH2OVapTentative * xyz_DelPress / Grav !!$ xy_H2OMassA = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_H2OMassA = xy_H2OMassA + xyz_DelH2OMass(:,:,k) !!$ end do !!$ ! modify cloud mass flux !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( xy_H2OMassA(i,j) < 0.0_DP ) then !!$ ! A safety factor ( 1.0_DP + 1.0d-5 ) is arbitrary. !!$ xy_CloudMassFluxBottom(i,j) = xy_CloudMassFluxBottom(i,j) & !!$ & * xy_H2OMassB(i,j) & !!$ & / ( ( xy_H2OMassB(i,j) - xy_H2OMassA(i,j) ) * ( 1.0_DP + 1.0d-5 ) ) !!$ end if !!$ end do !!$ end do ! Tendencies of specific temperature and humidity ! do k = 1, kmax xyz_DTempDtCumulus(:,:,k) = + xy_CldMassFluxBottom * xyz_GammaDSE(:,:,k) / CpDry xyz_DQVapDtCumulus(:,:,k) = + xy_CldMassFluxBottom * ( xyz_GammaMSE(:,:,k) - xyz_GammaDSE(:,:,k) ) / LatentHeat end do !!$ ! !!$ ! modification of tendency of temperature and water vapor in the mixed layer !!$ ! !!$ if ( FlagUniformMixedLayer ) then !!$ xy_SumTmp = 0.0_DP !!$ do k = 1, kmax !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( k <= xy_IndexMixLayTop(i,j) ) then !!$ xy_SumTmp(i,j) = xy_SumTmp(i,j) & !!$ & + xyz_DTempDtCumulus(i,j,k) & !!$ & * ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) ) !!$ end if !!$ end do !!$ end do !!$ end do !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ xy_SumTmp(i,j) = xy_SumTmp(i,j) & !!$ & / ( xyr_Press(i,j,0) - xyr_Press(i,j,xy_IndexMixLayTop(i,j)) ) !!$ end do !!$ end do !!$ do k = 1, kmax !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( k <= xy_IndexMixLayTop(i,j) ) then !!$ xyz_DTempDtCumulus(i,j,k) = xy_SumTmp(i,j) !!$ end if !!$ end do !!$ end do !!$ end do !!$ ! !!$ xy_SumTmp = 0.0_DP !!$ do k = 1, kmax !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( k <= xy_IndexMixLayTop(i,j) ) then !!$ xy_SumTmp(i,j) = xy_SumTmp(i,j) & !!$ & + xyz_DQVapDtCumulus(i,j,k) & !!$ & * ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) ) !!$ end if !!$ end do !!$ end do !!$ end do !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ xy_SumTmp(i,j) = xy_SumTmp(i,j) & !!$ & / ( xyr_Press(i,j,0) - xyr_Press(i,j,xy_IndexMixLayTop(i,j)) ) !!$ end do !!$ end do !!$ do k = 1, kmax !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( k <= xy_IndexMixLayTop(i,j) ) then !!$ xyz_DQVapDtCumulus(i,j,k) = xy_SumTmp(i,j) !!$ end if !!$ end do !!$ end do !!$ end do !!$ end if ! add tendencies to temperature and specific humidity ! xyz_Temp = xyz_Temp + xyz_DTempDtCumulus * 2.0_DP * DelTime xyz_QH2OVap = xyz_QH2OVap + xyz_DQVapDtCumulus * 2.0_DP * DelTime ! Precipitation rate at cloud top level ! unit is kg m-2 s-1 ! xyz_RainCumulus(:,:,l) = xy_CldMassFluxBottom * xyz_RainFactor(:,:,l) * xy_NormMassFluxCldTop * xy_CldQH2OLiqCldTop ! mass fix ! xyz_DelH2OMass = xyz_QH2OVap * xyz_DelPress / Grav ! total H2O mass in a vertical column xy_H2OMassB = 0.0_DP do k = kmax, 1, -1 xy_H2OMassB = xy_H2OMassB + xyz_DelH2OMass(:,:,k) end do do j = 1, jmax do i = 0, imax-1 if ( xy_H2OMassB(i,j) < 0.0_DP ) then call MessageNotify( 'E', module_name, 'Mass of water vapor in a column is negative (%d,%d), %f.', i = (/i,j/), d = (/xy_H2OMassB(i,j)/) ) end if end do end do ! negative mass is borrowed from above do k = 1, kmax-1 do j = 1, jmax do i = 0, imax-1 if ( xyz_DelH2OMass(i,j,k) < 0.0_DP ) then xyz_DelH2OMass(i,j,k+1) = xyz_DelH2OMass(i,j,k+1) + xyz_DelH2OMass(i,j,k) xyz_DelH2OMass(i,j,k ) = 0.0_DP end if end do end do end do k = kmax do j = 1, jmax do i = 0, imax-1 if ( xyz_DelH2OMass(i,j,k) < 0.0_DP ) then !!$ call MessageNotify( 'E', module_name, & !!$ & 'Mass of water vapor in the top layer is negative (%d,%d,%d), %f.', & !!$ & i = (/i,j,k/), d = (/xyz_DelH2OMass(i,j,k)/) ) !!$ !!$ xyz_RainCumulus(i,j,l) = xyz_RainCumulus(i,j,l) & !!$ & - xyz_DelH2OMass(i,j,k) / ( 2.0_DP * DelTime ) !!$ xyz_Temp (i,j,k) = xyz_Temp(i,j,k) & !!$ & - LatentHeat * xyz_DelH2OMass(i,j,k) / ( xyz_DelPress(i,j,k) / Grav )& !!$ & / CpDry xyz_DelH2OMass (i,j,k) = 0.0_DP end if end do end do !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( xyz_RainCumulus(i,j,l) < 0.0_DP ) then !!$ call MessageNotify( 'E', module_name, & !!$ & 'Mass of water vapor is insufficient at (%d,%d,%d), %f.', & !!$ & i = (/i,j,k/), d = (/xyz_RainCumulus(i,j,l)/) ) !!$ end if !!$ end do !!$ end do ! total H2O mass in a vertical column, again xy_H2OMassA = 0.0_DP do k = kmax, 1, -1 xy_H2OMassA = xy_H2OMassA + xyz_DelH2OMass(:,:,k) end do ! total mass in a vertical column is adjusted do j = 1, jmax do i = 0, imax-1 if ( xy_H2OMassA(i,j) > 0.0_DP ) then !!$ write( 6, * ) i, j, xy_H2OMassB(i,j), xy_H2OMassB(i,j) / xy_H2OMassA(i,j) do k = 1, kmax xyz_DelH2OMass(i,j,k) = xyz_DelH2OMass(i,j,k) * xy_H2OMassB(i,j) / xy_H2OMassA(i,j) end do else do k = 1, kmax xyz_DelH2OMass(i,j,k) = 0.0_DP end do end if end do end do xyz_QH2OVap = xyz_DelH2OMass / ( xyz_DelPress / Grav ) ! Detrainment mass tendency per unit mass (kg m-3 s-1 / ( kg m-3 ) = s-1). ! This corresponds to condensation rate (kg m-2 s-1) divided by layer thickness (m) ! and density (kg m-3), in other words. ! kg m-2 s-1 / ( Pa / ( m s-2 ) ) ! = kg m-2 s-1 Pa-1 m s-1 = kg m-2 (kg m s-2 m-2)-1 m s-2 ! = kg m-2 s-1 kg-1 m-1 s2 m2 m s-2 = s-1 if ( present( xyz_MoistConvDetTend ) ) then xyz_MoistConvDetTend(:,:,l) = xy_CldMassFluxBottom * xy_NormMassFluxCldTop / ( xyz_DelPress(:,:,l) / Grav ) end if if ( present( xyz_MoistConvSubsidMassFlux ) ) then ! Subsidence mass flux between the updrafts do k = 1, l-1 do j = 1, jmax do i = 0, imax-1 if ( k > xy_IndexMixLayTop(i,j) ) then DelNormMassFluxHalfLayer = - xy_EntParam(i,j) * xyz_BetaCldTop(i,j,k) * xyz_PotTemp(i,j,k) NormMassFlux = xyr_NormMassFlux(i,j,k-1) - DelNormMassFluxHalfLayer xyz_MoistConvSubsidMassFlux(i,j,k) = xyz_MoistConvSubsidMassFlux(i,j,k) + xy_CldMassFluxBottom(i,j) * NormMassFlux end if end do end do end do end if end do loop_cloud_top ! 温度変化率, 比湿変化率 ! Calculate specific humidity tendency and temperature tendency ! (In fact, temperature tendency does not need to calculate, here.) ! xyz_DTempDtCumulus = ( xyz_Temp - xyz_TempB ) / ( 2.0_DP * DelTime ) xyz_DQVapDtCumulus = ( xyz_QH2OVap - xyz_QH2OVapB ) / ( 2.0_DP * DelTime ) xyz_DTempDt = xyz_DTempDt + xyz_DTempDtCumulus xyz_DQVapDt = xyz_DQVapDt + xyz_DQVapDtCumulus ! Precipitation rate at the surface ! unit is kg m-2 s-1 ! !!$ xy_RainCumulus = 0.0d0 !!$ do k = kmax, 1, -1 !!$ xy_RainCumulus = xy_RainCumulus + xyz_RainCumulus(:,:,k) !!$ end do xyz_DQH2OLiqDt = xyz_RainCumulus / ( xyz_DelPress / Grav ) !!$ 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 !!$ !!$ xy_Rain = xy_Rain + xy_RainCumulus ! Calculation for debug ! check of conservation of water amount and internal energy ! !!$ xyz_DelVal = xyz_QH2OVapB * xyz_DelPress / Grav !!$ xy_SumValB = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValB = xy_SumValB + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xyz_DelVal = xyz_QH2OVap * xyz_DelPress / Grav !!$ xy_SumValA = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValA = xy_SumValA + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xy_SumValA = xy_SumValA + xy_RainCumulus * 2.0_DP * DelTime !!$ ! !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ Ratio = ( xy_SumValA(i,j) - xy_SumValB(i,j) ) & !!$ & / max( xy_SumValA(i,j), 1.0d-100 ) !!$ if ( abs( Ratio ) > 1.0d-14 ) then !!$ write( 6, * ) 'H2O: ', i, j, & !!$ & xy_SumValB(i,j), xy_SumValA(i,j), & !!$ & xy_RainCumulus(i,j) * 2.0_DP * DelTime, & !!$ & Ratio !!$ end if !!$ end do !!$ end do !!$ ! !!$ ! !!$ xyz_DelVal = CpDry * xyz_TempB * xyz_DelPress / Grav !!$ xy_SumValB = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValB = xy_SumValB + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xyz_DelVal = CpDry * xyz_Temp * xyz_DelPress / Grav !!$ xy_SumValA = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValA = xy_SumValA + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xy_SumValA = xy_SumValA - LatentHeat * xy_RainCumulus * 2.0_DP * DelTime !!$ ! !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ Ratio = ( xy_SumValA(i,j) - xy_SumValB(i,j) ) & !!$ & / max( xy_SumValA(i,j), 1.0d-100 ) !!$ if ( abs( Ratio ) > 1.0d-14 ) then !!$ write( 6, * ) 'CpT: ', i, j, & !!$ & xy_SumValB(i,j), xy_SumValA(i,j), & !!$ & - LatentHeat * xy_RainCumulus(i,j) * 2.0_DP * DelTime, & !!$ & Ratio !!$ end if !!$ end do !!$ end do ! calculation for output 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 ! ヒストリデータ出力 ! History data output ! call HistoryAutoPut( TimeN, 'RainCumulus' , xy_RainCumulus * LatentHeat ) call HistoryAutoPut( TimeN, 'DTempDtCumulus' , xyz_DTempDtCumulus ) call HistoryAutoPut( TimeN, 'DQVapDtCumulus' , xyz_DQVapDtCumulus ) call HistoryAutoPut( TimeN, 'RASMassFluxDistFunc', xyz_MassFluxDistFunc ) call HistoryAutoPut( TimeN, 'RASEntParam' , xyz_EntParam ) call HistoryAutoPut( TimeN, 'RASCWF' , xyz_CWF ) call HistoryAutoPut( TimeN, 'RASCWFCrtl' , xyz_CWFCrtl ) call HistoryAutoPut( TimeN, 'RASDCWFDtLS' , xyz_DCWFDtLS ) !!$ call HistoryAutoPut( TimeN, 'RASMixLayTopIndex' , real( xy_IndexMixLayTop ) ) !!$ if ( present( xyz_DQH2OLiqDt ) ) then !!$ !!$ ! unit is kg m-2 s-1 !!$ xyz_DDelLWDtCCPLV = xyz_RainCumulus !!$ !!$ ! 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 !!$ !!$ ! unit is s-1 !!$ xyz_DQH2OLiqDt = xyz_DDelLWDtCCPLV / ( xyz_DelPress / Grav ) !!$ !!$ end if ! 計算時間計測一時停止 ! Pause measurement of computation time ! call TimesetClockStop( module_name ) end subroutine RelaxedArakawaSchubert
Subroutine : | |||||
xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||||
xyz_Press(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||||
xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in )
| ||||
xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||||
xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in )
| ||||
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_QH2OVap(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_DQVapDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||||
xyz_DQH2OLiqDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out ) | ||||
xyz_MoistConvDetTend(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out ), optional | ||||
xyz_MoistConvSubsidMassFlux(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out ), optional |
relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化.
Change temperature and specific humidity by relaxed Arakawa-Schubert scheme
subroutine RelaxedArakawaSchubert1DWrapper( xy_SurfTemp, xyz_Press, xyr_Press, xyz_Exner, xyr_Exner, xyz_Temp, xyz_QH2OVap, xyz_DTempDt, xyz_DQVapDt, xyz_DQH2OLiqDt, xyz_MoistConvDetTend, xyz_MoistConvSubsidMassFlux ) ! ! relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化. ! ! Change temperature and specific humidity by relaxed Arakawa-Schubert scheme ! ! モジュール引用 ; 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, xyz_CalcDQVapSatDTemp ! Arakawa-Schubert scheme by Lord et al. (1982) ! Arakawa-Schubert scheme by Lord et al. (1982) ! use arakawa_schubert_L1982, only : ArakawaSchubertL1982CalcCWFCrtl ! 宣言文 ; Declaration statements ! real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax) ! Pressure 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 ) :: xyr_Exner (0:imax-1, 1:jmax, 0:kmax) ! Exner function real(DP), intent(inout) :: xyz_Temp (0:imax-1, 1:jmax, 1:kmax) ! Temperature real(DP), intent(inout) :: xyz_QH2OVap (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(inout) :: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax) ! 温度変化率. ! Temperature tendency real(DP), intent(inout) :: xyz_DQVapDt (0:imax-1, 1:jmax, 1:kmax) ! 比湿変化率. ! Specific humidity tendency real(DP), intent(out ) :: xyz_DQH2OLiqDt(0:imax-1, 1:jmax, 1:kmax) real(DP), intent(out ), optional :: xyz_MoistConvDetTend (0:imax-1, 1:jmax, 1:kmax) real(DP), intent(out ), optional :: xyz_MoistConvSubsidMassFlux(0:imax-1, 1:jmax, 1:kmax) ! 作業変数 ! Work variables ! real(DP) :: SurfTemp ! Pressure real(DP) :: z_Press (1:kmax) ! Pressure real(DP) :: r_Press (0:kmax) ! Pressure real(DP) :: z_Exner (1:kmax) ! Exner function real(DP) :: r_Exner (0:kmax) ! Exner function real(DP) :: z_Temp (1:kmax) ! Temperature real(DP) :: z_QH2OVap (1:kmax) ! $ q $ . 比湿. Specific humidity !!$ real(DP), intent(inout) :: xy_Rain (0:imax-1, 1:jmax) !!$ ! 降水量. !!$ ! Precipitation real(DP) :: z_DTempDt (1:kmax) ! 温度変化率. ! Temperature tendency real(DP) :: z_DQVapDt (1:kmax) ! 比湿変化率. ! Specific humidity tendency real(DP) :: z_DQH2OLiqDt(1:kmax) real(DP) :: z_MoistConvDetTend (1:kmax) real(DP) :: z_MoistConvSubsidMassFlux(1:kmax) 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_DelPress(0:imax-1, 1:jmax, 1:kmax) ! $ \Delta p $ ! real(DP) :: xyz_PotTemp (0:imax-1, 1:jmax, 1:kmax) ! Potential temperature ! real(DP) :: xyz_QH2OVapSat (0:imax-1, 1:jmax, 1:kmax) ! 飽和比湿. ! Saturation specific humidity. ! Dry and moist static energy in environment (Env) and cloud (Cld) ! real(DP) :: xyz_EnvDryStaticEne (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyr_EnvDryStaticEne (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyz_EnvMoistStaticEne (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyr_EnvMoistStaticEne (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyz_EnvMoistStaticEneSat(0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyr_EnvMoistStaticEneSat(0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyr_CldMoistStaticEne (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xy_Kernel (0:imax-1, 1:jmax) ! Tendency of cloud work function by cumulus convection, kernel real(DP) :: xy_CWF (0:imax-1, 1:jmax) ! Cloud work function real(DP) :: xyz_CWF (0:imax-1, 1:jmax, 1:kmax) ! Cloud work function ! (variable for output) real(DP) :: xy_DCWFDtLS (0:imax-1, 1:jmax) ! Tendency of cloud work function by large scale motion real(DP) :: xyz_DCWFDtLS (0:imax-1, 1:jmax, 1:kmax) ! Tendency of cloud work function by large scale motion ! (variable for output) real(DP) :: xy_CldMassFluxBottom (0:imax-1, 1:jmax) ! Cloud mass flux at cloud bottom real(DP) :: xyz_Beta (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_BetaCldTop (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_Gamma (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_GammaDSE (0:imax-1, 1:jmax, 1:kmax) ! Tendency of dry static energy per unit mass flux real(DP) :: xyz_GammaMSE (0:imax-1, 1:jmax, 1:kmax) ! Tendency of moist static energy per unit mass flux real(DP) :: xyz_Mu (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_Eps (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_PressCldBase (0:imax-1, 1:jmax) ! Pressure of cloud base real(DP) :: xyz_CWFCrtl (0:imax-1, 1:jmax, 1:kmax) ! "Critical value" of cloud work function real(DP) :: xyz_RainFactor (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_EntParam (0:imax-1, 1:jmax) ! Entrainment factor real(DP) :: xyz_EntParam (0:imax-1, 1:jmax, 1:kmax) ! Entrainment factor (variable for output) real(DP) :: xy_EntParamLL (0:imax-1, 1:jmax) ! Entrainment factor for a cloud with top at one layer ! higher level real(DP) :: xy_EntParamUL (0:imax-1, 1:jmax) ! Entrainment factor for a cloud with top at one layer ! lower level ! Difference of normalized mass flux between layer interface real(DP) :: xyz_DelNormMassFlux (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_DelNormMassFluxCldTop(0:imax-1, 1:jmax) ! Normalized mass flux at layer interface and cloud top real(DP) :: xyr_NormMassFlux (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xy_NormMassFluxCldTop (0:imax-1, 1:jmax) ! Liquid water at cloud top real(DP) :: xy_CldQH2OLiqCldTop (0:imax-1, 1:jmax) ! Mass flux distribution function real(DP) :: xyz_MassFluxDistFunc (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_DelH2OMass (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xy_H2OMassB (0:imax-1, 1:jmax) real(DP) :: xy_H2OMassA (0:imax-1, 1:jmax) real(DP) :: xyz_RainCumulus (0:imax-1, 1:jmax, 1:kmax) real(DP) :: xyz_DDelLWDtCCPLV(0:imax-1, 1:jmax, 1:kmax) logical :: xy_FlagCrossSatEquivPotTemp(0:imax-1, 1:jmax) ! ! Flag showing whether a parcel in cloud has moist static ! energy larger than environment's real(DP) :: xyr_QH2OVapSat (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xyr_TempAdiabAscent (0:imax-1, 1:jmax, 0:kmax) real(DP) :: xy_SurfPotTemp (0:imax-1, 1:jmax) !!$ real(DP) :: xyz_TempAdiabAscent (0:imax-1, 1:jmax, 1:kmax) ! Variables for looking for top of mixed layer ! logical :: xy_FlagMixLayTopFound (0:imax-1, 1:jmax) integer :: xy_IndexMixLayTop (0:imax-1, 1:jmax) ! Variables for modification of cloud mass flux ! real(DP) :: xyz_QH2OVapTentative (0:imax-1, 1:jmax, 1:kmax) real(DP) :: CldMassFluxCorFactor real(DP) :: xy_CldMassFluxCorFactor(0:imax-1, 1:jmax) real(DP) :: xyz_QH2OVapB(0:imax-1, 1:jmax, 1:kmax) ! 調節前の比湿. ! Specific humidity before adjustment ! Flags for modification of ! logical :: xy_FlagKernelNegative (0:imax-1, 1:jmax) logical :: xy_FlagNegH2OLiqCldTop(0:imax-1, 1:jmax) ! Variables for subsidence mass flux between updrafts ! real(DP) :: DelNormMassFluxHalfLayer real(DP) :: NormMassFlux ! Variables for debug ! !!$ real(DP) :: xyz_DelVal(0:imax-1, 1:jmax, 1:kmax) !!$ real(DP) :: xy_SumValB(0:imax-1, 1:jmax) !!$ real(DP) :: xy_SumValA(0:imax-1, 1:jmax) !!$ real(DP) :: Ratio real(DP) :: xy_SumTmp(0:imax-1, 1:jmax) 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 :: l integer :: m integer :: n ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. relaxed_arakawa_schubert_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 SurfTemp = xy_SurfTemp(i,j) do k = 1, kmax z_Press (k) = xyz_Press (i,j,k) z_Exner (k) = xyz_Exner (i,j,k) z_Temp (k) = xyz_Temp (i,j,k) z_QH2OVap(k) = xyz_QH2OVap(i,j,k) z_DTempDt(k) = xyz_DTempDt(i,j,k) z_DQVapDt(k) = xyz_DQVapDt(i,j,k) end do do k = 0, kmax r_Press (k) = xyr_Press (i,j,k) r_Exner (k) = xyr_Exner (i,j,k) end do call RelaxedArakawaSchubert1D( SurfTemp, z_Press, r_Press, z_Exner, r_Exner, z_Temp, z_QH2OVap, z_DTempDt, z_DQVapDt, z_DQH2OLiqDt, z_MoistConvDetTend, z_MoistConvSubsidMassFlux ) do k = 1, kmax xyz_Temp (i,j,k) = z_Temp (k) xyz_QH2OVap (i,j,k) = z_QH2OVap (k) xyz_DTempDt (i,j,k) = z_DTempDt (k) xyz_DQVapDt (i,j,k) = z_DQVapDt (k) xyz_DQH2OLiqDt(i,j,k) = z_DQH2OLiqDt(k) end do if ( present( xyz_MoistConvDetTend ) ) then do k = 1, kmax xyz_MoistConvDetTend(i,j,k) = z_MoistConvDetTend(k) end do end if if ( present( xyz_MoistConvSubsidMassFlux ) ) then do k = 1, kmax xyz_MoistConvSubsidMassFlux(i,j,k) = z_MoistConvSubsidMassFlux(k) end do end if end do end do ! calculation for output do k = 1, kmax xyz_DelPress(:,:,k) = xyr_Press(:,:,k-1) - xyr_Press(:,:,k) end do 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 ! ヒストリデータ出力 ! History data output ! call HistoryAutoPut( TimeN, 'RainCumulus' , xy_RainCumulus * LatentHeat ) call HistoryAutoPut( TimeN, 'DTempDtCumulus' , xyz_DTempDtCumulus ) call HistoryAutoPut( TimeN, 'DQVapDtCumulus' , xyz_DQVapDtCumulus ) call HistoryAutoPut( TimeN, 'RASMassFluxDistFunc', xyz_MassFluxDistFunc ) call HistoryAutoPut( TimeN, 'RASEntParam' , xyz_EntParam ) call HistoryAutoPut( TimeN, 'RASCWF' , xyz_CWF ) call HistoryAutoPut( TimeN, 'RASCWFCrtl' , xyz_CWFCrtl ) call HistoryAutoPut( TimeN, 'RASDCWFDtLS' , xyz_DCWFDtLS ) !!$ call HistoryAutoPut( TimeN, 'RASMixLayTopIndex' , real( xy_IndexMixLayTop ) ) !!$ if ( present( xyz_DQH2OLiqDt ) ) then !!$ !!$ ! unit is kg m-2 s-1 !!$ xyz_DDelLWDtCCPLV = xyz_RainCumulus !!$ !!$ ! 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 !!$ !!$ ! unit is s-1 !!$ xyz_DQH2OLiqDt = xyz_DDelLWDtCCPLV / ( xyz_DelPress / Grav ) !!$ !!$ end if ! 計算時間計測一時停止 ! Pause measurement of computation time ! call TimesetClockStop( module_name ) end subroutine RelaxedArakawaSchubert1DWrapper
Subroutine : |
moist_conv_adjust モジュールの初期化を行います. NAMELIST#moist_conv_adjust_nml の読み込みはこの手続きで行われます.
"moist_conv_adjust" module is initialized. "NAMELIST#moist_conv_adjust_nml" is loaded in this procedure.
This procedure input/output NAMELIST#relaxed_arakawa_schubert .
subroutine RelaxedArakawaSchubertInit ! ! moist_conv_adjust モジュールの初期化を行います. ! NAMELIST#moist_conv_adjust_nml の読み込みはこの手続きで行われます. ! ! "moist_conv_adjust" module is initialized. ! "NAMELIST#moist_conv_adjust_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 ! 文字列操作 ! Character handling ! use dc_string, only: StoA ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoAddVariable ! Arakawa-Schubert scheme by Lord et al. (1982) ! Arakawa-Schubert scheme by Lord et al. (1982) ! use arakawa_schubert_L1982, only : ArakawaSchubertL1982Init ! 宣言文 ; Declaration statements ! implicit none integer:: unit_nml ! NAMELIST ファイルオープン用装置番号. ! Unit number for NAMELIST file open integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT. ! IOSTAT of NAMELIST read character(STRING) :: VarName integer:: k ! NAMELIST 変数群 ! NAMELIST group name ! namelist /relaxed_arakawa_schubert/ AdjTimeConst, RainConversionFactor, FlagZeroCrtlCWF ! デフォルト値については初期化手続 "moist_conv_adjust#CumAdjInit" ! のソースコードを参照のこと. ! ! Refer to source codes in the initialization procedure ! "moist_conv_adjust#MoistConvAdjustInit" for the default values. ! ! 実行文 ; Executable statement ! if ( relaxed_arakawa_schubert_inited ) return ! デフォルト値の設定 ! Default values settings ! !!$ FlagUse = .true. !!$ FlagUniformMixedLayer = .false. AdjTimeConst = 7200.0_DP RainConversionFactor = -1.0_DP FlagZeroCrtlCWF = .false. ! NAMELIST の読み込み ! NAMELIST is input ! if ( trim(namelist_filename) /= '' ) then call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in) rewind( unit_nml ) read( unit_nml, nml = relaxed_arakawa_schubert, iostat = iostat_nml ) ! (out) close( unit_nml ) call NmlutilMsg( iostat_nml, module_name ) ! (in) end if ! Check values ! if ( RainConversionFactor > 1.0_DP ) then call MessageNotify( 'E', module_name, 'RainConversionFactor is %f, but it must be less than or equal to 1', d = (/ RainConversionFactor /) ) end if ! ヒストリデータ出力のためのへの変数登録 ! Register of variables for history data output ! call HistoryAutoAddVariable( 'RainCumulus', (/ 'lon ', 'lat ', 'time' /), 'precipitation by cumulus scheme, RAS', 'W m-2' ) call HistoryAutoAddVariable( 'DTempDtCumulus', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'cumulus condensation heating, RAS', 'K s-1' ) call HistoryAutoAddVariable( 'DQVapDtCumulus', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'cumulus condensation moistening, RAS', 'kg kg-1 s-1' ) call HistoryAutoAddVariable( 'RASMassFluxDistFunc', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'mass flux distribution function, RAS', 'kg m-2 s-1' ) call HistoryAutoAddVariable( 'RASEntParam', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'entrainment parameter', 'm-1' ) call HistoryAutoAddVariable( 'RASCWF', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'cloud work function', 'J kg-1' ) call HistoryAutoAddVariable( 'RASCWFCrtl', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'critical cloud work function', 'J kg-1' ) call HistoryAutoAddVariable( 'RASDCWFDtLS', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'time derivative of cloud work function by large scale', 'J kg-1 s-1' ) !!$ call HistoryAutoAddVariable( 'RASMixLayTopIndex', & !!$ & (/ 'lon ', 'lat ', 'time' /), & !!$ & 'index of top of mixed layer', '1' ) do k = 1, kmax write( VarName, '(a,i3.3)' ) 'RASCldTemp', k call HistoryAutoAddVariable( Varname, (/ 'lon ', 'lat ', 'sigm', 'time' /), 'temperature of cloud air', 'K' ) write( VarName, '(a,i3.3)' ) 'RASCldQH2OVap', k call HistoryAutoAddVariable( Varname, (/ 'lon ', 'lat ', 'sigm', 'time' /), 'mixing ratio of water vapor in cloud', '1' ) write( VarName, '(a,i3.3)' ) 'RASCldQH2OLiq', k call HistoryAutoAddVariable( Varname, (/ 'lon ', 'lat ', 'sigm', 'time' /), 'mixing ratio of liquid water in cloud', '1' ) end do ! Initialization of modules used in this module ! ! Arakawa-Schubert scheme by Lord et al. (1982) ! Arakawa-Schubert scheme by Lord et al. (1982) ! call ArakawaSchubertL1982Init ! 印字 ; Print ! call MessageNotify( 'M', module_name, '----- Initialization Messages -----' ) !!$ call MessageNotify( 'M', module_name, ' FlagUse = %b', l = (/ FlagUse /) ) !!$ call MessageNotify( 'M', module_name, ' FlagUniformMixedLayer = %b', l = (/ FlagUniformMixedLayer /) ) call MessageNotify( 'M', module_name, ' AdjTimeConst = %f', d = (/ AdjTimeConst /) ) call MessageNotify( 'M', module_name, ' RainConversionFactor = %f', d = (/ RainConversionFactor /) ) call MessageNotify( 'M', module_name, ' FlagZeroCrtlCWF = %b', l = (/ FlagZeroCrtlCWF /) ) call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) ) relaxed_arakawa_schubert_inited = .true. end subroutine RelaxedArakawaSchubertInit
Variable : | |||
FlagZeroCrtlCWF : | logical , save
|
Subroutine : | |||||
SurfTemp : | real(DP), intent(in )
| ||||
z_Press(1:kmax) : | real(DP), intent(in )
| ||||
r_Press(0:kmax) : | real(DP), intent(in )
| ||||
z_Exner(1:kmax) : | real(DP), intent(in )
| ||||
r_Exner(0:kmax) : | real(DP), intent(in )
| ||||
z_Temp(1:kmax) : | real(DP), intent(inout)
| ||||
z_QH2OVap(1:kmax) : | real(DP), intent(inout)
| ||||
z_DTempDt(1:kmax) : | real(DP), intent(inout)
| ||||
z_DQVapDt(1:kmax) : | real(DP), intent(inout)
| ||||
z_DQH2OLiqDt(1:kmax) : | real(DP), intent(out ) | ||||
z_MoistConvDetTend(1:kmax) : | real(DP), intent(out ), optional | ||||
z_MoistConvSubsidMassFlux(1:kmax) : | real(DP), intent(out ), optional | ||||
rz_CldTemp(0:kmax, 1:kmax) : | real(DP), intent(out ), optional | ||||
rz_CldQH2OVap(0:kmax, 1:kmax) : | real(DP), intent(out ), optional | ||||
rz_CldQH2OLiq(0:kmax, 1:kmax) : | real(DP), intent(out ), optional | ||||
rz_CldQH2OSol(0:kmax, 1:kmax) : | real(DP), intent(out ), optional |
relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化.
Change temperature and specific humidity by relaxed Arakawa-Schubert scheme
subroutine RAS1DTesting( SurfTemp, z_Press, r_Press, z_Exner, r_Exner, z_Temp, z_QH2OVap, z_DTempDt, z_DQVapDt, z_DQH2OLiqDt, z_MoistConvDetTend, z_MoistConvSubsidMassFlux, rz_CldTemp, rz_CldQH2OVap, rz_CldQH2OLiq, rz_CldQH2OSol ) ! ! relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化. ! ! Change temperature and specific humidity by relaxed Arakawa-Schubert scheme ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: Grav, GasRDry, CpDry, LatentHeat, LatentHeatFusion ! $ L $ [J kg-1] . ! 融解の潜熱. ! Latent heat of fusion ! 時刻管理 ! Time control ! use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! 飽和比湿の算出 ! Evaluate saturation specific humidity ! use saturate, only: a_CalcQVapSat, a_CalcDQVapSatDTemp ! Arakawa-Schubert scheme by Lord et al. (1982) ! Arakawa-Schubert scheme by Lord et al. (1982) ! use arakawa_schubert_L1982, only : ASL1982CalcCWFCrtl1D ! 雲関系ルーチン ! Cloud-related routines ! use cloud_utils, only : CloudUtilsWatFraction ! 宣言文 ; Declaration statements ! real(DP), intent(in ) :: SurfTemp ! Pressure real(DP), intent(in ) :: z_Press (1:kmax) ! Pressure real(DP), intent(in ) :: r_Press (0:kmax) ! Pressure real(DP), intent(in ) :: z_Exner (1:kmax) ! Exner function real(DP), intent(in ) :: r_Exner (0:kmax) ! Exner function real(DP), intent(inout) :: z_Temp (1:kmax) ! Temperature real(DP), intent(inout) :: z_QH2OVap (1:kmax) ! $ q $ . 比湿. Specific humidity !!$ real(DP), intent(inout) :: xy_Rain (0:imax-1, 1:jmax) !!$ ! 降水量. !!$ ! Precipitation real(DP), intent(inout) :: z_DTempDt (1:kmax) ! 温度変化率. ! Temperature tendency real(DP), intent(inout) :: z_DQVapDt (1:kmax) ! 比湿変化率. ! Specific humidity tendency real(DP), intent(out ) :: z_DQH2OLiqDt(1:kmax) real(DP), intent(out ), optional :: z_MoistConvDetTend (1:kmax) real(DP), intent(out ), optional :: z_MoistConvSubsidMassFlux(1:kmax) real(DP), intent(out ), optional :: rz_CldTemp (0:kmax, 1:kmax) real(DP), intent(out ), optional :: rz_CldQH2OVap(0:kmax, 1:kmax) real(DP), intent(out ), optional :: rz_CldQH2OLiq(0:kmax, 1:kmax) real(DP), intent(out ), optional :: rz_CldQH2OSol(0:kmax, 1:kmax) ! 作業変数 ! Work variables ! real(DP) :: z_Height (1:kmax) ! ! Height real(DP) :: r_Height (0:kmax) ! ! Height real(DP) :: RainCumulus ! 降水量. ! Precipitation real(DP) :: z_DTempDtCumulus (1:kmax) ! 温度変化率. ! Temperature tendency real(DP) :: z_DQVapDtCumulus (1:kmax) ! 比湿変化率. ! Specific humidity tendency real(DP) :: z_DelPress(1:kmax) ! $ \Delta p $ ! real(DP) :: z_PotTemp (1:kmax) ! Potential temperature ! real(DP) :: z_QH2OVapSat(1:kmax) ! 飽和比湿. ! Saturation specific humidity. ! Dry and moist static energy in environment (Env) and cloud (Cld) ! real(DP) :: z_EnvDryStaticEne (1:kmax) real(DP) :: r_EnvDryStaticEne (0:kmax) real(DP) :: z_EnvMoistStaticEne (1:kmax) real(DP) :: r_EnvMoistStaticEne (0:kmax) real(DP) :: z_EnvMoistStaticEneSat(1:kmax) real(DP) :: r_EnvMoistStaticEneSat(0:kmax) real(DP) :: z_EnvCondStaticEne (1:kmax) real(DP) :: r_CldMoistStaticEne (0:kmax) real(DP) :: r_CldCondStaticEne (0:kmax) real(DP) :: CldCondStaticEneCldTop real(DP) :: Kernel ! Tendency of cloud work function by cumulus convection, kernel real(DP) :: CWF ! Cloud work function real(DP) :: z_CWF(1:kmax) ! Cloud work function ! (variable for output) real(DP) :: DCWFDtLS ! Tendency of cloud work function by large scale motion real(DP) :: z_DCWFDtLS(1:kmax) ! Tendency of cloud work function by large scale motion ! (variable for output) real(DP) :: CldMassFluxBottom ! Cloud mass flux at cloud bottom real(DP) :: z_Beta (1:kmax) real(DP) :: z_BetaCldTop (1:kmax) real(DP) :: z_Gamma (1:kmax) real(DP) :: z_GammaDSE (1:kmax) ! Tendency of dry static energy per unit mass flux real(DP) :: z_GammaMSE (1:kmax) ! Tendency of moist static energy per unit mass flux real(DP) :: z_Mu (1:kmax) real(DP) :: z_Eps (1:kmax) real(DP) :: PressCldBase ! Pressure of cloud base real(DP) :: z_CWFCrtl (1:kmax) ! "Critical value" of cloud work function real(DP) :: z_RainFactor (1:kmax) real(DP) :: EntParam ! Entrainment factor real(DP) :: z_EntParam (1:kmax) ! Entrainment factor (variable for output) real(DP) :: EntParamLL ! Entrainment factor for a cloud with top at one layer ! higher level real(DP) :: EntParamUL ! Entrainment factor for a cloud with top at one layer ! lower level ! Difference of normalized mass flux between layer interface real(DP) :: z_DelNormMassFlux (1:kmax) real(DP) :: DelNormMassFluxCldTop ! Normalized mass flux at layer interface and cloud top real(DP) :: r_NormMassFlux (0:kmax) real(DP) :: NormMassFluxCldTop ! cloud total water real(DP) :: r_CldQH2OTot(0:kmax) ! cloud total water at cloud top real(DP) :: CldQH2OTotCldTop ! cloud condensate at cloud top real(DP) :: CldQH2OCondCldTop ! cloud water at cloud top real(DP) :: CldQH2OLiqCldTop ! cloud ice at cloud top real(DP) :: CldQH2OSolCldTop ! cloud ice at cloud top for save real(DP) :: CldQH2OSolCldTopB real(DP) :: WatFrac ! Mass flux distribution function real(DP) :: z_MassFluxDistFunc (1:kmax) real(DP) :: z_DelH2OMass (1:kmax) real(DP) :: H2OMassB real(DP) :: H2OMassA real(DP) :: z_RainCumulus (1:kmax) real(DP) :: NegDDelLWDt real(DP) :: z_DDelLWDtCCPLV(1:kmax) logical :: FlagCrossSatEquivPotTemp ! ! Flag showing whether a parcel in cloud has moist static ! energy larger than environment's real(DP) :: r_QH2OVapSat (0:kmax) real(DP) :: r_TempAdiabAscent (0:kmax) real(DP) :: SurfPotTemp !!$ real(DP) :: xyz_TempAdiabAscent (0:imax-1, 1:jmax, 1:kmax) ! Variables for looking for top of mixed layer ! logical :: FlagMixLayTopFound integer :: IndexMixLayTop ! Variables for modification of cloud mass flux ! real(DP) :: z_QH2OVapTentative (1:kmax) real(DP) :: CldMassFluxCorFactor real(DP) :: CldMassFluxCorFactorTentative real(DP) :: z_TempB (1:kmax) ! 調節前の温度. ! Temperature before adjustment real(DP) :: z_QH2OVapB(1:kmax) ! 調節前の比湿. ! Specific humidity before adjustment ! Flags for modification of ! logical :: FlagKernelNegative logical :: FlagNegH2OCondCldTop ! Variables for subsidence mass flux between updrafts ! real(DP) :: DelNormMassFluxHalfLayer real(DP) :: NormMassFlux ! Variables for debug ! !!$ real(DP) :: xyz_DelVal(0:imax-1, 1:jmax, 1:kmax) !!$ real(DP) :: xy_SumValB(0:imax-1, 1:jmax) !!$ real(DP) :: xy_SumValA(0:imax-1, 1:jmax) !!$ real(DP) :: Ratio real(DP) :: CldTempB real(DP) :: a_DQVapSatDTemp(1:1) real(DP) :: DelTemp real(DP) :: r_CldTemp (0:kmax) real(DP) :: r_CldQH2OVap(0:kmax) real(DP) :: r_CldQH2OLiq(0:kmax) real(DP) :: r_CldQH2OSol(0:kmax) real(DP) :: r_CldHeight (0:kmax) real(DP) :: SumTmp integer :: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer :: l integer :: m integer :: n ! Temporal real(DP) :: z_QH2OLiq(1:kmax) real(DP) :: z_QH2OSol(1:kmax) ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. relaxed_arakawa_schubert_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 計算時間計測開始 ! Start measurement of computation time ! !!$ call TimesetClockStart( module_name ) ! Temporal z_QH2OLiq = 0.0_DP z_QH2OSol = 0.0_DP ! 調節前 "Temp", "QH2OVap" の保存 ! Store "Temp", "QH2OVap" before adjustment ! z_TempB = z_Temp z_QH2OVapB = z_QH2OVap ! Preparation of variables ! ! ! Auxiliary variables ! Pressure difference between upper and lower interface of the layer do k = 1, kmax z_DelPress(k) = r_Press(k-1) - r_Press(k) end do ! beta do k = 1, kmax z_Beta(k) = CpDry / Grav * ( r_Exner(k-1) - r_Exner(k) ) end do do k = 1, kmax z_BetaCldTop(k) = CpDry / Grav * ( r_Exner(k-1) - z_Exner(k) ) end do ! ! Search for top of mixed layer (lifting condensation level) based on ! a description in p.684 of Arakawa and Shubert (1974). ! call RelaxedArakawaSchubertHeight1D( z_Temp, z_Exner, z_Beta, z_BetaCldTop, z_Height, r_Height ) ! !==================================== ! !!$ xyz_TempAdiabAscent(:,:,1) = xyz_Temp(:,:,1) !!$ do k = 2, kmax !!$ xyz_TempAdiabAscent(:,:,k) = & !!$ & xyz_Temp(:,:,1) - Grav / CpDry * ( xyz_Height(:,:,k) - xyz_Height(:,:,1) ) !!$ end do !!$ xyz_TempAdiabAscent = max( xyz_TempAdiabAscent, 1.0_DP ) !!$ xyz_QH2OVapSat = xyz_CalcQVapSat( xyz_TempAdiabAscent, xyz_Press ) !!$ xy_IndexMixLayTop = 1 !!$ xy_FlagMixLayTopFound = .false. !!$ do k = 2, kmax !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( ( xyz_QH2OVap(i,j,1) >= xyz_QH2OVapSat(i,j,k) ) .and. & !!$ & ( .not. xy_FlagMixLayTopFound(i,j) ) ) then !!$ xy_IndexMixLayTop (i,j) = k - 1 !!$ xy_FlagMixLayTopFound(i,j) = .true. !!$ end if !!$ end do !!$ end do !!$ end do ! !------------------------------------ ! !!$ xyr_TempAdiabAscent(:,:,0) = xy_SurfTemp !!$ do k = 1, kmax !!$ xyr_TempAdiabAscent(:,:,k) = & !!$ & xy_SurfTemp - Grav / CpDry * ( xyr_Height(:,:,k) - 0.0_DP ) !!$ end do !!$ xyr_TempAdiabAscent = max( xyr_TempAdiabAscent, 1.0_DP ) !!$ r_TempAdiabAscent(0) = SurfTemp SurfPotTemp = SurfTemp / r_Exner(0) do k = 1, kmax r_TempAdiabAscent(k) = SurfPotTemp * r_Exner(k) end do ! r_QH2OVapSat(0 ) = 1.0d100 r_QH2OVapSat(1:kmax-1) = a_CalcQVapSat( r_TempAdiabAscent(1:kmax-1), r_Press(1:kmax-1) ) r_QH2OVapSat(kmax ) = r_QH2OVapSat(kmax-1) ! IndexMixLayTop = 1 FlagMixLayTopFound = .false. do k = 2, kmax if ( ( z_QH2OVap(1) >= r_QH2OVapSat(k) ) .and. ( .not. FlagMixLayTopFound ) ) then IndexMixLayTop = k - 1 FlagMixLayTopFound = .true. end if end do ! !==================================== ! if ( .not. FlagMixLayTopFound ) then IndexMixLayTop = kmax - 1 end if ! ! Critical cloud work function ! if ( FlagZeroCrtlCWF ) then z_CWFCrtl = 0.0_DP else PressCldBase = r_Press(IndexMixLayTop) call ASL1982CalcCWFCrtl1D( PressCldBase, z_Press, z_CWFCrtl ) end if ! ! Rain conversion factor ! if ( RainConversionFactor < 0.0_DP ) then do k = 1, kmax if ( z_Press(k) < 500.0d2 ) then z_RainFactor(k) = 1.0_DP else if ( z_Press(k) < 800.0d2 ) then z_RainFactor(k) = 0.8_DP + ( 800.0d2 - z_Press(k) ) / 1500.0d2 else z_RainFactor(k) = 0.8_DP end if end do else z_RainFactor = RainConversionFactor end if z_RainCumulus (1) = 0.0_DP z_EntParam (1) = 0.0_DP z_CWF (1) = 0.0_DP z_DCWFDtLS (1) = 0.0_DP z_MassFluxDistFunc(1) = 0.0_DP if ( present( z_MoistConvDetTend ) ) then z_MoistConvDetTend(1) = 0.0_DP end if if ( present( z_MoistConvSubsidMassFlux ) ) then ! Subsidence mass flux between the updrafts ! Initialization ! z_MoistConvSubsidMassFlux = 0.0_DP end if r_CldTemp = 1.0d100 r_CldQH2OVap = 1.0d100 r_CldQH2OLiq = 1.0d100 r_CldQH2OSol = 1.0d100 l = 1 if ( present( rz_CldTemp ) ) rz_CldTemp (:,l) = r_CldTemp if ( present( rz_CldQH2OVap ) ) rz_CldQH2OVap(:,l) = r_CldQH2OVap if ( present( rz_CldQH2OLiq ) ) rz_CldQH2OLiq(:,l) = r_CldQH2OLiq if ( present( rz_CldQH2OSol ) ) rz_CldQH2OSol(:,l) = r_CldQH2OSol loop_cloud_top : do l = 2, kmax call RAS1DTestingCore01( l, z_Press, r_Press, z_Exner, r_Exner, z_Temp, z_QH2OVap, z_QH2OLiq, z_QH2OSol, IndexMixLayTop, z_DelPress, z_Beta, z_BetaCldTop, z_RainFactor, z_PotTemp, z_DelNormMassFlux, DelNormMassFluxCldTop, r_NormMassFlux, NormMassFluxCldTop, CldQH2OLiqCldTop, CWF, EntParam, z_Mu, z_Eps, z_Gamma, z_GammaDSE, z_GammaMSE, FlagNegH2OCondCldTop, rz_CldTemp, rz_CldQH2OVap, rz_CldQH2OLiq, rz_CldQH2OSol ) ! Time derivative of cloud work function by large scale motion ! DCWFDtLS = ( CWF - z_CWFCrtl(l) ) / ( 2.0_DP * DelTime ) ! for save z_CWF(l) = CWF ! for output z_EntParam(l) = EntParam ! for save z_DCWFDtLS(l) = DCWFDtLS ! Kernel, time derivative of cloud work function by cumulus convection per unit ! mass flux ! Kernel = z_Eps(IndexMixLayTop+1) * z_GammaMSE(IndexMixLayTop) - z_Eps(l) * r_NormMassFlux(l-1) * ( 1.0_DP + z_Gamma(l) ) * z_GammaDSE(l) do n = IndexMixLayTop+1, l-1 SumTmp = 0.0_DP do m = IndexMixLayTop+1, n SumTmp = SumTmp + z_DelNormMassFlux(m) * z_GammaMSE(m) end do Kernel = Kernel + ( z_Eps(n+1) + z_Mu(n) ) * ( z_GammaMSE(IndexMixLayTop) - SumTmp ) - ( z_Eps(n) * r_NormMassFlux(n-1) + z_Mu (n) * r_NormMassFlux(n ) ) * ( 1.0_DP + z_Gamma(n) ) * z_GammaDSE(n) end do ! Check whether kernel is positive or negative. ! if ( Kernel < 0.0_DP ) then FlagKernelNegative = .true. else FlagKernelNegative = .false. end if ! Load et al. (1982), p.108 Kernel = min( Kernel, -5.0d-3 ) ! Cloud mass flux at cloud bottom ! CldMassFluxBottom = - DCWFDtLS / Kernel ! ! mass flux has to be zero or positive CldMassFluxBottom = max( CldMassFluxBottom, 0.0_DP ) ! mass flux is zero if entrainment parameter is zero or negative if ( EntParam <= 0.0_DP ) then CldMassFluxBottom = 0.0_DP end if !!$ ! mass flux is zero if it is below lifting condensation level !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( .not. xy_FlagCrossSatEquivPotTemp(i,j) ) then !!$ xy_CloudMassFluxBottom(i,j) = 0.0_DP !!$ end if !!$ end do !!$ end do ! mass flux is zero if the LNB is unstable for updrafts ! (i.e., if the parcel is positively buoyant just above the LNB). ! See Lord et al. (1982), p.112, for more details. ! Strictly speaking, the process below is different from that ! proposed by Lord et al. (1982). Lord et al. (1982) compare ! entrainment parameters at 3 levels. But, entrainment ! parameters at 2 levels are compared below, because comparison ! of values between 2 levels seems to be sufficient. !!!$ if ( ( 3 <= l ) .and. ( l <= kmax-1 ) ) then !!!$ do j = 1, jmax !!!$ do i = 0, imax-1 !!!$ if ( ( xy_EntParamLL(i,j) < xy_EntParam (i,j) ) .and. & !!!$ & ( xy_EntParam (i,j) < xy_EntParamUL(i,j) ) ) then !!!$ if ( ( xy_EntParamLL(i,j) > 0.0_DP ) .and. & !!!$ & ( xy_EntParam (i,j) > 0.0_DP ) .and. & !!!$ & ( xy_EntParamUL(i,j) > 0.0_DP ) ) then !!!$ xy_CloudMassFluxBottom(i,j) = 0.0_DP !!!$ end if !!!$ end if !!!$ end do !!!$ end do !!!$ end if !!!$ if ( xy_IndexMixLayTop(i,j) == l ) then !!!$ if ( ( xy_EntParam (i,j) > 0.0_DP ) .and. & !!!$ & ( xy_EntParamUL(i,j) > 0.0_DP ) ) then !!!$ if ( xy_EntParam (i,j) < xy_EntParamUL(i,j) ) then !!!$ xy_CloudMassFluxBottom(i,j) = 0.0_DP !!!$ end if !!!$ end if !!!$ else if ( ( xy_IndexMixLayTop(i,j) < l ) .and. ( l <= kmax-1 ) ) then !!!$ if ( ( xy_EntParamLL(i,j) > 0.0_DP ) .and. & !!!$ & ( xy_EntParam (i,j) > 0.0_DP ) .and. & !!!$ & ( xy_EntParamUL(i,j) > 0.0_DP ) ) then !!!$ if ( ( xy_EntParamLL(i,j) < xy_EntParam (i,j) ) .and. & !!!$ & ( xy_EntParam (i,j) < xy_EntParamUL(i,j) ) ) then ! ! This was used in a version without ice. ! But, now, lines below are commented out, because EntParamUL is not ! set. (2014/02/02) ! !!$ if ( ( IndexMixLayTop <= l ) .and. ( l <= kmax-1 ) ) then !!$ if ( ( EntParam > 0.0_DP ) .and. & !!$ & ( EntParamUL > 0.0_DP ) ) then !!$ if ( EntParam < EntParamUL ) then !!$ CldMassFluxBottom = 0.0_DP !!$ end if !!$ end if !!$ end if ! ! mass flux is zero unless kernel is negative ! if ( .not. FlagKernelNegative ) then CldMassFluxBottom = 0.0_DP end if ! ! mass flux is zero if liquid water at a cloud top is negative ! if ( FlagNegH2OCondCldTop ) then CldMassFluxBottom = 0.0_DP end if ! ! multiply factor ! CldMassFluxBottom = CldMassFluxBottom * min( 2.0_DP * DelTime / AdjTimeConst, 1.0_DP ) ! ! for output z_MassFluxDistFunc(l) = CldMassFluxBottom ! Check values of cloud mass flux ! If water vapor amount transported by convection is larger than that in a ! column, cloud mass flux is reduced. ! ! tendency of specific humidity is calculated tentatively do k = 1, kmax z_DQVapDtCumulus(k) = + CldMassFluxBottom * ( z_GammaMSE(k) - z_GammaDSE(k) ) / LatentHeat end do ! total H2O mass in a vertical column after RAS z_QH2OVapTentative = z_QH2OVap + z_DQVapDtCumulus * 2.0_DP * DelTime CldMassFluxCorFactor = 1.0_DP do k = 1, kmax if ( z_QH2OVapTentative(k) < 0.0_DP ) then CldMassFluxCorFactorTentative = z_QH2OVap(k) / ( z_QH2OVap(k) - z_QH2OVapTentative(k) ) else CldMassFluxCorFactorTentative = 1.0_DP end if if ( CldMassFluxCorFactorTentative < CldMassFluxCorFactor ) then CldMassFluxCorFactor = CldMassFluxCorFactorTentative end if end do ! modify cloud mass flux CldMassFluxBottom = CldMassFluxCorFactor * CldMassFluxBottom !!$ do k = 1, kmax !!$ xyz_DQVapDtCumulus(:,:,k) = & !!$ & + xy_CloudMassFluxBottom * ( xyz_GammaMSE(:,:,k) - xyz_GammaDSE(:,:,k) ) & !!$ & / LatentHeat !!$ end do !!$ ! total H2O mass in a vertical column before RAS !!$ xyz_DelH2OMass = xyz_QH2OVap * xyz_DelPress / Grav !!$ xy_H2OMassB = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_H2OMassB = xy_H2OMassB + xyz_DelH2OMass(:,:,k) !!$ end do !!$ ! total H2O mass in a vertical column after RAS !!$ xyz_QH2OVapTentative = xyz_QH2OVap + xyz_DQVapDtCumulus * 2.0_DP * DelTime !!$ xyz_DelH2OMass = xyz_QH2OVapTentative * xyz_DelPress / Grav !!$ xy_H2OMassA = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_H2OMassA = xy_H2OMassA + xyz_DelH2OMass(:,:,k) !!$ end do !!$ ! modify cloud mass flux !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( xy_H2OMassA(i,j) < 0.0_DP ) then !!$ ! A safety factor ( 1.0_DP + 1.0d-5 ) is arbitrary. !!$ xy_CloudMassFluxBottom(i,j) = xy_CloudMassFluxBottom(i,j) & !!$ & * xy_H2OMassB(i,j) & !!$ & / ( ( xy_H2OMassB(i,j) - xy_H2OMassA(i,j) ) * ( 1.0_DP + 1.0d-5 ) ) !!$ end if !!$ end do !!$ end do call RAS1DTestingCore02( l, z_DelPress, z_GammaDSE, z_GammaMSE, CldMassFluxBottom, z_RainFactor, NormMassFluxCldTop, CldQH2OLiqCldTop, z_Temp, z_QH2OVap, z_DTempDtCumulus, z_DQVapDtCumulus, RainCumulus ) z_RainCumulus(l) = RainCumulus ! Detrainment mass tendency per unit mass (kg m-3 s-1 / ( kg m-3 ) = s-1). ! This corresponds to condensation rate (kg m-2 s-1) divided by layer thickness (m) ! and density (kg m-3), in other words. ! kg m-2 s-1 / ( Pa / ( m s-2 ) ) ! = kg m-2 s-1 Pa-1 m s-1 = kg m-2 (kg m s-2 m-2)-1 m s-2 ! = kg m-2 s-1 kg-1 m-1 s2 m2 m s-2 = s-1 if ( present( z_MoistConvDetTend ) ) then z_MoistConvDetTend(l) = CldMassFluxBottom * NormMassFluxCldTop / ( z_DelPress(l) / Grav ) end if if ( present( z_MoistConvSubsidMassFlux ) ) then ! Subsidence mass flux between the updrafts do k = 1, l-1 if ( k > IndexMixLayTop ) then DelNormMassFluxHalfLayer = - EntParam * z_BetaCldTop(k) * z_PotTemp(k) NormMassFlux = r_NormMassFlux(k-1) - DelNormMassFluxHalfLayer z_MoistConvSubsidMassFlux(k) = z_MoistConvSubsidMassFlux(k) + CldMassFluxBottom * NormMassFlux end if end do end if end do loop_cloud_top ! 温度変化率, 比湿変化率 ! Calculate specific humidity tendency and temperature tendency ! (In fact, temperature tendency does not need to calculate, here.) ! z_DTempDtCumulus = ( z_Temp - z_TempB ) / ( 2.0_DP * DelTime ) z_DQVapDtCumulus = ( z_QH2OVap - z_QH2OVapB ) / ( 2.0_DP * DelTime ) z_DTempDt = z_DTempDt + z_DTempDtCumulus z_DQVapDt = z_DQVapDt + z_DQVapDtCumulus ! Precipitation rate at the surface ! unit is kg m-2 s-1 ! !!$ xy_RainCumulus = 0.0d0 !!$ do k = kmax, 1, -1 !!$ xy_RainCumulus = xy_RainCumulus + xyz_RainCumulus(:,:,k) !!$ end do z_DQH2OLiqDt = z_RainCumulus / ( z_DelPress / Grav ) !!$ 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 !!$ !!$ xy_Rain = xy_Rain + xy_RainCumulus ! Calculation for debug ! check of conservation of water amount and internal energy ! !!$ xyz_DelVal = xyz_QH2OVapB * xyz_DelPress / Grav !!$ xy_SumValB = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValB = xy_SumValB + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xyz_DelVal = xyz_QH2OVap * xyz_DelPress / Grav !!$ xy_SumValA = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValA = xy_SumValA + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xy_SumValA = xy_SumValA + xy_RainCumulus * 2.0_DP * DelTime !!$ ! !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ Ratio = ( xy_SumValA(i,j) - xy_SumValB(i,j) ) & !!$ & / max( xy_SumValA(i,j), 1.0d-100 ) !!$ if ( abs( Ratio ) > 1.0d-14 ) then !!$ write( 6, * ) 'H2O: ', i, j, & !!$ & xy_SumValB(i,j), xy_SumValA(i,j), & !!$ & xy_RainCumulus(i,j) * 2.0_DP * DelTime, & !!$ & Ratio !!$ end if !!$ end do !!$ end do !!$ ! !!$ ! !!$ xyz_DelVal = CpDry * xyz_TempB * xyz_DelPress / Grav !!$ xy_SumValB = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValB = xy_SumValB + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xyz_DelVal = CpDry * xyz_Temp * xyz_DelPress / Grav !!$ xy_SumValA = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValA = xy_SumValA + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xy_SumValA = xy_SumValA - LatentHeat * xy_RainCumulus * 2.0_DP * DelTime !!$ ! !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ Ratio = ( xy_SumValA(i,j) - xy_SumValB(i,j) ) & !!$ & / max( xy_SumValA(i,j), 1.0d-100 ) !!$ if ( abs( Ratio ) > 1.0d-14 ) then !!$ write( 6, * ) 'CpT: ', i, j, & !!$ & xy_SumValB(i,j), xy_SumValA(i,j), & !!$ & - LatentHeat * xy_RainCumulus(i,j) * 2.0_DP * DelTime, & !!$ & Ratio !!$ end if !!$ end do !!$ end do ! calculation for output ! This calculation is meaningless because RainCumulus is not used below. z_RainCumulus = z_DQH2OLiqDt * ( z_DelPress / Grav ) RainCumulus = 0.0d0 do k = kmax, 1, -1 RainCumulus = RainCumulus + z_RainCumulus(k) end do !!$ if ( present( xyz_DQH2OLiqDt ) ) then !!$ !!$ ! unit is kg m-2 s-1 !!$ xyz_DDelLWDtCCPLV = xyz_RainCumulus !!$ !!$ ! 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 !!$ !!$ ! unit is s-1 !!$ xyz_DQH2OLiqDt = xyz_DDelLWDtCCPLV / ( xyz_DelPress / Grav ) !!$ !!$ end if ! 計算時間計測一時停止 ! Pause measurement of computation time ! !!$ call TimesetClockStop( module_name ) end subroutine RAS1DTesting
Subroutine : | |||
l : | integer , intent(in ) | ||
z_Press(1:kmax) : | real(DP), intent(in )
| ||
r_Press(0:kmax) : | real(DP), intent(in )
| ||
z_Exner(1:kmax) : | real(DP), intent(in )
| ||
r_Exner(0:kmax) : | real(DP), intent(in )
| ||
z_Temp(1:kmax) : | real(DP), intent(in )
| ||
z_QH2OVap(1:kmax) : | real(DP), intent(in )
| ||
z_QH2OLiq(1:kmax) : | real(DP), intent(in ) | ||
z_QH2OSol(1:kmax) : | real(DP), intent(in ) | ||
IndexMixLayTop : | integer , intent(in ) | ||
z_DelPress(1:kmax) : | real(DP), intent(in )
| ||
z_Beta(1:kmax) : | real(DP), intent(in ) | ||
z_BetaCldTop(1:kmax) : | real(DP), intent(in ) | ||
z_RainFactor(1:kmax) : | real(DP), intent(in ) | ||
z_PotTemp(1:kmax) : | real(DP), intent(out )
| ||
z_DelNormMassFlux(1:kmax) : | real(DP), intent(out ) | ||
DelNormMassFluxCldTop : | real(DP), intent(out )
| ||
r_NormMassFlux(0:kmax) : | real(DP), intent(out ) | ||
NormMassFluxCldTop : | real(DP), intent(out ) | ||
CldQH2OLiqCldTop : | real(DP), intent(out ) | ||
CWF : | real(DP), intent(out )
| ||
EntParam : | real(DP), intent(out )
| ||
z_Mu(1:kmax) : | real(DP), intent(out ) | ||
z_Eps(1:kmax) : | real(DP), intent(out ) | ||
z_Gamma(1:kmax) : | real(DP), intent(out ) | ||
z_GammaDSE(1:kmax) : | real(DP), intent(out )
| ||
z_GammaMSE(1:kmax) : | real(DP), intent(out )
| ||
FlagNegH2OCondCldTop : | logical , intent(out )
| ||
rz_CldTemp(0:kmax, 1:kmax) : | real(DP), intent(inout), optional | ||
rz_CldQH2OVap(0:kmax, 1:kmax) : | real(DP), intent(inout), optional | ||
rz_CldQH2OLiq(0:kmax, 1:kmax) : | real(DP), intent(inout), optional | ||
rz_CldQH2OSol(0:kmax, 1:kmax) : | real(DP), intent(inout), optional |
relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化.
Change temperature and specific humidity by relaxed Arakawa-Schubert scheme
subroutine RAS1DTestingCore01( l, z_Press, r_Press, z_Exner, r_Exner, z_Temp, z_QH2OVap, z_QH2OLiq, z_QH2OSol, IndexMixLayTop, z_DelPress, z_Beta, z_BetaCldTop, z_RainFactor, z_PotTemp, z_DelNormMassFlux, DelNormMassFluxCldTop, r_NormMassFlux, NormMassFluxCldTop, CldQH2OLiqCldTop, CWF, EntParam, z_Mu, z_Eps, z_Gamma, z_GammaDSE, z_GammaMSE, FlagNegH2OCondCldTop, rz_CldTemp, rz_CldQH2OVap, rz_CldQH2OLiq, rz_CldQH2OSol ) ! ! relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化. ! ! Change temperature and specific humidity by relaxed Arakawa-Schubert scheme ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: Grav, GasRDry, CpDry, LatentHeat, LatentHeatFusion ! $ L $ [J kg-1] . ! 融解の潜熱. ! Latent heat of fusion ! 時刻管理 ! Time control ! use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! 飽和比湿の算出 ! Evaluate saturation specific humidity ! use saturate, only: a_CalcQVapSat, a_CalcDQVapSatDTemp ! Arakawa-Schubert scheme by Lord et al. (1982) ! Arakawa-Schubert scheme by Lord et al. (1982) ! use arakawa_schubert_L1982, only : ASL1982CalcCWFCrtl1D ! 雲関系ルーチン ! Cloud-related routines ! use cloud_utils, only : CloudUtilsWatFraction ! 宣言文 ; Declaration statements ! integer , intent(in ) :: l real(DP), intent(in ) :: z_Press (1:kmax) ! Pressure real(DP), intent(in ) :: r_Press (0:kmax) ! Pressure real(DP), intent(in ) :: z_Exner (1:kmax) ! Exner function real(DP), intent(in ) :: r_Exner (0:kmax) ! Exner function real(DP), intent(in ) :: z_Temp (1:kmax) ! Temperature real(DP), intent(in ) :: z_QH2OVap (1:kmax) ! $ q $ . 比湿. Specific humidity real(DP), intent(in ) :: z_QH2OLiq(1:kmax) real(DP), intent(in ) :: z_QH2OSol(1:kmax) integer , intent(in ) :: IndexMixLayTop real(DP), intent(in ) :: z_DelPress(1:kmax) ! $ \Delta p $ ! real(DP), intent(in ) :: z_Beta (1:kmax) real(DP), intent(in ) :: z_BetaCldTop (1:kmax) real(DP), intent(in ) :: z_RainFactor (1:kmax) !!$ real(DP), intent(inout) :: xy_Rain (0:imax-1, 1:jmax) !!$ ! 降水量. !!$ ! Precipitation real(DP), intent(out ) :: z_PotTemp (1:kmax) ! Potential temperature ! ! Difference of normalized mass flux between layer interface real(DP), intent(out ) :: z_DelNormMassFlux (1:kmax) real(DP), intent(out ) :: DelNormMassFluxCldTop ! Normalized mass flux at layer interface and cloud top real(DP), intent(out ) :: r_NormMassFlux (0:kmax) real(DP), intent(out ) :: NormMassFluxCldTop ! cloud water at cloud top real(DP), intent(out ) :: CldQH2OLiqCldTop real(DP), intent(out ) :: CWF ! Cloud work function real(DP), intent(out ) :: EntParam ! Entrainment factor real(DP), intent(out ) :: z_Mu (1:kmax) real(DP), intent(out ) :: z_Eps (1:kmax) real(DP), intent(out ) :: z_Gamma (1:kmax) real(DP), intent(out ) :: z_GammaDSE (1:kmax) ! Tendency of dry static energy per unit mass flux real(DP), intent(out ) :: z_GammaMSE (1:kmax) ! Tendency of moist static energy per unit mass flux logical , intent(out ) :: FlagNegH2OCondCldTop ! Flags for modification of real(DP), intent(inout), optional :: rz_CldTemp (0:kmax, 1:kmax) real(DP), intent(inout), optional :: rz_CldQH2OVap(0:kmax, 1:kmax) real(DP), intent(inout), optional :: rz_CldQH2OLiq(0:kmax, 1:kmax) real(DP), intent(inout), optional :: rz_CldQH2OSol(0:kmax, 1:kmax) ! 作業変数 ! Work variables ! real(DP) :: z_Height (1:kmax) ! ! Height real(DP) :: r_Height (0:kmax) ! ! Height real(DP) :: z_QH2OVapSat(1:kmax) ! 飽和比湿. ! Saturation specific humidity. ! Dry and moist static energy in environment (Env) and cloud (Cld) ! real(DP) :: z_EnvDryStaticEne (1:kmax) real(DP) :: r_EnvDryStaticEne (0:kmax) real(DP) :: z_EnvMoistStaticEne (1:kmax) real(DP) :: r_EnvMoistStaticEne (0:kmax) real(DP) :: z_EnvMoistStaticEneSat(1:kmax) real(DP) :: r_EnvMoistStaticEneSat(0:kmax) real(DP) :: z_EnvCondStaticEne (1:kmax) real(DP) :: r_CldMoistStaticEne (0:kmax) real(DP) :: r_CldCondStaticEne (0:kmax) real(DP) :: CldCondStaticEneCldTop real(DP) :: z_CWF(1:kmax) ! Cloud work function ! (variable for output) real(DP) :: z_EntParam (1:kmax) ! Entrainment factor (variable for output) !!$ real(DP) :: EntParamLL !!$ ! Entrainment factor for a cloud with top at one layer !!$ ! higher level !!$ real(DP) :: EntParamUL !!$ ! Entrainment factor for a cloud with top at one layer !!$ ! lower level ! cloud total water real(DP) :: r_CldQH2OTot(0:kmax) ! cloud total water at cloud top real(DP) :: CldQH2OTotCldTop ! cloud condensate at cloud top real(DP) :: CldQH2OCondCldTop ! cloud ice at cloud top real(DP) :: CldQH2OSolCldTop ! cloud ice at cloud top for save real(DP) :: CldQH2OSolCldTopB real(DP) :: WatFrac ! Variables for debug ! !!$ real(DP) :: xyz_DelVal(0:imax-1, 1:jmax, 1:kmax) !!$ real(DP) :: xy_SumValB(0:imax-1, 1:jmax) !!$ real(DP) :: xy_SumValA(0:imax-1, 1:jmax) !!$ real(DP) :: Ratio real(DP) :: CldTempB real(DP) :: a_DQVapSatDTemp(1:1) real(DP) :: DelTemp real(DP) :: r_CldTemp (0:kmax) real(DP) :: r_CldQH2OVap(0:kmax) real(DP) :: r_CldQH2OLiq(0:kmax) real(DP) :: r_CldQH2OSol(0:kmax) real(DP) :: r_CldHeight (0:kmax) integer :: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer :: m integer :: n ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. relaxed_arakawa_schubert_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 計算時間計測開始 ! Start measurement of computation time ! !!$ call TimesetClockStart( module_name ) call RelaxedArakawaSchubertHeight1D( z_Temp, z_Exner, z_Beta, z_BetaCldTop, z_Height, r_Height ) ! Potential temperature ! z_PotTemp = z_Temp / z_Exner ! Saturation mixing ratio ! z_QH2OVapSat = a_CalcQVapSat( z_Temp, z_Press ) ! Calculation of dry and moist static energies ! z_EnvDryStaticEne = CpDry * z_Temp + Grav * z_Height z_EnvMoistStaticEne = z_EnvDryStaticEne + LatentHeat * z_QH2OVap ! k = 0 r_EnvDryStaticEne (k) = 1.0d100 r_EnvMoistStaticEne(k) = 1.0d100 do k = 1, kmax-1 r_EnvDryStaticEne (k) = ( z_EnvDryStaticEne (k) + z_EnvDryStaticEne (k+1) ) / 2.0_DP r_EnvMoistStaticEne(k) = ( z_EnvMoistStaticEne(k) + z_EnvMoistStaticEne(k+1) ) / 2.0_DP end do k = kmax r_EnvDryStaticEne (k) = z_EnvDryStaticEne (k) r_EnvMoistStaticEne(k) = z_EnvMoistStaticEne(k) ! Calculation of saturated moist static energy ! z_EnvMoistStaticEneSat = z_EnvDryStaticEne + LatentHeat * z_QH2OVapSat ! k = 0 r_EnvMoistStaticEneSat(k) = 1.0d100 do k = 1, kmax-1 r_EnvMoistStaticEneSat(k) = ( z_EnvMoistStaticEneSat(k) + z_EnvMoistStaticEneSat(k+1) ) / 2.0_DP end do k = kmax r_EnvMoistStaticEneSat(k) = z_EnvMoistStaticEneSat(k) ! Calculation of saturated moist static energy ! z_EnvCondStaticEne = z_EnvMoistStaticEne - LatentHeatFusion * z_QH2OSol ! Iteration for entrainment parameter determination ! ! Initialization ! CldQH2OSolCldTop = 0.0_DP ! loop_entparam : do m = 1, 100 ! save current value of cloud ice at cloud top CldQH2OSolCldTopB = CldQH2OSolCldTop ! cloud condensate static energy at cloud top CldCondStaticEneCldTop = z_EnvMoistStaticEneSat(l) - LatentHeatFusion * CldQH2OSolCldTop ! Entrainment parameter ! call RASEntParamWithIce1D( l, z_Beta, z_BetaCldTop, z_PotTemp, z_EnvMoistStaticEne, z_EnvCondStaticEne, CldCondStaticEneCldTop, IndexMixLayTop, EntParam ) ! subroutines below are commented out temporarily !!$ if ( l >= 3 ) then !!$ call RASEntParam1D( & !!$ & l-1, & ! (in) !!$ & z_Beta, z_BetaCldTop, z_PotTemp, & ! (in) !!$ & z_EnvMoistStaticEne, z_EnvMoistStaticEneSat, & ! (in) !!$ & IndexMixLayTop, & ! (in) !!$ & EntParamLL & ! (out) !!$ & ) !!$ else !!$ EntParamLL = 1.0d100 !!$ end if !!$ if ( l <= kmax-1 ) then !!$ call RASEntParam1D( & !!$ & l+1, & ! (in) !!$ & z_Beta, z_BetaCldTop, z_PotTemp, & ! (in) !!$ & z_EnvMoistStaticEne, z_EnvMoistStaticEneSat, & ! (in) !!$ & IndexMixLayTop, & ! (in) !!$ & EntParamUL & ! (out) !!$ & ) !!$ else !!$ EntParamUL = 1.0d100 !!$ end if ! for output z_EntParam(l) = EntParam ! Difference of normalized mass flux ! ! difference of normalized mass flux between layer bottom and top ! z_DelNormMassFlux(1) = 1.0d100 do k = 2, l-1 z_DelNormMassFlux(k) = - EntParam * z_Beta(k) * z_PotTemp(k) end do do k = l, kmax z_DelNormMassFlux(k) = 1.0d100 end do ! ! difference of normalized mass flux between layer bottom and mid-point ! DelNormMassFluxCldTop = - EntParam * z_BetaCldTop(l) * z_PotTemp(l) ! Normalized mass flux ! ! normalized mass flux at layer interface ! r_NormMassFlux(0) = 0.0_DP do k = 1, l-1 if ( k < IndexMixLayTop ) then r_NormMassFlux(k) = 0.0_DP else if ( k == IndexMixLayTop ) then r_NormMassFlux(k) = 1.0_DP else r_NormMassFlux(k) = r_NormMassFlux(k-1) - z_DelNormMassFlux(k) end if end do do k = l, kmax r_NormMassFlux(k) = 0.0_DP end do ! ! normalized mass flux at cloud top (at layer mid-point) ! NormMassFluxCldTop = r_NormMassFlux(l-1) - DelNormMassFluxCldTop ! Liquid water content at top of clouds ! If l is less than xy_IndexMixLayTop(i,j), i.e. the cloud top is below top of ! mixed layer, xy_SumTmp is zero, then, xy_CldQH2OLiqCldTop is also zero. ! if ( l > IndexMixLayTop ) then do k = 0, IndexMixLayTop-1 r_CldQH2OTot(k) = 1.0d100 end do k = IndexMixLayTop r_CldQH2OTot(k) = z_QH2OVap(IndexMixLayTop) do k = IndexMixLayTop+1, l-1 r_CldQH2OTot(k) = r_CldQH2OTot(k-1) - z_DelNormMassFlux(k) * ( z_QH2OVap(k) + z_QH2OLiq(k) + z_QH2OSol(k) ) end do CldQH2OTotCldTop = r_CldQH2OTot(l-1) - DelNormMassFluxCldTop * ( z_QH2OVap(l) + z_QH2OLiq(l) + z_QH2OSol(l) ) do k = l, kmax r_CldQH2OTot(k) = 1.0d100 end do else r_CldQH2OTot = 0.0_DP CldQH2OTotCldTop = 0.0_DP end if CldQH2OCondCldTop = CldQH2OTotCldTop / ( NormMassFluxCldTop + 1.0d-100 ) - z_QH2OVapSat(l) ! Check whether kernel is positive or negative. ! if ( CldQH2OCondCldTop < 0.0_DP ) then FlagNegH2OCondCldTop = .true. else FlagNegH2OCondCldTop = .false. end if ! avoid negative value CldQH2OCondCldTop = max( CldQH2OCondCldTop, 0.0_DP ) !!$ call CloudUtilsWatFraction( & !!$ & z_Temp(l), & ! (in ) !!$ & WatFrac & ! (out) !!$ & ) WatFrac = 1.0_DP CldQH2OLiqCldTop = CldQH2OCondCldTop * WatFrac CldQH2OSolCldTop = CldQH2OCondCldTop - CldQH2OLiqCldTop if ( abs( ( CldQH2OSolCldTop - CldQH2OSolCldTopB ) / ( CldQH2OSolCldTop + 1.0d-100 ) ) < 1.0d-10 ) exit loop_entparam end do loop_entparam if ( m >= 100 ) then call MessageNotify( 'E', module_name, 'Number of loop for entrainment parameter is too large, %d.', i = (/m/) ) end if ! Moist static energy in clouds ! r_CldCondStaticEne(0) = 1.0d100 do k = 1, l-1 if ( k < IndexMixLayTop ) then r_CldCondStaticEne(k) = 1.0d100 else if ( k == IndexMixLayTop ) then r_CldCondStaticEne(k) = z_EnvMoistStaticEne(IndexMixLayTop) else !!$ r_CldMoistStaticEne(k) = & !!$ & ( r_NormMassFlux(k-1) * r_CldMoistStaticEne(k-1) & !!$ & - z_DelNormMassFlux(k) * z_EnvMoistStaticEne(k) ) & !!$ & / r_NormMassFlux(k) r_CldCondStaticEne(k) = ( r_NormMassFlux(k-1) * r_CldCondStaticEne(k-1) - z_DelNormMassFlux(k) * z_EnvCondStaticEne(k) ) / r_NormMassFlux(k) end if end do do k = l, kmax r_CldCondStaticEne(k) = 1.0d100 end do r_CldMoistStaticEne = r_CldCondStaticEne if ( EntParam >= 0.0_DP ) then ! Calculation of cloud air temperature ! This value will not be used below. ! This is an attempt for next extention. ! do k = 0, IndexMixLayTop-1 r_CldTemp (k) = 1.0d100 r_CldQH2OVap(k) = 1.0d100 r_CldQH2OLiq(k) = 1.0d100 r_CldQH2OSol(k) = 1.0d100 r_CldHeight (k) = 1.0d100 end do k = IndexMixLayTop r_CldTemp (k) = z_Temp(k) r_CldQH2OVap(k) = z_QH2OVap(k) r_CldQH2OLiq(k) = 0.0_DP r_CldQH2OSol(k) = 0.0_DP r_CldHeight (k) = r_Height(k) do k = IndexMixLayTop+1, l-1 ! Iteration ! Initialization if ( k == IndexMixLayTop+1 ) then r_CldTemp(k) = z_Temp(k) else r_CldTemp(k) = r_CldTemp(k-1) end if ! loop_cloud_properties : do m = 1, 100 CldTempB = r_CldTemp(k) r_CldQH2OVap(k:k) = a_CalcQVapSat( r_CldTemp(k:k), r_Press(k:k) ) a_DQVapSatDTemp(1:1) = a_CalcDQVapSatDTemp( r_CldTemp(k:k), r_CldQH2OVap(k:k) ) r_CldHeight(k) = r_CldHeight(k-1) + z_Beta(k) * ( r_CldTemp(k-1) + r_CldTemp(k) ) / 2.0_DP / r_Exner(k) DelTemp = ( r_CldMoistStaticEne(k) - CpDry * r_CldTemp(k) - Grav * r_CldHeight(k) - LatentHeat * r_CldQH2OVap(k) ) / ( CpDry + LatentHeat * a_DQVapSatDTemp(1) ) r_CldTemp (k) = r_CldTemp (k) + DelTemp r_CldQH2OVap(k) = r_CldQH2OVap(k) + a_DQVapSatDTemp(1) * DelTemp !!$ write( 6, * ) EntParam, l, k, m, r_CldMoistStaticEne(k), Grav * r_CldHeight(k), r_CldTemp(k), r_CldQH2OVap(k) !!$ if ( abs( CldTempB - r_CldTemp(k) ) / CldTempB < 1.0d-3 ) & if ( abs( DelTemp ) < 1.0d-3 ) exit loop_cloud_properties end do loop_cloud_properties if ( m >= 100 ) then call MessageNotify( 'E', module_name, 'Number of loop for cloud properties is too large, %d.', i = (/m/) ) end if ! cloud water and cloud ice call CloudUtilsWatFraction( r_CldTemp(k), WatFrac ) WatFrac = 1.0_DP ! r_CldQH2OLiq(k) = ( r_CldQH2OTot(k) - r_CldQH2OVap(k) ) * WatFrac r_CldQH2OSol(k) = r_CldQH2OTot(k) - r_CldQH2OVap(k) - r_CldQH2OLiq(k) end do do k = l, kmax r_CldTemp (k) = 1.0d100 r_CldQH2OVap(k) = 1.0d100 r_CldQH2OLiq(k) = 1.0d100 r_CldQH2OSol(k) = 1.0d100 end do else r_CldTemp = 1.0d100 r_CldQH2OVap = 1.0d100 r_CldQH2OLiq = 1.0d100 r_CldQH2OSol = 1.0d100 end if if ( present( rz_CldTemp ) ) rz_CldTemp (:,l) = r_CldTemp if ( present( rz_CldQH2OVap ) ) rz_CldQH2OVap(:,l) = r_CldQH2OVap if ( present( rz_CldQH2OLiq ) ) rz_CldQH2OLiq(:,l) = r_CldQH2OLiq if ( present( rz_CldQH2OSol ) ) rz_CldQH2OSol(:,l) = r_CldQH2OSol !############################################### ! Check whether a parcel in cloud has moist static energy larger than environment's ! !!$ xy_FlagCrossSatEquivPotTemp = .false. !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ do k = xy_IndexMixLayTop(i,j), l-1 !!$ if ( xyr_EnvMoistStaticEneSat(i,j,k) < xyr_CldMoistStaticEne(i,j,k) ) then !!$ xy_FlagCrossSatEquivPotTemp(i,j) = .true. !!$ end if !!$ end do !!$ end do !!$ end do !############################################### ! Cloud work function ! ! Auxiliary variables ! z_Gamma = LatentHeat / CpDry * a_CalcDQVapSatDTemp( z_Temp, z_QH2OVapSat ) ! k = 1 z_Mu (k) = 1.0d100 z_Eps(k) = 1.0d100 do k = 2, kmax z_Mu (k) = ( z_Exner(k ) - r_Exner(k) ) / ( z_Exner(k) * ( 1.0_DP + z_Gamma(k) ) ) z_Eps(k) = ( r_Exner(k-1) - z_Exner(k) ) / ( z_Exner(k) * ( 1.0_DP + z_Gamma(k) ) ) end do ! ! Cloud work function ! CWF = 0.0_DP do k = 2, l-1 if ( k > IndexMixLayTop ) then CWF = CWF + z_Mu (k) * r_NormMassFlux(k ) * ( r_CldMoistStaticEne(k ) - z_EnvMoistStaticEneSat(k) ) + z_Eps(k) * r_NormMassFlux(k-1) * ( r_CldMoistStaticEne(k-1) - z_EnvMoistStaticEneSat(k) ) end if end do k = l if ( k > IndexMixLayTop ) then CWF = CWF + z_Eps(k) * r_NormMassFlux(k-1) * ( r_CldMoistStaticEne(k-1) - z_EnvMoistStaticEneSat(k) ) end if ! for save z_CWF(l) = CWF ! Tendency of dry static energy per unit mass flux ! do k = 1, l z_GammaDSE(k) = - Grav / z_DelPress(k) * ( r_NormMassFlux(k-1) * ( r_EnvDryStaticEne(k-1) - z_EnvDryStaticEne(k) ) + r_NormMassFlux(k ) * ( z_EnvDryStaticEne(k ) - r_EnvDryStaticEne(k) ) ) end do k = l z_GammaDSE(k) = z_GammaDSE(k) - Grav / z_DelPress(k) * LatentHeat * CldQH2OLiqCldTop * NormMassFluxCldTop * ( 1.0_DP - z_RainFactor(k) ) do k = l+1, kmax z_GammaDSE(k) = 0.0_DP end do ! Tendency of moist static energy per unit mass flux ! do k = 1, l z_GammaMSE(k) = - Grav / z_DelPress(k) * ( r_NormMassFlux(k-1) * ( r_EnvMoistStaticEne(k-1) - z_EnvMoistStaticEne(k) ) + r_NormMassFlux(k ) * ( z_EnvMoistStaticEne(k ) - r_EnvMoistStaticEne(k) ) ) end do k = l z_GammaMSE(k) = z_GammaMSE(k) + Grav / z_DelPress(k) * NormMassFluxCldTop * ( z_EnvMoistStaticEneSat(k) - z_EnvMoistStaticEne(k) ) do k = l+1, kmax z_GammaMSE(k) = 0.0_DP end do ! 計算時間計測一時停止 ! Pause measurement of computation time ! !!$ call TimesetClockStop( module_name ) end subroutine RAS1DTestingCore01
Subroutine : | |||||
l : | integer , intent(in ) | ||||
z_DelPress(1:kmax) : | real(DP), intent(in )
| ||||
z_GammaDSE(1:kmax) : | real(DP), intent(in )
| ||||
z_GammaMSE(1:kmax) : | real(DP), intent(in )
| ||||
CldMassFluxBottom : | real(DP), intent(in )
| ||||
z_RainFactor(1:kmax) : | real(DP), intent(in ) | ||||
NormMassFluxCldTop : | real(DP), intent(in ) | ||||
CldQH2OLiqCldTop : | real(DP), intent(in ) | ||||
z_Temp(1:kmax) : | real(DP), intent(inout)
| ||||
z_QH2OVap(1:kmax) : | real(DP), intent(inout)
| ||||
z_DTempDtCumulus(1:kmax) : | real(DP), intent(out )
| ||||
z_DQVapDtCumulus(1:kmax) : | real(DP), intent(out )
| ||||
RainCumulus : | real(DP), intent(inout) |
relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化.
Change temperature and specific humidity by relaxed Arakawa-Schubert scheme
subroutine RAS1DTestingCore02( l, z_DelPress, z_GammaDSE, z_GammaMSE, CldMassFluxBottom, z_RainFactor, NormMassFluxCldTop, CldQH2OLiqCldTop, z_Temp, z_QH2OVap, z_DTempDtCumulus, z_DQVapDtCumulus, RainCumulus ) ! ! relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化. ! ! Change temperature and specific humidity by relaxed Arakawa-Schubert scheme ! ! モジュール引用 ; USE statements ! ! 物理定数設定 ! Physical constants settings ! use constants, only: Grav, GasRDry, CpDry, LatentHeat, LatentHeatFusion ! $ L $ [J kg-1] . ! 融解の潜熱. ! Latent heat of fusion ! 時刻管理 ! Time control ! use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! 飽和比湿の算出 ! Evaluate saturation specific humidity ! use saturate, only: a_CalcQVapSat, a_CalcDQVapSatDTemp ! Arakawa-Schubert scheme by Lord et al. (1982) ! Arakawa-Schubert scheme by Lord et al. (1982) ! use arakawa_schubert_L1982, only : ASL1982CalcCWFCrtl1D ! 雲関系ルーチン ! Cloud-related routines ! use cloud_utils, only : CloudUtilsWatFraction ! 宣言文 ; Declaration statements ! integer , intent(in ) :: l real(DP), intent(in ) :: z_DelPress(1:kmax) ! $ \Delta p $ ! real(DP), intent(in ) :: z_GammaDSE (1:kmax) ! Tendency of dry static energy per unit mass flux real(DP), intent(in ) :: z_GammaMSE (1:kmax) ! Tendency of moist static energy per unit mass flux real(DP), intent(in ) :: CldMassFluxBottom ! Cloud mass flux at cloud bottom real(DP), intent(in ) :: z_RainFactor (1:kmax) real(DP), intent(in ) :: NormMassFluxCldTop real(DP), intent(in ) :: CldQH2OLiqCldTop real(DP), intent(inout) :: z_Temp (1:kmax) ! Temperature real(DP), intent(inout) :: z_QH2OVap (1:kmax) ! $ q $ . 比湿. Specific humidity !!$ real(DP), intent(inout) :: xy_Rain (0:imax-1, 1:jmax) !!$ ! 降水量. !!$ ! Precipitation real(DP), intent(out ) :: z_DTempDtCumulus (1:kmax) ! 温度変化率. ! Temperature tendency real(DP), intent(out ) :: z_DQVapDtCumulus (1:kmax) ! 比湿変化率. ! Specific humidity tendency real(DP), intent(inout) :: RainCumulus ! 作業変数 ! Work variables ! real(DP) :: z_DelH2OMass (1:kmax) real(DP) :: H2OMassB real(DP) :: H2OMassA integer :: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer :: m integer :: n ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. relaxed_arakawa_schubert_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 計算時間計測開始 ! Start measurement of computation time ! !!$ call TimesetClockStart( module_name ) ! Tendencies of specific temperature and humidity ! do k = 1, kmax z_DTempDtCumulus(k) = + CldMassFluxBottom * z_GammaDSE(k) / CpDry z_DQVapDtCumulus(k) = + CldMassFluxBottom * ( z_GammaMSE(k) - z_GammaDSE(k) ) / LatentHeat end do ! add tendencies to temperature and specific humidity ! z_Temp = z_Temp + z_DTempDtCumulus * 2.0_DP * DelTime z_QH2OVap = z_QH2OVap + z_DQVapDtCumulus * 2.0_DP * DelTime ! Precipitation rate at cloud top level ! unit is kg m-2 s-1 ! RainCumulus = CldMassFluxBottom * z_RainFactor(l) * NormMassFluxCldTop * CldQH2OLiqCldTop ! mass fix ! z_DelH2OMass = z_QH2OVap * z_DelPress / Grav ! total H2O mass in a vertical column H2OMassB = 0.0_DP do k = kmax, 1, -1 H2OMassB = H2OMassB + z_DelH2OMass(k) end do if ( H2OMassB < 0.0_DP ) then call MessageNotify( 'E', module_name, 'Mass of water vapor in a column is negative (%d,%d), %f.', i = (/0,0/), d = (/H2OMassB/) ) end if ! negative mass is borrowed from above do k = 1, kmax-1 if ( z_DelH2OMass(k) < 0.0_DP ) then z_DelH2OMass(k+1) = z_DelH2OMass(k+1) + z_DelH2OMass(k) z_DelH2OMass(k ) = 0.0_DP end if end do k = kmax if ( z_DelH2OMass(k) < 0.0_DP ) then z_DelH2OMass (k) = 0.0_DP end if ! total H2O mass in a vertical column, again H2OMassA = 0.0_DP do k = kmax, 1, -1 H2OMassA = H2OMassA + z_DelH2OMass(k) end do ! total mass in a vertical column is adjusted if ( H2OMassA > 0.0_DP ) then do k = 1, kmax z_DelH2OMass(k) = z_DelH2OMass(k) * H2OMassB / H2OMassA end do else do k = 1, kmax z_DelH2OMass(k) = 0.0_DP end do end if z_QH2OVap = z_DelH2OMass / ( z_DelPress / Grav ) ! 計算時間計測一時停止 ! Pause measurement of computation time ! !!$ call TimesetClockStop( module_name ) end subroutine RAS1DTestingCore02
Subroutine : | |
l : | integer , intent(in ) |
z_Beta(1:kmax) : | real(DP), intent(in ) |
z_BetaCldTop(1:kmax) : | real(DP), intent(in ) |
z_PotTemp(1:kmax) : | real(DP), intent(in ) |
z_EnvMoistStaticEne(1:kmax) : | real(DP), intent(in ) |
z_EnvMoistStaticEneSat(1:kmax) : | real(DP), intent(in ) |
IndexMixLayTop : | integer , intent(in ) |
EntParam : | real(DP), intent(out) |
エントレインメントパラメータの計算
Calculation of entrainment parameter
subroutine RASEntParam1D( l, z_Beta, z_BetaCldTop, z_PotTemp, z_EnvMoistStaticEne, z_EnvMoistStaticEneSat, IndexMixLayTop, EntParam ) ! ! エントレインメントパラメータの計算 ! ! Calculation of entrainment parameter ! ! モジュール引用 ; USE statements ! ! 宣言文 ; Declaration statements ! integer , intent(in ) :: l real(DP), intent(in ) :: z_Beta (1:kmax) real(DP), intent(in ) :: z_BetaCldTop (1:kmax) real(DP), intent(in ) :: z_PotTemp (1:kmax) real(DP), intent(in ) :: z_EnvMoistStaticEne (1:kmax) real(DP), intent(in ) :: z_EnvMoistStaticEneSat(1:kmax) integer , intent(in ) :: IndexMixLayTop real(DP), intent(out) :: EntParam ! 作業変数 ! Work variables ! integer :: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction ! 実行文 ; Executable statement ! ! Entrainment parameter ! EntParam = 0.0_DP do k = 2, l-1 if ( k > IndexMixLayTop ) then EntParam = EntParam + z_Beta(k) * z_PotTemp(k) * ( z_EnvMoistStaticEneSat(l) - z_EnvMoistStaticEne(k) ) end if end do if ( l > IndexMixLayTop ) then EntParam = EntParam + z_BetaCldTop(l) * z_PotTemp(l) * ( z_EnvMoistStaticEneSat(l) - z_EnvMoistStaticEne(l) ) EntParam = ( z_EnvMoistStaticEne(IndexMixLayTop) - z_EnvMoistStaticEneSat(l) ) / ( EntParam + 1.0d-100 ) end if end subroutine RASEntParam1D
Subroutine : | |
l : | integer , intent(in ) |
z_Beta(1:kmax) : | real(DP), intent(in ) |
z_BetaCldTop(1:kmax) : | real(DP), intent(in ) |
z_PotTemp(1:kmax) : | real(DP), intent(in ) |
z_EnvMoistStaticEne(1:kmax) : | real(DP), intent(in ) |
z_EnvCondStaticEne(1:kmax) : | real(DP), intent(in ) |
CldCondStaticEneCldTop : | real(DP), intent(in ) |
IndexMixLayTop : | integer , intent(in ) |
EntParam : | real(DP), intent(out) |
エントレインメントパラメータの計算
Calculation of entrainment parameter
subroutine RASEntParamWithIce1D( l, z_Beta, z_BetaCldTop, z_PotTemp, z_EnvMoistStaticEne, z_EnvCondStaticEne, CldCondStaticEneCldTop, IndexMixLayTop, EntParam ) ! ! エントレインメントパラメータの計算 ! ! Calculation of entrainment parameter ! ! モジュール引用 ; USE statements ! ! 宣言文 ; Declaration statements ! integer , intent(in ) :: l real(DP), intent(in ) :: z_Beta (1:kmax) real(DP), intent(in ) :: z_BetaCldTop (1:kmax) real(DP), intent(in ) :: z_PotTemp (1:kmax) real(DP), intent(in ) :: z_EnvMoistStaticEne (1:kmax) real(DP), intent(in ) :: z_EnvCondStaticEne (1:kmax) real(DP), intent(in ) :: CldCondStaticEneCldTop integer , intent(in ) :: IndexMixLayTop real(DP), intent(out) :: EntParam ! 作業変数 ! Work variables ! integer :: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction ! 実行文 ; Executable statement ! ! Entrainment parameter ! EntParam = 0.0_DP do k = 2, l-1 if ( k > IndexMixLayTop ) then EntParam = EntParam + z_Beta(k) * z_PotTemp(k) * ( CldCondStaticEneCldTop - z_EnvCondStaticEne(k) ) end if end do if ( l > IndexMixLayTop ) then EntParam = EntParam + z_BetaCldTop(l) * z_PotTemp(l) * ( CldCondStaticEneCldTop - z_EnvCondStaticEne(l) ) EntParam = ( z_EnvMoistStaticEne(IndexMixLayTop) - CldCondStaticEneCldTop ) / ( EntParam + 1.0d-100 ) end if end subroutine RASEntParamWithIce1D
Variable : | |||
RainConversionFactor : | real(DP), save
|
Subroutine : | |||||
SurfTemp : | real(DP), intent(in )
| ||||
z_Press(1:kmax) : | real(DP), intent(in )
| ||||
r_Press(0:kmax) : | real(DP), intent(in )
| ||||
z_Exner(1:kmax) : | real(DP), intent(in )
| ||||
r_Exner(0:kmax) : | real(DP), intent(in )
| ||||
z_Temp(1:kmax) : | real(DP), intent(inout)
| ||||
z_QH2OVap(1:kmax) : | real(DP), intent(inout)
| ||||
z_DTempDt(1:kmax) : | real(DP), intent(inout)
| ||||
z_DQVapDt(1:kmax) : | real(DP), intent(inout)
| ||||
z_DQH2OLiqDt(1:kmax) : | real(DP), intent(out ) | ||||
z_MoistConvDetTend(1:kmax) : | real(DP), intent(out ), optional | ||||
z_MoistConvSubsidMassFlux(1:kmax) : | real(DP), intent(out ), optional |
relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化.
Change temperature and specific humidity by relaxed Arakawa-Schubert scheme
subroutine RelaxedArakawaSchubert1D( SurfTemp, z_Press, r_Press, z_Exner, r_Exner, z_Temp, z_QH2OVap, z_DTempDt, z_DQVapDt, z_DQH2OLiqDt, z_MoistConvDetTend, z_MoistConvSubsidMassFlux ) ! ! relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化. ! ! Change temperature and specific humidity by relaxed Arakawa-Schubert scheme ! ! モジュール引用 ; 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: a_CalcQVapSat, a_CalcDQVapSatDTemp ! Arakawa-Schubert scheme by Lord et al. (1982) ! Arakawa-Schubert scheme by Lord et al. (1982) ! use arakawa_schubert_L1982, only : ASL1982CalcCWFCrtl1D ! 宣言文 ; Declaration statements ! real(DP), intent(in ) :: SurfTemp ! Pressure real(DP), intent(in ) :: z_Press (1:kmax) ! Pressure real(DP), intent(in ) :: r_Press (0:kmax) ! Pressure real(DP), intent(in ) :: z_Exner (1:kmax) ! Exner function real(DP), intent(in ) :: r_Exner (0:kmax) ! Exner function real(DP), intent(inout) :: z_Temp (1:kmax) ! Temperature real(DP), intent(inout) :: z_QH2OVap (1:kmax) ! $ q $ . 比湿. Specific humidity !!$ real(DP), intent(inout) :: xy_Rain (0:imax-1, 1:jmax) !!$ ! 降水量. !!$ ! Precipitation real(DP), intent(inout) :: z_DTempDt (1:kmax) ! 温度変化率. ! Temperature tendency real(DP), intent(inout) :: z_DQVapDt (1:kmax) ! 比湿変化率. ! Specific humidity tendency real(DP), intent(out ) :: z_DQH2OLiqDt(1:kmax) real(DP), intent(out ), optional :: z_MoistConvDetTend (1:kmax) real(DP), intent(out ), optional :: z_MoistConvSubsidMassFlux(1:kmax) ! 作業変数 ! Work variables ! real(DP) :: z_Height (1:kmax) ! ! Height real(DP) :: r_Height (0:kmax) ! ! Height real(DP) :: RainCumulus ! 降水量. ! Precipitation real(DP) :: z_DTempDtCumulus (1:kmax) ! 温度変化率. ! Temperature tendency real(DP) :: z_DQVapDtCumulus (1:kmax) ! 比湿変化率. ! Specific humidity tendency real(DP) :: z_DelPress(1:kmax) ! $ \Delta p $ ! real(DP) :: z_PotTemp (1:kmax) ! Potential temperature ! real(DP) :: z_QH2OVapSat(1:kmax) ! 飽和比湿. ! Saturation specific humidity. ! Dry and moist static energy in environment (Env) and cloud (Cld) ! real(DP) :: z_EnvDryStaticEne (1:kmax) real(DP) :: r_EnvDryStaticEne (0:kmax) real(DP) :: z_EnvMoistStaticEne (1:kmax) real(DP) :: r_EnvMoistStaticEne (0:kmax) real(DP) :: z_EnvMoistStaticEneSat(1:kmax) real(DP) :: r_EnvMoistStaticEneSat(0:kmax) real(DP) :: r_CldMoistStaticEne (0:kmax) real(DP) :: Kernel ! Tendency of cloud work function by cumulus convection, kernel real(DP) :: CWF ! Cloud work function real(DP) :: z_CWF(1:kmax) ! Cloud work function ! (variable for output) real(DP) :: DCWFDtLS ! Tendency of cloud work function by large scale motion real(DP) :: z_DCWFDtLS(1:kmax) ! Tendency of cloud work function by large scale motion ! (variable for output) real(DP) :: CldMassFluxBottom ! Cloud mass flux at cloud bottom real(DP) :: z_Beta (1:kmax) real(DP) :: z_BetaCldTop (1:kmax) real(DP) :: z_Gamma (1:kmax) real(DP) :: z_GammaDSE (1:kmax) ! Tendency of dry static energy per unit mass flux real(DP) :: z_GammaMSE (1:kmax) ! Tendency of moist static energy per unit mass flux real(DP) :: z_Mu (1:kmax) real(DP) :: z_Eps (1:kmax) real(DP) :: PressCldBase ! Pressure of cloud base real(DP) :: z_CWFCrtl (1:kmax) ! "Critical value" of cloud work function real(DP) :: z_RainFactor (1:kmax) real(DP) :: EntParam ! Entrainment factor real(DP) :: z_EntParam (1:kmax) ! Entrainment factor (variable for output) real(DP) :: EntParamLL ! Entrainment factor for a cloud with top at one layer ! higher level real(DP) :: EntParamUL ! Entrainment factor for a cloud with top at one layer ! lower level ! Difference of normalized mass flux between layer interface real(DP) :: z_DelNormMassFlux (1:kmax) real(DP) :: DelNormMassFluxCldTop ! Normalized mass flux at layer interface and cloud top real(DP) :: r_NormMassFlux (0:kmax) real(DP) :: NormMassFluxCldTop ! Liquid water at cloud top real(DP) :: CldQH2OLiqCldTop ! Mass flux distribution function real(DP) :: z_MassFluxDistFunc (1:kmax) real(DP) :: z_DelH2OMass (1:kmax) real(DP) :: H2OMassB real(DP) :: H2OMassA real(DP) :: z_RainCumulus (1:kmax) real(DP) :: NegDDelLWDt real(DP) :: z_DDelLWDtCCPLV(1:kmax) logical :: FlagCrossSatEquivPotTemp ! ! Flag showing whether a parcel in cloud has moist static ! energy larger than environment's real(DP) :: r_QH2OVapSat (0:kmax) real(DP) :: r_TempAdiabAscent (0:kmax) real(DP) :: SurfPotTemp !!$ real(DP) :: xyz_TempAdiabAscent (0:imax-1, 1:jmax, 1:kmax) ! Variables for looking for top of mixed layer ! logical :: FlagMixLayTopFound integer :: IndexMixLayTop ! Variables for modification of cloud mass flux ! real(DP) :: z_QH2OVapTentative (1:kmax) real(DP) :: CldMassFluxCorFactor real(DP) :: CldMassFluxCorFactorTentative real(DP) :: z_TempB (1:kmax) ! 調節前の温度. ! Temperature before adjustment real(DP) :: z_QH2OVapB(1:kmax) ! 調節前の比湿. ! Specific humidity before adjustment ! Flags for modification of ! logical :: FlagKernelNegative logical :: FlagNegH2OLiqCldTop ! Variables for subsidence mass flux between updrafts ! real(DP) :: DelNormMassFluxHalfLayer real(DP) :: NormMassFlux ! Variables for debug ! !!$ real(DP) :: xyz_DelVal(0:imax-1, 1:jmax, 1:kmax) !!$ real(DP) :: xy_SumValB(0:imax-1, 1:jmax) !!$ real(DP) :: xy_SumValA(0:imax-1, 1:jmax) !!$ real(DP) :: Ratio real(DP) :: r_CldTotWater(0:kmax) real(DP) :: CldTotWaterCldTop real(DP) :: SumTmp integer :: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction integer :: l integer :: m integer :: n ! 実行文 ; Executable statement ! ! 初期化確認 ! Initialization check ! if ( .not. relaxed_arakawa_schubert_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if ! 計算時間計測開始 ! Start measurement of computation time ! !!$ call TimesetClockStart( module_name ) ! 調節前 "Temp", "QH2OVap" の保存 ! Store "Temp", "QH2OVap" before adjustment ! z_TempB = z_Temp z_QH2OVapB = z_QH2OVap ! Preparation of variables ! ! ! Auxiliary variables ! Pressure difference between upper and lower interface of the layer do k = 1, kmax z_DelPress(k) = r_Press(k-1) - r_Press(k) end do ! beta do k = 1, kmax z_Beta(k) = CpDry / Grav * ( r_Exner(k-1) - r_Exner(k) ) end do do k = 1, kmax z_BetaCldTop(k) = CpDry / Grav * ( r_Exner(k-1) - z_Exner(k) ) end do ! ! Search for top of mixed layer (lifting condensation level) based on ! a description in p.684 of Arakawa and Shubert (1974). ! call RelaxedArakawaSchubertHeight1D( z_Temp, z_Exner, z_Beta, z_BetaCldTop, z_Height, r_Height ) ! !==================================== ! !!$ xyz_TempAdiabAscent(:,:,1) = xyz_Temp(:,:,1) !!$ do k = 2, kmax !!$ xyz_TempAdiabAscent(:,:,k) = & !!$ & xyz_Temp(:,:,1) - Grav / CpDry * ( xyz_Height(:,:,k) - xyz_Height(:,:,1) ) !!$ end do !!$ xyz_TempAdiabAscent = max( xyz_TempAdiabAscent, 1.0_DP ) !!$ xyz_QH2OVapSat = xyz_CalcQVapSat( xyz_TempAdiabAscent, xyz_Press ) !!$ xy_IndexMixLayTop = 1 !!$ xy_FlagMixLayTopFound = .false. !!$ do k = 2, kmax !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( ( xyz_QH2OVap(i,j,1) >= xyz_QH2OVapSat(i,j,k) ) .and. & !!$ & ( .not. xy_FlagMixLayTopFound(i,j) ) ) then !!$ xy_IndexMixLayTop (i,j) = k - 1 !!$ xy_FlagMixLayTopFound(i,j) = .true. !!$ end if !!$ end do !!$ end do !!$ end do ! !------------------------------------ ! !!$ xyr_TempAdiabAscent(:,:,0) = xy_SurfTemp !!$ do k = 1, kmax !!$ xyr_TempAdiabAscent(:,:,k) = & !!$ & xy_SurfTemp - Grav / CpDry * ( xyr_Height(:,:,k) - 0.0_DP ) !!$ end do !!$ xyr_TempAdiabAscent = max( xyr_TempAdiabAscent, 1.0_DP ) !!$ r_TempAdiabAscent(0) = SurfTemp SurfPotTemp = SurfTemp / r_Exner(0) do k = 1, kmax r_TempAdiabAscent(k) = SurfPotTemp * r_Exner(k) end do ! r_QH2OVapSat(0 ) = 1.0d100 r_QH2OVapSat(1:kmax-1) = a_CalcQVapSat( r_TempAdiabAscent(1:kmax-1), r_Press(1:kmax-1) ) r_QH2OVapSat(kmax ) = r_QH2OVapSat(kmax-1) ! IndexMixLayTop = 1 FlagMixLayTopFound = .false. do k = 2, kmax if ( ( z_QH2OVap(1) >= r_QH2OVapSat(k) ) .and. ( .not. FlagMixLayTopFound ) ) then IndexMixLayTop = k - 1 FlagMixLayTopFound = .true. end if end do ! !==================================== ! if ( .not. FlagMixLayTopFound ) then IndexMixLayTop = kmax - 1 end if ! ! Critical cloud work function ! if ( FlagZeroCrtlCWF ) then z_CWFCrtl = 0.0_DP else PressCldBase = r_Press(IndexMixLayTop) call ASL1982CalcCWFCrtl1D( PressCldBase, z_Press, z_CWFCrtl ) end if ! ! Rain conversion factor ! if ( RainConversionFactor < 0.0_DP ) then do k = 1, kmax if ( z_Press(k) < 500.0d2 ) then z_RainFactor(k) = 1.0_DP else if ( z_Press(k) < 800.0d2 ) then z_RainFactor(k) = 0.8_DP + ( 800.0d2 - z_Press(k) ) / 1500.0d2 else z_RainFactor(k) = 0.8_DP end if end do else z_RainFactor = RainConversionFactor end if z_RainCumulus (1) = 0.0_DP z_EntParam (1) = 0.0_DP z_CWF (1) = 0.0_DP z_DCWFDtLS (1) = 0.0_DP z_MassFluxDistFunc(1) = 0.0_DP if ( present( z_MoistConvDetTend ) ) then z_MoistConvDetTend(1) = 0.0_DP end if if ( present( z_MoistConvSubsidMassFlux ) ) then ! Subsidence mass flux between the updrafts ! Initialization ! z_MoistConvSubsidMassFlux = 0.0_DP end if loop_cloud_top : do l = 2, kmax call RelaxedArakawaSchubertHeight1D( z_Temp, z_Exner, z_Beta, z_BetaCldTop, z_Height, r_Height ) ! Potential temperature ! z_PotTemp = z_Temp / z_Exner ! Saturation mixing ratio ! z_QH2OVapSat = a_CalcQVapSat( z_Temp, z_Press ) ! Calculation of dry and moist static energies ! z_EnvDryStaticEne = CpDry * z_Temp + Grav * z_Height z_EnvMoistStaticEne = z_EnvDryStaticEne + LatentHeat * z_QH2OVap ! k = 0 r_EnvDryStaticEne (k) = 1.0d100 r_EnvMoistStaticEne(k) = 1.0d100 do k = 1, kmax-1 r_EnvDryStaticEne (k) = ( z_EnvDryStaticEne (k) + z_EnvDryStaticEne (k+1) ) / 2.0_DP r_EnvMoistStaticEne(k) = ( z_EnvMoistStaticEne(k) + z_EnvMoistStaticEne(k+1) ) / 2.0_DP end do k = kmax r_EnvDryStaticEne (k) = z_EnvDryStaticEne (k) r_EnvMoistStaticEne(k) = z_EnvMoistStaticEne(k) ! Calculation of saturated moist static energy ! z_EnvMoistStaticEneSat = z_EnvDryStaticEne + LatentHeat * z_QH2OVapSat ! k = 0 r_EnvMoistStaticEneSat(k) = 1.0d100 do k = 1, kmax-1 r_EnvMoistStaticEneSat(k) = ( z_EnvMoistStaticEneSat(k) + z_EnvMoistStaticEneSat(k+1) ) / 2.0_DP end do k = kmax r_EnvMoistStaticEneSat(k) = z_EnvMoistStaticEneSat(k) ! Entrainment parameter ! call RASEntParam1D( l, z_Beta, z_BetaCldTop, z_PotTemp, z_EnvMoistStaticEne, z_EnvMoistStaticEneSat, IndexMixLayTop, EntParam ) if ( l >= 3 ) then call RASEntParam1D( l-1, z_Beta, z_BetaCldTop, z_PotTemp, z_EnvMoistStaticEne, z_EnvMoistStaticEneSat, IndexMixLayTop, EntParamLL ) else EntParamLL = 1.0d100 end if if ( l <= kmax-1 ) then call RASEntParam1D( l+1, z_Beta, z_BetaCldTop, z_PotTemp, z_EnvMoistStaticEne, z_EnvMoistStaticEneSat, IndexMixLayTop, EntParamUL ) else EntParamUL = 1.0d100 end if ! for output z_EntParam(l) = EntParam ! Difference of normalized mass flux ! ! difference of normalized mass flux between layer bottom and top ! z_DelNormMassFlux(1) = 1.0d100 do k = 2, l-1 z_DelNormMassFlux(k) = - EntParam * z_Beta(k) * z_PotTemp(k) end do do k = l, kmax z_DelNormMassFlux(k) = 1.0d100 end do ! ! difference of normalized mass flux between layer bottom and mid-point ! DelNormMassFluxCldTop = - EntParam * z_BetaCldTop(l) * z_PotTemp(l) ! Normalized mass flux ! ! normalized mass flux at layer interface ! r_NormMassFlux(0) = 0.0_DP do k = 1, l-1 if ( k < IndexMixLayTop ) then r_NormMassFlux(k) = 0.0_DP else if ( k == IndexMixLayTop ) then r_NormMassFlux(k) = 1.0_DP else r_NormMassFlux(k) = r_NormMassFlux(k-1) - z_DelNormMassFlux(k) end if end do do k = l, kmax r_NormMassFlux(k) = 0.0_DP end do ! ! normalized mass flux at cloud top (at layer mid-point) ! NormMassFluxCldTop = r_NormMassFlux(l-1) - DelNormMassFluxCldTop ! Liquid water content at top of clouds ! If l is less than xy_IndexMixLayTop(i,j), i.e. the cloud top is below top of ! mixed layer, xy_SumTmp is zero, then, xy_CldQH2OLiqCldTop is also zero. ! if ( l > IndexMixLayTop ) then !!$ SumTmp = z_QH2OVap(IndexMixLayTop) !!$ do k = IndexMixLayTop+1, l-1 !!$ SumTmp = SumTmp & !!$ & - z_DelNormMassFlux(k) * z_QH2OVap(k) !!$ end do !!$ SumTmp = SumTmp & !!$ & - DelNormMassFluxCldTop * z_QH2OVap(l) do k = 0, IndexMixLayTop-1 r_CldTotWater(k) = 0.0_DP end do k = IndexMixLayTop r_CldTotWater(k) = z_QH2OVap(IndexMixLayTop) do k = IndexMixLayTop+1, l-1 r_CldTotWater(k) = r_CldTotWater(k-1) - z_DelNormMassFlux(k) * z_QH2OVap(k) end do CldTotWaterCldTop = r_CldTotWater(l-1) - DelNormMassFluxCldTop * z_QH2OVap(l) do k = l, kmax r_CldTotWater(k) = 0.0_DP end do else r_CldTotWater = 0.0_DP CldTotWaterCldTop = 0.0_DP end if CldQH2OLiqCldTop = CldTotWaterCldTop / ( NormMassFluxCldTop + 1.0d-100 ) - z_QH2OVapSat(l) ! Check whether kernel is positive or negative. ! if ( CldQH2OLiqCldTop < 0.0_DP ) then FlagNegH2OLiqCldTop = .true. else FlagNegH2OLiqCldTop = .false. end if ! avoid negative value CldQH2OLiqCldTop = max( CldQH2OLiqCldTop, 0.0_DP ) ! Moist static energy in clouds ! r_CldMoistStaticEne(0) = 1.0d100 do k = 1, l-1 if ( k < IndexMixLayTop ) then r_CldMoistStaticEne(k) = 1.0d100 else if ( k == IndexMixLayTop ) then r_CldMoistStaticEne(k) = z_EnvMoistStaticEne(IndexMixLayTop) else r_CldMoistStaticEne(k) = ( r_NormMassFlux(k-1) * r_CldMoistStaticEne(k-1) - z_DelNormMassFlux(k) * z_EnvMoistStaticEne(k) ) / r_NormMassFlux(k) end if end do do k = l, kmax r_CldMoistStaticEne(k) = 1.0d100 end do !############################################### ! Check whether a parcel in cloud has moist static energy larger than environment's ! !!$ xy_FlagCrossSatEquivPotTemp = .false. !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ do k = xy_IndexMixLayTop(i,j), l-1 !!$ if ( xyr_EnvMoistStaticEneSat(i,j,k) < xyr_CldMoistStaticEne(i,j,k) ) then !!$ xy_FlagCrossSatEquivPotTemp(i,j) = .true. !!$ end if !!$ end do !!$ end do !!$ end do !############################################### ! Cloud work function ! ! Auxiliary variables ! z_Gamma = LatentHeat / CpDry * a_CalcDQVapSatDTemp( z_Temp, z_QH2OVapSat ) ! k = 1 z_Mu (k) = 1.0d100 z_Eps(k) = 1.0d100 do k = 2, kmax z_Mu (k) = ( z_Exner(k ) - r_Exner(k) ) / ( z_Exner(k) * ( 1.0_DP + z_Gamma(k) ) ) z_Eps(k) = ( r_Exner(k-1) - z_Exner(k) ) / ( z_Exner(k) * ( 1.0_DP + z_Gamma(k) ) ) end do ! ! Cloud work function ! CWF = 0.0_DP do k = 2, l-1 if ( k > IndexMixLayTop ) then CWF = CWF + z_Mu (k) * r_NormMassFlux(k ) * ( r_CldMoistStaticEne(k ) - z_EnvMoistStaticEneSat(k) ) + z_Eps(k) * r_NormMassFlux(k-1) * ( r_CldMoistStaticEne(k-1) - z_EnvMoistStaticEneSat(k) ) end if end do k = l if ( k > IndexMixLayTop ) then CWF = CWF + z_Eps(k) * r_NormMassFlux(k-1) * ( r_CldMoistStaticEne(k-1) - z_EnvMoistStaticEneSat(k) ) end if ! for save z_CWF(l) = CWF ! Time derivative of cloud work function by large scale motion ! DCWFDtLS = ( CWF - z_CWFCrtl(l) ) / ( 2.0_DP * DelTime ) ! for save z_DCWFDtLS(l) = DCWFDtLS ! Tendency of dry static energy per unit mass flux ! do k = 1, l z_GammaDSE(k) = - Grav / z_DelPress(k) * ( r_NormMassFlux(k-1) * ( r_EnvDryStaticEne(k-1) - z_EnvDryStaticEne(k) ) + r_NormMassFlux(k ) * ( z_EnvDryStaticEne(k ) - r_EnvDryStaticEne(k) ) ) end do k = l z_GammaDSE(k) = z_GammaDSE(k) - Grav / z_DelPress(k) * LatentHeat * CldQH2OLiqCldTop * NormMassFluxCldTop * ( 1.0_DP - z_RainFactor(k) ) do k = l+1, kmax z_GammaDSE(k) = 0.0_DP end do ! Tendency of moist static energy per unit mass flux ! do k = 1, l z_GammaMSE(k) = - Grav / z_DelPress(k) * ( r_NormMassFlux(k-1) * ( r_EnvMoistStaticEne(k-1) - z_EnvMoistStaticEne(k) ) + r_NormMassFlux(k ) * ( z_EnvMoistStaticEne(k ) - r_EnvMoistStaticEne(k) ) ) end do k = l z_GammaMSE(k) = z_GammaMSE(k) + Grav / z_DelPress(k) * NormMassFluxCldTop * ( z_EnvMoistStaticEneSat(k) - z_EnvMoistStaticEne(k) ) do k = l+1, kmax z_GammaMSE(k) = 0.0_DP end do ! Kernel, time derivative of cloud work function by cumulus convection per unit ! mass flux ! Kernel = z_Eps(IndexMixLayTop+1) * z_GammaMSE(IndexMixLayTop) - z_Eps(l) * r_NormMassFlux(l-1) * ( 1.0_DP + z_Gamma(l) ) * z_GammaDSE(l) do n = IndexMixLayTop+1, l-1 SumTmp = 0.0_DP do m = IndexMixLayTop+1, n SumTmp = SumTmp + z_DelNormMassFlux(m) * z_GammaMSE(m) end do Kernel = Kernel + ( z_Eps(n+1) + z_Mu(n) ) * ( z_GammaMSE(IndexMixLayTop) - SumTmp ) - ( z_Eps(n) * r_NormMassFlux(n-1) + z_Mu (n) * r_NormMassFlux(n ) ) * ( 1.0_DP + z_Gamma(n) ) * z_GammaDSE(n) end do ! Check whether kernel is positive or negative. ! if ( Kernel < 0.0_DP ) then FlagKernelNegative = .true. else FlagKernelNegative = .false. end if ! Load et al. (1982), p.108 Kernel = min( Kernel, -5.0d-3 ) ! Cloud mass flux at cloud bottom ! CldMassFluxBottom = - DCWFDtLS / Kernel ! ! mass flux has to be zero or positive CldMassFluxBottom = max( CldMassFluxBottom, 0.0_DP ) ! mass flux is zero if entrainment parameter is zero or negative if ( EntParam <= 0.0_DP ) then CldMassFluxBottom = 0.0_DP end if !!$ ! mass flux is zero if it is below lifting condensation level !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( .not. xy_FlagCrossSatEquivPotTemp(i,j) ) then !!$ xy_CloudMassFluxBottom(i,j) = 0.0_DP !!$ end if !!$ end do !!$ end do ! mass flux is zero if the LNB is unstable for updrafts ! (i.e., if the parcel is positively buoyant just above the LNB). ! See Lord et al. (1982), p.112, for more details. ! Strictly speaking, the process below is different from that ! proposed by Lord et al. (1982). Lord et al. (1982) compare ! entrainment parameters at 3 levels. But, entrainment ! parameters at 2 levels are compared below, because comparison ! of values between 2 levels seems to be sufficient. !!$ if ( ( 3 <= l ) .and. ( l <= kmax-1 ) ) then !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( ( xy_EntParamLL(i,j) < xy_EntParam (i,j) ) .and. & !!$ & ( xy_EntParam (i,j) < xy_EntParamUL(i,j) ) ) then !!$ if ( ( xy_EntParamLL(i,j) > 0.0_DP ) .and. & !!$ & ( xy_EntParam (i,j) > 0.0_DP ) .and. & !!$ & ( xy_EntParamUL(i,j) > 0.0_DP ) ) then !!$ xy_CloudMassFluxBottom(i,j) = 0.0_DP !!$ end if !!$ end if !!$ end do !!$ end do !!$ end if !!$ if ( xy_IndexMixLayTop(i,j) == l ) then !!$ if ( ( xy_EntParam (i,j) > 0.0_DP ) .and. & !!$ & ( xy_EntParamUL(i,j) > 0.0_DP ) ) then !!$ if ( xy_EntParam (i,j) < xy_EntParamUL(i,j) ) then !!$ xy_CloudMassFluxBottom(i,j) = 0.0_DP !!$ end if !!$ end if !!$ else if ( ( xy_IndexMixLayTop(i,j) < l ) .and. ( l <= kmax-1 ) ) then !!$ if ( ( xy_EntParamLL(i,j) > 0.0_DP ) .and. & !!$ & ( xy_EntParam (i,j) > 0.0_DP ) .and. & !!$ & ( xy_EntParamUL(i,j) > 0.0_DP ) ) then !!$ if ( ( xy_EntParamLL(i,j) < xy_EntParam (i,j) ) .and. & !!$ & ( xy_EntParam (i,j) < xy_EntParamUL(i,j) ) ) then if ( ( IndexMixLayTop <= l ) .and. ( l <= kmax-1 ) ) then if ( ( EntParam > 0.0_DP ) .and. ( EntParamUL > 0.0_DP ) ) then if ( EntParam < EntParamUL ) then CldMassFluxBottom = 0.0_DP end if end if end if ! ! mass flux is zero unless kernel is negative ! if ( .not. FlagKernelNegative ) then CldMassFluxBottom = 0.0_DP end if ! ! mass flux is zero if liquid water at a cloud top is negative ! if ( FlagNegH2OLiqCldTop ) then CldMassFluxBottom = 0.0_DP end if ! ! multiply factor ! CldMassFluxBottom = CldMassFluxBottom * min( 2.0_DP * DelTime / AdjTimeConst, 1.0_DP ) ! ! for output z_MassFluxDistFunc(l) = CldMassFluxBottom ! Check values of cloud mass flux ! If water vapor amount transported by convection is larger than that in a ! column, cloud mass flux is reduced. ! ! tendency of specific humidity is calculated tentatively do k = 1, kmax z_DQVapDtCumulus(k) = + CldMassFluxBottom * ( z_GammaMSE(k) - z_GammaDSE(k) ) / LatentHeat end do ! total H2O mass in a vertical column after RAS z_QH2OVapTentative = z_QH2OVap + z_DQVapDtCumulus * 2.0_DP * DelTime CldMassFluxCorFactor = 1.0_DP do k = 1, kmax if ( z_QH2OVapTentative(k) < 0.0_DP ) then CldMassFluxCorFactorTentative = z_QH2OVap(k) / ( z_QH2OVap(k) - z_QH2OVapTentative(k) ) else CldMassFluxCorFactorTentative = 1.0_DP end if if ( CldMassFluxCorFactorTentative < CldMassFluxCorFactor ) then CldMassFluxCorFactor = CldMassFluxCorFactorTentative end if end do ! modify cloud mass flux CldMassFluxBottom = CldMassFluxCorFactor * CldMassFluxBottom !!$ do k = 1, kmax !!$ xyz_DQVapDtCumulus(:,:,k) = & !!$ & + xy_CloudMassFluxBottom * ( xyz_GammaMSE(:,:,k) - xyz_GammaDSE(:,:,k) ) & !!$ & / LatentHeat !!$ end do !!$ ! total H2O mass in a vertical column before RAS !!$ xyz_DelH2OMass = xyz_QH2OVap * xyz_DelPress / Grav !!$ xy_H2OMassB = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_H2OMassB = xy_H2OMassB + xyz_DelH2OMass(:,:,k) !!$ end do !!$ ! total H2O mass in a vertical column after RAS !!$ xyz_QH2OVapTentative = xyz_QH2OVap + xyz_DQVapDtCumulus * 2.0_DP * DelTime !!$ xyz_DelH2OMass = xyz_QH2OVapTentative * xyz_DelPress / Grav !!$ xy_H2OMassA = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_H2OMassA = xy_H2OMassA + xyz_DelH2OMass(:,:,k) !!$ end do !!$ ! modify cloud mass flux !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( xy_H2OMassA(i,j) < 0.0_DP ) then !!$ ! A safety factor ( 1.0_DP + 1.0d-5 ) is arbitrary. !!$ xy_CloudMassFluxBottom(i,j) = xy_CloudMassFluxBottom(i,j) & !!$ & * xy_H2OMassB(i,j) & !!$ & / ( ( xy_H2OMassB(i,j) - xy_H2OMassA(i,j) ) * ( 1.0_DP + 1.0d-5 ) ) !!$ end if !!$ end do !!$ end do ! Tendencies of specific temperature and humidity ! do k = 1, kmax z_DTempDtCumulus(k) = + CldMassFluxBottom * z_GammaDSE(k) / CpDry z_DQVapDtCumulus(k) = + CldMassFluxBottom * ( z_GammaMSE(k) - z_GammaDSE(k) ) / LatentHeat end do !!$ ! !!$ ! modification of tendency of temperature and water vapor in the mixed layer !!$ ! !!$ if ( FlagUniformMixedLayer ) then !!$ xy_SumTmp = 0.0_DP !!$ do k = 1, kmax !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( k <= xy_IndexMixLayTop(i,j) ) then !!$ xy_SumTmp(i,j) = xy_SumTmp(i,j) & !!$ & + xyz_DTempDtCumulus(i,j,k) & !!$ & * ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) ) !!$ end if !!$ end do !!$ end do !!$ end do !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ xy_SumTmp(i,j) = xy_SumTmp(i,j) & !!$ & / ( xyr_Press(i,j,0) - xyr_Press(i,j,xy_IndexMixLayTop(i,j)) ) !!$ end do !!$ end do !!$ do k = 1, kmax !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( k <= xy_IndexMixLayTop(i,j) ) then !!$ xyz_DTempDtCumulus(i,j,k) = xy_SumTmp(i,j) !!$ end if !!$ end do !!$ end do !!$ end do !!$ ! !!$ xy_SumTmp = 0.0_DP !!$ do k = 1, kmax !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( k <= xy_IndexMixLayTop(i,j) ) then !!$ xy_SumTmp(i,j) = xy_SumTmp(i,j) & !!$ & + xyz_DQVapDtCumulus(i,j,k) & !!$ & * ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) ) !!$ end if !!$ end do !!$ end do !!$ end do !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ xy_SumTmp(i,j) = xy_SumTmp(i,j) & !!$ & / ( xyr_Press(i,j,0) - xyr_Press(i,j,xy_IndexMixLayTop(i,j)) ) !!$ end do !!$ end do !!$ do k = 1, kmax !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( k <= xy_IndexMixLayTop(i,j) ) then !!$ xyz_DQVapDtCumulus(i,j,k) = xy_SumTmp(i,j) !!$ end if !!$ end do !!$ end do !!$ end do !!$ end if ! add tendencies to temperature and specific humidity ! z_Temp = z_Temp + z_DTempDtCumulus * 2.0_DP * DelTime z_QH2OVap = z_QH2OVap + z_DQVapDtCumulus * 2.0_DP * DelTime ! Precipitation rate at cloud top level ! unit is kg m-2 s-1 ! z_RainCumulus(l) = CldMassFluxBottom * z_RainFactor(l) * NormMassFluxCldTop * CldQH2OLiqCldTop ! mass fix ! z_DelH2OMass = z_QH2OVap * z_DelPress / Grav ! total H2O mass in a vertical column H2OMassB = 0.0_DP do k = kmax, 1, -1 H2OMassB = H2OMassB + z_DelH2OMass(k) end do if ( H2OMassB < 0.0_DP ) then call MessageNotify( 'E', module_name, 'Mass of water vapor in a column is negative (%d,%d), %f.', i = (/0,0/), d = (/H2OMassB/) ) end if ! negative mass is borrowed from above do k = 1, kmax-1 if ( z_DelH2OMass(k) < 0.0_DP ) then z_DelH2OMass(k+1) = z_DelH2OMass(k+1) + z_DelH2OMass(k) z_DelH2OMass(k ) = 0.0_DP end if end do k = kmax if ( z_DelH2OMass(k) < 0.0_DP ) then !!$ call MessageNotify( 'E', module_name, & !!$ & 'Mass of water vapor in the top layer is negative (%d,%d,%d), %f.', & !!$ & i = (/i,j,k/), d = (/xyz_DelH2OMass(i,j,k)/) ) !!$ !!$ xyz_RainCumulus(i,j,l) = xyz_RainCumulus(i,j,l) & !!$ & - xyz_DelH2OMass(i,j,k) / ( 2.0_DP * DelTime ) !!$ xyz_Temp (i,j,k) = xyz_Temp(i,j,k) & !!$ & - LatentHeat * xyz_DelH2OMass(i,j,k) / ( xyz_DelPress(i,j,k) / Grav )& !!$ & / CpDry z_DelH2OMass (k) = 0.0_DP end if !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ if ( xyz_RainCumulus(i,j,l) < 0.0_DP ) then !!$ call MessageNotify( 'E', module_name, & !!$ & 'Mass of water vapor is insufficient at (%d,%d,%d), %f.', & !!$ & i = (/i,j,k/), d = (/xyz_RainCumulus(i,j,l)/) ) !!$ end if !!$ end do !!$ end do ! total H2O mass in a vertical column, again H2OMassA = 0.0_DP do k = kmax, 1, -1 H2OMassA = H2OMassA + z_DelH2OMass(k) end do ! total mass in a vertical column is adjusted if ( H2OMassA > 0.0_DP ) then !!$ write( 6, * ) i, j, xy_H2OMassB(i,j), xy_H2OMassB(i,j) / xy_H2OMassA(i,j) do k = 1, kmax z_DelH2OMass(k) = z_DelH2OMass(k) * H2OMassB / H2OMassA end do else do k = 1, kmax z_DelH2OMass(k) = 0.0_DP end do end if z_QH2OVap = z_DelH2OMass / ( z_DelPress / Grav ) ! Detrainment mass tendency per unit mass (kg m-3 s-1 / ( kg m-3 ) = s-1). ! This corresponds to condensation rate (kg m-2 s-1) divided by layer thickness (m) ! and density (kg m-3), in other words. ! kg m-2 s-1 / ( Pa / ( m s-2 ) ) ! = kg m-2 s-1 Pa-1 m s-1 = kg m-2 (kg m s-2 m-2)-1 m s-2 ! = kg m-2 s-1 kg-1 m-1 s2 m2 m s-2 = s-1 if ( present( z_MoistConvDetTend ) ) then z_MoistConvDetTend(l) = CldMassFluxBottom * NormMassFluxCldTop / ( z_DelPress(l) / Grav ) end if if ( present( z_MoistConvSubsidMassFlux ) ) then ! Subsidence mass flux between the updrafts do k = 1, l-1 if ( k > IndexMixLayTop ) then DelNormMassFluxHalfLayer = - EntParam * z_BetaCldTop(k) * z_PotTemp(k) NormMassFlux = r_NormMassFlux(k-1) - DelNormMassFluxHalfLayer z_MoistConvSubsidMassFlux(k) = z_MoistConvSubsidMassFlux(k) + CldMassFluxBottom * NormMassFlux end if end do end if end do loop_cloud_top ! 温度変化率, 比湿変化率 ! Calculate specific humidity tendency and temperature tendency ! (In fact, temperature tendency does not need to calculate, here.) ! z_DTempDtCumulus = ( z_Temp - z_TempB ) / ( 2.0_DP * DelTime ) z_DQVapDtCumulus = ( z_QH2OVap - z_QH2OVapB ) / ( 2.0_DP * DelTime ) z_DTempDt = z_DTempDt + z_DTempDtCumulus z_DQVapDt = z_DQVapDt + z_DQVapDtCumulus ! Precipitation rate at the surface ! unit is kg m-2 s-1 ! !!$ xy_RainCumulus = 0.0d0 !!$ do k = kmax, 1, -1 !!$ xy_RainCumulus = xy_RainCumulus + xyz_RainCumulus(:,:,k) !!$ end do z_DQH2OLiqDt = z_RainCumulus / ( z_DelPress / Grav ) !!$ 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 !!$ !!$ xy_Rain = xy_Rain + xy_RainCumulus ! Calculation for debug ! check of conservation of water amount and internal energy ! !!$ xyz_DelVal = xyz_QH2OVapB * xyz_DelPress / Grav !!$ xy_SumValB = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValB = xy_SumValB + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xyz_DelVal = xyz_QH2OVap * xyz_DelPress / Grav !!$ xy_SumValA = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValA = xy_SumValA + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xy_SumValA = xy_SumValA + xy_RainCumulus * 2.0_DP * DelTime !!$ ! !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ Ratio = ( xy_SumValA(i,j) - xy_SumValB(i,j) ) & !!$ & / max( xy_SumValA(i,j), 1.0d-100 ) !!$ if ( abs( Ratio ) > 1.0d-14 ) then !!$ write( 6, * ) 'H2O: ', i, j, & !!$ & xy_SumValB(i,j), xy_SumValA(i,j), & !!$ & xy_RainCumulus(i,j) * 2.0_DP * DelTime, & !!$ & Ratio !!$ end if !!$ end do !!$ end do !!$ ! !!$ ! !!$ xyz_DelVal = CpDry * xyz_TempB * xyz_DelPress / Grav !!$ xy_SumValB = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValB = xy_SumValB + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xyz_DelVal = CpDry * xyz_Temp * xyz_DelPress / Grav !!$ xy_SumValA = 0.0_DP !!$ do k = kmax, 1, -1 !!$ xy_SumValA = xy_SumValA + xyz_DelVal(:,:,k) !!$ end do !!$ ! !!$ xy_SumValA = xy_SumValA - LatentHeat * xy_RainCumulus * 2.0_DP * DelTime !!$ ! !!$ do j = 1, jmax !!$ do i = 0, imax-1 !!$ Ratio = ( xy_SumValA(i,j) - xy_SumValB(i,j) ) & !!$ & / max( xy_SumValA(i,j), 1.0d-100 ) !!$ if ( abs( Ratio ) > 1.0d-14 ) then !!$ write( 6, * ) 'CpT: ', i, j, & !!$ & xy_SumValB(i,j), xy_SumValA(i,j), & !!$ & - LatentHeat * xy_RainCumulus(i,j) * 2.0_DP * DelTime, & !!$ & Ratio !!$ end if !!$ end do !!$ end do ! calculation for output ! This calculation is meaningless because RainCumulus is not used below. z_RainCumulus = z_DQH2OLiqDt * ( z_DelPress / Grav ) RainCumulus = 0.0d0 do k = kmax, 1, -1 RainCumulus = RainCumulus + z_RainCumulus(k) end do !!$ if ( present( xyz_DQH2OLiqDt ) ) then !!$ !!$ ! unit is kg m-2 s-1 !!$ xyz_DDelLWDtCCPLV = xyz_RainCumulus !!$ !!$ ! 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 !!$ !!$ ! unit is s-1 !!$ xyz_DQH2OLiqDt = xyz_DDelLWDtCCPLV / ( xyz_DelPress / Grav ) !!$ !!$ end if ! 計算時間計測一時停止 ! Pause measurement of computation time ! !!$ call TimesetClockStop( module_name ) end subroutine RelaxedArakawaSchubert1D
Subroutine : | |
l : | integer , intent(in ) |
xyz_Beta(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
xyz_BetaCldTop(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
xyz_PotTemp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
xyz_EnvMoistStaticEne(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
xyz_EnvMoistStaticEneSat(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
xy_IndexMixLayTop(0:imax-1, 1:jmax) : | integer , intent(in ) |
xy_EntParam(0:imax-1, 1:jmax) : | real(DP), intent(out) |
エントレインメントパラメータの計算
Calculation of entrainment parameter
subroutine RelaxedArakawaSchubertEntParam( l, xyz_Beta, xyz_BetaCldTop, xyz_PotTemp, xyz_EnvMoistStaticEne, xyz_EnvMoistStaticEneSat, xy_IndexMixLayTop, xy_EntParam ) ! ! エントレインメントパラメータの計算 ! ! Calculation of entrainment parameter ! ! モジュール引用 ; USE statements ! ! 宣言文 ; Declaration statements ! integer , intent(in ) :: l real(DP), intent(in ) :: xyz_Beta (0:imax-1, 1:jmax, 1:kmax) real(DP), intent(in ) :: xyz_BetaCldTop (0:imax-1, 1:jmax, 1:kmax) real(DP), intent(in ) :: xyz_PotTemp (0:imax-1, 1:jmax, 1:kmax) real(DP), intent(in ) :: xyz_EnvMoistStaticEne (0:imax-1, 1:jmax, 1:kmax) real(DP), intent(in ) :: xyz_EnvMoistStaticEneSat(0:imax-1, 1:jmax, 1:kmax) integer , intent(in ) :: xy_IndexMixLayTop (0:imax-1, 1:jmax) real(DP), intent(out) :: xy_EntParam (0:imax-1, 1:jmax) ! 作業変数 ! Work variables ! 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 ! ! Entrainment parameter ! xy_EntParam = 0.0_DP do k = 2, l-1 do j = 1, jmax do i = 0, imax-1 if ( k > xy_IndexMixLayTop(i,j) ) then xy_EntParam(i,j) = xy_EntParam(i,j) + xyz_Beta(i,j,k) * xyz_PotTemp(i,j,k) * ( xyz_EnvMoistStaticEneSat(i,j,l) - xyz_EnvMoistStaticEne(i,j,k) ) end if end do end do end do do j = 1, jmax do i = 0, imax-1 if ( l > xy_IndexMixLayTop(i,j) ) then xy_EntParam(i,j) = xy_EntParam(i,j) + xyz_BetaCldTop(i,j,l) * xyz_PotTemp(i,j,l) * ( xyz_EnvMoistStaticEneSat(i,j,l) - xyz_EnvMoistStaticEne(i,j,l) ) xy_EntParam(i,j) = ( xyz_EnvMoistStaticEne(i,j,xy_IndexMixLayTop(i,j)) - xyz_EnvMoistStaticEneSat(i,j,l) ) / ( xy_EntParam(i,j) + 1.0d-100 ) end if end do end do end subroutine RelaxedArakawaSchubertEntParam
Subroutine : | |
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
xyz_Beta(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
xyz_BetaCldTop(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
xyz_Height(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out) |
xyr_Height(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) |
高度の計算
Calculation of height
subroutine RelaxedArakawaSchubertHeight( xyz_Temp, xyz_Exner, xyz_Beta, xyz_BetaCldTop, xyz_Height, xyr_Height ) ! ! 高度の計算 ! ! Calculation of height ! ! モジュール引用 ; USE statements ! ! 宣言文 ; Declaration statements ! real(DP), intent(in ) :: xyz_Temp (0:imax-1, 1:jmax, 1:kmax) real(DP), intent(in ) :: xyz_Exner (0:imax-1, 1:jmax, 1:kmax) real(DP), intent(in ) :: xyz_Beta (0:imax-1, 1:jmax, 1:kmax) real(DP), intent(in ) :: xyz_BetaCldTop(0:imax-1, 1:jmax, 1:kmax) real(DP), intent(out) :: xyz_Height (0:imax-1, 1:jmax, 1:kmax) real(DP), intent(out) :: xyr_Height (0:imax-1, 1:jmax, 0:kmax) ! 作業変数 ! Work variables ! real(DP) :: xyz_PotTemp(0:imax-1, 1:jmax, 1:kmax) integer :: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction ! 実行文 ; Executable statement ! xyz_PotTemp = xyz_Temp / xyz_Exner xyr_Height(:,:,0) = 0.0_DP do k = 1, kmax xyz_Height(:,:,k) = xyr_Height(:,:,k-1) + xyz_BetaCldTop(:,:,k) * xyz_PotTemp(:,:,k) xyr_Height(:,:,k) = xyr_Height(:,:,k-1) + xyz_Beta (:,:,k) * xyz_PotTemp(:,:,k) end do end subroutine RelaxedArakawaSchubertHeight
Subroutine : | |
z_Temp(1:kmax) : | real(DP), intent(in ) |
z_Exner(1:kmax) : | real(DP), intent(in ) |
z_Beta(1:kmax) : | real(DP), intent(in ) |
z_BetaCldTop(1:kmax) : | real(DP), intent(in ) |
z_Height(1:kmax) : | real(DP), intent(out) |
r_Height(0:kmax) : | real(DP), intent(out) |
高度の計算
Calculation of height
subroutine RelaxedArakawaSchubertHeight1D( z_Temp, z_Exner, z_Beta, z_BetaCldTop, z_Height, r_Height ) ! ! 高度の計算 ! ! Calculation of height ! ! モジュール引用 ; USE statements ! ! 宣言文 ; Declaration statements ! real(DP), intent(in ) :: z_Temp (1:kmax) real(DP), intent(in ) :: z_Exner (1:kmax) real(DP), intent(in ) :: z_Beta (1:kmax) real(DP), intent(in ) :: z_BetaCldTop(1:kmax) real(DP), intent(out) :: z_Height (1:kmax) real(DP), intent(out) :: r_Height (0:kmax) ! 作業変数 ! Work variables ! real(DP) :: z_PotTemp(1:kmax) character(STRING) :: VarName integer :: k ! 鉛直方向に回る DO ループ用作業変数 ! Work variables for DO loop in vertical direction ! 実行文 ; Executable statement ! z_PotTemp = z_Temp / z_Exner r_Height(0) = 0.0_DP do k = 1, kmax z_Height(k) = r_Height(k-1) + z_BetaCldTop(k) * z_PotTemp(k) r_Height(k) = r_Height(k-1) + z_Beta (k) * z_PotTemp(k) end do end subroutine RelaxedArakawaSchubertHeight1D
Constant : | |||
module_name = ‘relaxed_arakawa_schubert‘ : | character(*), parameter
|
Variable : | |||
relaxed_arakawa_schubert_inited = .false. : | logical, save
|
Constant : | |||
version = ’$Name: dcpam5-20140204-5 $’ // ’$Id: relaxed_arakawa_schubert.f90,v 1.8 2014-02-04 10:24:42 yot Exp $’ : | character(*), parameter
|