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

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

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

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

  ! 格子点設定
  ! grid point setting
  !
  use gridset, only: GridsetSet, &
    &                imax,       &
    &                jmax,       &
    &                kmax

  ! 2方向近似放射モデル
  ! two stream approximation
  !
  use rad_rte_two_stream_app,                &
    &          only: RadRTETwoStreamAppInit, &
    &                RadRTETwoStreamAppLW,   &
    &                RadRTETwoStreamAppSW

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

  ! 物理定数設定2
  ! physical constants setteing 2
  !
  use planetconst, only: Grav, &
    &                    AvogadroNum

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

  ! opacity 計算
  ! opacity calculation
  !
  use mod_LinebyLine, only: KtoTauSSA


  ! 宣言文 ; Declaration statements
  !
  implicit none

  ! 設定変数, 他
  ! setteing parameter, values
  !
  integer :: Flag_InputFormat
                          ! 入力物理量の指定フラグ
                          ! flag of input value 
                          ! 1: optical depth and single scattering albedo
                          ! 2: absorption, scattering coefficient
  character(99) :: InFileName
                          ! 入力 ファイル名
                          ! input file name
  character(99) :: OutFileFlux
                          ! 出力ファイル名
                          ! output file name
  real(8) :: SolTemp
                          ! 中心星表面温度 [K]
                          ! surface temperature of central star [K]
  real(8) :: SolConst
                          ! 太陽定数 [W m-2]
                          ! solar constant [W m-2]
  real(8) :: InAngle
                          ! sec (太陽光入射角) [1]
                          ! sec (angle of incidence) [1]
  real(8) :: SurfAlbedo
                          ! 地表面アルベド [1]
                          ! surface albedo [1]
  real(8), allocatable :: r_Press(:)
                          ! 気圧 (半整数レベル) [Pa]
                          ! pressure (half level) [Pa]
  real(8), allocatable :: r_Temp(:)
                          ! 気温 (半整数レベル) [K]
                          ! temperature (half level) [K]
  real(8), allocatable :: rm_MixRatio(:,:)
                          ! 体積混合比 (半整数レベル) [1]
                          ! volume mixing ratio (half level) [1]
  real(8), allocatable :: m_MolWt(:)
                          ! 分子量 [kg mol-1]
                          ! molecular weight [kg mol-1]
  real(8), allocatable :: w_WaveNum(:)
                          ! 波数 [m-1]
                          ! wavenumber [m-1]
  real(8), allocatable :: rw_OptDepTOA(:,:)
                          ! 大気上端からの光学的厚さ (半整数レベル) [1]
                          ! optical depth from top of atmosphere (half level) [1]
  real(8), allocatable :: zw_SingleScatA(:,:)
                          ! 1次散乱アルベド (full level) [1]
                          ! single scattering albedo (full level) [1]
  real(8), allocatable :: rmw_LAbsorpK(:,:,:)
                          ! 線吸収の吸収係数 [m2 kg-1]
                          ! absorption coefficient for line absorption [m2 kg-1]
  real(8), allocatable :: rmw_CAbsorpK(:,:,:)
                          ! 連続吸収の吸収係数 [m2 kg-1]
                          ! absorption coefficient for continuum absotption [m2 kg-1]
  real(8), allocatable :: rmw_RayScatK(:,:,:)
                          ! Rayleigh 散乱の吸収係数 [m2 kg-1]
                          ! scattering coefficient for Rayleigh scattering [m2 kg-1]

  ! 計算される物理量
  ! calculated values
  !
  real(8), allocatable :: rw_RadUwFlux(:,:)
                          ! 惑星大気の射出に由来する上向きフラックス [W m-2]
                          ! upward flux originated from emittance of planetary atmospheire [W m-2]
  real(8), allocatable :: rw_RadDwFlux(:,:)
                          ! 惑星大気の射出に由来する下向きフラックス [W m-2]
                          ! downward flux originated from emittance of planetary atmospheire [W m-2]
  real(8), allocatable :: rw_SolUwFlux(:,:)
                          ! 中心星入射に由来する上向きフラックス [W m-2]
                          ! upward flux originated from emittance of central star [W m-2]
  real(8), allocatable :: rw_SolDwFlux(:,:)
                          ! 中心星入射に由来する下向きフラックス [W m-2]
                          ! downward flux originated from emittance of central star [W m-2]

  integer :: fnml = 12
                          ! 装置番号(name list file)
                          ! device number : namelist file
  integer :: wmax
                          ! 波数レベル最大値
                          ! wavenumber level maximum
  integer :: mmax
                          ! 分子数最大値
                          ! molecular number maxim
  integer :: k_w, k_r, k_m
                          ! 作業変数
                          ! work variables
  real(8) :: RatioSolarConst
                          ! 中心星の黒体放射に対する太陽定数の割合
                          ! ratio of solar constant to black body radiation
  real(8), pointer :: a_Axis(:) => null()
                          ! 配列のサイズを調べるための配列
                          ! dummy array

  ! for dcpam5 radiation code
  !
  !integer :: imax
  !integer :: jmax
  !integer :: kmax
  real(8), allocatable :: xyz_SSA          (:, :, :) 
                          ! single scattering albedo
  real(8), allocatable :: xyz_AF           (:, :, :) 
                          ! asynmetory factor
  real(8), allocatable :: xyr_OptDep       (:, :, :) 
                          ! optical depth from top of atmosphere
  real(8), allocatable :: xy_SurfAlbedo    (:, :) 
                          ! surface albedo
  real(8), allocatable :: xyr_PFInted      (:, :, :) 
                          ! planck function integral
  real(8), allocatable :: xy_SurfPFInted   (:, :) 
                          ! planck function in land surface
  real(8), allocatable :: xy_SurfDPFDTInted(:, :) 
                          ! planck function の温度微分
  real(8), allocatable :: xyr_RadUwFlux    (:, :, :) 
                          ! upward flux
  real(8), allocatable :: xyr_RadDwFlux    (:, :, :) 
                          ! downward flux
  real(8), allocatable :: xyra_DelRadUwFlux(:, :, :, :) 
                          ! 
  real(8), allocatable :: xyra_DelRadDwFlux(:, :, :, :)
                          !
  real(8)              :: SolarFluxTOA
                          ! solar flux of top of atmosphere
  real(8), allocatable :: xy_InAngle       (:, :)
                          ! sec (angle of incidence)


  ! namelist file の設定
  ! setting of namelist file
  !
  namelist /flux_nml/ &
    & Flag_InputFormat,&
    & InFileName,      &
    & OutFileFlux,     &
    & SurfAlbedo,      & 
    & SolTemp,         & 
    & SolConst,        & 
    & InAngle


  ! namelist file の読み込み
  ! read namelist file
  !
  open(fnml, file='dcrtm.nml')
  read(fnml, nml=flux_nml)
  close(fnml)
  print *, SurfAlbedo, SolTemp, SolConst, InAngle

  ! 配列サイズを調べる
  ! 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 )

  !print *, wmax
  !print *, ResWN
  !wmax = abs((MaxWN - MinWN)/ResWN)


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


  call RadRTETwoStreamAppInit

  !== set1: grid
  imax = 1
  jmax = 1
  !kmax = 32

  call GridsetSet(imax, jmax, kmax)


  allocate (xyz_SSA          (0:imax-1, 1:jmax, 1:kmax) )
  allocate (xyz_AF           (0:imax-1, 1:jmax, 1:kmax) )
  allocate (xyr_OptDep       (0:imax-1, 1:jmax, 0:kmax) )
  allocate (xy_SurfAlbedo    (0:imax-1, 1:jmax)         )
  allocate (xyr_PFInted      (0:imax-1, 1:jmax, 0:kmax) )
  allocate (xy_SurfPFInted   (0:imax-1, 1:jmax)         )
  allocate (xy_SurfDPFDTInted(0:imax-1, 1:jmax)         )
  allocate (xyr_RadUwFlux    (0:imax-1, 1:jmax, 0:kmax) )
  allocate (xyr_RadDwFlux    (0:imax-1, 1:jmax, 0:kmax) )
  allocate (xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) )
  allocate (xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) )
  allocate (xy_InAngle       (0:imax-1, 1:jmax)         ) 


  allocate( r_Press       (0:kmax        ) )
  allocate( r_Temp        (0:kmax        ) )
  allocate( w_WaveNum     (        0:wmax) )
  allocate( rw_RadUwFlux  (0:kmax, 0:wmax) )
  allocate( rw_RadDwFlux  (0:kmax, 0:wmax) )
  allocate( rw_SolUwFlux  (0:kmax, 0:wmax) )
  allocate( rw_SolDwFlux  (0:kmax, 0:wmax) )
  allocate( rw_OptDepTOA  (0:kmax, 0:wmax) )
  allocate( zw_SingleScatA(1:kmax, 0:wmax) )


  ! 鉛直プロファイルの読み込み
  ! read profile data
  !
  if(Flag_InputFormat == 1) then

    call HistoryGet(''//InFileName//'', 'r_Press', r_Press)
    call HistoryGet(''//InFileName//'', 'w_WaveNum', w_WaveNum)
    call HistoryGet(''//InFileName//'', 'r_Temp', r_Temp)
    call HistoryGet(''//InFileName//'', 'rw_OptDep', rw_OptDepTOA)
    call HistoryGet(''//InFileName//'', 'zw_SingleScatA', zw_SingleScatA)

  elseif(Flag_InputFormat == 2) then

    ! 配列サイズを調べる
    ! check array size
    !
    call HistoryGetPointer(''//InFileName//'', 'm_MolWt', a_Axis)

    mmax = size( a_Axis ) - 1

    Nullify( a_Axis )

    !print *, mmax

    allocate( m_MolWt       (        0:mmax        ) )
    allocate( rm_MixRatio   (0:kmax, 0:mmax        ) )
    allocate( rmw_LAbsorpK  (0:kmax, 0:mmax, 0:wmax) )
    allocate( rmw_CAbsorpK  (0:kmax, 0:mmax, 0:wmax) )
    allocate( rmw_RayScatK  (0:kmax, 0:mmax, 0:wmax) )

    call HistoryGet(''//InFileName//'', 'r_Press', r_Press)
    call HistoryGet(''//InFileName//'', 'r_Temp', r_Temp)
    call HistoryGet(''//InFileName//'', 'm_MolWt', m_MolWt)
    call HistoryGet(''//InFileName//'', 'rm_MixRatio', rm_MixRatio)
    call HistoryGet(''//InFileName//'', 'w_WaveNum', w_WaveNum)
    call HistoryGet(''//InFileName//'', 'rmw_LAbsorpK', rmw_LAbsorpK)
    call HistoryGet(''//InFileName//'', 'rmw_CAbsorpK', rmw_CAbsorpK)
    call HistoryGet(''//InFileName//'', 'rmw_RayScatK', rmw_RayScatK)

    call KtoTauSSA (                                          & ! 2016/08/05
      & kmax, mmax, wmax,                                     &
      & r_Press, r_Temp, rm_MixRatio(0:kmax,0:mmax), m_MolWt, &
      & rmw_RayScatK, rmw_LAbsorpK, rmw_CAbsorpK,             &
      & rw_OptDepTOA, zw_SingleScatA                          )

  else

    print *, 'invalid value (Flag_InputFormat):', Flag_InputFormat
    print *, 'stop'
    stop
 
  end if


  ! フラックスの計算
  ! calculate flux
  !

  ! dcpam5
  !xyz_SSA          (:, :, :) 
  xyz_AF           (:, :, :) = 0.0_8
  !xyr_OptDep       (:, :, :) 
  xy_SurfAlbedo    (:, :)    = SurfAlbedo
  !xyr_PFInted      (:, :, :)
  !xy_SurfPFInted   (:, :)
  xy_SurfDPFDTInted(:, :)    = 0.0_8
  xy_InAngle       (:, :)    = InAngle


  ! dcpam5 radiation code
  !
  do k_w = 0, wmax
    if( w_WaveNum(k_w) == 0.0_8 ) then
      rw_SolUwFlux(0:kmax,k_w) = 0.0_8
      rw_SolDwFlux(0:kmax,k_w) = 0.0_8
      rw_RadUwFlux(0:kmax,k_w) = 0.0_8
      rw_RadDwFlux(0:kmax,k_w) = 0.0_8
    else    

!      if( w_WaveNum(k_w) <= 300000.0_8 ) then
!        xy_SurfAlbedo(:,:) = 0.0_8
!      else
      xy_SurfAlbedo(:,:) = SurfAlbedo
!      endif

      do k_r = 1, kmax
        xyz_SSA(:,:,k_r) = zw_SingleScatA(k_r,k_w)

        if(xyz_SSA(0,1,k_r) >= 1.0_8) then ! 2014/04/21
          xyz_SSA(:,:,k_r) = 1.0_8 - 10.0**(-16) !2014/04/22 !0.999999_8
        end if

      end do

      do k_r = 0, kmax
        xyr_OptDep(:,:,k_r) = rw_OptDepTOA(k_r, k_w)
      end do

      xy_SurfPFInted(:, :) = PI * integral_planck(w_WaveNum(k_w) - 0.5_8, w_WaveNum(k_w) + 0.5_8, r_Temp(0))

      do k_r = 0, kmax
        xyr_PFInted(:, :,k_r) = PI * integral_planck(w_WaveNum(k_w) - 0.5_8, w_WaveNum(k_w) + 0.5_8, r_Temp(k_r))
      end do

      SolarFluxTOA = PI * RatioSolarConst * integral_planck(w_WaveNum(k_w) - 0.5_8, w_WaveNum(k_w) + 0.5_8, SolTemp)

      ! dcpam5 radiation code

      call RadRTETwoStreamAppSW(            &
        & xyz_SSA, xyz_AF,                  & ! (in)
        & xyr_OptDep,                       & ! (in)
        & xy_SurfAlbedo,                    & ! (in)
        & SolarFluxTOA, xy_InAngle,         & ! (in)
        & xyr_RadUwFlux, xyr_RadDwFlux      & ! (out)
        & )

      rw_SolUwFlux(0:kmax,k_w) = xyr_RadUwFlux(0,1,0:kmax)
      rw_SolDwFlux(0:kmax,k_w) = xyr_RadDwFlux(0,1,0:kmax)

      call RadRTETwoStreamAppLW(                           &
        & xyz_SSA, xyz_AF,                                 & ! (in)
        & xyr_OptDep,                                      & ! (in)
        & xy_SurfAlbedo,                                   & ! (in)
        & xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted,  & ! (in)
        & xyr_RadUwFlux, xyr_RadDwFlux,                    & ! (out)
        & xyra_DelRadUwFlux, xyra_DelRadDwFlux             & ! (out)
        & )

      rw_RadUwFlux(0:kmax,k_w) = xyr_RadUwFlux(0,1,0:kmax)
      rw_RadDwFlux(0:kmax,k_w) = xyr_RadDwFlux(0,1,0:kmax)

    end if
  end do

  ! 出力ファイルの設定
  ! setting of output file
  !
  call HistoryCreate(                                    &
    & file=''//OutFileFlux//'', &
    & title='Flux; '//OutFileFlux//'', &
    & source='main_flux.f90;'//InFileName//'', &
    & institution='DCRTM',       &
    & dims=(/'r_Press  ', 'w_WaveNum'/), dimsizes=(/kmax+1,wmax+1/),      &
    & longnames=(/'pressure  ', 'wavenumber'/),       &
    & units=(/'Pa ', 'm-1'/)                    )

  call HistoryPut('r_Press  ',r_Press)
  call HistoryPut('w_WaveNum',w_WaveNum)

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

  call HistoryAddVariable( &
    & varname='rw_RadFluxUp', dims=(/'r_Press  ', 'w_WaveNum'/), &
    & longname='upward flux(planet)', units='W m-2 m', xtype='double')
  call HistoryPut('rw_RadFluxUp',rw_RadUwFlux)

  call HistoryAddVariable( &
    & varname='rw_RadFluxDn', dims=(/'r_Press  ', 'w_WaveNum'/), &
    & longname='downward flux(planet)', units='W m-2 m', xtype='double')
  call HistoryPut('rw_RadFluxDn',rw_RadDwFlux)

  call HistoryAddVariable( &
    & varname='rw_SolFluxUp', dims=(/'r_Press  ', 'w_WaveNum'/), &
    & longname='upward flux(solar)', units='W m-2 m', xtype='double')
  call HistoryPut('rw_SolFluxUp',rw_SolUwFlux)

  call HistoryAddVariable( &
    & varname='rw_SolFluxDn', dims=(/'r_Press  ', 'w_WaveNum'/), &
    & longname='downward flux(solar)', units='W m-2 m', xtype='double')
  call HistoryPut('rw_SolFluxDn',rw_SolDwFlux)

  call HistoryClose


  deallocate( r_Press           )
  deallocate( r_Temp            )
  deallocate( w_WaveNum         )
  deallocate( rw_RadUwFlux      )
  deallocate( rw_RadDwFlux      )
  deallocate( rw_SolUwFlux      )
  deallocate( rw_SolDwFlux      )
  deallocate( rw_OptDepTOA      )
  deallocate( zw_SingleScatA    )

  deallocate (xyz_SSA           )
  deallocate (xyz_AF            )
  deallocate (xyr_OptDep        )
  deallocate (xy_SurfAlbedo     )
  deallocate (xyr_PFInted       )
  deallocate (xy_SurfPFInted    )
  deallocate (xy_SurfDPFDTInted )
  deallocate (xyr_RadUwFlux     )
  deallocate (xyr_RadDwFlux     )
  deallocate (xyra_DelRadUwFlux )
  deallocate (xyra_DelRadDwFlux )
  deallocate (xy_InAngle        )

end program main_flux
