!2013/02/14
!2013/07/12
!2014/05/15
!2014/06/04

program main_flux

  use gtool_history
  use mod_flux
  implicit none

  real(8) :: StBo = 5.670373e-8_8 ! CODATA 2010
                              ! $ \sigma_{SB} $ [W m-2 K-4].
                              ! ステファンボルツマン定数.
                              ! Stefan-Boltzmann constant  
  real(8) :: r_TempIn 
  real(8) :: r_TempOut
  real(8) :: r_FluxIn
  real(8) :: r_FluxOut
  real(8) :: eps !0.1_8 ! 0.0001_8 2013/07/10
  integer :: k
  real(8) :: r_FluxTop
  real(8) :: r_FluxSurface
  real(8) :: r_TempSurface
  integer :: LayNum
  integer :: kmax
  integer :: NumMinWaveNum
  integer :: NumRangeWaveNum
  character( 4) :: RunCode
  integer :: fo = 10
  integer :: mfc = 11
  character( 3) :: charEPS
  character(10) :: DataBase
  character( 3) :: charData

  real(8), allocatable :: r_Press(:) !in
  real(8), allocatable :: r_Temp (:) !in
  real(8), allocatable :: WaveNum(:) !in
  real(8), allocatable :: z_DelOptDep(:,:) !
  real(8), allocatable :: r_RadUwFlux(:,:) !out
  real(8), allocatable :: r_RadDwFlux(:,:) !out
  real(8), allocatable :: r_OptDepTOA(:,:) !in

  !==read condition file
  namelist /flux_icrccm_nml/ &
    & RunCode,         &
    & DataBase,        &
    & kmax,            &
    & NumMinWaveNum,   &
    & NumRangeWaveNum, &
    & eps,             &
    & charEPS

  !DataBase = 'HITRAN1996'
  !kmax = 32
  !NumMinWaveNum = 1
  !NumRangeWaveNum = 300000
  !RunCode = 'C35x'
  !eps = 0.001
  !charEPS = 'e-3'

  open(mfc, file='Flux_ICRCCM_onishi.nml')
  read(mfc, nml=flux_icrccm_nml)
  close(mfc)
  print *, RunCode, DataBase, kmax, NumMinWaveNum, NumRangeWaveNum, eps, charEPS

  write(charData, '(A1, A2)') 'H', DataBase(9:10)

  allocate( r_Press(0:kmax) )                        !in
  allocate( r_Temp (0:kmax) )                        !in
  allocate( WaveNum(1:NumRangeWaveNum)            )  !in
  allocate( z_DelOptDep(1:kmax, 1:NumRangeWaveNum) ) !
  allocate( r_RadUwFlux(0:kmax, 1:NumRangeWaveNum) ) !out
  allocate( r_RadDwFlux(0:kmax, 1:NumRangeWaveNum) ) !out
  allocate( r_OptDepTOA(0:kmax, 1:NumRangeWaveNum) ) !in
  !== read vertical profile & optical depth
  !
  call HistoryGet('./output_OptDepTOA/OptDepTOA_ICRCCM_'//charData//'_'//RunCode//'.nc', 'p', r_Press)
  call HistoryGet('./output_OptDepTOA/OptDepTOA_ICRCCM_'//charData//'_'//RunCode//'.nc', 'wn', WaveNum, range='wn=1:300000')
  call HistoryGet('./output_OptDepTOA/OptDepTOA_ICRCCM_'//charData//'_'//RunCode//'.nc', 'temp', r_Temp)
  call HistoryGet('./output_OptDepTOA/OptDepTOA_ICRCCM_'//charData//'_'//RunCode//'.nc', 'OptDep', r_OptDepTOA, range='wn=1:300000')

  do k = 1, kmax
    z_DelOptDep(k, :) = r_OptDepTOA(k-1, :) - r_OptDepTOA(k, :)
  end do

  !==setting
  r_TempSurface = r_Temp(0)
  r_RadDwFlux(kmax, :) = 0.0_8

  r_FluxTop = 0.0_8

  do k = 1, NumRangeWaveNum !WaveNumber loop
    ! print *, "**** wavenumber ****", k

    ! Upward Flux
    do LayNum = 0, kmax-1 ! layers loop

      if( LayNum == 0) then
        r_FluxIn = 3.1415926535897932_8 * integral_planck(WaveNum(k) - 0.5_8, WaveNum(k) + 0.5_8, r_TempSurface)
        r_RadUwFlux(0, k) = r_FluxIn
        ! print *, ":::: calculation of upward flux (non gray) ::::"
        ! print *, "==== land surface ====", "TEMP:", r_TempSurface
      else
        r_FluxIn = r_FluxOut
        ! print *, "----------------------", "TEMP:", Temp_US(LayNum)
      end if

      call calc_Flux(                                                  &
        & eps, WaveNum(k) - 0.5_8, WaveNum(k) + 0.5_8,                       & !(in ) 今回は幅1[m-1]
        & z_DelOptDep(LayNum+1,k), r_Temp(LayNum), r_Temp(LayNum + 1), r_FluxIn, & !(in )
        & r_FluxOut                                                    & !(out)
        &                                                              )
      ! print *, "↓", r_FluxOut
      r_RadUwFlux(LayNum+1, k) = r_FluxOut
    end do ! LayNum in k
    
    ! print *, "==== top of Atmosphere ===="

    ! print *, "==================downword flux=================="

    ! Downward Flux
    ! print *, "**** wavenumber ****", k
    do LayNum = kmax, 1, -1

      !!!! read(fileNum(LayNum), *) WaveNum, z_DelOptDep
      ! print *, WaveNum, z_DelOptDep(LayNum)

      if( LayNum == kmax) then
        r_FluxIn = r_FluxTop
        ! print *, ":::: calculation of downward flux (non gray) ::::"
        ! print *, "==== top of Atmosphere ====", "TEMP:", Temp_US(LayNum + 1)
      else
        r_FluxIn = r_FluxOut
        ! print *, "----------------------", "TEMP:", Temp_US(LayNum + 1)
      end if

      call calc_Flux(                                                  &
        & eps, WaveNum(k) - 0.5_8, WaveNum(k) + 0.5_8,                       & !(in ) 今回は幅1[m-1]
        & z_DelOptDep(LayNum,k), r_Temp(LayNum), r_Temp(LayNum-1), r_FluxIn, & !(in )
        & r_FluxOut                                                    & !(out)
        &                                                              )
      ! print *, "↓", r_FluxOut
      r_RadDwFlux(LayNum-1, k) = r_FluxOut
    end do ! LayNum in k

  end do ! k

  open(fo, file='Flux_ICRCCM_'//charData//'_'//RunCode//'_onishi_'//charEPS//'')
  write(fo, *) '# Flux; ICRCCM; '//DataBase//','//RunCode//', '//charEPS//''
  write(fo, *) '# wavenumber[m-1], FluxUp(r_p=0, r_p=32), FluxDn(r_p=0, r_p=32)'

  do k = NumMinWaveNum, NumRangeWaveNum
    write(fo, '(67E22.12e3)') WaveNum(k), r_RadUwFlux(0:32,k), r_RadDwFlux(0:32,k)
  end do

  close(fo)

  call HistoryCreate(                                    &
    & file='Flux_ICRCCM_'//charData//'_'//RunCode//'_onishi_'//charEPS//'.nc', &
    & title='Flux; ICRCCM;'//DataBase//','//RunCode//'', &
    & source='main_flux.f90;'//DataBase//','//RunCode//'', &
    & institution='DCRTM',       &
    & dims=(/'p ', 'wn'/), dimsizes=(/     33, 300000/),      &
    & longnames=(/'pressure  ', 'wavenumber'/),       &
    & units=(/'Pa ', 'm-1'/)                    )

  call HistoryPut('p ',r_Press)
  call HistoryPut('wn',WaveNum)

  call HistoryAddVariable( &
    & varname='temp', dims=(/'p'/), &
    & longname='temperature', units='K', xtype='double')
  call HistoryPut('temp',r_Temp)

  call HistoryAddVariable( &
    & varname='FluxUp', dims=(/'p ', 'wn'/), &
    & longname='upward flux', units='W m-2 m', xtype='double')
  call HistoryPut('FluxUp',r_RadUwFlux)

  call HistoryAddVariable( &
    & varname='FluxDn', dims=(/'p ', 'wn'/), &
    & longname='downward flux', units='W m-2 m', xtype='double')
  call HistoryPut('FluxDn',r_RadDwFlux)


  call HistoryClose

end program main_flux
