!c Description:
!c   断熱曲線の計算プログラム
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Histry: 
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2003-03-31    新規作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved


program obradbt
  use gtool_history
  use prm_common, only: obrinit, struct, property, emol_ini, &
       & smax, emax, step, pres_ini, pres_end, temp_ini, cnv, pseudosw 
  use interface_mingibbs
  use interface_mingibbs_chk
  use interface_entropy
  use interface_mol2emol
  use interface_rmcnd
  use interface_molfraction
  use interface_exec_time
  
  implicit none
  real(8), parameter          :: temp_chg = -0.5d0
  real(8), allocatable        :: pres_map(:)
  real(8), allocatable        :: entr_map(:)
  real(8), allocatable        :: temp_map(:)
  real(8), allocatable        :: emol_map(:,:)
  real(8), allocatable        :: mol_map(:,:)
  real(8), allocatable        :: molfr_map(:,:)
  real(8), allocatable        :: mol(:)
  real(8), allocatable        :: molfr(:)
  real(8), allocatable        :: emol(:)
  real(8)                     :: entr1
  real(8)                     :: entr2
  real(8)                     :: entr3
  real(8)                     :: entr_gas
  real(8)                     :: temp1
  real(8)                     :: temp2
  real(8)                     :: temp3
  integer                  :: i, l
  integer                  :: time1(8)
  integer                  :: time2(8)
  integer                  :: time(2)
  logical                  :: esw
  
!!!
!!!実行開始時間を計測
!!!
  call date_and_time(values = time1)


!!!
!!!データ呼び出し
!!!
  call obrinit


!!!
!!!配列の割り当て
!!!
  allocate(pres_map(step), &
       & entr_map(step), &
       & temp_map(step), &
       & mol_map(step, smax), &
       & molfr_map(step, smax), &
       & emol_map(step, emax), &
       & mol(smax))

!!!
!!!NetCDF ファイルの用意
!!!
  call output_gtool4_init


!!!
!!!計算に用いる圧力の設定
!!!
  esw = .false.
  pres_map = pres_ini
  do i = 2, step
!     pres_map(i) = pres_ini &
!          & * 10.0d0**((i - 1) * dlog10(pres_end / pres_ini) / (step - 1))
     pres_map(i) = pres_ini + (i - 1) * (pres_end - pres_ini) / (step - 1)
  end do
  call HistoryPut('p', real(pres_map, 4))


!!!
!!!初期温度, 初期圧力を別途計算
!!!
  l = 1
  temp2 = temp_ini
  call mingibbs_chk(temp2, pres_map(l), emol_ini, struct, property, mol)
  call entropy(temp2, pres_map(l), mol, &
       & property, entr_all=entr2, entr_gas=entr_gas)
  call store_data(l, esw)
  
  

!!!
!!!温度を変化させながら計算
!!!
  CHANGE_PRES: do l = 2, step
!     write(*,*) "START CHANGE_PRES"
     
     !前の圧力での収束温度でのエントロピーだけ別途計算
     temp1 = temp_map(l - 1) + 2.0d0
     call mingibbs_chk(temp1, pres_map(l), &
          & emol_map(l-1,:), struct, property, mol, mol_map(l-1,:))
     call entropy(temp1, pres_map(l), mol, property, entr_all = entr1)

     !念のためのチェック
     if (entr1 < entr_map(l-1)) then
        write(*,*) "the entropy is lower than the pervious entropy"
        write(*,*) "temp ", temp1
        write(*,*) "entr ", entr1, entr_map(l-1)
        stop
     end if

     !ステップ1: エントロピーの値を挟む温度を求める
     i = 1
     STEP1: do 
        !温度定義
        temp3 = temp1 + i * temp_chg
        
        !エントロピーの計算
        call mingibbs_chk(temp3, pres_map(l), &
             & emol_map(l-1,:), struct, property, mol, mol_map(l-1,:))
        call entropy(temp3, pres_map(l), mol, property, entr_all = entr3)

        !ループ終了条件
        if (entr3 < entr_map(l-1)) exit STEP1
        
        !ループを返すための操作
        i = i + 1
        entr1 = entr3
        temp1 = temp3
     end do STEP1
     
     !ステップ 2: 温度の中点でのエントロピーを求め, 1 ステップ前での
     !            エントロピーと同じになるまで収束させる.
     STEP2: do 
        !中間点での値を計算する
        temp2 = (temp1 + temp3) / 2.0d0
        
        call mingibbs_chk(temp2, pres_map(l), &
             & emol_map(l-1,:), struct, property, mol, mol_map(l-1,:))
        call entropy(temp2, pres_map(l), &
             & mol, property, entr_all=entr2, entr_gas=entr_gas)
        
        !ループ終了条件. 
        if (abs((entr2 - entr_map(l-1))/entr_map(l-1)) <= cnv) then
           temp_map(l)   = temp2
           call store_data(l, esw)
           write(*,*) "(", l, ") ", "Entropy is now converged" 
           exit STEP2
        end if
        
        !数値的に温度が等しい場合にはメッセージを出力. 
        if (temp1 == temp2 .OR. temp2 == temp3) then
           write(*,*) "Conserved Entropy is slightly changed", &
                & entr_map(l-1), " -> ", entr2

!           ! temp3 で計算しなおし
!           call mingibbs_chk(temp3, pres_map(l), &
!                & emol_map(l-1,:), struct, property, mol, mol_map(l-1,:))
!           call entropy(temp3, pres_map(l), &
!                & mol, property, entr_all=entr2, entr_gas=entr_gas)           

           esw = .true. 
           call store_data(l, esw)
           esw = .false.

           write(*,*) "(", l, ") ", "Entropy is now converged" 
           exit STEP2
        end if
        
        !ステップを進めるための処理
        if (entr3 < entr_map(l-1) .and. entr_map(l-1) < entr2) then
           temp1 = temp2
           entr1 = entr2
        elseif (entr2 <= entr_map(l-1) .and. entr_map(l-1) < entr1) then
           temp3 = temp2
           entr3 = entr2
        end if
     end do STEP2
  end do CHANGE_PRES


!!!
!!!NetCDF ファイルを閉じる
!!!
  call HistoryClose


!!!
!!!実行開始時間を計測
!!!
  call date_and_time(values = time2)
  call exec_time(time1, time2, time)
  write(*,*) "計算実行時間: ", time(1), "分", time(2), "秒"

contains

  subroutine output_gtool4_init
    use prm_common, only: smax, emax, planet, AdbtFile 
    use gtool_history
    
    implicit none
    real(8)                 :: spc(smax)
    real(8)                 :: elm(emax)
    integer              :: num  
    integer              :: i 
    character(1000)         :: dataset
    
    
    !ファイルに書き出すべき初期値を変数に入れる
    write(dataset, *) "pseudosw= ", pseudosw, "  ", &
         & planet, property%name
    num = len_trim(dataset)

    !大域属性, 座標軸の設定
    call HistoryCreate(&
         & file=AdbtFile, &
         & title="Equilibrium Calculation", &
         & source=dataset(:num+1), &
         & institution="GFD_Dennou_Club OBORO Project", &
         & dims=(/'spc ', 'elm ', 'pres'/), &
         & dimsizes=(/int(smax, 4), int(emax, 4), 0/), &
         & longnames=(/"chemical species", "element         ", "Pressure        "/), &
         & units=(/"1 ", "1 ", "Pa"/), &
         & origin=0.0,  &
         & interval=0.0)
    
    !次元の定義
    spc = (/(i, i = 1, smax)/)
    call HistoryPut('spc', real(spc, 4))
    elm = (/(i, i = 1, emax)/)
    call HistoryPut('elm', real(elm, 4))
    
    !変数定義, エントロピー
    call HistoryAddVariable( &
         & varname="entr", &
         & dims=(/'pres'/), &
         & longname="Entropy", &
         & units='J', &
         & xtype="double" )
    
    !変数定義, 温度
    call HistoryAddVariable( &
         & varname="temp", &
         & dims=(/'pres'/), &
         & longname="Temperature", &
         & units='K', &
         & xtype="double" )
    
    !変数定義, モル
    call HistoryAddVariable( &
         & varname="mol", &           
         & dims=(/'spc ', 'pres'/), &
         & longname="Chemical Species Abundance", &
         & units='mol', &
         & xtype="double" )

    !変数定義, モル比
    call HistoryAddVariable( &
         & varname="molfr", &           
         & dims=(/'spc ', 'pres'/), &
         & longname="Mol Fraction", &
         & units='1', &
         & xtype="double" )

    !変数定義, 元素量
    call HistoryAddVariable( &
         & varname="emol", &           
         & dims=(/'elm ', 'pres'/), &
         & longname="Element Abundance", &
         & units='mol', &
         & xtype="double" )

  end subroutine output_gtool4_init


  subroutine store_data(l, sw)
    use prm_common, only: ErrCode
    implicit none
    integer           :: l
    logical           :: sw

    ! エントロピーが収束しない場合の処理. 欠損値として扱う
    if (sw) then 
       temp_map(l)    = temp2
       mol_map(l, :)  = ErrCode 
       molfr_map(l,:) = ErrCode 
       entr_map(l)    = entr_map(l-1)
       emol_map(l,:) = emol_map(l-1,:)

    else

       temp_map(l) = temp2
       mol_map(l, :) = mol
       call molfraction(mol, property%phase, molfr)
       molfr_map(l,:) = molfr(:)
       if (pseudosw) then 
          entr_map(l)   = entr_gas
          write(*,*) entr_gas
          call rmcnd(mol, property)
          call mol2emol(mol, struct, emol)
          emol_map(l, :) = emol(:)
       else
          entr_map(l)   = entr2
          emol_map(l, :) = emol_ini   !元素量は保存
       end if
    end if

    call HistoryPut('pres', pres_map(l))
    call HistoryPut('temp', temp_map(l))
    call HistoryPut('mol', mol_map(l,:))
    call HistoryPut('molfr', molfr_map(l,:))
    call HistoryPut('entr', entr_map(l))
    call HistoryPut('emol', emol_map(l,:))
    
  end subroutine store_data
  
  
end program obradbt
