subroutine RadDennouAGCMFinalize
    !
    ! ���¹ã�¿ã�¼ã�����¡ã�¤ã�������­ã�¼ã�ºã��, 
    ! �¢ã�¸ã�¥ã�¼ã����������°ã���²ã��ä»���解é�¤ã��è¡����¾ã��. 
    !
    ! Close a restart file, and 
    ! deallocate variables in this module. 
    !

    ! �¢ã�¸ã�¥ã�¼ã����� ; USE statements
    !

    ! ���¹ã�¿ã�¼ã�����¼ã�¿å�¥å�ºå��
    ! Restart data input/output
    !
    use gtool_history, only: HistoryClose

    ! 宣�� ; Declaration statements
    !
    implicit none

    ! ���� ; Executable statement
    !

    if ( .not. rad_DennouAGCM_inited ) return


    ! �����������¤ã�¸æ�»ã��
    ! Return to default values
    !
    Old_Flux_saved = .false.

    ! ���¹ã�¿ã�¼ã�����¡ã�¤ã�������­ã�¼ã��
    ! close a restart file
    !
    call HistoryClose( history = gthst_rst ) ! (inout)

    ! �²ã��ä»���解é��
    ! Deallocation
    !
    if ( allocated( xy_IncomRadSFlux       ) ) deallocate( xy_IncomRadSFlux       )
    if ( allocated( xy_InAngle             ) ) deallocate( xy_InAngle             )
    if ( allocated( xy_TempSave            ) ) deallocate( xy_TempSave            )
    if ( allocated( xyr_RadSUwFluxSave     ) ) deallocate( xyr_RadSUwFluxSave     )
    if ( allocated( xyr_RadSDwFluxSave     ) ) deallocate( xyr_RadSDwFluxSave     )
    if ( allocated( xyr_RadLUwFluxSave     ) ) deallocate( xyr_RadLUwFluxSave     )
    if ( allocated( xyr_RadLDwFluxSave     ) ) deallocate( xyr_RadLDwFluxSave     )
    if ( allocated( xyra_DelRadLUwFluxSave ) ) deallocate( xyra_DelRadLUwFluxSave )
    if ( allocated( xyra_DelRadLDwFluxSave ) ) deallocate( xyra_DelRadLDwFluxSave )

    rad_DennouAGCM_inited = .false.


  end subroutine RadDennouAGCMFinalize

  subroutine RadDennouAGCMFlux( xyz_Temp, xyz_QVap, xyr_Press, xy_SurfTemp, xy_SurfAlbedo, xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )
    !
    ! æ¸�º¦, æ¯�æ¹�, æ°��§ã����, �¾å����������¹ã��è¨�ç®����¾ã��. 
    !
    ! Calculate radiation flux from temperature, specific humidity, and 
    ! air pressure. 
    !

    ! �¢ã�¸ã�¥ã�¼ã����� ; USE statements
    !

    ! �­æ³¢�¥å� (å¤��½å�¥å�)
    ! Short wave (insolation) incoming
    !
    use rad_short_income, only: RadShortIncome

    ! ����å®��°è¨­å®�
    ! Physical constants settings
    !
    use constants, only: Grav ! $ g $ [m s-2]. 
                              ! ���������. 
                              ! Gravitational acceleration

    ! ���»ç���
    ! Time control
    !
    use timeset, only: TimeN, EndTime, TimesetClockStart, TimesetClockStop

    ! ���¹ã�¿ã�¼ã�����¼ã�¿å�ºå��
    ! Restart data output
    !
    use gtool_history, only: HistoryPut, HistorySetTime

    ! �������°ç�����¼ã���£ã������
    ! Utilities for debug
    !
    use dc_trace, only: DbgMessage, BeginSub, EndSub

    ! 宣�� ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyz_Temp  (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     æ¸�º¦. Temperature
    real(DP), intent(in):: xyz_QVap  (0:imax-1, 1:jmax, 1:kmax)
                              ! $ q $ .     ��. Specific humidity
    real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . æ°��� (���´æ�°ã������). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! �°è¡¨�¢æ¸©åº�. 
                              ! Surface temperature
    real(DP), intent(in):: xy_SurfAlbedo (0:imax-1, 1:jmax)
                              ! �°è¡¨�¢ã������. 
                              ! Surface albedo

    real(DP), intent(out):: xyr_RadSUwFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! �­æ³¢ (�¥å�) ����������. 
                              ! Upward shortwave (insolation) flux
    real(DP), intent(out):: xyr_RadSDwFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! �­æ³¢ (�¥å�) ����������. 
                              ! Downward shortwave (insolation) flux

    real(DP), intent(out):: xyr_RadLUwFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! �·æ³¢����������. 
                              ! Upward longwave flux
    real(DP), intent(out):: xyr_RadLDwFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! �·æ³¢����������. 
                              ! Downward longwave flux
    real(DP), intent(out):: xyra_DelRadLUwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! �·æ³¢�°è¡¨æ¸�º¦å¤���. 
                              ! 
    real(DP), intent(out):: xyra_DelRadLDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! �·æ³¢�°è¡¨æ¸�º¦å¤���. 
                              ! 



    ! �業��
    ! Work variables
    !
    real(DP):: xyr_ColDenQVap (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho q \, dz $ . 
                              ! ���´å±¤ k ����ä¸�空ã��°´�¸æ���������å¯�åº�. 
                              ! Column density of water vapor above vertical level k. 
    real(DP):: xyr_ColDenDryAir (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho \, dz $ . 
                              ! ���´å±¤ k ����ä¸�空ã���æ°���������å¯�åº�. 
                              ! Column density of air above vertical level k. 
    integer:: k               ! ���´æ�¹å�������� DO ���¼ã�����業å���
                              ! Work variables for DO loop in vertical direction
    logical:: flag_rst_output
                              ! ���¹ã�¿ã�¼ã�����¡ã�¤ã���ºå����������. 
                              ! Flag for output of a restart file

    real(DP) :: xy_CosSZA(0:imax-1, 1:jmax)
    real(DP) :: DistFromStarScld
    real(DP) :: DiurnalMeanFactor

    integer:: i
    integer:: j


    real(DP) :: MaxError


    ! ���� ; Executable statement
    !

    ! ������確è�
    ! Initialization check
    !
    if ( .not. rad_DennouAGCM_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! �����������
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! ���´å±¤ k ����ä¸�空ã��������å¯�åº����ç®�
    ! Calculate column density above vertical level k
    !
    xyr_ColDenQVap  (:,:,kmax) = 0.
    xyr_ColDenDryAir(:,:,kmax) = 0.

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

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

    ! �·æ³¢���������¹ã�����
    ! Calculate long wave flux
    !
    if ( ( TimeN - PrevTimeLong >= IntTimeLong ) .or. ( .not. Old_Flux_saved ) ) then

      if ( .not. Old_Flux_saved ) then
        PrevTimeLong = TimeN
      else
        PrevTimeLong = PrevTimeLong + IntTimeLong
      end if


      call LongFlux( xyz_Temp, xy_SurfTemp, xyr_ColDenQVap, xyr_ColDenDryAir, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )


    ! �������¤ã������
    ! Use values in last time
    !
    else

      xyr_RadLUwFlux     = xyr_RadLUwFluxSave
      xyr_RadLDwFlux     = xyr_RadLDwFluxSave
      xyra_DelRadLUwFlux = xyra_DelRadLUwFluxSave
      xyra_DelRadLDwFlux = xyra_DelRadLDwFluxSave

      if ( .not. flag_rst_input ) then
        do k = 0, kmax
          xyr_RadLUwFlux(:,:,k) = xyr_RadLUwFlux(:,:,k) + xyra_DelRadLUwFlux(:,:,k,1) * ( xyz_Temp(:,:,1) - xy_TempSave )
          xyr_RadLDwFlux(:,:,k) = xyr_RadLDwFlux(:,:,k) + xyra_DelRadLDwFlux(:,:,k,1) * ( xyz_Temp(:,:,1) - xy_TempSave )

          xyra_DelRadLUwFlux(:,:,k,1) = xyra_DelRadLUwFlux(:,:,k,1) / ( xy_TempSave**3 ) * ( xyz_Temp(:,:,1)**3 )
          xyra_DelRadLDwFlux(:,:,k,1) = xyra_DelRadLDwFlux(:,:,k,1) / ( xy_TempSave**3 ) * ( xyz_Temp(:,:,1)**3 )
        end do
      else
        flag_rst_input = .false.
        call DbgMessage( '%c: restart data is used. ', c1 = module_name )
      end if

    end if



    ! �­æ³¢ (�¥å�) ���������¹ã�����
    ! Calculate short wave (insolation)
    !
    if ( ( TimeN - PrevTimeShort >= IntTimeShort ) .or. ( .not. Old_Flux_saved ) ) then

      if ( .not. Old_Flux_saved ) then
        PrevTimeShort = TimeN
      else
        PrevTimeShort = PrevTimeShort + IntTimeShort
      end if

      ! �­æ³¢�¥å����ç®�
      ! Calculate short wave (insolation) incoming radiation
      !
      call RadShortIncome( xy_CosZet         = xy_CosSZA, xy_InAngle        = xy_InAngle, DistFromStarScld  = DistFromStarScld, DiurnalMeanFactor = DiurnalMeanFactor )


      ! �­æ³¢���������¹ã���ç®�
      ! Calculate short wave (insolation) flux
      !

      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_CosSZA(i,j) > 0.0_DP ) then
            xy_IncomRadSFlux(i,j) = SolarConst / DistFromStarScld**2 * xy_CosSZA(i,j)
          else
            xy_IncomRadSFlux(i,j) = 0.0_DP
          end if
        end do
      end do
      ! Correction of incoming solar flux when diurnal mean insolation is assumed.
      ! If diurnal mean insolation is assumed, DiurnalMeanFactor = 1 / PI.
      ! If diurnal mean insolation is not assumed, DiurnalMeanFactor = 1.
      xy_IncomRadSFlux = xy_IncomRadSFlux * DiurnalMeanFactor


      ! 大æ��¢ã������������
      ! Taking atmospheric albedo into consideration
      !
      xy_IncomRadSFlux = xy_IncomRadSFlux * ( 1.0d0 - ShortAtmosAlbedo )


      call ShortFlux( xyr_ColDenQVap, xyr_ColDenDryAir, xy_SurfAlbedo, xy_IncomRadSFlux, xyr_RadSUwFlux, xyr_RadSDwFlux )


    else
      ! �������¤ã������
      ! Use values in last time

      xyr_RadSUwFlux = xyr_RadSUwFluxSave
      xyr_RadSDwFlux = xyr_RadSDwFluxSave
    end if



    ! ä»���è¨�ç®������¤ã��ä¿�å­�
    ! Save calculated values in this time
    !
    xy_TempSave            = xyz_Temp (:,:,1)

    xyr_RadSUwFluxSave     = xyr_RadSUwFlux
    xyr_RadSDwFluxSave     = xyr_RadSDwFlux

    xyr_RadLUwFluxSave     = xyr_RadLUwFlux
    xyr_RadLDwFluxSave     = xyr_RadLDwFlux
    xyra_DelRadLUwFluxSave = xyra_DelRadLUwFlux
    xyra_DelRadLDwFluxSave = xyra_DelRadLDwFlux

    if ( .not. Old_Flux_saved ) Old_Flux_saved = .true.

    ! ���¹ã�¿ã�¼ã�����¡ã�¤ã�����ºå���¿ã�¤ã���³ã�°ã�����§ã����
    ! Check output timing of a restart file
    !
!!
!! old code to be deleted
!! 
!    flag_rst_output = ( TimeN - PrevRstOutputTime >= RstFileIntTime )
!    if ( TimeN >= EndTime .and. .not. flag_rst_output_end ) then
!      flag_rst_output = .true.
!      flag_rst_output_end = .true.
!    end if
!    flag_rst_output = ( .not. TimeN == PrevRstOutputTime ) .and. flag_rst_output


    if ( TimeN - PrevRstOutputTime >= RstFileIntTime ) then
       flag_rst_output = .true.
    else
       flag_rst_output = .false.
    end if
    if ( TimeN >= EndTime .and. .not. flag_rst_output_end ) then
      flag_rst_output = .true.
      flag_rst_output_end = .true.
    end if
    if ( ( .not. TimeN == PrevRstOutputTime ) .and. flag_rst_output ) then
       flag_rst_output = .true.
    else
       flag_rst_output = .false.
    end if


    if ( flag_rst_output ) then
      ! 次å������, ä»������ºå�� (å¸���) ���� ��ä¿�å­�
      ! Save output time (expected) in this time, for next time
      !
      PrevRstOutputTime = PrevRstOutputTime + RstFileIntTime

      ! ���»ã��¨­å®�
      ! Set time
      !
      call HistorySetTime( timed = TimeN, history = gthst_rst )

      ! ���¼ã�¿å�ºå��
      ! Data output
      !
      call HistoryPut( 'PrevTimeLong', PrevTimeLong, history = gthst_rst ) ! (in)
      call HistoryPut( 'PrevTimeShort', PrevTimeShort, history = gthst_rst ) ! (in)
      call HistoryPut( 'SurfTemp', xy_TempSave, history = gthst_rst ) ! (in)
      call HistoryPut( 'RadLUwFlux', xyr_RadLUwFluxSave, history = gthst_rst ) ! (in)
      call HistoryPut( 'RadLDwFlux', xyr_RadLDwFluxSave, history = gthst_rst ) ! (in)
      call HistoryPut( 'RadSUwFlux', xyr_RadSUwFluxSave, history = gthst_rst ) ! (in)
      call HistoryPut( 'RadSDwFlux', xyr_RadSDwFluxSave, history = gthst_rst ) ! (in)
      call HistoryPut( 'DelRadLUwFlux', xyra_DelRadLUwFluxSave, history = gthst_rst ) ! (in)
      call HistoryPut( 'DelRadLDwFlux', xyra_DelRadLDwFluxSave, history = gthst_rst ) ! (in)

    end if

    ! ��������������
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine RadDennouAGCMFlux

  subroutine RadDennouAGCMInit( flag_rst )
    !
    ! rad_DennouAGCM �¢ã�¸ã�¥ã�¼ã������������è¡����¾ã��. 
    ! NAMELIST#rad_DennouAGCM_nml ����¿è¾¼�¿ã��������ç¶����§è�����¾ã��. 
    !
    ! "rad_DennouAGCM" module is initialized. 
    ! "NAMELIST#rad_DennouAGCM_nml" is loaded in this procedure. 
    !

    ! �¢ã�¸ã�¥ã�¼ã����� ; USE statements
    !

    ! �ºå�����¡ã�¤ã�����ºæ������
    ! Basic information for output files
    ! 
    use fileset, only: FileTitle, FileSource, FileInstitution
                              ! ���¼ã�¿ã���¡ã�¤ã������çµ�������´ã����çµ�ç¹�/��äº�. 
                              ! Institution or person that changes data files for the last time

    ! �����»æ�°å­¦å®��°è¨­å®�
    ! Physical and mathematical constants settings
    !
    use constants0, only: PI                    ! $ \pi $ .
                              ! �����.  Circular constant

    ! 座æ����¼ã�¿è¨­å®�
    ! Axes data settings
    !
    use axesset, only: x_Lon, x_Lon_Weight, y_Lat, y_Lat_Weight, z_Sigma, r_Sigma, z_DelSigma
                              ! $ \Delta \sigma $ (�´æ��). 
                              ! $ \Delta \sigma $ (Full)

    ! ���»ç���
    ! Time control
    !
    use timeset, only: RestartTime           ! ���¹ã�¿ã�¼ã����å§�����. 
                              ! Retart time of calculation


    ! NAMELIST ���¡ã�¤ã���¥å�����¢ã�������¼ã���£ã������
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! �����¥æ��������±ã��
    ! Calendar and Date handler
    !
    use dc_calendar, only: DCCalConvertByUnit

    ! ���¡ã�¤ã���¥å�ºå��è£���
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! ç¨��¥å�������¡ã��
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! æ¨�æº��ºå�����ç½����. Unit number of standard output

    ! çµ��¿è¾¼�¿é�¢æ�� PRESENT ���¡å¼µ���¢æ��
    ! Extended functions of intrinsic function "PRESENT"
    !
    use dc_present, only: present_and_true

    ! ��������
    ! Character handling
    !
    use dc_string, only: toChar

    ! ���¹ã�¿ã�¼ã�����¼ã�¿å�¥å�ºå��
    ! Restart data input/output
    !
    use gtool_history, only: HistoryCreate, HistoryAddAttr, HistoryAddVariable, HistoryPut, HistoryGet, HistoryGetAttr

    ! �­æ³¢�¥å� (å¤��½å�¥å�)
    ! Short wave (insolation) incoming
    !
    use rad_short_income, only : RadShortIncomeInit

    ! �£ä¹±���¡è������¾å�ä¼����¹ç�å¼�
    ! Radiative transfer equation without considering scattering
    !
    use rad_rte_nonscat, only : RadRTENonScatInit

    ! 宣�� ; Declaration statements
    !
    implicit none
    logical, intent(in), optional:: flag_rst
                              ! ���¹ã�¿ã�¼ã���§ã����������示ã��������. 
                              ! .true. ��ä¸����������´å��, 
                              ! �·æ³¢�¾å�, �­æ³¢�¾å����¢ã�������¹ã�¿ã�¼ã��
                              ! ���¡ã�¤ã����å¿�è¦��������¾ã��. 
                              ! ���¹ã�¿ã�¼ã�����¡ã�¤ã�����¢ã�������±ã��
                              ! NAMELIST#rad_DennouAGCM_nml
                              ! �§æ��å®������¾ã��. 
                              ! ������������ .false. �§ã��. 
                              ! 
                              ! Flag for restart. 
                              ! If .true. is given, 
                              ! a restart file for long radiation 
                              ! and short radiation is needed. 
                              ! Information about the restart file 
                              ! is specified by "NAMELIST#rad_DennouAGCM_nml".
                              ! Default value is .false.
                              ! 

    character(STRING):: RstInputFile
                              ! �¥å���������¹ã�¿ã�¼ã�����¼ã�¿ã�����¡ã�¤ã����
                              ! Filename of input restart data

    character(STRING):: RstOutputFile
                              ! �ºå���������¹ã�¿ã�¼ã�����¼ã�¿ã�����¡ã�¤ã����
                              ! Filename of output restart data
    character(STRING):: time_range
                              ! ���»ã����å®�. 
                              ! Specification of time

    character(TOKEN):: dummy_str
                              ! �¥å�����§ã�������������¼å���
                              ! Dummy variable for check of input
    logical:: get_err
                              ! �¥å�����������¼ã������. 
                              ! Error flag for input

    integer:: unit_nml        ! NAMELIST ���¡ã�¤ã�����¼ã���³ç���ç½����. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読ã�¿è¾¼�¿æ���� IOSTAT. 
                              ! IOSTAT of NAMELIST read

    character(STRING):: title_msg
                              ! è¡������������¡ã���»ã�¼ã��. 
                              ! Message added to title

    real(DP):: origin_time
                              ! ���������. 
                              ! Start time of calculation
    character(12):: time_unit
                              ! �¥æ������ä½�. Units of date and time

    logical:: flag_mpi_init

    real(DP):: BandWeightSum  ! ���³ã�����§ã�¤ã������
                              ! Sum of band weights
    integer:: bn              ! æ³¢é�·ã���¤ã�������� DO ���¼ã�����業å���
                              ! Work variables for DO loop in wavenumber bands


    ! NAMELIST å¤��°ç¾¤
    ! NAMELIST group name
    !
    namelist /rad_DennouAGCM_nml/ SolarConst, DelTimeLongValue, DelTimeLongUnit, DelTimeShortValue, DelTimeShortUnit, LongBandNum, LongAbsorpCoefQVap, LongAbsorpCoefDryAir, LongBandWeight, LongPathLengthFact, ShortBandNum, ShortAbsorpCoefQVap, ShortAbsorpCoefDryAir, ShortBandWeight, ShortSecScat, ShortAtmosAlbedo, RstInputFile, RstOutputFile
          !
          ! �����������¤ã���¤ã��������������ç¶� "rad_DennouAGCM#RadInit" 
          ! ���½ã�¼ã�¹ã�³ã�¼ã�������§ã������. 
          !
          ! Refer to source codes in the initialization procedure
          ! "rad_DennouAGCM#RadInit" for the default values. 
          !

    ! ���� ; Executable statement
    !

    if ( rad_DennouAGCM_inited ) return


!    flag_mpi_init = .false.
    flag_mpi_init = .true.

    ! �����������¤ã��¨­å®�
    ! Default values settings
    !

    ! �·æ³¢���������¹ç������
    ! Information for long wave flux
    !
    PrevTimeLong = RestartTime

    DelTimeLongValue = 3.0_DP
    DelTimeLongUnit  = 'hrs.'

    LongBandNum      = 4
    LongAbsorpCoefQVap     = -999.9_DP
    LongAbsorpCoefDryAir   = -999.9_DP
    LongBandWeight         = -999.9_DP
    LongAbsorpCoefQVap    (1:LongBandNum) = (/ 8.0_DP, 1.0_DP, 0.1_DP, 0.0_DP /)
    LongAbsorpCoefDryAir  (1:LongBandNum) = (/ 0.0_DP, 0.0_DP, 0.0_DP, 5.0e-5_DP /)
    LongBandWeight        (1:LongBandNum) = (/ 0.2_DP, 0.1_DP, 0.1_DP, 0.6_DP /)
    LongPathLengthFact = 1.5_DP

    ! �­æ³¢���������¹ç������
    ! Information for short wave flux
    !
    SolarConst             = 1380.0_DP

    PrevTimeShort = RestartTime

    DelTimeShortValue = 1.0_DP
    DelTimeShortUnit  = 'hrs.'

    ShortBandNum = 1
    ShortAbsorpCoefQVap    = -999.9_DP
    ShortAbsorpCoefDryAir  = -999.9_DP
    ShortBandWeight        = -999.9_DP
    ShortAbsorpCoefQVap   (1:ShortBandNum) = (/ 0.002_DP /)
    ShortAbsorpCoefDryAir (1:ShortBandNum) = (/ 0.0_DP /)
    ShortBandWeight       (1:ShortBandNum) = (/ 1.0_DP /)
    ShortSecScat           = 1.66_DP
    ShortAtmosAlbedo       = 0.2_DP

    ! ���¹ã�¿ã�¼ã�����¡ã�¤ã������
    ! Information about a restart file
    !
    RstInputFile  = ''
    RstOutputFile = 'rst_rad.nc'

    ! 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 = rad_DennouAGCM_nml, iostat = iostat_nml )             ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if


    if ( LongBandNum <= 0 ) then
      call MessageNotify( 'E', module_name, 'LongBandNum has to be positive.' )
    end if
    if ( ShortBandNum <= 0 ) then
      call MessageNotify( 'E', module_name, 'ShortBandNum has to be positive.' )
    end if


    ! ���³ã�����§ã�¤ã����¨­å®�
    ! Configure band weight
    !
    BandWeightSum = 0.
    do bn = 1, LongBandNum
      BandWeightSum =  BandWeightSum + LongBandWeight(bn)
    end do
    do bn = 1, LongBandNum
      LongBandWeight(bn) = LongBandWeight(bn) / BandWeightSum
    end do
    !
    BandWeightSum = 0.
    do bn = 1, ShortBandNum
      BandWeightSum = BandWeightSum + ShortBandWeight(bn)
    end do
    do bn = 1, ShortBandNum
      ShortBandWeight(bn) = ShortBandWeight(bn) / BandWeightSum
    end do



    ! ��������������
    ! Handle interval time
    !
    IntTimeLong  = DCCalConvertByUnit( DelTimeLongValue,  DelTimeLongUnit,  'sec' ) ! (in)
    IntTimeShort = DCCalConvertByUnit( DelTimeShortValue, DelTimeShortUnit, 'sec' ) ! (in)

    ! ���³ã����, �¸å�ä¿���, ���³ã�����§ã�¤ã�������§ã����
    ! Check number of band, absorption coefficients, band weight
    !
    call NmlutilAryValid( module_name, LongAbsorpCoefQVap, 'LongAbsorpCoefQVap', LongBandNum,        'LongBandNum' )         ! (in)

    call NmlutilAryValid( module_name, LongAbsorpCoefDryAir, 'LongAbsorpCoefDryAir', LongBandNum,          'LongBandNum' )           ! (in)

    call NmlutilAryValid( module_name, LongBandWeight, 'LongBandWeight', LongBandNum,    'LongBandNum' )      ! (in)

    call NmlutilAryValid( module_name, ShortAbsorpCoefQVap, 'ShortAbsorpCoefQVap', ShortBandNum,        'ShortBandNum' )         ! (in)

    call NmlutilAryValid( module_name, ShortAbsorpCoefDryAir, 'ShortAbsorpCoefDryAir', ShortBandNum,          'ShortBandNum' )           ! (in)

    call NmlutilAryValid( module_name, ShortBandWeight, 'ShortBandWeight', ShortBandNum,    'ShortBandNum' )      ! (in)

    ! �­æ³¢�¥å�����°ã���²ä�
    ! Allocate variables for short wave (insolation) incoming radiation
    !
    allocate( xy_IncomRadSFlux (0:imax-1, 1:jmax) )
    allocate( xy_InAngle (0:imax-1, 1:jmax) )

    ! ä¿�å­�������°ã���²ã��ä»���
    ! Allocate variables for saving
    !
    allocate( xy_TempSave            (0:imax-1, 1:jmax) )
    allocate( xyr_RadSUwFluxSave     (0:imax-1, 1:jmax, 0:kmax) )
    allocate( xyr_RadSDwFluxSave     (0:imax-1, 1:jmax, 0:kmax) )
    allocate( xyr_RadLUwFluxSave     (0:imax-1, 1:jmax, 0:kmax) )
    allocate( xyr_RadLDwFluxSave     (0:imax-1, 1:jmax, 0:kmax) )
    allocate( xyra_DelRadLUwFluxSave (0:imax-1, 1:jmax, 0:kmax, 0:1) )
    allocate( xyra_DelRadLDwFluxSave (0:imax-1, 1:jmax, 0:kmax, 0:1) )

    ! ���¹ã�¿ã�¼ã�����¡ã�¤ã�����¥å��
    ! Input restart file
    !
    if ( present_and_true( flag_rst ) ) then

      if ( trim(RstInputFile) == '' ) then
        call MessageNotify( 'E', module_name, 'a restart file is needed. ' // 'Specify the restart file to "RstInputFile" in NAMELIST "rad_DennouAGCM_nml"' )
      end if

      ! ���»æ���±ã���å¾�
      ! Get time information
      !
      time_range = 'time=' // toChar( RestartTime )

      ! ���¡ã�¤ã�������¡ã��確è�
      ! Conform an existence of an input file
      ! 
      call HistoryGetAttr( RstInputFile, 'lon', 'units', dummy_str, flag_mpi_split = flag_mpi_init, err = get_err )                                  ! (out)

      if ( get_err ) then
        call MessageNotify( 'E', module_name, 'restart/initial data file "%c" is not found.', c1 = trim(RstInputFile) )
      end if

      ! �¥å��
      ! Input
      !
      call HistoryGet( RstInputFile, 'SurfTemp', xy_TempSave, range = time_range, flag_mpi_split = flag_mpi_init ) ! (in) optional
      call HistoryGet( RstInputFile, 'RadSUwFlux', xyr_RadSUwFluxSave, range = time_range, flag_mpi_split = flag_mpi_init ) ! (in) optional
      call HistoryGet( RstInputFile, 'RadSDwFlux', xyr_RadSDwFluxSave, range = time_range, flag_mpi_split = flag_mpi_init ) ! (in) optional
      call HistoryGet( RstInputFile, 'RadLUwFlux', xyr_RadLUwFluxSave, range = time_range, flag_mpi_split = flag_mpi_init ) ! (in) optional
      call HistoryGet( RstInputFile, 'RadLDwFlux', xyr_RadLDwFluxSave, range = time_range, flag_mpi_split = flag_mpi_init ) ! (in) optional
      call HistoryGet( RstInputFile, 'DelRadLUwFlux', xyra_DelRadLUwFluxSave, range = time_range, flag_mpi_split = flag_mpi_init ) ! (in) optional
      call HistoryGet( RstInputFile, 'DelRadLDwFlux', xyra_DelRadLDwFluxSave, range = time_range, flag_mpi_split = flag_mpi_init ) ! (in) optional

      call HistoryGet( RstInputFile, 'PrevTimeLong', PrevTimeLong, range = time_range, flag_mpi_split = flag_mpi_init ) ! (in) optional

      call HistoryGet( RstInputFile, 'PrevTimeShort', PrevTimeShort, range = time_range, flag_mpi_split = flag_mpi_init ) ! (in) optional

      flag_rst_input = .true.
      Old_Flux_saved = .true.
    else
      RstInputFile = ''
      flag_rst_input = .false.
      Old_Flux_saved = .false.
    end if

    ! �ºå������������¨­å®�
    ! Configure time interval of output
    !
    PrevRstOutputTime = RestartTime

    ! ���¹ã�¿ã�¼ã�����¡ã�¤ã�������
    ! Create a restart file
    !
    title_msg = ' restart data for "' // module_name // '" module'
!!$    origin_time = EvalByUnit( RestartTime, RstFileIntUnit )
!!$    time_unit = RstFileIntUnit
    time_unit = 'sec'

    call HistoryCreate( file = RstOutputFile, title = trim(FileTitle) // trim(title_msg), source = FileSource, institution = FileInstitution, dims = (/ 'lon  ', 'lat  ', 'sig  ', 'sigm ', 'sorbl', 'time ' /), dimsizes = (/ imax, jmax, kmax, kmax + 1, 2, 0 /), longnames = (/ 'longitude                             ', 'latitude                              ', 'sigma at layer midpoints              ', 'sigma at layer end-points (half level)', 'surface or bottom layer               ', 'time                                  ' /), units = (/ 'degree_east ', 'degree_north', '1           ', '1           ', '1           ', time_unit /), xtypes = (/'double', 'double', 'double', 'double', 'int   ', 'double'/), origind = RestartTime, intervald = RstFileIntValue, flag_mpi_split = flag_mpi_init, history = gthst_rst )            ! (out) optional

    ! 座æ����¼ã�¿ã��¨­å®�
    ! Axes data settings
    !
    call HistoryAddAttr( 'lon', attrname = 'standard_name', value = 'longitude', history = gthst_rst )                    ! (inout)
    call HistoryAddAttr( 'lat', attrname = 'standard_name', value = 'latitude', history = gthst_rst )                    ! (inout)
    call HistoryAddAttr( 'sig', attrname = 'standard_name', value = 'atmosphere_sigma_coordinate', history = gthst_rst )                    ! (inout)
    call HistoryAddAttr( 'sigm', attrname = 'standard_name', value = 'atmosphere_sigma_coordinate', history = gthst_rst )                    ! (inout)
    call HistoryAddAttr( 'time', attrname = 'standard_name', value = 'time', history = gthst_rst )                    ! (inout)
    call HistoryAddAttr( 'sig', attrname = 'positive', value = 'down', history = gthst_rst )                    ! (inout)
    call HistoryAddAttr( 'sigm', attrname = 'positive', value = 'down', history = gthst_rst )                    ! (inout)

    call HistoryPut( 'lon', x_Lon / PI * 180.0_DP, history = gthst_rst )           ! (inout)
    call HistoryPut( 'lat', y_Lat / PI * 180.0_DP, history = gthst_rst )           ! (inout)
    call HistoryPut( 'sig', z_Sigma, history = gthst_rst )           ! (inout)
    call HistoryPut( 'sigm', r_Sigma, history = gthst_rst )           ! (inout)
    call HistoryPut( 'sorbl', (/ 0, 1 /), history = gthst_rst )           ! (inout)

    ! 座æ����¿ã��¨­å®�
    ! Axes weights settings
    !
    call HistoryAddVariable( 'lon_weight', (/'lon'/), 'weight for integration in longitude', 'radian', xtype = 'double', history = gthst_rst )                              ! (inout)
    call HistoryAddAttr( 'lon', attrname = 'gt_calc_weight', value = 'lon_weight', history = gthst_rst )                     ! (inout)
    call HistoryPut( 'lon_weight', x_Lon_Weight, history = gthst_rst )                     ! (inout)

    call HistoryAddVariable( 'lat_weight', (/'lat'/), 'weight for integration in latitude', units = 'radian', xtype = 'double', history = gthst_rst )                                     ! (inout)
    call HistoryAddAttr( 'lat', attrname = 'gt_calc_weight', value = 'lat_weight', history = gthst_rst )                     ! (inout)
    call HistoryPut( 'lat_weight', y_Lat_Weight, history = gthst_rst )                     ! (inout)

    call HistoryAddVariable( 'sig_weight', (/'sig'/), 'weight for integration in sigma', '1', xtype = 'double', history = gthst_rst )                     ! (inout)
    call HistoryAddAttr( 'sig', attrname = 'gt_calc_weight', value = 'sig_weight', history = gthst_rst )                     ! (inout)
    call HistoryPut( 'sig_weight', z_DelSigma, history = gthst_rst )                     ! (inout)

    call HistoryAddVariable( 'PrevTimeLong', (/ 'time' /), 'previous time at which longwave flux is calculated', 'sec', xtype = 'double', history = gthst_rst )                       ! (inout)
    call HistoryAddVariable( 'PrevTimeShort', (/ 'time' /), 'previous time at which shortwave flux is calculated', 'sec', xtype = 'double', history = gthst_rst )                       ! (inout)

    call HistoryAddVariable( 'SurfTemp', (/ 'lon ', 'lat ', 'time' /), 'surface temperature', 'K', xtype = 'double', history = gthst_rst )                       ! (inout)
    call HistoryAddVariable( 'RadSUwFlux', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'upward shortwave flux', 'W m-2', xtype = 'double', history = gthst_rst )                       ! (inout)
    call HistoryAddVariable( 'RadSDwFlux', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'downward shortwave flux', 'W m-2', xtype = 'double', history = gthst_rst )                       ! (inout)
    call HistoryAddVariable( 'RadLUwFlux', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'upward longwave flux', 'W m-2', xtype = 'double', history = gthst_rst )                       ! (inout)
    call HistoryAddVariable( 'RadLDwFlux', (/ 'lon ', 'lat ', 'sigm', 'time' /), 'downward longwave flux', 'W m-2', xtype = 'double', history = gthst_rst )                       ! (inout)
    call HistoryAddVariable( 'DelRadLUwFlux', (/ 'lon  ', 'lat  ', 'sigm ', 'sorbl', 'time ' /), 'longwave flux tendency with surface temperature', 'W m-2 K-1', xtype = 'double', history = gthst_rst )                                             ! (inout)
    call HistoryAddVariable( 'DelRadLDwFlux', (/ 'lon  ', 'lat  ', 'sigm ', 'sorbl', 'time ' /), 'longwave flux tendency with surface temperature', 'W m-2 K-1', xtype = 'double', history = gthst_rst )                                             ! (inout)



    ! Initialization of modules used in this module
    !

    ! �­æ³¢�¥å� (å¤��½å�¥å�)
    ! Short wave (insolation) incoming
    !
    call RadShortIncomeInit

    ! �£ä¹±���¡è������¾å�ä¼����¹ç�å¼�
    ! Radiative transfer equation without considering scattering
    !
    call RadRTENonScatInit


    ! �°å� ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'Restart:' )

    if ( trim(RstInputFile) == '' ) then
      call MessageNotify( 'M', module_name, '  InputFile  = <no input>', c1 = trim( RstInputFile ) )
    else
      call MessageNotify( 'M', module_name, '  InputFile       = %c', c1 = trim( RstInputFile ) )
      call MessageNotify( 'M', module_name, '    PrevTimeLong  = %f [%c]', d = (/ PrevTimeLong /), c1 = 'sec' )
      call MessageNotify( 'M', module_name, '    PrevTimeShort = %f [%c]', d = (/ PrevTimeShort /), c1 = 'sec' )
    end if

    call MessageNotify( 'M', module_name, '  OutputFile = %c', c1 = trim( RstOutputFile ) )
    call MessageNotify( 'M', module_name, '  IntTime    = %f [%c] (same as IntTime in "restart_file_io" module)', d = (/ RstFileIntValue /), c1 = trim( RstFileIntUnit ) )
!
    call MessageNotify( 'M', module_name, '  SolarConst   = %f', d = (/ SolarConst /) )
    call MessageNotify( 'M', module_name, 'DelTime:' )
    call MessageNotify( 'M', module_name, '  DelTimeLong  = %f [%c]', d = (/ DelTimeLongValue /), c1 = trim( DelTimeLongUnit ) )
    call MessageNotify( 'M', module_name, '  DelTimeShort = %f [%c]', d = (/ DelTimeShortValue /), c1 = trim( DelTimeShortUnit ) )
!
    call MessageNotify( 'M', module_name, 'LongFlux:' )
    call MessageNotify( 'M', module_name, '  LongBandNum            = %d', i = (/ LongBandNum /) )
    call MessageNotify( 'M', module_name, '  LongAbsorpCoefQVap     = (/ %*r /)', r = real( LongAbsorpCoefQVap(1:LongBandNum) ), n = (/ LongBandNum /) )
    call MessageNotify( 'M', module_name, '  LongAbsorpCoefDryAir   = (/ %*r /)', r = real( LongAbsorpCoefDryAir(1:LongBandNum) ), n = (/ LongBandNum /) )
    call MessageNotify( 'M', module_name, '  LongBandWeight         = (/ %*r /)', r = real( LongBandWeight(1:LongBandNum) ), n = (/ LongBandNum /) )
    call MessageNotify( 'M', module_name, '  LongPathLengthFact     = %f', d = (/ LongPathLengthFact /) )
!
    call MessageNotify( 'M', module_name, 'ShortFlux:' )
    call MessageNotify( 'M', module_name, '  ShortBandNum           = %d', i = (/ ShortBandNum /) )
    call MessageNotify( 'M', module_name, '  ShortAbsorpCoefQVap    = (/ %*r /)', r = real( ShortAbsorpCoefQVap(1:ShortBandNum) ), n = (/ ShortBandNum /) )
    call MessageNotify( 'M', module_name, '  ShortAbsorpCoefDryAir  = (/ %*r /)', r = real( ShortAbsorpCoefDryAir(1:ShortBandNum) ), n = (/ ShortBandNum /) )
    call MessageNotify( 'M', module_name, '  ShortBandWeight        = (/ %*r /)', r = real( ShortBandWeight(1:ShortBandNum) ), n = (/ ShortBandNum /) )
    call MessageNotify( 'M', module_name, '  ShortSecScat           = %f', d = (/ ShortSecScat /) )
    call MessageNotify( 'M', module_name, '  ShortAtmosAlbedo       = %f', d = (/ ShortAtmosAlbedo /) )

    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    rad_DennouAGCM_inited = .true.

  end subroutine RadDennouAGCMInit

  subroutine LongFlux( xyz_Temp, xy_SurfTemp, xyr_ColDenQVap, xyr_ColDenDryAir, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )
    !
    ! �·æ³¢���������¹ã���ç®�
    !
    ! Calculate long wave flux
    !

    ! �¢ã�¸ã�¥ã�¼ã����� ; USE statements
    !

    ! �����»æ�°å­¦å®��°è¨­å®�
    ! Physical and mathematical constants settings
    !
    use constants0, only: StB                   ! $ \sigma_{SB} $ . 
                              ! �¹ã�����¡ã�³ã���������³å���. 
                              ! Stefan-Boltzmann constant

    ! �£ä¹±���¡è������¾å�ä¼����¹ç�å¼�
    ! Radiative transfer equation without considering scattering
    !
    use rad_rte_nonscat, only : RadRTENonScat


    ! 宣�� ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     æ¸�º¦. Temperature

    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! �°è¡¨�¢æ¸©åº�. 
                              ! Surface temperature
    real(DP), intent(in):: xyr_ColDenQVap (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho q \, dz $ . 
                              ! ���´å±¤ k ����ä¸�空ã��°´�¸æ���������å¯�åº�. 
                              ! Column density of water vapor above vertical level k. 
    real(DP), intent(in):: xyr_ColDenDryAir (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho \, dz $ . 
                              ! ���´å±¤ k ����ä¸�空ã���æ°���������å¯�åº�. 
                              ! Column density of air above vertical level k. 
    real(DP), intent(out):: xyr_RadLUwFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! �·æ³¢����������. 
                              ! Upward longwave flux
    real(DP), intent(out):: xyr_RadLDwFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! �·æ³¢����������. 
                              ! Downward longwave flux
    real(DP), intent(out):: xyra_DelRadLUwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! �·æ³¢�°è¡¨æ¸�º¦å¤���. 
                              ! 
    real(DP), intent(out):: xyra_DelRadLDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! �·æ³¢�°è¡¨æ¸�º¦å¤���. 
                              ! 

    ! �業��
    ! Work variables
    !
    real(DP):: xyrr_Trans(0:imax-1, 1:jmax, 0:kmax, 0:kmax)
                              ! ������. 
                              ! Transmission coefficient
    real(DP):: xyz_PiB (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \pi B = \sigma T^{4} $
    real(DP):: xy_SurfPiB (0:imax-1, 1:jmax)
                              ! �°è¡¨�¢ã�� $ \pi B $ . 
                              ! $ \pi B $ on surface
    real(DP):: xy_PiDBDT1    (0:imax-1, 1:jmax)
                              ! $ \pi DBDT = 4 \sigma T^{3} at lowest level$
    real(DP):: xy_SurfPiDBDT (0:imax-1, 1:jmax)
                              ! �°è¡¨�¢ã�� $ \pi DBDT $ . 
                              ! $ \pi B $ on surface

    real(DP):: BandWeightSum  ! ���³ã�����§ã�¤ã������
                              ! Sum of band weights

    integer:: k, kk           ! ���´æ�¹å�������� DO ���¼ã�����業å���
                              ! Work variables for DO loop in vertical direction
    integer:: bn              ! æ³¢é�·ã���¤ã�������� DO ���¼ã�����業å���
                              ! Work variables for DO loop in wavenumber bands

    ! ���� ; Executable statement
    !


    ! ����¢æ�°è�ç®�
    ! Calculate transmission functions
    !
    !   Initialization
    !
    xyrr_Trans = 0.
    !
    do bn = 1, LongBandNum

      do k = 0, kmax
        do kk = k, kmax
          xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,k,kk) + LongBandWeight(bn) * exp( - LongPathLengthFact * (   LongAbsorpCoefQVap(bn) * abs(   xyr_ColDenQVap(:,:,kk) - xyr_ColDenQVap(:,:,k)  ) + LongAbsorpCoefDryAir(bn) * abs(   xyr_ColDenDryAir(:,:,kk) - xyr_ColDenDryAir(:,:,k)  ) ) )
        end do
      end do

    end do

    do k = 0, kmax
      do kk = 0, k-1
        xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,kk,k)
      end do
    end do


    ! $ \pi B $, $ \pi DBDT $ ����
    ! Calculate $ \pi B $ and $ \pi DBDT $
    !
    xyz_PiB       = StB * ( xyz_Temp**4 )
    xy_SurfPiB    = StB * ( xy_SurfTemp**4 )
    xy_PiDBDT1    = 4.0_DP * xyz_PiB(:,:,1) / xyz_Temp(:,:,1)
    xy_SurfPiDBDT = 4.0_DP * xy_SurfPiB / xy_SurfTemp



    call RadRTENonScat( xyz_PiB, xy_SurfPiB, xy_PiDBDT1, xy_SurfPiDBDT, xyrr_Trans, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )


  end subroutine LongFlux

  subroutine ShortFlux( xyr_ColDenQVap, xyr_ColDenDryAir, xy_SurfAlbedo, xy_TOARadSDwFlux, xyr_RadSUwFlux, xyr_RadSDwFlux )
    !
    ! �­æ³¢���������¹ã��è¨�ç®����¾ã��.
    !
    ! Calculate short wave flux. 
    !

    ! �¢ã�¸ã�¥ã�¼ã����� ; USE statements
    !

    ! 宣�� ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyr_ColDenQVap (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho q \, dz $ . 
                              ! ���´å±¤ k ����ä¸�空ã��°´�¸æ���������å¯�åº�. 
                              ! Column density of water vapor above vertical level k. 
    real(DP), intent(in):: xyr_ColDenDryAir (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho \, dz $ . 
                              ! ���´å±¤ k ����ä¸�空ã���æ°���������å¯�åº�. 
                              ! Column density of air above vertical level k. 
    real(DP), intent(in):: xy_SurfAlbedo (0:imax-1, 1:jmax)
                              ! �°è¡¨�¢ã������. 
                              ! Surface albedo
    real(DP), intent(in):: xy_TOARadSDwFlux (0:imax-1, 1:jmax)
                              ! �­æ³¢ (�¥å�) ����������. 
                              ! Shortwave (insolation) flux at the top of the atmosphere
    real(DP), intent(out):: xyr_RadSUwFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! �­æ³¢ (�¥å�) ����������. 
                              ! Upward shortwave (insolation) flux
    real(DP), intent(out):: xyr_RadSDwFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! �­æ³¢ (�¥å�) ����������. 
                              ! Downward shortwave (insolation) flux

    ! �業��
    ! Work variables
    !
    real(DP) :: xyr_TransForDwFlux(0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Transmittance from level k to top of the atmosphere
    real(DP) :: xyr_TransForUwFlux(0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Transmittance from surface to level k
    integer:: k               ! ���´æ�¹å�������� DO ���¼ã�����業å���
                              ! Work variables for DO loop in vertical direction
    integer:: bn              ! æ³¢é�·ã���¤ã�������� DO ���¼ã�����業å���
                              ! Work variables for DO loop in wavenumber bands

    ! ���� ; Executable statement
    !

    xyr_TransForDwFlux = 0.0_DP
    do k = 0, kmax
      do bn = 1, ShortBandNum
        xyr_TransForDwFlux(:,:,k) = xyr_TransForDwFlux(:,:,k) + ShortBandWeight(bn) * exp( - xy_InAngle * (   ShortAbsorpCoefQVap  (bn) * xyr_ColDenQVap  (:,:,k) + ShortAbsorpCoefDryAir(bn) * xyr_ColDenDryAir(:,:,k) ) )
      end do
    end do
    !
    xyr_TransForUwFlux = 0.0_DP
    do k = 0, kmax
      do bn = 1, ShortBandNum
        xyr_TransForUwFlux(:,:,k) = xyr_TransForUwFlux(:,:,k) + ShortBandWeight(bn) * exp( - xy_InAngle * (   ShortAbsorpCoefQVap  (bn) * xyr_ColDenQVap  (:,:,0) + ShortAbsorpCoefDryAir(bn) * xyr_ColDenDryAir(:,:,0) ) ) * xy_SurfAlbedo * exp( - ShortSecScat * (   ShortAbsorpCoefQVap(bn) * ( xyr_ColDenQVap  (:,:,0) - xyr_ColDenQVap  (:,:,k) ) + ShortAbsorpCoefDryAir(bn) * ( xyr_ColDenDryAir(:,:,0) - xyr_ColDenDryAir(:,:,k) ) ) )
      end do
    end do

    do k = 0, kmax
      xyr_RadSDwFlux(:,:,k) = xy_TOARadSDwFlux * xyr_TransForDwFlux(:,:,k)
      xyr_RadSUwFlux(:,:,k) = xy_TOARadSDwFlux * xyr_TransForUwFlux(:,:,k)
    end do


  end subroutine ShortFlux