!= 雲の光学特性計算プログラム
!
!= optical properties of Venus cloud (1993) calculation program
!
! Authors::   Masanori Onishi
! Version::   $Id: main_dcrtm.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20170209 $
! Copyright:: Copyright (C) GFD Dennou Club, 2016. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]

!2017/02/07
!2017/02/09

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

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

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


  ! Mie 理論計算プログラム
  ! Mie theory calulcation program
  !
  use mod_clscat6, only: miescat

  ! 内挿計算プログラム
  ! interpolation program
  !
  use interpol, only: PointInterPol

  ! 宣言文 ; Declaration statements
  !
  implicit none

  integer :: Flag_OUTPUT
                          ! 出力フォーマット指定フラグ
                          ! flag of output data file format
                          ! 1: netcdf file, 2: text file
  integer :: fo   = 10
                          ! 装置番号(output text file)
                          ! device number : output text file
  integer :: fi   = 11
                          ! 装置番号(input text file)
                          ! device number : input text file
  integer :: fnml = 12
                          ! 装置番号(name list file)
                          ! device number : namelist file
  character(99) :: InFileName
                          ! 鉛直プロファイル入力ファイル名
                          ! input file name (vertical profile)
  character(99) :: WaveNumData
                          ! 計算波数入力ファイル名
                          ! input file name (wavenumber)
  character(99) :: OutFileParticle = 'Particle_dcrtm.nc'
                          ! 出力ファイル名
                          ! output file name
  integer :: kmax, k_r
                          ! 鉛直グリッド用作業変数
                          ! work variables for vertical grid
  integer :: mmax, k_m
                          ! 粒子種類グリッド用作業変数
                          ! work variables for particle species grid
  integer :: wmax, k_w
                          ! 波数グリッド用作業変数
                          ! work variables for wavenumber grid
  integer :: vmax, k_v
                          ! 複素屈折率データファイル用作業変数
                          ! work variables for refractive index file
  integer, allocatable :: z_LayNum       (:)
                          ! 鉛直層の番号
                          ! number of vertical layer
  integer, allocatable :: m_ParNum       (:)
                          ! 粒子種類の番号
                          ! number of particle species
  real(8), allocatable :: r_Press        (:)
                          ! 圧力 
                          ! pressure [Pa]
  real(8), allocatable :: w_WaveNum      (:)
                          ! 波数
                          ! wavenumber [m-1]
  real(8), allocatable :: m_Rmean        (:)
                          ! 平均粒径
                          ! size distribution parameter (RM)
                          ! [ log normal: v=exp(-(ln(R/RM)/ln(RS))**2 / 2) ]
  real(8), allocatable :: m_Rsigma       (:)
                          ! 粒径の広がりパラメタ
                          ! size distribution parameter (RS)
                          ! [ log normal: v=exp(-(ln(R/RM)/ln(RS))**2 / 2) ]
  real(8), allocatable :: w_Rmed         (:)
                          ! 媒質の屈折率 (実部)
                          !  (real) refractive index of surrounding medium
  real(8) :: Rre = 1.429D+0
                          ! 粒子の屈折率 (実部)
                          !  (real) refractive index of particle
  real(8) :: Rim = 1.420D-8
                          ! 粒子の屈折率 (実部)
                          !  (imaginary) refractive index of particle
  character(99), allocatable :: m_CharRefIndexFile(:)
                          ! 粒子の複素屈折率データファイル名
                          ! file name of refractive index of particle
  real(8), allocatable :: v_RefWaveNum   (:)
                          ! 複素屈折率データファイル内の波数
                          ! wavenumber [m-1] in refractive index file
  real(8), allocatable :: v_RefReal      (:)
                          ! 複素屈折率データファイル内の屈折率の実部
                          ! real part of refractive index in refractive index file
  real(8), allocatable :: v_RefImag      (:)
                          ! 複素屈折率データファイル内の屈折率の虚部
                          ! imaginary part of refractive index in refractive index file
  real(8), allocatable :: m_WaveNumRef   (:)
                          ! reference wavenumber [m-1]
  real(8), allocatable :: m_RmedRef      (:)
                          ! reference wavenumber での媒質の屈折率実部
                          ! real part of refractive index in reference wavenumber [m-1]
  real(8), allocatable :: zm_DelOptDepRef(:,:)
                          ! reference wavenumber での光学的厚さの鉛直分布 (各鉛直大気層)
                          ! optical depth in reference wavenumber (each layer)
  real(8), allocatable :: mw_ExtCrossSec (:,:)
                          ! 消散断面積
                          ! extinction cross section
  real(8), allocatable :: mw_SingleScatA (:,:)
                          ! 一次散乱アルベド
                          ! single scattering albedo
  real(8), allocatable :: mw_AsymFact    (:,:)
                          ! 非対称因子
                          ! asymmetry factor
  real(8), pointer :: a_Axis(:) => null()
                          ! 配列のサイズを調べるための配列
                          ! dummy array
  ! 出力変数
  ! output variables
  real(8), allocatable :: zmw_DelOptDep  (:,:,:)
                          ! 光学的厚さ (each layer)
                          ! optical depth (each layer)
  real(8), allocatable :: zmw_SingleScatA(:,:,:)
                          ! 一次散乱アルベド
                          ! single scattering albedo
  real(8), allocatable :: zmw_AsymFact   (:,:,:)
                          ! 非対称因子
                          ! asymmetry factor

  ! namelist file の設定
  ! setting of namelist file
  !
  namelist /cloud_nml/ &
    & InFileName,      &
    & WaveNumData,     &
    & Flag_OUTPUT,     &
    & OutFileParticle

  ! namelist file の読み込み
  ! read namelist file
  !
  open(fnml, file='dcrtm.nml')
  read(fnml, nml=cloud_nml)
  close(fnml)

  ! 設定値のチェック
  ! check setting values
  !
  print *, '== Setting Values =='
  print *, 'InFileName:', InFileName
  print *, 'WaveNumData:', WaveNumData

  print *, 'Flag_OUTPUT:', Flag_OUTPUT
    if(Flag_OUTPUT == 1) then
      print *, '* output file is NETCDF format'
    elseif(Flag_OUTPUT == 2) then
      print *, '* output file is TEXT format'
    else
      print *, '* invalid value (Flag_OUTPUT = 1 or 2): stop'
      stop
    end if

  print *, 'OutFileParticle:', OutFileParticle

  ! 鉛直プロファイルの読み込み
  ! read vertical profile
  !
  open(fi, file=InFileName)
  read(fi, *) !-- number of species and boundary
  read(fi, *) mmax, kmax
  kmax = kmax -1
  print *, '-- number of species and boundary'
  print *, mmax, kmax

  allocate ( m_ParNum          (       1:mmax) )
  allocate ( m_CharRefIndexFile(       1:mmax) )
  allocate ( m_WaveNumRef      (       1:mmax) )
  allocate ( m_RmedRef         (       1:mmax) )
  allocate ( m_Rmean           (       1:mmax) )
  allocate ( m_Rsigma          (       1:mmax) )
  allocate ( zm_DelOptDepRef   (1:kmax,1:mmax) )
  allocate ( r_Press           (0:kmax       ) )
  allocate ( z_LayNum          (1:kmax       ) )

  read(fi, *) !-- size distribution parameter (RM,RS)
  read(fi, *) ! [ log normal: v=exp(-(ln(R/RM)/ln(RS))**2 / 2) ]

  print *, '-- size distribution parameter (RM,RS)'
  do k_m = 1, mmax
    read(fi, *) m_Rmean(k_m), m_Rsigma(k_m)
    print *, m_Rmean(k_m), m_Rsigma(k_m)
  end do ! k_m

  read(fi, *) !-- refractive index data file
  print *, '-- refractive index data file'
  do k_m = 1, mmax
    !read(fi, *) m_CharRefIndexFile(k_m)
    read(fi, '(a)') m_CharRefIndexFile(k_m)
    print *, m_CharRefIndexFile(k_m)
  end do ! k_m

  read(fi, *) !-- reference wavenumber [m-1]
  print *, '-- reference wavenumber [m-1]'
  do k_m = 1, mmax
    read(fi, *) m_WaveNumRef(k_m)
    print *, m_WaveNumRef(k_m)
  end do ! k_m

  read(fi, *) !-- (real) refractive index of surrounding medium
  print *, '-- (real) refractive index of surrounding medium'
  do k_m = 1, mmax
    read(fi, *) m_RmedRef(k_m)
    print *, m_RmedRef(k_m)
  end do ! k_m

  read(fi, *) !-- pressure [Pa], temperature [K]
  read(fi, *) !-- optical depth

  print *, '-- boundary pressure [Pa]'
  print *, '-- optical depth'

  read(fi, *) r_Press(0)
  print *, r_Press(0)

  do k_r = 1, kmax
    read(fi, *) zm_DelOptDepRef(k_r,:)
    print *, zm_DelOptDepRef(k_r,:)
    read(fi, *) r_Press(k_r)
    print *, r_Press(k_r)
  end do

  close(fi)

  ! 波数データの読み込み
  ! read wavenumber data
  !
  open(fi, file=WaveNumData)
  read(fi, *)
  read(fi, *) wmax

  wmax = wmax !-1
  print *, '-- number of wavenumbers (>=2)'
  if(wmax < 2) then 
    print *, 'invalid number (wmax >= 2):', wmax
    stop
  else
    print *, wmax
  end if

  allocate( w_WaveNum      (              1:wmax) )
  allocate( w_Rmed         (              1:wmax) )
  allocate( mw_ExtCrossSec (       1:mmax,0:wmax) )
  allocate( mw_SingleScatA (       1:mmax,0:wmax) )
  allocate( mw_AsymFact    (       1:mmax,0:wmax) )
  allocate( zmw_DelOptDep  (1:kmax,1:mmax,1:wmax) )
  allocate( zmw_SingleScatA(1:kmax,1:mmax,1:wmax) )
  allocate( zmw_AsymFact   (1:kmax,1:mmax,1:wmax) )

  read(fi, *) !'wavenumber [m-1]'
  print *, '-- target wavenumber [m-1]'
  do k_w = 1, wmax
    read(fi, *) w_WaveNum(k_w)
    print *, w_WaveNum(k_w)
  end do

  read(fi, *) !'-- (real) refractive index of surrounding medium'
  print *, '-- refractive index of surrounding medium'
  do k_w = 1, wmax
    read(fi, *) w_Rmed(k_w)
    print *, w_Rmed(k_w)
  end do

  close(fi)



  do k_m = 1, mmax

    ! 複素屈折率データの読み込み
    ! read refractive index
    !
    call HistoryGetPointer(''//m_CharRefIndexFile(k_m)//'', 'w_WaveNum', a_Axis )
    vmax = size( a_Axis )
    Nullify( a_Axis )
    print *, vmax

    allocate( v_RefWaveNum (1:vmax) )
    allocate( v_RefReal    (1:vmax) ) 
    allocate( v_RefImag    (1:vmax) )

    call HistoryGet(''//m_CharRefIndexFile(k_m)//'', 'w_WaveNum', v_RefWaveNum)
    call HistoryGet(''//m_CharRefIndexFile(k_m)//'', 'w_RefReal', v_RefReal)
    call HistoryGet(''//m_CharRefIndexFile(k_m)//'', 'w_RefImag', v_RefImag)

    if( (m_WaveNumRef(k_m) > v_RefWaveNum(vmax)) .or. &
      & (m_WaveNumRef(k_m) < v_RefWaveNum(   1)) ) then
      print *, 'invalid reference wavenumber, wavenumber:', m_WaveNumRef(k_m)
      print *, v_RefWaveNum(1), 'to', v_RefWaveNum(vmax), '[m-1] is OK.'
      stop
    end if

    do k_w = 1, wmax !in k_m
      if( (w_WaveNum(k_w) > v_RefWaveNum(vmax)) .or. &
        & (w_WaveNum(k_w) < v_RefWaveNum(   1)) ) then
        print *, 'invalid wavenumber, wavenumber:', w_WaveNum(k_w)
        print *, v_RefWaveNum(1), 'to', v_RefWaveNum(vmax), '[m-1] is OK.'
        stop
      end if
    end do ! k_w in k_m


    ! 光学特性の計算
    ! calculate optical properties
    !
    do k_w = 1, wmax ! in k_m

      print *, 'wavenumber:', w_WaveNum(k_w)

      call PointInterPol(                   & ! 2017/01/26 modify
        & 1, vmax, v_RefWaveNum, v_RefReal, & ! (in )
        & w_WaveNum(k_w),                   & ! (in )
        & Rre                               ) ! (out)

      print *, 'real part of refractive index:', Rre

      call PointInterPol(                   & ! 2017/01/26 modify
        & 1, vmax, v_RefWaveNum, v_RefImag, & ! (in )
        & w_WaveNum(k_w),                   & ! (in )
        & Rim                               ) ! (out)

      print *, 'imaginary part of refractive index', Rim

      call miescat(                    & ! 2017/01/23 onishi
        & w_WaveNum(k_w),              & ! (in )
        & m_Rmean(k_m), m_Rsigma(k_m), & ! (in )
        & w_Rmed(k_w), Rre, Rim,       & ! (in )
        & mw_ExtCrossSec(k_m,k_w),     & ! (out)
        & mw_SingleScatA(k_m,k_w),     & ! (out)
        & mw_AsymFact(k_m,k_w)         ) ! (out)

      print *, 'ExtCrossSec', mw_ExtCrossSec(k_m,k_w)
      print *, 'SingleScatA', mw_SingleScatA(k_m,k_w)
      print *, 'AsymFact', mw_AsymFact(k_m,k_w)

    end do ! k_w, in k_m


    ! reference wavenumber での ExtCorssSec の計算
    ! calculation of ExCrossSec for reference wavenumber
    !
    do k_w = 1, wmax ! k_m
      if(m_WaveNumRef(k_m) .eq. w_WaveNum(k_w)) then
        mw_ExtCrossSec(k_m,0) = mw_ExtCrossSec(k_m,k_w)
        print *, 'ExtCrossSecRef', mw_ExtCrossSec(k_m,0)
        exit
      end if

      if(k_w .eq. wmax) then

        print *, 'wavenumber(ref):', m_WaveNumRef(k_m)
  
        call PointInterPol(                   & ! 2017/01/26 modify
          & 1, vmax, v_RefWaveNum, v_RefReal, & ! (in )
          & m_WaveNumRef(k_m),                & ! (in )
          & Rre                               ) ! (out)

        print *, 'real part of refractive index:', Rre

        call PointInterPol(                   & ! 2017/01/26 modify
          & 1, vmax, v_RefWaveNum, v_RefImag, & ! (in )
          & m_WaveNumRef(k_m),                & ! (in )
          & Rim                               ) ! (out)

        print *, 'imaginary part of refractive index', Rim

        call miescat(                    & ! 2017/01/23 onishi
          & m_WaveNumRef(k_m),           & ! (in )
          & m_Rmean(k_m), m_Rsigma(k_m), & ! (in )
          & m_RmedRef(k_w), Rre, Rim,    & ! (in )
          & mw_ExtCrossSec(k_m,0),       & ! (out)
          & mw_SingleScatA(k_m,0),       & ! (out)
          & mw_AsymFact(k_m,0)           ) ! (out)

        print *, 'ExtCrossSecRef', mw_ExtCrossSec(k_m,0)
        print *, 'SingleScatARef', mw_SingleScatA(k_m,0)
        print *, 'AsymFactRef', mw_AsymFact(k_m,0)

      end if
    end do !k_w in k_m

    if(mw_ExtCrossSec(k_m,0) .eq. 0.0_8) then
      print *, 'ExtCrossSecRef = 0', mw_ExtCrossSec(k_m,0)
      print *, 'k_m:', k_m
      stop
    end if

    ! DelOptDep, SingleScatA, AsymFact の計算
    !
    do k_w = 1, wmax ! in k_m
      do k_r = 1, kmax ! in k_w, k_m
        if(zm_DelOptDepRef(k_r,k_m) .eq. 0.0_8) then
          zmw_DelOptDep  (k_r,k_m,k_w) = 0.0_8
          zmw_SingleScatA(k_r,k_m,k_w) = 0.0_8
          zmw_AsymFact   (k_r,k_m,k_w) = 0.0_8
        else
          zmw_DelOptDep  (k_r,k_m,k_w) = zm_DelOptDepRef(k_r,k_m)* &
            & mw_ExtCrossSec(k_m,k_w)/mw_ExtCrossSec(k_m,0) 
      
          zmw_SingleScatA(k_r,k_m,k_w) = mw_SingleScatA(k_m,k_w)
          zmw_AsymFact   (k_r,k_m,k_w) = mw_AsymFact   (k_m,k_w)
        end if
      end do ! k_r in k_w
    end do ! k_w 

    deallocate( v_RefWaveNum )
    deallocate( v_RefReal    ) 
    deallocate( v_RefImag    )

  end do ! k_m


  ! データ出力
  ! data output
  !
  if(Flag_OUTPUT == 1) then

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

    do k_m = 1, mmax
      m_ParNum(k_m) = k_m
    end do

    call HistoryCreate(                      &
      & file=''//OutFileParticle//'',        &
      & title=''//OutFileParticle//'',       &
      & source='initial_pollackcloud.f90',       &
      & institution='DCRTM',                 &
      & dims=(/'r_Press  ','z_LayNum ','p_ParNum ','w_WaveNum'/), &
      & dimsizes=(/kmax+1,kmax,mmax,wmax/),   &
      & 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 ',m_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',zmw_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',zmw_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',zmw_AsymFact)

    call HistoryClose

  elseif(Flag_OUTPUT == 2) then

    open(fo, file=''//OutFileParticle//'')

    write(fo, *) '-- number of species, wavenumbers (>=2), layers (boundary-1)'
    write(fo, *) mmax, wmax, kmax !2  3  32
    write(fo, *) '-- boundary pressure [Pa] (layers + 1)'
    do k_r = 0, kmax
      write(fo, *) r_Press(k_r)
    end do

    do k_m = 1, mmax
      write(fo, *) '== species', k_m

      do k_w = 1, wmax
        write(fo, *) '-- wavenumber [m-1]', k_w
        write(fo, *) '-- optical depth, single scattering albedo, asymmetry factor'

        write(fo, *) w_WaveNum(k_w)
        do k_r = 1, kmax
          write(fo, *) zmw_DelOptDep(k_r,k_m,k_w),   &
            &          zmw_SingleScatA(k_r,k_m,k_w), &
            &          zmw_AsymFact(k_r,k_m,k_w)
        end do ! k_r in k_w, k_m

      end do ! k_w in k_m

    end do ! k_m

    close(fo)

  else
    print *, 'invalid Flag_OUTPUT (1 or 2):', Flag_OUTPUT
    stop
  end if


end program main_cloud
