subroutine PhyRadFluxClose( phy_radflx, err )
    !
    ! PHYRADFLX 型の変数の終了処理を行います.
    ! なお, 与えられた *phy_radflx* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Deconstructor of "PHYRADFLX".
    ! Note that if *phy_radflx* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYRADFLX), intent(inout):: phy_radflx
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyRadFluxClose'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_radflx % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYRADFLX'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  "PHYRADFLX" の設定の消去
    !  Clear the settings for "PHYRADFLX"
    !-----------------------------------------------------------------
    deallocate( phy_radflx % y_Lat )

    deallocate( phy_radflx % xy_IncomRadSFlux )
    deallocate( phy_radflx % xy_InAngle )

    deallocate( phy_radflx % xy_SurfTemp )
    deallocate( phy_radflx % xy_SurfAlbedo )

    deallocate( phy_radflx % xy_TempSave  )
    deallocate( phy_radflx % xyr_RadLFluxSave )
    deallocate( phy_radflx % xyr_RadSFluxSave )
    deallocate( phy_radflx % xyra_DelRadLFluxSave )

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
    phy_radflx % initialized = .false.
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyRadFluxClose

  subroutine PhyRadFluxCreate( phy_radflx, imax, jmax, kmax, y_Lat, Grav, Cp, StB, xy_SurfTemp, xy_SurfAlbedo, current_time_value, current_time_unit, delta_time_value, delta_time_unit, delta_time_value_RadLong, delta_time_unit_RadLong, delta_time_value_RadShort, delta_time_unit_RadShort, nmlfile, err )
    !
    ! PHYRADFLX 型の変数の初期設定を行います.
    ! 他のサブルーチンを使用する前に必ずこのサブルーチンによって
    ! PHYRADFLX 型の変数を初期設定してください.
    !
    ! なお, 与えられた *phy_radflx* が既に初期設定されている場合,
    ! プログラムはエラーを発生させます.
    !
    ! NAMELIST を利用する場合には引数 *nmlfile* に NAMELIST ファイル名
    ! を与えてください. NAMELIST 変数群の詳細に関しては 
    ! NAMELIST#phy_radiation_flux_nml を参照してください. 
    !
    ! Constructor of "PHYRADFLX".
    ! Initialize *phy_radflx* by this subroutine, 
    ! before other procedures are used, 
    !
    ! Note that if *phy_radflx* is already initialized 
    ! by this procedure, error is occurred.
    !
    ! In order to use NAMELIST, specify a NAMELIST filename to 
    ! argument *nmlfile*. See "NAMELIST#phy_radiation_flux_nml"
    ! for details about a NAMELIST group.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_present, only: present_and_not_empty, present_and_true
    use dc_message, only: MessageNotify
    use dc_error, only: StoreError, DC_NOERR, DC_EALREADYINIT, DC_EARGLACK, DC_ENEGATIVE, DC_ENOFILEREAD
    use dc_date, only: Create
    implicit none
    type(PHYRADFLX), intent(inout):: phy_radflx
    integer, intent(in):: imax ! 経度格子点数. 
                   ! Number of grid points in longitude
    integer, intent(in):: jmax ! 緯度格子点数. 
                   ! Number of grid points in latitude
    integer, intent(in):: kmax ! 鉛直層数. 
                   ! Number of vertical level
    real(DP), intent(in):: Grav      ! $ g $ .      重力加速度.     Gravitational acceleration
    real(DP), intent(in):: Cp        ! $ C_p $ .    大気定圧比熱.   Specific heat of air at constant pressure
    real(DP), intent(in):: StB       ! $ \sigma_{SB} $ .       ステファンボルツマン定数. Stefan-Boltzmann constant
    real(DP), intent(in):: y_Lat (0:jmax-1)
                              ! 緯度. Latitude
    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 0:jmax-1)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in):: xy_SurfAlbedo (0:imax-1, 0:jmax-1)
                              ! 地表アルベド. 
                              ! Surface albedo
    real(DP), intent(in), optional:: current_time_value
                              ! 現在時刻の値. Value of current time
    character(*), intent(in), optional:: current_time_unit
                              ! 現在時刻の単位. Unit of current time
    real(DP), intent(in), optional:: delta_time_value
                              ! $ \Delta t $ . タイムステップの値. Value of time step
    character(*), intent(in), optional:: delta_time_unit
                              ! タイムステップの単位. Unit of time step
    real(DP), intent(in), optional:: delta_time_value_RadLong
                              ! 長波フラックスを計算するステップ間隔の値. 
                              ! Value of step interval of long wave flux calculation
    character(*), intent(in), optional:: delta_time_unit_RadLong
                              ! 長波フラックスを計算するステップ間隔の単位. 
                              ! Unit of step interval of long wave flux calculation
    real(DP), intent(in), optional:: delta_time_value_RadShort
                              ! 短波 (日射) フラックスを計算するステップ間隔の値. 
                              ! Value of step interval of short wave (insolation) flux calculation
    character(*), intent(in), optional:: delta_time_unit_RadShort
                              ! 短波 (日射) フラックスを計算するステップ間隔の単位. 
                              ! Unit of step interval of short wave (insolation) flux calculation
    character(*), intent(in), optional:: nmlfile
                              ! NAMELIST ファイルの名称. 
                              ! この引数に空文字以外を与えた場合, 
                              ! 指定されたファイルから 
                              ! NAMELIST 変数群を読み込みます. 
                              ! ファイルを読み込めない場合にはエラーを
                              ! 生じます.
                              !
                              ! NAMELIST 変数群の詳細に関しては 
                              ! NAMELIST#phy_radiation_flux_nml 
                              ! を参照してください. 
                              !
                              ! NAMELIST file name. 
                              ! If nonnull character is specified to
                              ! this argument, 
                              ! NAMELIST group name is loaded from the 
                              ! file. 
                              ! If the file can not be read, 
                              ! an error occurs.
                              ! 
                              ! See "NAMELIST#phy_radiation_flux_nml" 
                              ! for details about a NAMELIST group.
                              ! 
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyRadFluxCreate'
  continue
    call BeginSub( subname, version )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( phy_radflx % initialized ) then
      stat = DC_EALREADYINIT
      cause_c = 'PHYRADFLX'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  引数の正当性のチェック
    !  Validate arguments
    !-----------------------------------------------------------------
    if (imax < 1) then
      stat = DC_ENEGATIVE
      cause_c = 'imax'
      goto 999
    end if
    if (jmax < 1) then
      stat = DC_ENEGATIVE
      cause_c = 'jmax'
      goto 999
    end if
    if (kmax < 1) then
      stat = DC_ENEGATIVE
      cause_c = 'kmax'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  "PHYRADFLX" の設定
    !  Configure the settings for "PHYRADFLX"
    !-----------------------------------------------------------------
    phy_radflx % imax    = imax   
    phy_radflx % jmax    = jmax   
    phy_radflx % kmax    = kmax   
    phy_radflx % Grav    = Grav   
    phy_radflx % Cp      = Cp   
    phy_radflx % StB     = StB    

    allocate( phy_radflx % y_Lat (0:phy_radflx%jmax-1) )
    phy_radflx % y_Lat   = y_Lat    

    allocate( phy_radflx % xy_SurfTemp (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1) )
    allocate( phy_radflx % xy_SurfAlbedo (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1) )
    phy_radflx % xy_SurfTemp   = xy_SurfTemp   
    phy_radflx % xy_SurfAlbedo = xy_SurfAlbedo 

    !-----------------------------------------------------------------
    !  時刻の設定
    !  Configure time
    !-----------------------------------------------------------------
    if ( present(current_time_value) .and. present(current_time_unit) ) then
      call Create( diff = phy_radflx % current_time, value = current_time_value, unit = current_time_unit )            ! (in)
    else
      call Create( diff = phy_radflx % current_time, value = 0.0_DP, unit = 'sec' )            ! (in)
    end if

    if ( present(delta_time_value) .and. present(delta_time_unit) ) then
      call Create( diff = phy_radflx % delta_time, value = delta_time_value, unit = delta_time_unit )            ! (in)
    else
      call Create( diff = phy_radflx % delta_time, value = 1.0_DP, unit = 'sec' )             ! (in)
    end if

    if ( present(delta_time_value_RadLong) .and. present(delta_time_unit_RadLong) ) then
      call Create( diff = phy_radflx % delta_time_RadLong, value = delta_time_value_RadLong, unit = delta_time_unit_RadLong )            ! (in)
    else
      call Create( diff = phy_radflx % delta_time_RadLong, value = -1.0_DP, unit = 'sec' )             ! (in)
    end if

    if ( present(delta_time_value_RadShort) .and. present(delta_time_unit_RadShort) ) then
      call Create( diff = phy_radflx % delta_time_RadShort, value = delta_time_value_RadShort, unit = delta_time_unit_RadShort )            ! (in)
    else
      call Create( diff = phy_radflx % delta_time_RadShort, value = -1.0_DP, unit = 'sec' )             ! (in)
    end if

    !-----------------------------------------------------------------
    !  保存用変数の割付
    !  Allocate variables for save
    !-----------------------------------------------------------------
    phy_radflx % save_flag = .false.
    allocate( phy_radflx % xy_TempSave (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1) )
    allocate( phy_radflx % xyr_RadLFluxSave (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax) )
    allocate( phy_radflx % xyr_RadSFluxSave (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax) )
    allocate( phy_radflx % xyra_DelRadLFluxSave (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax, 0:1) )

    phy_radflx % xy_TempSave = 0.0_DP
    phy_radflx % xyr_RadLFluxSave = 0.0_DP
    phy_radflx % xyr_RadSFluxSave = 0.0_DP
    phy_radflx % xyra_DelRadLFluxSave = 0.0_DP

    !-----------------------------------------------------------------
    !  NAMELIST の取得
    !  Get NAMELIST
    !-----------------------------------------------------------------

    !-------------------------
    !  デフォルト値
    !  Default values
!!$    phy_radflx % param_r = 0.0_DP
!!$    phy_radflx % param_c = 'hogehoge'

    !-------------------------
    !  オプショナル引数からの値
    !  Values from optional arguments
!!$    phy_radflx % param_i = param_i
!!$    if ( present(param_r) )  phy_radflx % param_r = param_r
!!$    if ( present(param_c) )  phy_radflx % param_c = param_c

    !-------------------------
    !  NAMELIST からの値
    !  Values from NAMELIST

!!$    if ( present_and_not_empty(nmlfile) ) then
!!$      call MessageNotify( 'M', subname, &
!!$        & 'Loading NAMELIST file "%c" ...', &
!!$        & c1 = trim(nmlfile) )
!!$      call NmlRead ( nmlfile = nmlfile, &      ! (in)
!!$        & param_i = phy_radflx % param_i, &   ! (inout)
!!$        & param_r = phy_radflx % param_r, &   ! (inout)
!!$        & param_c_ = phy_radflx % param_c, &  ! (inout)
!!$        & err = err )                          ! (out)
!!$      if ( present_and_true(err) ) then
!!$        call MessageNotify( 'W', subname, &
!!$          & '"%c" can not be read.', &
!!$          & c1 = trim(nmlfile) )
!!$        stat = DC_ENOFILEREAD
!!$        cause_c = nmlfile
!!$        goto 999
!!$      end if
!!$    end if

    !-----------------------------------------------------------------
    !  設定値の正当性のチェック
    !  Validate setting values
    !-----------------------------------------------------------------
!!$    if ( phy_radflx % param_i < 0 ) then
!!$      stat = DC_ENEGATIVE
!!$      cause_c = 'param_i'
!!$      goto 999
!!$    end if


    !-----------------------------------------------------------------
    !  初期化フラグの設定
    !  Set initialization flag 
    !-----------------------------------------------------------------
    phy_radflx % initialized = .true.

    !-----------------------------------------------------------------
    !  短波入射の計算
    !  Calculate short wave (insolation) incoming radiation
    !-----------------------------------------------------------------
    allocate( phy_radflx % xy_IncomRadSFlux (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1) )
    allocate( phy_radflx % xy_InAngle (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1) )

    call ShortIncoming( phy_radflx = phy_radflx, xy_IncomRadSFlux = phy_radflx % xy_IncomRadSFlux, xy_InAngle = phy_radflx % xy_InAngle )              ! (out)

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyRadFluxCreate

  logical function PhyRadFluxInitialized( phy_radflx ) result(result)
    !
    ! *phy_radflx* が初期設定されている場合には .true. が,
    ! 初期設定されていない場合には .false. が返ります.
    !
    ! If *phy_radflx* is initialized, .true. is returned.
    ! If *phy_radflx* is not initialized, .false. is returned.
    !
    implicit none
    type(PHYRADFLX), intent(in):: phy_radflx
  continue
    result = phy_radflx % initialized
  end function PhyRadFluxInitialized

  subroutine PhyRadFluxLongFlux( phy_radflx, xyr_RadLFlux, xyra_DelRadLFlux, xyz_Temp, xy_SurfTemp, xyr_TauQVap, xyr_TauDryAir, err )
    !
    ! 長波フラックスの計算
    !
    ! なお, 与えられた *phy_radflx* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Calculate long wave flux
    !
    ! If *phy_radflx* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYRADFLX), intent(inout):: phy_radflx
    real(DP), intent(out):: xyr_RadLFlux (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 長波フラックス. 
                              ! Long wave flux
    real(DP), intent(out):: xyra_DelRadLFlux (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with long wave
    real(DP), intent(in):: xyz_Temp (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax-1)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(in):: xy_SurfTemp (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in):: xyr_TauQVap (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 光学的厚さ：水
                              ! Optical depth of water
    real(DP), intent(in):: xyr_TauDryAir (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 光学的厚さ：空気
                              ! Optical depth of air
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 
    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: kmax ! 鉛直層数. 
                   ! Number of vertical level

    real(DP):: StB       ! $ \sigma_{SB} $ .       ステファンボルツマン定数. Stefan-Boltzmann constant

    real(DP):: xyr_Trans (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 透過係数. 
                              ! Transmission coefficient
    real(DP):: xyr_Trans1 (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 1/2 レベルからの透過係数. 
                              ! Transmission coefficient above 1/2 level
    real(DP):: xyr_Trans2 (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 3/2 レベルからの透過係数. 
                              ! Transmission coefficient above 3/2 level
    real(DP):: xyz_PiB (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax-1)
                              ! $ \pi B = \sigma T^{4} $
    real(DP):: xy_SurfPiB (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1)
                              ! 地表面の $ \pi B $ . 
                              ! $ \pi B $ on surface

    integer:: kk , bn
    real(DP):: BandWeightSum

    integer, parameter :: BandNumber = 4
                              ! 長波バンド数. 
                              ! Number of long wave band

    real(DP):: AbsorpCoeffQVap(BandNumber)
                              ! 水の吸収係数. 
                              ! Absorption coefficient of water
    real(DP):: AbsorpCoeffDryAir(BandNumber)
                              ! 空気の吸収係数. 
                              ! Absorption coefficient of air
    real(DP):: BandWeight(BandNumber)
                              ! バンドウェイト. 
                              ! Band weight
    real(DP):: PathLengthFact
                              ! 光路長のファクタ. 
                              ! Factor of optical length

    integer:: k               ! DO ループ用作業変数
                              ! Work variables for DO loop

    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyRadFluxLongFlux'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_radflx % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYRADFLX'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  *phy_radflx* に格納されている設定値の取り出し
    !  Fetch setting values stored in *phy_radflx*
    !-----------------------------------------------------------------
    kmax    = phy_radflx % kmax   
    StB     = phy_radflx % StB    

    !-----------------------------------------------------------------
    !  吸収係数などの設定
    !  Configure absorption coefficients etc.
    !-----------------------------------------------------------------
    AbsorpCoeffQVap = (/ 8.0_DP, 1.0_DP, 0.1_DP, 0.0_DP /)
    AbsorpCoeffDryAir = (/ 0.0_DP, 0.0_DP, 0.0_DP, 5.0e-5_DP /)
    BandWeight = (/ 0.2_DP, 0.1_DP, 0.1_DP, 0.6_DP /)
    PathLengthFact = 1.5_DP


    !-----------------------------------------------------------------
    !  バンドウェイトの設定
    !  Configure band weight
    !-----------------------------------------------------------------
    BandWeightSum = 0.0_DP

    do bn = 1, BandNumber
      BandWeightSum =  BandWeightSum +  BandWeight(bn)
    end do

    do bn = 1, BandNumber
      BandWeight(bn) = BandWeight(bn) / BandWeightSum
    end do

    !-----------------------------------------------------------------
    !  $ \pi B $ の計算
    !  Calculate $ \pi B $
    !-----------------------------------------------------------------
    xyz_PiB    = StB * ( xyz_Temp**4 )
    xy_SurfPiB = StB * ( xy_SurfTemp**4 )

    do k = 0, kmax

      !---------------------------------------------------------------
      !  透過関数計算
      !  Calculate transmission functions
      !---------------------------------------------------------------
      xyr_Trans = 0.0_DP
      
      do bn = 1, BandNumber
        do kk = 0, kmax
          xyr_Trans(:,:,kk) = xyr_Trans(:,:,kk) + BandWeight(bn) * exp( - PathLengthFact * (   AbsorpCoeffQVap(bn) * abs(   xyr_TauQVap(:,:,kk) - xyr_TauQVap(:,:,k)  ) + AbsorpCoeffDryAir(bn) * abs(   xyr_TauDryAir(:,:,kk) - xyr_TauDryAir(:,:,k)  ) ) )
        end do
      end do
      

      !---------------------------------------------------------------
      !  放射フラックス計算
      !  Calculate radiation flux
      !---------------------------------------------------------------
      xyr_RadLFlux(:,:,k) = xy_SurfPiB * xyr_Trans(:,:,0)
 
      do kk = 0, kmax-1
        xyr_RadLFlux(:,:,k) = xyr_RadLFlux(:,:,k) - xyz_PiB(:,:,kk) * ( xyr_Trans(:,:,kk) - xyr_Trans(:,:,kk+1) )
      end do

      !---------------------------------------------------------------
      !  補正項計算用透過関数計算
      !  Calculate transmission functions for correction terms
      !---------------------------------------------------------------
      xyr_Trans1(:,:,k) = xyr_Trans(:,:,0)
      xyr_Trans2(:,:,k) = xyr_Trans(:,:,1)

    end do

    !-----------------------------------------------------------------
    !  長波地表温度変化の計算
    !  Calclate surface temperature tendency with long wave
    !-----------------------------------------------------------------
    do k = 0, kmax
      xyra_DelRadLFlux(:,:,k,0) = 4.0_DP * xy_SurfPiB / xy_SurfTemp * xyr_Trans1(:,:,k)

      xyra_DelRadLFlux(:,:,k,1) = 4.0_DP * xyz_PiB(:,:,0) / xyz_Temp(:,:,0) * ( xyr_Trans2(:,:,k) - xyr_Trans1(:,:,k) )
    end do

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyRadFluxLongFlux

  subroutine PhyRadFluxNmlRead( nmlfile, err )
    !
    ! NAMELIST ファイル *nmlfile* から値を入力するための
    ! 内部サブルーチンです. Create 内で呼び出されることを
    ! 想定しています.
    !
    ! 値が NAMELIST ファイル内で指定されていない場合には,
    ! 入力された値がそのまま返ります.
    !
    ! なお, *nmlfile* に空文字が与えられた場合, または
    ! 与えられた *nmlfile* を読み込むことができない場合, 
    ! プログラムはエラーを発生させます.
    !
    ! This is an internal subroutine to input values from 
    ! NAMELIST file *nmlfile*. This subroutine is expected to be
    ! called by "Create".
    !
    ! A value not specified in NAMELIST file is returned
    ! without change.
    !
    ! If *nmlfile* is empty, or *nmlfile* can not be read, 
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_iounit, only: FileOpen
    use dc_message, only: MessageNotify
    use dc_present, only: present_and_true
    use dc_error, only: StoreError, DC_NOERR, DC_ENOFILEREAD
    implicit none
    character(*), intent(in):: nmlfile
                              ! NAMELIST ファイルの名称. 
                              ! NAMELIST file name
!!$    integer, intent(inout):: param_i
!!$    real(DP), intent(inout):: param_r
!!$    character(*), intent(inout):: param_c_
!!$    character(TOKEN):: param_c
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

!!$    namelist /phy_radiation_flux_nml/ &
!!$      & param_i, param_r, param_c
                              ! phy_radiation_flux モジュール用
                              ! NAMELIST 変数群名.
                              !
                              ! phy_radiation_flux#Create を使用する際に, 
                              ! オプショナル引数 *nmlfile* へ NAMELIST 
                              ! ファイル名を指定することで, そのファイルから
                              ! この NAMELIST 変数群を読み込みます.
                              !
                              ! NAMELIST group name for 
                              ! "phy_radiation_flux" module.
                              ! 
                              ! If a NAMELIST filename is specified to 
                              ! an optional argument *nmlfile* 
                              ! when "phy_radiation_flux#Create" is used, 
                              ! this NAMELIST group is loaded from 
                              ! the file.

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
!!$    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
!!$                              ! IOSTAT of NAMELIST read
    character(*), parameter:: subname = 'PhyRadFluxNmlRead'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''



    !-----------------------------------------------------------------
    !  文字型引数を NAMELIST 変数群へ代入
    !  Substitute character arguments to NAMELIST group
    !-----------------------------------------------------------------
!!$    param_c = param_c_

    !----------------------------------------------------------------
    !  NAMELIST ファイルのオープン
    !  Open NAMELIST file
    !----------------------------------------------------------------
    call FileOpen( unit = unit_nml, file = nmlfile, mode = 'r', err = err )                   ! (out)
    if ( present_and_true(err) ) then
      stat = DC_ENOFILEREAD
      cause_c = nmlfile
      goto 999
    end if


    !-----------------------------------------------------------------
    !  NAMELIST 変数群の取得
    !  Get NAMELIST group
    !-----------------------------------------------------------------
!!$    read( unit = unit_nml, & ! (in)
!!$      & nml = phy_radiation_flux_nml, iostat = iostat_nml ) ! (out)
!!$    if ( iostat_nml == 0 ) then
!!$      call MessageNotify( 'M', subname, &
!!$        & 'NAMELIST group "%c" is loaded from "%c".', &
!!$        & c1 = 'phy_radiation_flux_nml', c2 = trim(nmlfile) )
!!$      write(STDOUT, nml = phy_radiation_flux_nml)
!!$    else
!!$      call MessageNotify( 'W', subname, &
!!$        & 'NAMELIST group "%c" is not found in "%c" (iostat=%d).', &
!!$        & c1 = 'phy_radiation_flux_nml', c2 = trim(nmlfile), &
!!$        & i = (/iostat_nml/) )
!!$    end if

    close( unit_nml )

    !-----------------------------------------------------------------
    !  NAMELIST 変数群を文字型引数へ代入
    !  Substitute NAMELIST group to character arguments
    !-----------------------------------------------------------------
!!$    param_c_ = param_c

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyRadFluxNmlRead

  subroutine PhyRadFluxPutLine( phy_radflx, unit, indent, err )
    !
    ! 引数 *phy_radflx* に設定されている情報を印字します.
    ! デフォルトではメッセージは標準出力に出力されます. 
    ! *unit* に装置番号を指定することで, 出力先を変更することが可能です.
    !
    ! Print information of *phy_radflx*.
    ! By default messages are output to standard output.
    ! Unit number for output can be changed by *unit* argument.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    use dc_date, only: toChar
    implicit none
    type(PHYRADFLX), intent(in):: phy_radflx
    integer, intent(in), optional:: unit
                              ! 出力先の装置番号.
                              ! デフォルトの出力先は標準出力.
                              !
                              ! Unit number for output.
                              ! Default value is standard output.
    character(*), intent(in), optional:: indent
                              ! 表示されるメッセージの字下げ.
                              !
                              ! Indent of displayed messages.
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    integer:: out_unit
    integer:: indent_len
    character(STRING):: indent_str
    character(*), parameter:: subname = 'PhyRadFluxPutLine'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( present(unit) ) then
      out_unit = unit
    else
      out_unit = STDOUT
    end if

    indent_len = 0
    indent_str = ''
    if ( present(indent) ) then
      if ( len(indent) /= 0 ) then
        indent_len = len(indent)
        indent_str(1:indent_len) = indent
      end if
    end if


    !-----------------------------------------------------------------
    !  "PHYRADFLX" の設定の印字
    !  Print the settings for "PHYRADFLX"
    !-----------------------------------------------------------------
    if ( phy_radflx % initialized ) then
      call Printf( out_unit, indent_str(1:indent_len) // '#<PHYRADFLX:: @initialized=%y', l = (/phy_radflx % initialized/) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @imax=%d @jmax=%d @kmax=%d', i = (/phy_radflx % imax, phy_radflx % jmax, phy_radflx % kmax/) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @Grav=%f @Cp=%f @StB=%f', d = (/ phy_radflx % Grav, phy_radflx % Cp, phy_radflx % StB/) )

      call PutLine( phy_radflx % y_Lat, unit = out_unit, lbounds = lbound(phy_radflx % y_Lat), ubounds = ubound(phy_radflx % y_Lat), indent = indent_str(1:indent_len) // ' @y_Lat=' )

      call PutLine( phy_radflx % xy_SurfTemp, unit = out_unit, lbounds = lbound(phy_radflx % xy_SurfTemp), ubounds = ubound(phy_radflx % xy_SurfTemp), indent = indent_str(1:indent_len) // ' @xy_SurfTemp=' )

      call PutLine( phy_radflx % xy_SurfAlbedo, unit = out_unit, lbounds = lbound(phy_radflx % xy_SurfAlbedo), ubounds = ubound(phy_radflx % xy_SurfAlbedo), indent = indent_str(1:indent_len) // ' @xy_SurfAlbedo=' )

      call PutLine( phy_radflx % xy_IncomRadSFlux, unit = out_unit, lbounds = lbound(phy_radflx % xy_IncomRadSFlux), ubounds = ubound(phy_radflx % xy_IncomRadSFlux), indent = indent_str(1:indent_len) // ' @xy_IncomRadSFlux=' )

      call PutLine( phy_radflx % xy_InAngle, unit = out_unit, lbounds = lbound(phy_radflx % xy_InAngle), ubounds = ubound(phy_radflx % xy_InAngle), indent = indent_str(1:indent_len) // ' @xy_InAngle=' )

      call Printf( out_unit, indent_str(1:indent_len) // ' @current_time=%c', c1 = trim( toChar(phy_radflx % current_time) ) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @delta_time=%c', c1 = trim( toChar(phy_radflx % delta_time) ) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @delta_time_RadLong=%c', c1 = trim( toChar(phy_radflx % delta_time_RadLong) ) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @delta_time_RadShort=%c', c1 = trim( toChar(phy_radflx % delta_time_RadShort) ) )

      call Printf( out_unit, indent_str(1:indent_len) // '>' )
    else
      call Printf( out_unit, indent_str(1:indent_len) // '#<PHYRADFLX:: @initialized=%y>', l = (/phy_radflx % initialized/) )
    end if

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyRadFluxPutLine

  subroutine PhyRadFluxRadiationDelTemp( phy_radflx, xyr_RadLFlux, xyr_RadSFlux, xyr_Press, xyz_DRadLTempDt, xyz_DRadSTempDt, err )
    !
    ! 放射による温度変化率を計算します. 
    !
    ! なお, 与えられた *phy_radflx* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Calculate temperature tendency with radiation.
    !
    ! If *phy_radflx* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYRADFLX), intent(inout):: phy_radflx
    real(DP), intent(in):: xyr_RadLFlux (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 長波フラックス. 
                              ! Long wave flux
    real(DP), intent(in):: xyr_RadSFlux (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 短波 (日射) フラックス. 
                              ! Short wave (insolation) flux
    real(DP), intent(in):: xyr_Press (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! $ P_s $ . 地表面気圧 (半整数レベル). 
                              ! Surface pressure (half level)
    real(DP), intent(out):: xyz_DRadLTempDt (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax-1)
                              ! 長波加熱率. 
                              ! Temperature tendency with long wave
    real(DP), intent(out):: xyz_DRadSTempDt (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax-1)
                              ! 短波加熱率. 
                              ! Temperature tendency with short wave

    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: kmax ! 鉛直層数. 
                   ! Number of vertical level

    real(DP):: Grav      ! $ g $ .      重力加速度.     Gravitational acceleration
    real(DP):: Cp        ! $ C_p $ .    大気定圧比熱.   Specific heat of air at constant pressure

    integer:: k               ! DO ループ用作業変数
                              ! Work variables for DO loop

    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyRadFluxRadiationDelTemp'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_radflx % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYRADFLX'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  *phy_radflx* に格納されている設定値の取り出し
    !  Fetch setting values stored in *phy_radflx*
    !-----------------------------------------------------------------
    kmax        = phy_radflx % kmax   
    Grav        = phy_radflx % Grav   
    Cp          = phy_radflx % Cp   

    !----------------------------------------------------------------
    !  放射冷却率の演算
    !  Calculate radiation cooling rate
    !----------------------------------------------------------------
    do k = 0, kmax - 1
      xyz_DRadLTempDt(:,:,k) = ( xyr_RadLFlux(:,:,k) - xyr_RadLFlux(:,:,k+1) ) / ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) / Cp * Grav
      xyz_DRadSTempDt(:,:,k) = ( xyr_RadSFlux(:,:,k) - xyr_RadSFlux(:,:,k+1) ) / ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) / Cp * Grav
    end do

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyRadFluxRadiationDelTemp

  subroutine PhyRadFluxRadiationFlux( phy_radflx, xyz_Temp, xyz_QVap, xyr_Press, xyr_RadLFlux, xya_SurfRadLMatrix, xyra_DelRadLFlux, xyr_RadSFlux, err )
    !
    ! 温度, 比湿, 気圧から, 放射フラックスを計算します. 
    !
    ! なお, 与えられた *phy_radflx* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Calculate radiation flux from temperature, specific humidity, and 
    ! air pressure. 
    !
    ! If *phy_radflx* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub, DbgMessage
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    use dc_date, only: EvalSec, EvalByUnit, mod, operator(*), operator(==), operator(<), operator(/), operator(+), operator(-)
    implicit none
    type(PHYRADFLX), intent(inout):: phy_radflx
    real(DP), intent(in):: xyz_Temp (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax-1)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(in):: xyz_QVap (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax-1)
                              ! $ q $ .     比湿. Specific humidity
    real(DP), intent(in):: xyr_Press (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! $ P_s $ . 地表面気圧 (半整数レベル). 
                              ! Surface pressure (half level)
    real(DP), intent(out):: xyr_RadLFlux (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 長波フラックス. 
                              ! Long wave flux
    real(DP), intent(out):: xya_SurfRadLMatrix (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, -1:1)
                              ! $ T $ 陰解行列: 地表. 
                              ! $ T $ implicit matrix: surface
    real(DP), intent(out):: xyra_DelRadLFlux (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with long wave
    real(DP), intent(out):: xyr_RadSFlux (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 短波 (日射) フラックス. 
                              ! Short wave (insolation) flux

    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: kmax ! 鉛直層数. 
                   ! Number of vertical level
    real(DP):: Grav      ! $ g $ .      重力加速度.     Gravitational acceleration
    real(DP):: xy_SurfTemp (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1)
                              ! 地表面温度. 
                              ! Surface temperature

    real(DP):: xy_TempSave (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1)
                              ! $ T $ .     温度 (保存用). Temperature (for save)
    real(DP):: xyr_RadLFluxSave (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 長波フラックス (保存用). 
                              ! Long wave flux (for save)
    real(DP):: xyr_RadSFluxSave (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 短波 (日射) フラックス (保存用). 
                              ! Short wave (insolation) flux (for save)
    real(DP):: xyra_DelRadLFluxSave (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax, 0:1)
                              ! 長波地表温度変化 (保存用). 
                              ! Surface temperature tendency with long wave (for save)

    real(DP):: xyr_TauQVap (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 光学的厚さ：水
                              ! Optical depth of water
    real(DP):: xyr_TauDryAir (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 光学的厚さ：空気
                              ! Optical depth of air

    integer:: k               ! DO ループ用作業変数
                              ! Work variables for DO loop


    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyRadFluxRadiationFlux'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_radflx % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYRADFLX'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  *phy_radflx* に格納されている設定値の取り出し
    !  Fetch setting values stored in *phy_radflx*
    !-----------------------------------------------------------------
    kmax    = phy_radflx % kmax   
    Grav    = phy_radflx % Grav   

    xy_SurfTemp   = phy_radflx % xy_SurfTemp   

    xyr_RadSFlux = 0.0_DP
    xyr_RadSFlux(:,:,kmax) = phy_radflx % xy_IncomRadSFlux

    xy_TempSave          = phy_radflx % xy_TempSave         
    xyr_RadLFluxSave     = phy_radflx % xyr_RadLFluxSave     
    xyr_RadSFluxSave     = phy_radflx % xyr_RadSFluxSave     
    xyra_DelRadLFluxSave = phy_radflx % xyra_DelRadLFluxSave 

    !-----------------------------------------------------------------
    !  光学的厚さの計算
    !  Calculate optical depth
    !-----------------------------------------------------------------
    xyr_TauQVap   = 0.0_DP
    xyr_TauDryAir = 0.0_DP

    do k = kmax-1 , 0, -1
      xyr_TauQVap(:,:,k) = xyr_TauQVap(:,:,k+1) + xyz_QVap(:,:,k) * ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) / Grav

      xyr_TauDryAir(:,:,k) = xyr_TauDryAir(:,:,k+1) + ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) / Grav
    end do


    !-----------------------------------------------------------------
    !  長波フラックスの算出
    !  Calculate long wave flux
    !-----------------------------------------------------------------

    !---------------------------------------------
    !  current_time / delta_time_RadLong == 0, または
    !  delta_time_RadLong 未設定, または
    !  まだ一度も計算していない場合には計算を行う. 
    !
    !  If current_time / delta_time_RadLong == 0, or
    !  "delta_time_RadLong" is not configured, or
    !  calculation was never performed, 
    !  calculation is performed. 
    if ( mod( phy_radflx % current_time, phy_radflx % delta_time_RadLong) == 0 .or. EvalSec( phy_radflx % delta_time_RadLong ) < 0.0_DP .or. .not. phy_radflx % save_flag                               ) then

      call DbgMessage( 'LongFlux is calculated' )

      call LongFlux( phy_radflx = phy_radflx, xyr_RadLFlux = xyr_RadLFlux, xyra_DelRadLFlux = xyra_DelRadLFlux, xyz_Temp = xyz_Temp, xy_SurfTemp = xy_SurfTemp, xyr_TauQVap = xyr_TauQVap, xyr_TauDryAir = xyr_TauDryAir )        ! (in)
       
    !-----------------------------------
    !  前回の値を利用
    !  Use values in last time
    else

      call DbgMessage( 'LongFlux is not calculated. Save values are used' )

      xyr_RadLFlux = xyr_RadLFluxSave
      xyra_DelRadLFlux = xyra_DelRadLFluxSave
       
      do k = 0, kmax
        xyr_RadLFlux(:,:,k) = xyr_RadLFlux(:,:,k) + xyra_DelRadLFlux(:,:,k,1) * ( xyz_Temp(:,:,0) - xy_TempSave )

        xyra_DelRadLFlux(:,:,k,1) = xyra_DelRadLFlux(:,:,k,1) / ( xy_TempSave**3 ) * ( xyz_Temp(:,:,0)**3 )
      end do
    end if


    !-----------------------------------------------------------------
    !  長波陰解用行列の算出
    !  Calculate long wave implicit matrix
    !-----------------------------------------------------------------
    xya_SurfRadLMatrix(:,:,0)  = xyra_DelRadLFlux(:,:,0,0)
    xya_SurfRadLMatrix(:,:,1)  = xyra_DelRadLFlux(:,:,0,1)
    xya_SurfRadLMatrix(:,:,-1) = 0.0_DP


    !-----------------------------------------------------------------
    !  短波 (日射) フラックスの算出
    !  Calculate short wave (insolation)
    !-----------------------------------------------------------------

    !---------------------------------------------
    !  current_time / delta_time_RadShort == 0, または
    !  delta_time_RadShort 未設定, または
    !  まだ一度も計算していない場合には計算を行う. 
    !
    !  If current_time / delta_time_RadShort == 0, or
    !  "delta_time_RadShort" is not configured, or
    !  calculation was never performed, 
    !  calculation is performed. 
    if ( mod( phy_radflx % current_time, phy_radflx % delta_time_RadShort) == 0 .or. EvalSec( phy_radflx % delta_time_RadShort ) < 0.0_DP .or. .not. phy_radflx % save_flag                               ) then

      call DbgMessage( 'ShortFlux is calculated' )

      !---------------------------------
      !  短波フラックスの計算
      !  Calculate short wave (insolation) flux
      call ShortFlux( phy_radflx = phy_radflx, xyr_RadSFlux = xyr_RadSFlux, xyr_TauQVap = xyr_TauQVap, xyr_TauDryAir = xyr_TauDryAir )        ! (in)


    !-----------------------------------
    !  前回の値を利用
    !  Use values in last time
    else

      call DbgMessage( 'ShortFlux is not calculated. Save values are used' )

      xyr_RadSFlux = xyr_RadSFluxSave
    end if


    !-----------------------------------------------------------------
    !  今回計算した値を保存
    !  Save calculated values in this time
    !-----------------------------------------------------------------
    phy_radflx % xy_TempSave          = xyz_Temp (:,:,0)
    phy_radflx % xyr_RadLFluxSave     = xyr_RadLFlux
    phy_radflx % xyr_RadSFluxSave     = xyr_RadSFlux
    phy_radflx % xyra_DelRadLFluxSave = xyra_DelRadLFlux

    phy_radflx % save_flag = .true.

    !-----------------------------------------------------------------
    !  時刻の更新
    !  Update time
    !-----------------------------------------------------------------
    phy_radflx % current_time = phy_radflx % current_time + phy_radflx % delta_time

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyRadFluxRadiationFlux

  subroutine PhyRadFluxShortFlux( phy_radflx, xyr_RadSFlux, xyr_TauQVap, xyr_TauDryAir, err )
    !
    ! 長波フラックスの計算
    !
    ! なお, 与えられた *phy_radflx* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Calculate long wave flux
    !
    ! If *phy_radflx* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYRADFLX), intent(inout):: phy_radflx
    real(DP), intent(inout):: xyr_RadSFlux (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 短波 (日射) フラックス. 
                              ! Short wave (insolation) flux

    real(DP), intent(in):: xyr_TauQVap (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 光学的厚さ：水
                              ! Optical depth of water
    real(DP), intent(in):: xyr_TauDryAir (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1, 0:phy_radflx%kmax)
                              ! 光学的厚さ：空気
                              ! Optical depth of air
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 
    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: kmax ! 鉛直層数. 
                   ! Number of vertical level
    real(DP):: xy_InAngle (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1)
                              ! sec (入射角). 
                              ! sec (angle of incidence)

    real(DP):: xy_SurfAlbedo (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1)
                              ! 地表アルベド. 
                              ! Surface albedo

    integer:: bn
    real(DP):: BandWeightSum

    integer, parameter :: BandNumber = 1
                              ! 短波バンド数. 
                              ! Number of short wave band

    real(DP):: AbsorpCoeffQVap(BandNumber)
                              ! 水の吸収係数. 
                              ! Absorption coefficient of water
    real(DP):: AbsorpCoeffDryAir(BandNumber)
                              ! 空気の吸収係数. 
                              ! Absorption coefficient of air
    real(DP):: BandWeight(BandNumber)
                              ! バンドウェイト. 
                              ! Band weight
    real(DP):: SECSCT
                              ! 散乱の $ sec \zeta $
                              ! $ sec \zeta $ of scattering

    integer:: k               ! DO ループ用作業変数
                              ! Work variables for DO loop

    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyRadFluxShortFlux'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_radflx % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYRADFLX'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  *phy_radflx* に格納されている設定値の取り出し
    !  Fetch setting values stored in *phy_radflx*
    !-----------------------------------------------------------------
    kmax          = phy_radflx % kmax   
    xy_InAngle    = phy_radflx % xy_InAngle
    xy_SurfAlbedo = phy_radflx % xy_SurfAlbedo

    !-----------------------------------------------------------------
    !  吸収係数などの設定
    !  Configure absorption coefficients etc.
    !-----------------------------------------------------------------
    AbsorpCoeffQVap = (/ 0.002_DP /)
    AbsorpCoeffDryAir = (/ 0.0_DP /)
    BandWeight = (/ 1.0_DP /)
    SECSCT = 1.66_DP

    !-----------------------------------------------------------------
    !  バンドウェイトの設定
    !  Configure band weight
    !-----------------------------------------------------------------
    BandWeightSum = 0.0_DP

    do bn = 1, BandNumber
      BandWeightSum =  BandWeightSum +  BandWeight(bn)
    end do

    do bn = 1, BandNumber
      BandWeight(bn) = BandWeight(bn) / BandWeightSum
    end do

    do bn = 1, BandNumber
      do k = 0, kmax

        !-------------------------------------------------------------
        !  各レベルでの下向き透過
        !  Downward transmission on each level
        !-------------------------------------------------------------
        if ( k /= kmax ) then
          xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) + BandWeight(bn) * xyr_RadSFlux(:,:,kmax) * exp( - xy_InAngle(:,:) * (   AbsorpCoeffQVap(bn) * xyr_TauQVap(:,:,k) + AbsorpCoeffDryAir(bn) * xyr_TauDryAir(:,:,k) ) )
        end if
        
        !-------------------------------------------------------------
        !  各レベルでの上向き透過
        !  Upward transmission on each level
        !-------------------------------------------------------------
        xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) - BandWeight(bn) * xyr_RadSFlux(:,:,kmax) * exp( - xy_InAngle(:,:) * (   AbsorpCoeffQVap(bn) * xyr_TauQVap(:,:,0) + AbsorpCoeffDryAir(bn) * xyr_TauDryAir(:,:,0) ) ) * xy_SurfAlbedo * exp( - SECSCT * (   AbsorpCoeffQVap(bn) * ( xyr_TauQVap(:,:,0) - xyr_TauQVap(:,:,k) ) + AbsorpCoeffDryAir(bn) * ( xyr_TauDryAir(:,:,0) - xyr_TauDryAir(:,:,k) ) ) )
      end do
    end do

    !-----------------------------------------------------------------
    !  吸収なしの場合
    !  In the case of no absorption
    !-----------------------------------------------------------------
    if ( BandNumber == 0 ) then
      do k = 0, kmax
        xyr_RadSFlux(:,:,k) = ( 1.0_DP - xy_SurfAlbedo ) * xyr_RadSFlux(:,:,kmax)
      end do
    end if

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyRadFluxShortFlux

  subroutine PhyRadFluxShortIncoming( phy_radflx, xy_IncomRadSFlux, xy_InAngle, err )
    !
    ! 短波入射を計算します.
    !
    ! なお, 与えられた *phy_radflx* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Calculate short wave (insolation) incoming radiation. 
    !
    ! If *phy_radflx* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYRADFLX), intent(inout):: phy_radflx
    real(DP), intent(out):: xy_IncomRadSFlux (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1)
                              ! 短波 (日射) フラックス. 
                              ! Short wave (insolation) flux
    real(DP), intent(out):: xy_InAngle (0:phy_radflx%imax-1, 0:phy_radflx%jmax-1)
                              ! sec (入射角). 
                              ! sec (angle of incidence)
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: imax ! 経度格子点数. 
                   ! Number of grid points in longitude
    integer:: jmax ! 緯度格子点数. 
                   ! Number of grid points in latitude

    real(DP):: y_Lat (0:phy_radflx%jmax-1)
                              ! 緯度 (ラジアン). Latitude (radian)

    real(DP), parameter:: SolarCoeff = 1380.0_DP
                              ! 太陽定数. 
                              ! Solar constant
    real(DP), parameter:: AtmosAlbedo = 0.2_DP
                              ! 大気アルベド. 
                              ! Albedo of air

!!$    real(DP), parameter:: EPSORB = 23.5_DP
!!$                              ! 赤道傾斜角. 
!!$                              ! Inclination of equator to orbit
!!$    real(DP), parameter:: EQNORB = 110.0_DP
!!$                              ! 昇降点黄経. 
!!$                              ! 
    real(DP), parameter:: AINS   = 0.127_DP
                              ! 年平均入射の係数
                              ! Coefficient of annual mean incoming radiation
    real(DP), parameter:: BINS   = 0.183_DP
                              ! 年平均入射の係数
                              ! Coefficient of annual mean incoming radiation
    real(DP), parameter:: AZET   = 0.410_DP
                              ! 年平均入射角の係数
                              ! Coefficient of annual mean incoming radiation
    real(DP), parameter:: BZET   = 0.590_DP
                              ! 年平均入射角の係数
                              ! Coefficient of annual mean incoming radiation

    integer:: i, j            ! DO ループ用作業変数
                              ! Work variables for DO loop

    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyRadFluxShortIncoming'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_radflx % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYRADFLX'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  *phy_radflx* に格納されている設定値の取り出し
    !  Fetch setting values stored in *phy_radflx*
    !-----------------------------------------------------------------
    imax    = phy_radflx % imax   
    jmax    = phy_radflx % jmax   
    y_Lat   = phy_radflx % y_Lat   

    !-----------------------------------------------------------------
    !  年, 日平均日射の計算
    !  Calculate annual mean, daily mean insolation
    !-----------------------------------------------------------------
    do i = 0, imax - 1
      do j = 0, jmax - 1
        xy_IncomRadSFlux(i,j) = - SolarCoeff * (1.0_DP - AtmosAlbedo ) * ( AINS + BINS * cos( y_Lat(j) )**2 )

        if ( xy_IncomRadSFlux(i,j) < 0.0_DP ) then
          xy_InAngle(i,j) = 1.0_DP / ( AZET + BZET * cos( y_Lat(j) )**2 )
        else
          xy_IncomRadSFlux(i,j) = 0.0_DP
          xy_InAngle(i,j) = 0.0_DP
        end if
      end do
    end do

    !-----------------------------------------------------------------
    !  季節変化, 日変化がある場合の計算
    !  Calculate with seasonal change and diurnal change
    !-----------------------------------------------------------------
    ! 近く AGCM5 から移植予定. 
    ! Importation from AGCM5 is planed within days


    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyRadFluxShortIncoming