Class major_comp_phase_change
In: saturate/major_comp_phase_change.f90

主成分相変化

Phase change of atmospheric major component

Note that Japanese and English are described in parallel.

Procedures List

!$ ! DryConvAdjust :乾燥対流調節
!$ ! ———— :————
!$ ! DryConvAdjust :Dry convective adjustment

NAMELIST

NAMELIST#major_comp_phase_change_nml

Methods

Included Modules

gridset dc_types namelist_util dc_message timeset gtool_historyauto constants saturate_major_comp composition auxiliary mass_fixer dc_iounit dc_string

Public Instance methods

Subroutine :
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in )
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (half level)
xyz_Press(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
: $ p $ . 気圧 (整数レベル). Air pressure (full level)
xy_Ps(0:imax-1, 1:jmax) :real(DP), intent(inout)
: $ T $ . 温度. Temperature
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(inout)
: $ T $ . 温度. Temperature
xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) :real(DP), intent(inout)
xyz_U(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(inout)
xyz_V(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(inout)
xy_SurfMajCompIce(0:imax-1, 1:jmax) :real(DP), intent(inout)
: Surface major component ice amount

CO2 相変化

CO2 phase change

[Source]

  subroutine MajorCompPhaseChangeInAtm( xyr_Press, xyz_Press, xy_Ps, xyz_Temp, xyzf_QMix, xyz_U, xyz_V, xy_SurfMajCompIce )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

    ! モジュール引用 ; USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only : ncmax, IndexTKE

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    use auxiliary, only: AuxVars

    ! 主成分相変化
    ! Phase change of atmospheric major component
    !
    use saturate_major_comp, only : SaturateMajorCompCondTemp, SaturateMajorCompInqLatentHeat

    ! 質量の補正
    ! Mass fixer
    !
    use mass_fixer, only: MassFixerColumn

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in   ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in   ):: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
                              ! $ p $ . 気圧 (整数レベル). 
                              ! Air pressure (full level)
    real(DP), intent(inout):: xy_Ps            (0:imax-1, 1:jmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyz_Temp         (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyzf_QMix        (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP), intent(inout):: xyz_U            (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout):: xyz_V            (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout):: xy_SurfMajCompIce(0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount

    ! 作業変数
    ! Work variables
    !
    real(DP):: LatentHeatMajCompSubl

    real(DP):: xy_PsB              (0:imax-1, 1:jmax)
    real(DP):: xy_PsA              (0:imax-1, 1:jmax)
    real(DP):: xyr_PressB          (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_PressA          (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyz_TempB           (0:imax-1, 1:jmax, 1:kmax)
                              ! 調節前の温度. 
                              ! Temperature before adjustment
    real(DP):: xyzf_QMixB          (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)

    real(DP):: xyz_DelAtmMass      (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xy_TempTmp          (0:imax-1, 1:jmax)
    real(DP):: xy_DTempDtSubl      (0:imax-1, 1:jmax)
    real(DP):: xy_DTempDtCond      (0:imax-1, 1:jmax)
    real(DP):: xyr_AtmMassFallFlux (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyz_TempTmp         (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_QH2OVapTmp      (0:imax-1, 1:jmax, 1:kmax)

    real(DP):: xyz_DTempDt         (0:imax-1, 1:jmax, 1:kmax)
                              ! 温度変化率. 
                              ! Temperature tendency
    real(DP):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
                              ! 
                              ! Surface major component ice tendency
    real(DP):: xy_DPsDt            (0:imax-1, 1:jmax)

    real(DP):: xy_DSurfMajCompIceDtOneLayer(0:imax-1, 1:jmax)
    real(DP):: xyr_DPressDt                (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyz_TempCond   (0:imax-1, 1:jmax, 1:kmax)

    integer :: mmax
    real(DP):: xyza_Array       (0:imax-1, 1:jmax, 1:kmax, 1:ncmax+1+1)
    logical :: a_FlagSurfaceSink                          (1:ncmax+1+1)
    real(DP):: xyra_MassFlow    (0:imax-1, 1:jmax, 0:kmax, 1:ncmax+1+1)

    real(DP):: xyrf_MassFlow    (0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
    real(DP):: xyr_MomXFlow     (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_MomYFlow     (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyz_DelAtmMassB   (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_DelAtmMassA   (0:imax-1, 1:jmax, 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:: m
    integer:: n

    logical :: FlagCheckPs


    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. major_comp_phase_change_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    if ( .not. FlagMajCompPhaseChange ) return


    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! Set latent heat
    LatentHeatMajCompSubl = SaturateMajorCompInqLatentHeat()


    FlagCheckPs = .false.
    do j = 1, jmax
      do i = 0, imax-1
        if ( xyr_Press(i,j,0) > 1.0e4_DP ) then
          FlagCheckPs = .true.
        end if
      end do
    end do
    if ( FlagCheckPs ) then
      call MessageNotify( 'W', module_name, 'Surface pressure is greater than 10000 Pa.' )
    end if


    ! Store variables
    !
    xyz_TempB  = xyz_Temp
    xyzf_QMixB = xyzf_QMix


    call SaturateMajorCompCondTemp( xyz_Press, xyz_TempCond )


    do k = 1, kmax
      xyz_DelAtmMass(:,:,k) = ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do

    xyr_AtmMassFallFlux(:,:,kmax) = 0.0_DP
    do k = kmax, 1, -1
      ! sublimation of falling condensate
!!$      xy_DTempDtSubl = &
!!$        & - LatentHeatMajCompSubl * xyr_AtmMassFallFlux(:,:,k) &
!!$        &   / ( CpDry * xyz_DelAtmMass(:,:,k) )
!!$      xyz_Temp(:,:,k) = xyz_Temp(:,:,k) + xy_DTempDtSubl * ( 2.0_DP * DelTime )
!!$      xyr_AtmMassFallFlux(:,:,k) = 0.0_DP

      ! condensation
      xy_TempTmp = xyz_Temp(:,:,k)
      xyz_Temp(:,:,k) = max( xyz_TempCond(:,:,k), xyz_Temp(:,:,k) )
      xy_DTempDtCond = ( xyz_Temp(:,:,k) - xy_TempTmp ) / ( 2.0_DP * DelTime )
      xy_DSurfMajCompIceDtOneLayer = CpDry * xyz_DelAtmMass(:,:,k) * xy_DTempDtCond / LatentHeatMajCompSubl
      xyr_AtmMassFallFlux(:,:,k-1) = xyr_AtmMassFallFlux(:,:,k) + xy_DSurfMajCompIceDtOneLayer
    end do

    xyr_DPressDt = - xyr_AtmMassFallFlux * Grav


    xy_DPsDt             = xyr_DPressDt(:,:,0)
    xy_DSurfMajCompIceDt = - xy_DPsDt / Grav



    ! packing
    if ( FlagModMom ) then
      mmax = ncmax + 1 + 1
    else
      mmax = ncmax
    end if
    do m = 1, ncmax
      n = m
      xyza_Array(:,:,:,m) = xyzf_QMix(:,:,:,n)
      if ( n == IndexTKE ) then
        a_FlagSurfaceSink(m) = .true.
      else
        a_FlagSurfaceSink(m) = .false.
      end if
    end do
    if ( FlagModMom ) then
      m = ncmax
      m = m + 1
      xyza_Array(:,:,:,m) = xyz_U
      a_FlagSurfaceSink(m) = .true.
      m = m + 1
      xyza_Array(:,:,:,m) = xyz_V
      a_FlagSurfaceSink(m) = .true.
    end if

    call MajorCompPhaseChangeCalcFlow( xyr_Press, xyr_DPressDt, mmax, a_FlagSurfaceSink(1:mmax), xyza_Array(:,:,:,1:mmax), xyra_MassFlow(:,:,:,1:mmax) )

    ! unpacking
    do m = 1, ncmax
      xyrf_MassFlow(:,:,:,m) = xyra_MassFlow(:,:,:,m)
    end do
    if ( FlagModMom ) then
      m = ncmax
      m = m + 1
      xyr_MomXFlow = xyra_MassFlow(:,:,:,m)
      m = m + 1
      xyr_MomYFlow = xyra_MassFlow(:,:,:,m)
    end if


    ! Adjustment
    !   preparation
    xy_PsB = xy_Ps
    xy_PsA = xy_PsB + xy_DPsDt * ( 2.0_DP * DelTime )

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    xyz_TempTmp    = 300.0_DP
    xyz_QH2OVapTmp =   0.0_DP
    call AuxVars( xy_PsB, xyz_TempTmp, xyz_QH2OVapTmp, xyr_Press = xyr_PressB )
    call AuxVars( xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp, xyr_Press = xyr_PressA )
    do k = 1, kmax
      xyz_DelAtmMassB(:,:,k) = ( xyr_PressB(:,:,k-1) - xyr_PressB(:,:,k) ) / Grav
      xyz_DelAtmMassA(:,:,k) = ( xyr_PressA(:,:,k-1) - xyr_PressA(:,:,k) ) / Grav
    end do
    !   Atmospheric composition
    do n = 1, ncmax
      do k = 1, kmax
        xyzf_QMix(:,:,k,n) = (   xyz_DelAtmMassB(:,:,k) * xyzf_QMix(:,:,k,n) - ( xyrf_MassFlow(:,:,k,n) - xyrf_MassFlow(:,:,k-1,n) ) ) / xyz_DelAtmMassA(:,:,k)
      end do
    end do
    if ( FlagModMom ) then
      do k = 1, kmax
        ! Zonal wind
        xyz_U(:,:,k) = (   xyz_DelAtmMassB(:,:,k) * xyz_U(:,:,k) - ( xyr_MomXFlow(:,:,k) - xyr_MomXFlow(:,:,k-1) ) ) / xyz_DelAtmMassA(:,:,k)
        ! Meridional wind
        xyz_V(:,:,k) = (   xyz_DelAtmMassB(:,:,k) * xyz_V(:,:,k) - ( xyr_MomYFlow(:,:,k) - xyr_MomYFlow(:,:,k-1) ) ) / xyz_DelAtmMassA(:,:,k)
      end do
    end if


    ! Surface major component ice adjustment
    xy_SurfMajCompIce = xy_SurfMajCompIce + xy_DSurfMajCompIceDt * ( 2.0_DP * DelTime )
    ! Surface pressure adjustment
    xy_Ps = xy_PsA


    call AuxVars( xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp, xyr_Press = xyr_PressA )
    ! 成分の質量の補正
    ! Fix masses of constituents
    !
    call MassFixerColumn( xyr_PressA, xyzf_QMix )

    ! Check
    call MajorCompPhaseChangeConsChk( a_FlagSurfaceSink, xyz_DelAtmMassB, xyzf_QMixB, xyz_DelAtmMassA, xyzf_QMix )

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DTempDtMajCompPhaseChange', xyz_DTempDt )


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine MajorCompPhaseChangeInAtm
Subroutine :
ArgFlagMajCompPhaseChange :logical , intent(in)
CondMajCompName :character(*), intent(in)

major_comp_phase_change モジュールの初期化を行います. NAMELIST#major_comp_phase_change_nml の読み込みはこの手続きで行われます.

"major_comp_phase_change" module is initialized. "NAMELIST#major_comp_phase_change_nml" is loaded in this procedure.

This procedure input/output NAMELIST#major_comp_phase_change_nml .

[Source]

  subroutine MajorCompPhaseChangeInit( ArgFlagMajCompPhaseChange, CondMajCompName )
    !
    ! major_comp_phase_change モジュールの初期化を行います. 
    ! NAMELIST#major_comp_phase_change_nml の読み込みはこの手続きで行われます. 
    !
    ! "major_comp_phase_change" module is initialized. 
    ! "NAMELIST#major_comp_phase_change_nml" is loaded in this procedure. 
    !

    ! モジュール引用 ; USE statements
    !

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: StoA

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! 補助的な変数を計算するサブルーチン・関数群
    ! Subroutines and functions for calculating auxiliary variables
    !
    use auxiliary, only : AuxVarsInit

    ! 主成分相変化
    ! Phase change of atmospheric major component
    !
    use saturate_major_comp, only : SaturateMajorCompInit

    ! 質量の補正
    ! Mass fixer
    !
    use mass_fixer, only : MassFixerInit


    ! 宣言文 ; Declaration statements
    !
    implicit none

    logical     , intent(in) :: ArgFlagMajCompPhaseChange
    character(*), intent(in) :: CondMajCompName


    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /major_comp_phase_change_nml/ FlagModMom

          ! デフォルト値については初期化手続 "major_comp_phase_change#MajorCompPhaseChangeInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "major_comp_phase_change#MajorCompPhaseChangeInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( major_comp_phase_change_inited ) return


    FlagMajCompPhaseChange = ArgFlagMajCompPhaseChange


    ! デフォルト値の設定
    ! Default values settings
    !
    FlagModMom = .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 = major_comp_phase_change_nml, iostat = iostat_nml )                 ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
!      if ( iostat_nml == 0 ) write( STDOUT, nml = cumulus_adjust_nml )
    end if


    if ( FlagMajCompPhaseChange ) then
      ! 主成分相変化
      ! Phase change of atmospheric major component
      !
      call SaturateMajorCompInit( CondMajCompName )
    end if

    ! 補助的な変数を計算するサブルーチン・関数群
    ! Subroutines and functions for calculating auxiliary variables
    !
    call AuxVarsInit

    ! 質量の補正
    ! Mass fixer
    !
    call MassFixerInit


    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'DSurfTempDtMajCompPhaseChange', (/ 'lon ', 'lat ', 'time' /), 'heating by major component phase change', 'K s-1' )
    call HistoryAutoAddVariable( 'DTempDtMajCompPhaseChange', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'heating by major component phase change', 'K s-1' )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  FlagModMom = %b', l = (/ FlagModMom /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    major_comp_phase_change_inited = .true.

  end subroutine MajorCompPhaseChangeInit
Subroutine :
xy_DPsDt(0:imax-1, 1:jmax) :real(DP), intent(in )
xy_DSurfMajCompIceDt(0:imax-1, 1:jmax) :real(DP), intent(in )
xy_Ps(0:imax-1, 1:jmax) :real(DP), intent(inout)
xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) :real(DP), intent(inout)
xyz_U(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(inout)
xyz_V(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(inout)
xy_SurfMajCompIce(0:imax-1, 1:jmax) :real(DP), intent(inout)
: Surface major component ice amount

CO2 相変化

CO2 phase change

[Source]

  subroutine MajorCompPhaseChangeOnGround( xy_DPsDt, xy_DSurfMajCompIceDt, xy_Ps, xyzf_QMix, xyz_U, xyz_V, xy_SurfMajCompIce )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

    ! モジュール引用 ; USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only : ncmax, IndexTKE

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    use auxiliary, only: AuxVars

    ! 質量の補正
    ! Mass fixer
    !
    use mass_fixer, only: MassFixerColumn


    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in   ):: xy_DPsDt            (0:imax-1, 1:jmax)
    real(DP), intent(in   ):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
    real(DP), intent(inout):: xy_Ps            (0:imax-1, 1:jmax)
    real(DP), intent(inout):: xyzf_QMix        (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP), intent(inout):: xyz_U            (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout):: xyz_V            (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout):: xy_SurfMajCompIce   (0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount

    ! 作業変数
    ! Work variables
    !
    real(DP):: xy_PsB              (0:imax-1, 1:jmax)
    real(DP):: xy_PsA              (0:imax-1, 1:jmax)
    real(DP):: xyr_PressB          (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_PressA          (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyzf_QMixB          (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)

    real(DP):: xyz_TempTmp         (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_QH2OVapTmp      (0:imax-1, 1:jmax, 1:kmax)

    real(DP):: xyr_DPressDt        (0:imax-1, 1:jmax, 0:kmax)

    integer :: mmax
    real(DP):: xyza_Array       (0:imax-1, 1:jmax, 1:kmax, 1:ncmax+1+1)
    logical :: a_FlagSurfaceSink                          (1:ncmax+1+1)
    real(DP):: xyra_MassFlow    (0:imax-1, 1:jmax, 0:kmax, 1:ncmax+1+1)

    real(DP):: xyrf_MassFlow    (0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
    real(DP):: xyr_MomXFlow     (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_MomYFlow     (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyz_DelAtmMassB   (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_DelAtmMassA   (0:imax-1, 1:jmax, 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:: m
    integer:: n


    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. major_comp_phase_change_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    if ( .not. FlagMajCompPhaseChange ) return


    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    xy_PsB = xy_Ps
    xy_PsA = xy_PsB + xy_DPsDt * ( 2.0_DP * DelTime )

    xyzf_QMixB = xyzf_QMix

    xyz_TempTmp    = 300.0_DP
    xyz_QH2OVapTmp =   0.0_DP

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    call AuxVars( xy_PsB, xyz_TempTmp, xyz_QH2OVapTmp, xyr_Press = xyr_PressB )
    call AuxVars( xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp, xyr_Press = xyr_PressA )


    xyr_DPressDt = ( xyr_PressA - xyr_PressB ) / ( 2.0_DP * Deltime )


    ! packing
    if ( FlagModMom ) then
      mmax = ncmax + 1 + 1
    else
      mmax = ncmax
    end if
    do m = 1, ncmax
      n = m
      xyza_Array(:,:,:,m) = xyzf_QMix(:,:,:,n)
      if ( n == IndexTKE ) then
        a_FlagSurfaceSink(m) = .true.
      else
        a_FlagSurfaceSink(m) = .false.
      end if
    end do
    if ( FlagModMom ) then
      m = ncmax
      m = m + 1
      xyza_Array(:,:,:,m) = xyz_U
      a_FlagSurfaceSink(m) = .true.
      m = m + 1
      xyza_Array(:,:,:,m) = xyz_V
      a_FlagSurfaceSink(m) = .true.
    end if

    call MajorCompPhaseChangeCalcFlow( xyr_PressB, xyr_DPressDt, mmax, a_FlagSurfaceSink(1:mmax), xyza_Array(:,:,:,1:mmax), xyra_MassFlow(:,:,:,1:mmax) )

    ! unpacking
    do m = 1, ncmax
      xyrf_MassFlow(:,:,:,m) = xyra_MassFlow(:,:,:,m)
    end do
    if ( FlagModMom ) then
      m = ncmax
      m = m + 1
      xyr_MomXFlow = xyra_MassFlow(:,:,:,m)
      m = m + 1
      xyr_MomYFlow = xyra_MassFlow(:,:,:,m)
    end if


    ! Adjustment
    !   preparation
    do k = 1, kmax
      xyz_DelAtmMassB(:,:,k) = ( xyr_PressB(:,:,k-1) - xyr_PressB(:,:,k) ) / Grav
      xyz_DelAtmMassA(:,:,k) = ( xyr_PressA(:,:,k-1) - xyr_PressA(:,:,k) ) / Grav
    end do
    !   Atmospheric composition
    do n = 1, ncmax
      do k = 1, kmax
        xyzf_QMix(:,:,k,n) = (   xyz_DelAtmMassB(:,:,k) * xyzf_QMix(:,:,k,n) - ( xyrf_MassFlow(:,:,k,n) - xyrf_MassFlow(:,:,k-1,n) ) ) / xyz_DelAtmMassA(:,:,k)
      end do
    end do
    if ( FlagModMom ) then
      do k = 1, kmax
        ! Zonal wind
        xyz_U(:,:,k) = (   xyz_DelAtmMassB(:,:,k) * xyz_U(:,:,k) - ( xyr_MomXFlow(:,:,k) - xyr_MomXFlow(:,:,k-1) ) ) / xyz_DelAtmMassA(:,:,k)
        ! Meridional wind
        xyz_V(:,:,k) = (   xyz_DelAtmMassB(:,:,k) * xyz_V(:,:,k) - ( xyr_MomYFlow(:,:,k) - xyr_MomYFlow(:,:,k-1) ) ) / xyz_DelAtmMassA(:,:,k)
      end do
    end if

    !   Surface major component ice
    xy_SurfMajCompIce = xy_SurfMajCompIce + xy_DSurfMajCompIceDt * ( 2.0_DP * DelTime )
    !   Surface pressure
    xy_Ps = xy_PsA


    call AuxVars( xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp, xyr_Press = xyr_PressA )
    ! 成分の質量の補正
    ! Fix masses of constituents
    !
    call MassFixerColumn( xyr_PressA, xyzf_QMix )


    ! Check
    call MajorCompPhaseChangeConsChk( a_FlagSurfaceSink, xyz_DelAtmMassB, xyzf_QMixB, xyz_DelAtmMassA, xyzf_QMix )


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )


  end subroutine MajorCompPhaseChangeOnGround
major_comp_phase_change_inited
Variable :
major_comp_phase_change_inited = .false. :logical, save, public
: 初期設定フラグ. Initialization flag

Private Instance methods

FlagMajCompPhaseChange
Variable :
FlagMajCompPhaseChange :logical, save
FlagModMom
Variable :
FlagModMom :logical, save
Subroutine :
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in )
: pressure
xyr_DPressDt(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in )
mmax :integer , intent(in )
a_FlagSurfaceSink(1:mmax) :logical , intent(in )
xyza_Array(0:imax-1, 1:jmax, 1:kmax, 1:mmax) :real(DP), intent(in )
xyra_MassFlow(0:imax-1, 1:jmax, 0:kmax, 1:mmax) :real(DP), intent(out)

CO2 相変化

CO2 phase change

[Source]

  subroutine MajorCompPhaseChangeCalcFlow( xyr_Press, xyr_DPressDt, mmax, a_FlagSurfaceSink, xyza_Array, xyra_MassFlow )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

    ! モジュール引用 ; USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime               ! $ \Delta t $

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav                  ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in ):: xyr_Press    (0:imax-1, 1:jmax, 0:kmax)
                              ! pressure
    real(DP), intent(in ):: xyr_DPressDt (0:imax-1, 1:jmax, 0:kmax)
    integer , intent(in ):: mmax
    logical , intent(in ):: a_FlagSurfaceSink(1:mmax)
    real(DP), intent(in ):: xyza_Array   (0:imax-1, 1:jmax, 1:kmax, 1:mmax)
    real(DP), intent(out):: xyra_MassFlow(0:imax-1, 1:jmax, 0:kmax, 1:mmax)

    ! 作業変数
    ! Work variables
    !
    real(DP):: xyr_DPPress  (0:imax-1, 1:jmax, 0:kmax)
                              ! pressure at departure point
    real(DP):: DelAtmMass
    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:: k2              ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: m


    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. major_comp_phase_change_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    if ( .not. FlagMajCompPhaseChange ) return


    xyr_DPPress = xyr_Press + xyr_DPressDt * ( 2.0_DP * DelTime )

    ! check
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyr_DPPress(i,j,k-1) < xyr_DPPress(i,j,k) ) then
            call MessageNotify( 'E', module_name, 'Order of departure points are inappropriate, P(k=%d)=%f < P(k=%d)=%f.', i = (/ k-1, k /), d = (/ xyr_DPPress(i,j,k-1), xyr_DPPress(i,j,k) /) )
          end if
        end do
      end do
    end do

    xyra_MassFlow = 0.0_DP
    do k = 0, kmax-1
      do j = 1, jmax
        do i = 0, imax-1

!!$          if ( xyr_DPressDt(i,j,k) >= 0.0_DP ) then
          if ( xyr_DPPress(i,j,k) > xyr_Press(i,j,k) ) then

            sum_upward_mass_transport : do k2 = k, 1, -1
              if ( xyr_DPPress(i,j,k) > xyr_Press(i,j,k2-1) ) then
                DelAtmMass = ( xyr_Press(i,j,k2-1) - xyr_Press(i,j,k2) ) / Grav
                do m = 1, mmax
                  xyra_MassFlow(i,j,k,m) = xyra_MassFlow(i,j,k,m) + xyza_Array(i,j,k2,m) * DelAtmMass
                end do
              else
                DelAtmMass = ( xyr_DPPress(i,j,k) - xyr_Press(i,j,k2) ) / Grav
                do m = 1, mmax
                  xyra_MassFlow(i,j,k,m) = xyra_MassFlow(i,j,k,m) + xyza_Array(i,j,k2,m) * DelAtmMass
                end do
                exit sum_upward_mass_transport
              end if
            end do sum_upward_mass_transport

          else

            sum_downward_mass_transport : do k2 = k+1, kmax
              if ( xyr_DPPress(i,j,k) < xyr_Press(i,j,k2  ) ) then
                DelAtmMass = ( xyr_Press(i,j,k2-1) - xyr_Press(i,j,k2) ) / Grav
                do m = 1, mmax
                  xyra_MassFlow(i,j,k,m) = xyra_MassFlow(i,j,k,m) - xyza_Array(i,j,k2,m) * DelAtmMass
                end do
              else
                DelAtmMass = ( xyr_Press(i,j,k2-1) - xyr_DPPress(i,j,k) ) / Grav
                do m = 1, mmax
                  xyra_MassFlow(i,j,k,m) = xyra_MassFlow(i,j,k,m) - xyza_Array(i,j,k2,m) * DelAtmMass
                end do
                exit sum_downward_mass_transport
              end if
            end do sum_downward_mass_transport

          end if

        end do
      end do

    end do


    do k = 0, 0
      do j = 1, jmax
        do i = 0, imax-1

          ! not surface sink
          if ( xyr_DPressDt(i,j,k) <= 0.0_DP ) then
            do m = 1, mmax
              if ( .not. a_FlagSurfaceSink(m) ) then
                xyra_MassFlow(i,j,k,m) = 0.0_DP
              end if
            end do
          end if

        end do
      end do

    end do



  end subroutine MajorCompPhaseChangeCalcFlow
Subroutine :
a_FlagSurfaceSink(1:ncmax) :logical , intent(in)
xyz_DelAtmMassB(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, ncmax) :real(DP), intent(in)
xyz_DelAtmMassA(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, ncmax) :real(DP), intent(in)

[Source]

  subroutine MajorCompPhaseChangeConsChk( a_FlagSurfaceSink, xyz_DelAtmMassB, xyzf_QMixB, xyz_DelAtmMassA, xyzf_QMixA )

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only : ncmax

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure

    logical , intent(in) :: a_FlagSurfaceSink(1:ncmax)
    real(DP), intent(in) :: xyz_DelAtmMassB(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyzf_QMixB     (0:imax-1, 1:jmax, 1:kmax, ncmax)
    real(DP), intent(in) :: xyz_DelAtmMassA(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyzf_QMixA     (0:imax-1, 1:jmax, 1:kmax, ncmax)

    ! Local variables
    !
    real(DP) :: ValB
    real(DP) :: ValA
    real(DP) :: xyf_SumB(0:imax-1, 1:jmax, 1:ncmax)
    real(DP) :: xyf_SumA(0:imax-1, 1:jmax, 1:ncmax)
    real(DP) :: Ratio
    integer  :: i
    integer  :: j
    integer  :: k
    integer  :: n


    xyf_SumB = 0.0_DP
    xyf_SumA = 0.0_DP
    do n = 1, ncmax
      do k = kmax, 1, -1
        xyf_SumB(:,:,n) = xyf_SumB(:,:,n) + xyz_DelAtmMassB(:,:,k) * xyzf_QMixB(:,:,k,n)
        xyf_SumA(:,:,n) = xyf_SumA(:,:,n) + xyz_DelAtmMassA(:,:,k) * xyzf_QMixA(:,:,k,n)
      end do
    end do
    do n = 1, ncmax
      if ( .not. a_FlagSurfaceSink(n) ) then
        do j = 1, jmax
          do i = 0, imax-1
            ValB = xyf_SumB(i,j,n)
            ValA = xyf_SumA(i,j,n)

            Ratio = ( ValA - ValB ) / ( ValA + 1.0d-100 )
            if ( abs( Ratio ) > 1.0d-10 ) then
              if ( ( ValB < 0.0_DP ) .and. ( abs( ValB ) < 1.0e-20_DP ) ) then
                ! Do nothing
              else
                call MessageNotify( 'M', module_name, 'Mass No. %d is not conserved, %f.', i = (/ n /), d = (/ Ratio /) )
              end if
            end if
          end do
        end do
      end if
    end do


  end subroutine MajorCompPhaseChangeConsChk
Subroutine :
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in )
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (half level)
xyz_Press(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
: $ p $ . 気圧 (整数レベル). Air pressure (full level)
xy_Ps(0:imax-1, 1:jmax) :real(DP), intent(inout)
: $ T $ . 温度. Temperature
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(inout)
: $ T $ . 温度. Temperature
xy_SurfMajCompIce(0:imax-1, 1:jmax) :real(DP), intent(inout)
: Surface major component ice amount

CO2 相変化

CO2 phase change

[Source]

  subroutine MajorCompPhaseChangeInAtmBK( xyr_Press, xyz_Press, xy_Ps, xyz_Temp, xy_SurfMajCompIce )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

    ! モジュール引用 ; USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure

    ! 主成分相変化
    ! Phase change of atmospheric major component
    !
    use saturate_major_comp, only : SaturateMajorCompCondTemp, SaturateMajorCompInqLatentHeat


    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in   ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in   ):: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
                              ! $ p $ . 気圧 (整数レベル). 
                              ! Air pressure (full level)
    real(DP), intent(inout):: xy_Ps            (0:imax-1, 1:jmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyz_Temp         (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xy_SurfMajCompIce(0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount

    ! 作業変数
    ! Work variables
    !
    real(DP):: xyz_TempB           (0:imax-1, 1:jmax, 1:kmax)
                              ! 調節前の温度. 
                              ! Temperature before adjustment
    real(DP):: xyz_DTempDt         (0:imax-1, 1:jmax, 1:kmax)
                              ! 温度変化率. 
                              ! Temperature tendency
    real(DP):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
                              ! 
                              ! Surface major component ice tendency

    real(DP):: xyz_TempCond   (0:imax-1, 1:jmax, 1:kmax)

    real(DP):: LatentHeatMajCompSubl

    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

    logical :: FlagCheckPs


    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. major_comp_phase_change_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    if ( .not. FlagMajCompPhaseChange ) return


    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! Set latent heat
    LatentHeatMajCompSubl = SaturateMajorCompInqLatentHeat()


    FlagCheckPs = .false.
    do j = 1, jmax
      do i = 0, imax-1
        if ( xyr_Press(i,j,0) > 1.0e4_DP ) then
          FlagCheckPs = .true.
        end if
      end do
    end do
    if ( FlagCheckPs ) then
      call MessageNotify( 'W', module_name, 'Surface pressure is greater than 10000 Pa.' )
    end if


    ! 調節前 "Temp" の保存
    ! Store "Temp" before adjustment
    !
    xyz_TempB    = xyz_Temp


    call SaturateMajorCompCondTemp( xyz_Press, xyz_TempCond )


    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_Temp(i,j,k) < xyz_TempCond(i,j,k) ) then
            xyz_Temp(i,j,k) = xyz_TempCond(i,j,k)
          end if
        end do
      end do
    end do

    ! 温度変化率
    ! Calculate temperature tendency
    !
    xyz_DTempDt = ( xyz_Temp - xyz_TempB ) / ( 2.0_DP * DelTime )


    !
    ! Surface major component ice adjustment
    !
    xy_DSurfMajCompIceDt = 0.0_DP
    do k = kmax, 1, -1
      xy_DSurfMajCompIceDt = xy_DSurfMajCompIceDt + CpDry * xyz_DTempDt(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav / LatentHeatMajCompSubl
    end do
    xy_SurfMajCompIce = xy_SurfMajCompIce + xy_DSurfMajCompIceDt * ( 2.0_DP * DelTime )


    !
    ! Surface pressure adjustment
    !
    xy_Ps = xy_Ps - xy_DSurfMajCompIceDt * Grav * ( 2.0_DP * DelTime )


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DTempDtMajCompPhaseChange', xyz_DTempDt )


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine MajorCompPhaseChangeInAtmBK
Subroutine :
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in )
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (half level)
xyz_Press(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
: $ p $ . 気圧 (整数レベル). Air pressure (full level)
xy_Ps(0:imax-1, 1:jmax) :real(DP), intent(inout)
: $ T $ . 温度. Temperature
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(inout)
: $ T $ . 温度. Temperature
xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) :real(DP), intent(inout)
xy_SurfMajCompIce(0:imax-1, 1:jmax) :real(DP), intent(inout)
: Surface major component ice amount

CO2 相変化

CO2 phase change

[Source]

  subroutine MajorCompPhaseChangeInAtmBK2( xyr_Press, xyz_Press, xy_Ps, xyz_Temp, xyzf_QMix, xy_SurfMajCompIce )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

    ! モジュール引用 ; USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only : ncmax, IndexTKE

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    use auxiliary, only: AuxVars

    ! 主成分相変化
    ! Phase change of atmospheric major component
    !
    use saturate_major_comp, only : SaturateMajorCompCondTemp, SaturateMajorCompInqLatentHeat

    ! 質量の補正
    ! Mass fixer
    !
    use mass_fixer, only: MassFixerColumn

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in   ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in   ):: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
                              ! $ p $ . 気圧 (整数レベル). 
                              ! Air pressure (full level)
    real(DP), intent(inout):: xy_Ps            (0:imax-1, 1:jmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyz_Temp         (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyzf_QMix        (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP), intent(inout):: xy_SurfMajCompIce(0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount

    ! 作業変数
    ! Work variables
    !
    real(DP):: LatentHeatMajCompSubl

    real(DP):: xy_PsB              (0:imax-1, 1:jmax)
    real(DP):: xy_PsA              (0:imax-1, 1:jmax)
    real(DP):: xyr_PressB          (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_PressA          (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyz_TempB           (0:imax-1, 1:jmax, 1:kmax)
                              ! 調節前の温度. 
                              ! Temperature before adjustment
    real(DP):: xyzf_QMixB          (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)


    real(DP):: xyz_TempTmp         (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_QH2OVapTmp      (0:imax-1, 1:jmax, 1:kmax)

    real(DP):: xyz_DTempDt         (0:imax-1, 1:jmax, 1:kmax)
                              ! 温度変化率. 
                              ! Temperature tendency
    real(DP):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
                              ! 
                              ! Surface major component ice tendency
    real(DP):: xy_DPsDt            (0:imax-1, 1:jmax)

    real(DP):: xy_DSurfMajCompIceDtOneLayer(0:imax-1, 1:jmax)
    real(DP):: xyr_DPressDt                (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyz_TempCond   (0:imax-1, 1:jmax, 1:kmax)

    integer :: mmax
    real(DP):: xyza_Array       (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    logical :: a_FlagSurfaceSink(1:ncmax)
    real(DP):: xyra_MassFlow    (0:imax-1, 1:jmax, 0:kmax, 1:ncmax)

    real(DP):: xyrf_MassFlow    (0:imax-1, 1:jmax, 0:kmax, 1:ncmax)

    real(DP):: xyz_DelAtmMassB   (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_DelAtmMassA   (0:imax-1, 1:jmax, 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:: m
    integer:: n

    logical :: FlagCheckPs


    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. major_comp_phase_change_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    if ( .not. FlagMajCompPhaseChange ) return


    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! Set latent heat
    LatentHeatMajCompSubl = SaturateMajorCompInqLatentHeat()


    FlagCheckPs = .false.
    do j = 1, jmax
      do i = 0, imax-1
        if ( xyr_Press(i,j,0) > 1.0e4_DP ) then
          FlagCheckPs = .true.
        end if
      end do
    end do
    if ( FlagCheckPs ) then
      call MessageNotify( 'W', module_name, 'Surface pressure is greater than 10000 Pa.' )
    end if


    ! Store variables
    !
    xyz_TempB  = xyz_Temp
    xyzf_QMixB = xyzf_QMix


    call SaturateMajorCompCondTemp( xyz_Press, xyz_TempCond )


    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_Temp(i,j,k) < xyz_TempCond(i,j,k) ) then
            xyz_Temp(i,j,k) = xyz_TempCond(i,j,k)
          end if
        end do
      end do
    end do

    ! 温度変化率
    ! Calculate temperature tendency
    !
    xyz_DTempDt = ( xyz_Temp - xyz_TempB ) / ( 2.0_DP * DelTime )


    xyr_DPressDt(:,:,kmax) = 0.0_DP
    do k = kmax, 1, -1
      xy_DSurfMajCompIceDtOneLayer = CpDry * xyz_DTempDt(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav / LatentHeatMajCompSubl
      xyr_DPressDt(:,:,k-1) = xyr_DPressDt(:,:,k) - xy_DSurfMajCompIceDtOneLayer * Grav
    end do

    xy_DPsDt             = xyr_DPressDt(:,:,0)
    xy_DSurfMajCompIceDt = - xy_DPsDt / Grav



    ! packing
    mmax = ncmax
    do m = 1, mmax
      n = m
      xyza_Array(:,:,:,m) = xyzf_QMix(:,:,:,n)
      if ( n == IndexTKE ) then
        a_FlagSurfaceSink(m) = .true.
      else
        a_FlagSurfaceSink(m) = .false.
      end if
    end do

    call MajorCompPhaseChangeCalcFlow( xyr_Press, xyr_DPressDt, mmax, a_FlagSurfaceSink, xyza_Array, xyra_MassFlow )

    ! unpacking
    do m = 1, mmax
      xyrf_MassFlow(:,:,:,m) = xyra_MassFlow(:,:,:,m)
    end do


    ! Adjustment
    !   preparation
    xy_PsB = xy_Ps
    xy_PsA = xy_PsB + xy_DPsDt * ( 2.0_DP * DelTime )

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    xyz_TempTmp    = 300.0_DP
    xyz_QH2OVapTmp =   0.0_DP
    call AuxVars( xy_PsB, xyz_TempTmp, xyz_QH2OVapTmp, xyr_Press = xyr_PressB )
    call AuxVars( xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp, xyr_Press = xyr_PressA )
    do k = 1, kmax
      xyz_DelAtmMassB(:,:,k) = ( xyr_PressB(:,:,k-1) - xyr_PressB(:,:,k) ) / Grav
      xyz_DelAtmMassA(:,:,k) = ( xyr_PressA(:,:,k-1) - xyr_PressA(:,:,k) ) / Grav
    end do
    !   Atmospheric composition
    do n = 1, ncmax
      do k = 1, kmax
        xyzf_QMix(:,:,k,n) = (   xyz_DelAtmMassB(:,:,k) * xyzf_QMix(:,:,k,n) - ( xyrf_MassFlow(:,:,k,n) - xyrf_MassFlow(:,:,k-1,n) ) ) / xyz_DelAtmMassA(:,:,k)
      end do
    end do


    ! Surface major component ice adjustment
    xy_SurfMajCompIce = xy_SurfMajCompIce + xy_DSurfMajCompIceDt * ( 2.0_DP * DelTime )
    ! Surface pressure adjustment
    xy_Ps = xy_PsA


    call AuxVars( xy_PsA, xyz_TempTmp, xyz_QH2OVapTmp, xyr_Press = xyr_PressA )
    ! 成分の質量の補正
    ! Fix masses of constituents
    !
    call MassFixerColumn( xyr_PressA, xyzf_QMix )

    ! Check
    call MajorCompPhaseChangeConsChk( a_FlagSurfaceSink, xyz_DelAtmMassB, xyzf_QMixB, xyz_DelAtmMassA, xyzf_QMix )

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DTempDtMajCompPhaseChange', xyz_DTempDt )


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine MajorCompPhaseChangeInAtmBK2
Subroutine :
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in )
: $ hat{p} $ . 気圧 (半整数レベル). Air pressure (half level)
xyz_Press(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
: $ p $ . 気圧 (整数レベル). Air pressure (full level)
xyz_Height(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xy_Ps(0:imax-1, 1:jmax) :real(DP), intent(inout)
: $ T $ . 温度. Temperature
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(inout)
: $ T $ . 温度. Temperature
xy_SurfMajCompIce(0:imax-1, 1:jmax) :real(DP), intent(inout)
: Surface major component ice amount

CO2 相変化

CO2 phase change

[Source]

  subroutine MajorCompPhaseChangeInAtmTest( xyr_Press, xyz_Press, xyz_Height, xy_Ps, xyz_Temp, xy_SurfMajCompIce )
    !
    ! CO2 相変化
    !
    ! CO2 phase change
    !

    ! モジュール引用 ; USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure

    ! 主成分相変化
    ! Phase change of atmospheric major component
    !
    use saturate_major_comp, only : SaturateMajorCompCondTemp, SaturateMajorCompInqLatentHeat


    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in   ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in   ):: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
                              ! $ p $ . 気圧 (整数レベル). 
                              ! Air pressure (full level)
    real(DP), intent(in   ):: xyz_Height(0:imax-1, 1:jmax, 1:kmax)
                              ! 
                              ! 
    real(DP), intent(inout):: xy_Ps            (0:imax-1, 1:jmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xyz_Temp         (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(inout):: xy_SurfMajCompIce(0:imax-1, 1:jmax)
                              !
                              ! Surface major component ice amount

    ! 作業変数
    ! Work variables
    !
    real(DP):: xyz_TempB           (0:imax-1, 1:jmax, 1:kmax)
                              ! 調節前の温度. 
                              ! Temperature before adjustment
    real(DP):: xyz_DelAtmMass      (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Atmospheric mass in a layer
    real(DP):: xy_FallingIce       (0:imax-1, 1:jmax)
                              !
                              !
    real(DP):: xyz_DelMajCompIce   (0:imax-1, 1:jmax, 1:kmax)
                              !
                              !
    real(DP):: xy_DelSurfMajCompIce(0:imax-1, 1:jmax)
                              !
                              !
    real(DP):: xyz_DTempDt         (0:imax-1, 1:jmax, 1:kmax)
                              ! 温度変化率. 
                              ! Temperature tendency
    real(DP):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
                              ! 
                              ! Surface major component ice tendency

    real(DP):: xyz_TempCond   (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xy_SurfTempCond(0:imax-1, 1:jmax)
    real(DP):: SpecHeatCO2Ice

    real(DP):: LatentHeatMajCompSubl

    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

    logical :: FlagCheckPs


    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. major_comp_phase_change_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    if ( .not. FlagMajCompPhaseChange ) return


    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! Set latent heat
    LatentHeatMajCompSubl = SaturateMajorCompInqLatentHeat()


    FlagCheckPs = .false.
    do j = 1, jmax
      do i = 0, imax-1
        if ( xyr_Press(i,j,0) > 1.0e4_DP ) then
          FlagCheckPs = .true.
        end if
      end do
    end do
    if ( FlagCheckPs ) then
      call MessageNotify( 'W', module_name, 'Surface pressure is greater than 10000 Pa.' )
    end if


    ! 調節前 "Temp" の保存
    ! Store "Temp" before adjustment
    !
    xyz_TempB    = xyz_Temp


    call SaturateMajorCompCondTemp( xyz_Press, xyz_TempCond )


    do k = 1, kmax
      xyz_DelAtmMass(:,:,k) = ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) ) / Grav
    end do


    k = kmax
    do j = 1, jmax
      do i = 0, imax-1
        if ( xyz_Temp(i,j,k) < xyz_TempCond(i,j,k) ) then
          xyz_DelMajCompIce(i,j,k) = CpDry * ( xyz_TempCond(i,j,k) - xyz_Temp(i,j,k) ) * xyz_DelAtmMass(i,j,k) / LatentHeatMajCompSubl
          xyz_Temp(i,j,k) = xyz_TempCond(i,j,k)
        else
          xyz_DelMajCompIce(i,j,k) = 0.0_DP
        end if
      end do
    end do
    !
    xy_FallingIce = 0.0_DP
    do k = kmax-1, 1, -1
      xy_FallingIce = xy_FallingIce + xyz_DelMajCompIce(:,:,k+1)
      do j = 1, jmax
        do i = 0, imax-1
          SpecHeatCO2Ice = 349.0_DP + 4.8_DP * xyz_TempCond(i,j,k)
          !                                            Forget et al. (1998)
          xyz_DelMajCompIce(i,j,k) = CpDry * ( xyz_TempCond(i,j,k) - xyz_Temp(i,j,k) ) * xyz_DelAtmMass(i,j,k) / LatentHeatMajCompSubl - (   Grav * ( xyz_Height(i,j,k+1) - xyz_Height(i,j,k) ) + SpecHeatCO2Ice * ( xyz_TempCond(i,j,k+1) - xyz_TempCond(i,j,k) ) ) / LatentHeatMajCompSubl * xy_FallingIce(i,j)
          if ( ( xy_FallingIce(i,j) + xyz_DelMajCompIce(i,j,k) ) >= 0.0_DP ) then
            xyz_Temp(i,j,k) = xyz_TempCond(i,j,k)
          else
            xyz_DelMajCompIce(i,j,k) = - xy_FallingIce(i,j)
            xyz_Temp(i,j,k) = xyz_Temp(i,j,k) + ( - LatentHeatMajCompSubl + Grav * ( xyz_Height(i,j,k+1) - xyz_Height(i,j,k) ) + SpecHeatCO2Ice * ( xyz_TempCond(i,j,k+1) - xyz_TempCond(i,j,k) ) ) / ( CpDry * xyz_DelAtmMass(i,j,k) ) * xy_FallingIce(i,j)
          end if
        end do
      end do
    end do


    ! Ice falling on the surface
    !   This may result in supersaturation in the lowest level.
    !
    xy_FallingIce = xy_FallingIce + xyz_DelMajCompIce(:,:,1)
    k = 1
    do j = 1, jmax
      do i = 0, imax-1
        xy_DelSurfMajCompIce(i,j) = - (   Grav * ( xyz_Height(i,j,1) - 0.0_DP ) ) / LatentHeatMajCompSubl * xy_FallingIce(i,j)



          SpecHeatCO2Ice = 349.0_DP + 4.8_DP * xy_SurfTempCond(i,j)
          !                                            Forget et al. (1998)
          xy_DelSurfMajCompIce(i,j) = - (   Grav * ( xyz_Height(i,j,1) - 0.0_DP ) + SpecHeatCO2Ice * ( xyz_TempCond(i,j,1) - xy_SurfTempCond(i,j) ) ) / LatentHeatMajCompSubl * xy_FallingIce(i,j)


        if ( ( xy_FallingIce(i,j) + xy_DelSurfMajCompIce(i,j) ) >= 0.0_DP ) then
          ! Part of ice sublimes.
          ! NOTE: In this case, temperature in the lowest layer should be 
          ! condensation temperature. So, actually, the set of temperature is 
          ! meaningless.
          xyz_Temp(i,j,k) = xyz_TempCond(i,j,k)
        else
          ! All falling ice sublimes.
          ! NOTE: The formulation below is different from that by Forget et al.
          ! (1998). The latent heat and heat by potential energy release and 
          ! heating ice is distributed in the lowest layer in this model, not
          ! to the soil.
          xy_DelSurfMajCompIce(i,j) = - xy_FallingIce(i,j)
          xyz_Temp(i,j,k) = xyz_Temp(i,j,k) + ( - LatentHeatMajCompSubl + Grav * ( xyz_Height(i,j,1) - 0.0_DP ) + SpecHeatCO2Ice * ( xyz_TempCond(i,j,1) - xy_SurfTempCond(i,j) ) ) / ( CpDry * xyz_DelAtmMass(i,j,k) ) * xy_FallingIce(i,j)
       end if
       end do
    end do


    ! 温度変化率
    ! Calculate temperature tendency
    !
    xyz_DTempDt = ( xyz_Temp - xyz_TempB ) / ( 2.0_DP * DelTime )


    !
    ! Surface major component ice adjustment
    !
    xy_DSurfMajCompIceDt = 0.0_DP
    do k = kmax, 1, -1
      xy_DSurfMajCompIceDt = xy_DSurfMajCompIceDt + xyz_DelMajCompIce(:,:,k)
    end do
    xy_DSurfMajCompIceDt = xy_DSurfMajCompIceDt + xy_DelSurfMajCompIce(i,j)
    !
    xy_SurfMajCompIce = xy_SurfMajCompIce + xy_DSurfMajCompIceDt * ( 2.0_DP * DelTime )


    !
    ! Surface pressure adjustment
    !
    xy_Ps = xy_Ps - xy_DSurfMajCompIceDt * Grav * ( 2.0_DP * DelTime )


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DTempDtMajCompPhaseChange', xyz_DTempDt )


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine MajorCompPhaseChangeInAtmTest
module_name
Constant :
module_name = ‘major_comp_phase_change :character(*), parameter
: モジュールの名称. Module name
version
Constant :
version = ’$Name: $’ // ’$Id: major_comp_phase_change.f90,v 1.3 2014/05/07 09:39:21 murashin Exp $’ :character(*), parameter
: モジュールのバージョン Module version