program main

  use gtool_history                                 ! モジュール指定

  implicit none

  integer                :: kmax                    ! グリッド数
  integer                :: NTemp                   ! 
  integer                :: NVMRH2O
  real(8), allocatable   :: r_Press(:)              ! Pressure
  real(8), allocatable   :: r_Temp(:)               ! Temperature

  integer                :: NMol
  integer, allocatable   :: m_MolNum(:)             ! Molecular number
  real(8), allocatable   :: m_MolWt(:)              ! Molecular weight

  real(8), allocatable   :: rm_VMR(:,:)             ! Volume Mixing ratio

  real(8), allocatable   :: z_Press(:)              ! Pressure
  real(8)                :: Ps
  real(8)                :: SurfTemp
  real(8), allocatable   :: z_Temp(:)               ! Temperature
  real(8), allocatable   :: zm_MMR(:,:)             ! Mass Mixing ratio

  logical :: FlagIncludeDry
  !
  integer :: FNVMRH2OStartIndex
  !
  real(8) :: PressMin
  real(8) :: PressMax
  integer :: NPressPer10

  real(8) :: TempMin
  real(8) :: TempMax
  real(8) :: TempInc

  real(8) :: VMRH2OMin
  real(8) :: VMRH2OMax
  integer :: NVMRH2OPer10
  !
  real(8) :: VMRCO2
  real(8) :: VMRO3

  real(8) :: Press
  real(8) :: Temp

  character(128) :: ncfn
  character(128) :: ctlfn = 'mkprofile.conf'

  integer :: FileUnit
  logical :: FlagEx
  logical :: FlagOpened
  integer :: iostat

  integer :: k
  integer :: m
  integer :: iVMRH2O
  integer :: iVMRH2Os
  integer :: iTemp

  namelist /range_nml/ &
    & FlagIncludeDry, &
    & FNVMRH2OStartIndex, &
    & PressMax, PressMin, NPressPer10, &
    & TempMin, TempMax, TempInc, &
    & VMRH2OMin, VMRH2OMax, NVMRH2OPer10, &
    & VMRCO2, VMRO3

  FlagIncludeDry = .true.

  FNVMRH2OStartIndex = 1

  PressMax    = 1.0d5
  PressMin    = 1.0d2
  NPressPer10 = 10

  TempMin = 150.0d0
  TempMax = 350.0d0
  TempInc =  50.0d0

  VMRH2OMin    = 0.5d-2
  VMRH2OMax    = 0.5d-2
  NVMRH2OPer10 = 1.0
  !
  VMRCO2       = 369.0d-6
  VMRO3        = 1.0d-6


  inquire( file = ctlfn, opened = FlagOpened, exist = FlagEx )
  if ( .not. FlagEx ) then
    stop 'Configure file does not exist.'
  end if
  if ( FlagOpened ) then
    stop 'Configure file is opened.'
  end if
  search_fileunit : do FileUnit = 50, 99
    inquire( unit = FileUnit, opened = FlagOpened, iostat = iostat )
    if ( .not. FlagOpened ) exit search_fileunit
  end do search_fileunit

  open( FileUnit, file=ctlfn, status='unknown' )
  read( FileUnit, nml = range_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist is something wrong.'
  end if
  write( 6, nml = range_nml )
  close( FileUnit )


  kmax   = int( log10( PressMax / PressMin ) * NPressPer10 ) + 1 + 1
  NTemp  = int( ( TempMax - TempMin ) / TempInc ) + 1
  NVMRH2O = int( log10( VMRH2OMax / VMRH2OMin ) * NVMRH2OPer10 ) + 1


  NMol  = 3
  allocate( m_MolNum(NMol) )
  allocate( m_MolWt (NMol) )

  do m = 1, NMol
    m_MolNum(m) = m
  end do
  m_MolWt (1) = 0.0180152636021376
  m_MolWt (2) = 0.044009767472744
  m_MolWt (3) = 0.0479970946907997


  allocate( r_Press(kmax) )
  allocate( r_Temp (kmax) )

  allocate( rm_VMR(kmax, NMol) )


  do k = 1, kmax
!!$    Press = 10.0d0**( log10( PressMax ) - dble( k - 1 ) / dble(NPressPer10) )
    if ( k == 1 ) then
      Press = 0.0d0
    else
      Press = 10.0d0**( log10( PressMin ) + dble( k-1-1 ) / dble(NPressPer10) )
    end if
    r_Press(k) = Press
  end do

  do m = 1+1, NMol
    do k = 1, kmax
      if ( NMol == 2 ) then
        rm_VMR(k,m) = VMRCO2
      end if
      if ( NMol == 3 ) then
        rm_VMR(k,m) = VMRO3
      end if
      if ( m > 3 ) then
        rm_VMR(k,m) = 0.0d0
      end if
    end do
  end do

  if ( FlagIncludeDry ) then
    iVMRH2Os = 0
  else
    iVMRH2Os = 1
  end if


  allocate( z_Press(kmax-1) )
  allocate( z_Temp (kmax-1) )
  allocate( zm_MMR (kmax-1,NMol) )


  do iVMRH2O = iVMRH2Os, NVMRH2O

    ! mixing ratio of H2O
    if ( iVMRH2O == 0 ) then
      do k = 1, kmax
        rm_VMR(k,1) = 0.0d0
      end do
    else
      do k = 1, kmax
        rm_VMR(k,1) = &
          & 10.0d0**( log10( VMRH2OMin ) + dble( iVMRH2O - 1 ) / dble(NVMRH2OPer10) )
      end do
    end if

    do iTemp = 1, NTemp

      Temp = TempMin + TempInc * ( iTemp - 1 )

      r_Temp = Temp


      ! In practice, this will not be used. 
      do k = 1, kmax-1
        z_Press(k) = ( r_Press(k) - r_Press(k+1 ) ) / 2.0d0
      end do
      Ps       = r_Press(1)
      SurfTemp = Temp
      z_Temp   = Temp
      zm_MMR   = 0.0d0


      if ( iVMRH2O == 0 ) then
        write( ncfn, '(a,i4.4,a,i4.4,a)' ) &
          & 'data_profile/Profile_VMRH2O', iVMRH2O, &
          & '_T', int( Temp ), 'K.nc'
      else
        write( ncfn, '(a,i4.4,a,i4.4,a)' ) &
          & 'data_profile/Profile_VMRH2O', iVMRH2O + FNVMRH2OStartIndex - 1, &
          & '_T', int( Temp ), 'K.nc'
      end if
!!$        & 'data_profile/Profile_VMRH2O', iVMRH2O, '_T', int( Temp ), 'K.nc'
!!$        & 'data_profile/Profile_VMRH2O', iVMRH2O, '_T', iTemp, '.nc'

      call HistoryCreate( &                             ! ヒストリー作成
        & file=ncfn, title='Atmospheric profile', &
        & source='Program for making profiles',   &
        & institution='GFD_Dennou Club dcmodel project',       &
        & dims=(/'r_Press ', 'z_Press ', 'm_MolNum'/), &
        & dimsizes=(/kmax, kmax-1, NMol/), &
        & xtypes=(/'double', 'double', 'double'/), &
        & longnames=(/'Pressure        ', 'Pressure        ', &
        & 'molecular number' /),       &
        & units=(/'Pa', 'Pa', '1 '/) )

      call HistoryPut('r_Press',r_Press)                     ! 次元変数出力
      call HistoryPut('m_MolNum',m_MolNum)                   ! 次元変数出力

      call HistoryAddVariable( &                        ! 変数定義
        & varname='r_Temp', dims=(/'r_Press'/), &
        & longname='air_temperature', units='K', xtype='double')
      call HistoryPut('r_Temp', r_Temp)                      ! 変数出力

      call HistoryAddVariable( &                        ! 変数定義
        & varname='m_MolWt', dims=(/'m_MolNum'/), &
        & longname='molecular weight', units='kg mol-1', xtype='double')
      call HistoryPut('m_MolWt', m_MolWt)                      ! 変数出力

      call HistoryAddVariable( &                        ! 変数定義
        & varname='rm_MixRatio', dims=(/'r_Press ', 'm_MolNum'/), &
        & longname='volume mixing ratio', units='1', xtype='double')
      call HistoryPut('rm_MixRatio', rm_VMR)              ! 変数出力


      !----------------------------------------------------

      call HistoryPut('z_Press', z_Press)                        ! 次元変数出力

      call HistoryAddVariable( &                                 ! 変数定義
        & varname='SurfTemp', dims=(/''/), &
        & longname='surface_temperature', units='K', xtype='double')
      call HistoryPut('SurfTemp', SurfTemp)                      ! 変数出力

      call HistoryAddVariable( &                        ! 変数定義
        & varname='Ps', dims=(/''/), &
        & longname='surface_pressure', units='Pa', xtype='double')
      call HistoryPut('Ps', Ps)                      ! 変数出力

      call HistoryAddVariable( &                        ! 変数定義
        & varname='z_Temp', dims=(/'z_Press'/), &
        & longname='air_temperature', units='K', xtype='double')
      call HistoryPut('z_Temp', z_Temp)                      ! 変数出力

      call HistoryAddVariable( &                        ! 変数定義
        & varname='zm_MassMixRatio', dims=(/'z_Press ', 'm_MolNum'/), &
        & longname='mass mixing ratio', units='1', xtype='double')
      call HistoryPut('zm_MassMixRatio', zm_MMR)              ! 変数出力

      call HistoryClose

    end do
  end do

  stop
end program main
