!= ブランク関数計算プログラム
!
!= planck function calculation program
!
! Authors::   Masanori Onishi
! Version::   $Id: gallery_blackbody.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20161126 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2016. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]

program gallery_blackbody
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !

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

  ! netcdf データ入出力
  ! netcdf data input or output
  !
  use gtool_history

  ! 物理定数設定1
  ! physical constants setteing 1
  !
  use constants0, only: PI, &
    &                   StB

  ! プランク関数計算
  ! planck function calculation
  !
  use planckf, only: fBk

  implicit none

  character(99) :: OutFileName
  integer :: k_w, wmax
  real(8) :: Temp
  real(8) :: SolConst
  real(8) :: RatioSolarConst
  real(8), allocatable :: w_WaveNum(:)
  real(8), allocatable :: w_Bk(:)
  real(8) :: MinWN
  real(8) :: MaxWN
  real(8) :: ResWN
  character( 6) :: CharTemp
  character( 6) :: CharSolConst

  OutFileName = 'blackbody.nc'

  Temp = 5800.0_8      ! 太陽の輝度温度
  SolConst = 1366.0_8  ! 太陽定数

  RatioSolarConst = SolConst/(StB*Temp*Temp*Temp*Temp)

  write(CharTemp, '(I6)') int(Temp)
  !print *, CharTemp
  write(CharSolConst, '(I6)') int(SolCOnst)
  !print *, CharSolConst

  MinWN = 1d+3
  MaxWN = 2d+7
  ResWN = 1d+2

  wmax = int( (MaxWN - MinWN)/ResWN )

  allocate( w_WaveNum(0:wmax) )
  allocate( w_Bk     (0:wmax) )

  !function fBk(WaveNum, Temp) result(a)
  do k_w = 0, wmax
    w_WaveNum(k_w) = MinWN + ResWN*k_w
    w_Bk(k_w) = PI * RatioSolarConst * fBk(w_WaveNum(k_w), Temp)
    !print *, w_Bk(k_w)
  end do

  ! 出力ファイルの設定
  ! setting of output file
  !
  call HistoryCreate(                                     &
    & file=''//OutFileName//'',                           &
    & title='black body spectrum; '//OutFileName//'',     &
    & source='gallery_blackbody.f90;Temp[K]='//CharTemp// &
    & ',SolConst[W m-2]='//CharSolConst//'',              &
    & institution='DCRTM',                                &
    & dims=(/'w_WaveNum'/), dimsizes=(/wmax+1/),          &
    & longnames=(/'wavenumber'/),                         &
    & units=(/'m-1'/)                                     )

  call HistoryPut('w_WaveNum',w_WaveNum)

  call HistoryAddVariable(                                &
    & varname='w_Bk', dims=(/'w_WaveNum'/),               &
    & longname='planck function', units='W m-2 m', xtype='double')
  call HistoryPut('w_Bk',w_Bk)

  call HistoryClose

end program gallery_blackbody
