!= dcrtm 鉛直フラックス 計算 プログラム
!
!= dcrtm vertical flux calculation main program
!
! Authors::   Masanori Onishi
! Version::   $Id: main_dcrtm.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20160705 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2015. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]

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

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

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

  ! 物理定数設定2
  ! physical constants setteing 2
  !
  use planetconst, only: LightSpeed, &
    &                    PlanckConst,&
    &                    BoltzmannConst

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

  ! 宣言文 ; Declaration statements
  !
  implicit none

  ! 設定変数, 他
  ! setteing parameter, values
  !
  character(99) :: InFileName
                          ! 入力 ファイル名
                          ! input file name
  character(99) :: OutFileName
                          ! 出力ファイル名
                          ! output file name
  real(8) :: ResWN
                          ! 波数分解能 [m-1]
                          ! resolution of wavenumber [m-1]
  real(8), allocatable :: r_Press(:)
                          ! 気圧 (半整数レベル) [Pa]
                          ! pressure (half level) [Pa]
  real(8), allocatable :: w_WaveNum(:)
                          ! 波数 [m-1]
                          ! wavenumber [m-1]
  real(8), allocatable :: rw_RadFluxUp(:,:)
                          ! 惑星大気の射出に由来する上向きフラックス [W m-2 m]
                          ! upward flux originated from emittance of planetary atmospheire [W m-2 m]
  real(8), allocatable :: rw_RadFluxDn(:,:)
                          ! 惑星大気の射出に由来する下向きフラックス [W m-2 m]
                          ! downward flux originated from emittance of planetary atmospheire [W m-2 m]
  real(8), allocatable :: rw_SolFluxUp(:,:)
                          ! 中心星入射に由来する上向きフラックス [W m-2 m]
                          ! upward flux originated from emittance of central star [W m-2 m]
  real(8), allocatable :: rw_SolFluxDn(:,:)
                          ! 中心星入射に由来する下向きフラックス [W m-2 m]
                          ! downward flux originated from emittance of central star [W m-2 m]

  integer :: fnml = 12
                          ! 装置番号(name list file)
                          ! device number : namelist file
  integer :: wmax
                          ! 波数レベル最大値
                          ! wavenumber level maximum
  integer :: kmax
                          ! 鉛直レベル最大値
                          ! vertical level maximum
  integer :: k_w, k_r
                          ! 作業変数
                          ! work variables
  real(8), pointer :: a_Axis(:) => null()
                          ! 配列のサイズを調べるための配列
                          ! dummy array

  ! 計算される物理量
  ! calculated values
  !
  real(8), allocatable :: w_BrightT(:)
                          ! 輝度温度
                          ! brightness temperature [K]

  ! namelist file の設定
  ! setting of namelist file
  !
  namelist /brightt_nml/ &
    & InFileName,      &
    & OutFileName

  ! namelist file の読み込み
  ! read namelist file
  !
  open(fnml, file='dcrtm.nml')
  read(fnml, nml=brightt_nml)
  close(fnml)
  !print *, MinWN, MaxWN, ResWN

  ! 配列サイズを調べる
  ! check array size
  !
  call HistoryGetPointer(''//InFileName//'', 'r_Press', a_Axis )
  kmax = size( a_Axis ) - 1
  Nullify( a_Axis )
  !print *, kmax

  call HistoryGetPointer(''//InFileName//'', 'w_WaveNum', a_Axis )
  wmax = size( a_Axis ) - 1
  ResWN = a_Axis(2) - a_Axis(1)
  Nullify( a_Axis )


  allocate( r_Press       (0:kmax        ) )
  allocate( w_WaveNum     (        0:wmax) )
  allocate( rw_RadFluxUp  (0:kmax, 0:wmax) )
  allocate( w_BrightT     (        0:wmax) )

  ! 鉛直プロファイルの読み込み
  ! read profile data
  !
  call HistoryGet(''//InFileName//'', 'r_Press', r_Press)
  call HistoryGet(''//InFileName//'', 'w_WaveNum', w_WaveNum)
  call HistoryGet(''//InFileName//'', 'rw_RadFluxUp', rw_RadFluxUp)

  ! 輝度温度の計算
  ! brightness temperature calculation
  !
  k_r = kmax
  do k_w = 0, wmax
    if(w_WaveNum(k_w) .ne. 0.0_8) then

      w_BrightT(k_w) = LightSpeed*PlanckConst*w_WaveNum(k_w)/(BoltzmannConst * log( &
        & 2.0_8 * PI * PlanckConst * LightSpeed * LightSpeed*w_WaveNum(k_w)*w_WaveNum(k_w)*w_WaveNum(k_w) &
        & /rw_RadFluxUp(k_r,k_w) + 1.0_8) )

    else
      w_BrightT(k_w) = 0.0_8    
    end if

  end do

  ! 出力ファイルの設定
  ! setting of output file
  !
  call HistoryCreate(                                    &
    & file=''//OutFileName//'', &
    & title='Brightness Temperature; '//OutFileName//'', &
    & source='calc_brightnesstemp.f90;'//InFileName//'', &
    & institution='DCRTM',       &
    & dims=(/'w_WaveNum'/), dimsizes=(/wmax+1/),      &
    & longnames=(/'wavenumber'/),       &
    & units=(/'m-1'/)                    )

  call HistoryPut('w_WaveNum',w_WaveNum)

  call HistoryAddVariable( &
    & varname='w_BrightT', dims=(/'w_WaveNum'/), &
    & longname='brightness temperature', units='K', xtype='double')
  call HistoryPut('w_BrightT',w_BrightT)

  call HistoryClose

end program calc_brightnesstemp
