!= dcrtm 粒子の光学特性の netcdf data 作成プログラム
!
!= dcrtm particle optical property netcdf data making program
!
! Authors::   Masanori Onishi
! Version::   $Id: initialdata_particle.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20161110 $
! Copyright:: Copyright (C) GFD Dennou Club, 2016. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]

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

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

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

  ! 宣言文 ; Declaration statements
  !
  implicit none

  ! 設定変数, 他
  ! setteing parameter, values
  !
  integer :: fi = 11
  integer :: fnml = 12
  character(99) :: InFileParticle
  character(99) :: OutFileParticle  = 'OptDep_Particle.nc'
  integer :: kmax, k_r
  integer :: pmax, k_p
  integer :: wmax, k_w
  integer, allocatable :: p_ParNum(:)
  real(8), allocatable :: w_WaveNum(:)
  real(8), allocatable :: r_Press(:)
  integer, allocatable :: z_LayNum(:)
  real(8), allocatable :: zpw_DelOptDep(:,:,:)
  real(8), allocatable :: zpw_SingleScatA(:,:,:)
  real(8), allocatable :: zpw_AsymFact(:,:,:)
  integer :: LayNum, k_l, k_lm
  character(99) :: charDMMY

  ! namelist file の設定
  ! setting of namelist file
  !
  namelist /initial_particle_nml/ &
    & InFileParticle,     &
    & OutFileParticle

  ! namelist file の読み込み
  ! read namelist file
  !
  open(fnml, file='dcrtm.nml')
  read(fnml, nml=initial_particle_nml)
  close(fnml)
  print *, InFileParticle, OutFileParticle

  open(fi, file=InFIleParticle)
  read(fi, *) charDMMY !-- number of species, wavenumbers(>=2), layers
  read(fi, *) pmax, wmax, kmax

  wmax = wmax - 1

  print *, pmax, wmax, kmax

  allocate( r_Press        (0:kmax                ) )
  allocate( w_WaveNum      (                0:wmax) )
  allocate( p_ParNum       (        1:pmax        ) )
  allocate( z_LayNum       (1:kmax                ) )
  allocate( zpw_DelOptDep  (1:kmax, 1:pmax, 0:wmax) )
  allocate( zpw_SingleScatA(1:kmax, 1:pmax, 0:wmax) )
  allocate( zpw_AsymFact   (1:kmax, 1:pmax, 0:wmax) )

  do k_r = 1, kmax
    z_LayNum(k_r) = k_r
  end do

  read(fi, *) charDMMY !-- boundary pressure [Pa] (layer + 1)
  do k_r = 0, kmax
    read(fi, *) r_Press(k_r)
  end do

  do k_p = 1, pmax

    read(fi, *) charDMMY !== species 1
    print *, '== species:', k_p
    do k_w = 0, wmax

      read(fi, *) charDMMY !-- wavenumber 1 [m-1]
      read(fi, *) charDMMY !-- optical depth, single scattering albedo, asymmetry factor
      read(fi, *) w_WaveNum(k_w)

      print *, 'WN', w_WaveNum(k_w)
    
      do k_r = 1, kmax
        read(fi, *) zpw_DelOptDep  (k_r,k_p,k_w), &
          &         zpw_SingleScatA(k_r,k_p,k_w), &
          &         zpw_AsymFact   (k_r,k_p,k_w)
        print *, zpw_DelOptDep(k_r,k_p,k_w), zpw_SingleScatA(k_r,k_p,k_w), zpw_AsymFact(k_r,k_p,k_w)
      end do

    end do

  end do

  do k_p = 1, pmax
    p_ParNum(k_p) = k_p
  end do

  call HistoryCreate(                      &
    & file=''//OutFileParticle//'',        &
    & title=''//OutFileParticle//'',       &
    & source=''//InFileParticle//'',       &
    & institution='DCRTM',                 &
    & dims=(/'r_Press  ','z_LayNum ','p_ParNum ','w_WaveNum'/), &
    & dimsizes=(/kmax+1,kmax,pmax,wmax+1/),   &
    & longnames=(/'pressure(boundary)',   &
    &             'layer number      ',   &
    &             'particle number   ',   &
    &             'wavenumber        '/), &
    & units=(/'Pa ', '1  ', '1  ', 'm-1'/)    )

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

  call HistoryAddVariable( &
    & varname='zpw_DelOptDep', dims=(/'z_LayNum ','p_ParNum ','w_WaveNum'/), &
    & longname='optical depth of layer (particle)', units='1', xtype='double')
  call HistoryPut('zpw_DelOptDep',zpw_DelOptDep)

  call HistoryAddVariable( &
    & varname='zpw_SingleScatA', dims=(/'z_LayNum ','p_ParNum ','w_WaveNum'/), &
    & longname='single scattering albedo (particle)', units='1', xtype='double')
  call HistoryPut('zpw_SingleScatA',zpw_SingleScatA)

  call HistoryAddVariable( &
    & varname='zpw_AsymFact', dims=(/'z_LayNum ','p_ParNum ','w_WaveNum'/), &
    & longname='asymmetry factor (particle)', units='1', xtype='double')
  call HistoryPut('zpw_AsymFact',zpw_AsymFact)

  call HistoryClose


end program initialdata_parcitle
