!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-06-02    新規作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

program OboroAtmos
  use netcdf
  use gt4_history
  use prm_common, only: obrinit, property, AdbtFile, AtmosFile, &
       & grav, smax, emax
  use prm_phys, only: gasr
  use interface_netcdf_open
  use interface_netcdf_def 
  use interface_netcdf_redef
  use interface_netcdf_get
  use interface_netcdf_put
  use interface_netcdf_close
  use interface_molweight
  use interface_cp_noblegas
  use interface_cp_spl
  use interface_lapserate
  use interface_clouddensity

  implicit none
  real(8), allocatable        :: pres(:)       !圧力
  real(8), allocatable        :: entr(:)       !エントロピー
  real(8), allocatable        :: temp(:)       !温度
  real(8), allocatable        :: vtemp(:)      !仮温度
  real(8), allocatable        :: mol(:,:)      !化学種のモル数
  real(8), allocatable        :: emol(:,:)     !元素のモル数
  real(8), allocatable        :: molfr(:,:)    !化学種のモル比 
  real(8), allocatable        :: molwt(:)      !分子量
  real(8), allocatable        :: molwt_dry(:)  !乾燥大気の分子量
  real(8), allocatable        :: dtdz(:)       !温度減率
  real(8), allocatable        :: dtdz_low(:)   !  凝結成分少ない近似
  real(8), allocatable        :: dtdz_high(:)  !  凝結成分多い近似
  real(8), allocatable        :: dtdz_dry(:)   !  乾燥大気
  real(8), allocatable        :: vdtdz(:)      !仮温度減率
  real(8), allocatable        :: vdtdz_low(:)  !  凝結成分少ない近似
  real(8), allocatable        :: vdtdz_high(:) !  凝結成分多い近似
  real(8), allocatable        :: vdtdz_dry(:)  !  乾燥大気
  real(8), allocatable        :: stab(:)       !静的安定度
  real(8), allocatable        :: vstab(:)      !仮温度での静的安定度
  real(8), allocatable        :: vstab_low(:)  !  凝結成分少ない近似
  real(8), allocatable        :: vstab_high(:) !  凝結成分多い近似
  real(8), allocatable        :: cp(:)         !比熱
  real(8), allocatable        :: cp_dry(:)     !乾燥大気の比熱
  real(8), allocatable        :: cp_tmp(:,:)   !比熱(作業配列)
  real(8), allocatable        :: cldns(:,:)    !雲密度
  real(8), allocatable        :: mol_tmp(:,:)
  real(8), allocatable        :: molfr_tmp(:,:)
  real(8), allocatable        :: emol_tmp(:,:)
  real(8), allocatable        :: A(:)
  real(8), allocatable        :: B(:)
  integer                     :: smax_gas      !気相の化学種の数
  integer                     :: num           !番号
  integer                     :: snum          !化学種番号
  integer                     :: step          !計算ステップ数
  integer                     :: ncid
  character(30)               :: long_name     !名前
  character(10)               :: name          !名前  
  character(10)               :: dim_name1
  character(10)               :: dim_name2
  character(10)               :: units         !名前

  real(8), parameter          :: LatentHeat = 40.66d3
  integer, parameter          :: cond = 3


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


!!!
!!!NetCDF ファイルのオープン
!!!次元は pres, spc, elm であることを知っている
!!!
  !ファイルのオープン
  call netcdf_open(trim(AdbtFile), ncid)

  !変数の取得. 取り出す変数は pres, spc, elm, temp, mol, molfr 
  name = 'pres'
  call netcdf_get(ncid, name, pres)
  name = 'temp'
  call netcdf_get(ncid, name, temp)
  name = 'entr'
  call netcdf_get(ncid, name, entr)
  name = 'mol'
  call netcdf_get(ncid, name, mol_tmp)
  name = 'molfr'
  call netcdf_get(ncid, name, molfr_tmp)
  name = 'emol'
  call netcdf_get(ncid, name, emol_tmp) 

  !配列の割り当て
  step = size(pres)
  allocate(mol(step, smax), &
       & molfr(step, smax), &       
       & emol(step, emax))
  mol   = transpose(mol_tmp)
  molfr = transpose(molfr_tmp)
  emol  = transpose(emol_tmp)

  !ファイルを閉じる
  call netcdf_close(ncid)
  
  
!!!
!!!配列の割り当て
!!!
  smax_gas = count(property%phase == 1)
  allocate(molwt(step), &
       & molwt_dry(step), &
       & vtemp(step), &
       & dtdz(step), &
       & dtdz_low(step), &
       & dtdz_high(step), &
       & dtdz_dry(step), &
       & vdtdz(step), &
       & vdtdz_low(step), &
       & vdtdz_high(step), &
       & vdtdz_dry(step), &
       & stab(step), &
       & vstab(step), &
       & vstab_low(step), &
       & vstab_high(step), &
       & cp(step), &
       & cp_dry(step), &
       & cp_tmp(step, smax_gas), &
       & A(step), B(step))

  
!!!
!!!分子量の計算. 乾燥大気の分子量は大気最上層での分子量とする. 
!!!
  !気塊の分子量
  do num = 1, step 
     call molweight(molfr(num, :), property, molwt(num))
  end do
  
  !乾燥大気の分子量
  molwt_dry = molwt(step)
  
  
!!!
!!!比熱の計算. 乾燥大気の比熱は大気最上層での分子量とする. 
!!!
  !気塊の比熱
  do num = 1, step
     do snum = 1, smax
        if (property(snum)%method == 1) then 
           call cp_noblegas(temp(num), cp=cp_tmp(num, snum))
        elseif (property(snum)%method == 2) then
           call cp_spl(temp(num), property(snum), cp=cp_tmp(num, snum))
        end if
     end do
     
     cp(num) = dot_product(molfr(num, 1: smax_gas), cp_tmp(num,:))
  end do
  
  !乾燥大気の比熱. 比熱は温度(圧力)依存だが, モル比は一定. 
  do num = 1, step   
     cp_dry(num) = &
          & dot_product(molfr(step, 1: smax_gas), cp_tmp(num, 1: smax_gas))
  end do
  
  
!!!
!!!仮温度の計算
!!!
  vtemp = temp * (molwt_dry / molwt)
  
  
!!!
!!!温度減率の計算
!!!
  !温度の減率
  call lapserate(temp, pres, grav, molwt, dtdz)
  
  !仮温度の減率
  call lapserate(vtemp, pres, grav, molwt_dry, vdtdz)

  !係数. 温度 Temp で評価
  A(:) = LatentHeat * Molfr(:, cond) / ( GasR * Temp(:) )
  B(:) = ( LatentHeat ** 2.0d0 ) * MolFr(:, cond) &
       & / ( cp(:) * GasR * ( Temp(:) ** 2.0d0 ) )

  write(*,*) Molfr(:, cond)
  write(*,*) ""
  write(*,*) A
  write(*,*) ""
  write(*,*) B
  
  !温度減率
!  dtdz     =  Grav * MolWt * ( 1 + A ) / ( Cp * ( 1 + B ) )
  dtdz_low  =  Grav * MolWt * ( 1 + A - B ) / Cp
  dtdz_high =  Grav * MolWt * A / ( Cp * B )
 
 
  !仮温度減率
  vdtdz_low  =  dtdz_low  * MolWt_dry / MolWt
  vdtdz_high =  dtdz_high * MolWt_dry / MolWt
  

!!!
!!!乾燥大気の断熱温度減率
!!!
  !温度減率
  dtdz_dry = grav * molwt / cp
  
  !仮温度減率
  vdtdz_dry = grav * molwt_dry / cp_dry 
  
  
!!!
!!!静的安定度
!!!
  !静的安定度
  stab = - (grav * (dtdz - dtdz_dry) * 1.0D-3) / temp
  
  !仮温度での静的安定度
  vstab = - (grav * (vdtdz - vdtdz_dry) * 1.0D-3) / vtemp
  
  !仮温度での静的安定度
  vstab_low = - (grav * (vdtdz_low - vdtdz_dry) * 1.0D-3) / vtemp
  
  !仮温度での静的安定度
  vstab_high = - (grav * (vdtdz_high - vdtdz_dry) * 1.0D-3) / vtemp
  
!!!
!!!雲量
!!!
  call clouddensity(temp, pres, property, molfr, cldns)
  
  
!!!
!!!NetCDF ファイルの用意
!!!
  call output_gtool4_init


!!!
!!!物理量の書き出し
!!!
  call HistoryPut('entr', entr)
  call HistoryPut('temp', temp)
  call HistoryPut('vtemp', vtemp)
  call HistoryPut('mol', mol)
  call HistoryPut('emol', emol)
  call HistoryPut('molfr', molfr)
  call HistoryPut('molwt', molwt)
  call HistoryPut('molwt_dry', molwt_dry)
  call HistoryPut('dtdz',      dtdz)
  call HistoryPut('dtdz_low',  dtdz_low)
  call HistoryPut('dtdz_high', dtdz_high)
  call HistoryPut('dtdz_dry',  dtdz_dry)
  call HistoryPut('vdtdz',      vdtdz)
  call HistoryPut('vdtdz_low',  vdtdz_low)
  call HistoryPut('vdtdz_high', vdtdz_high)
  call HistoryPut('vdtdz_dry',  vdtdz_dry)
  call HistoryPut('stab', stab)
  call HistoryPut('vstab',      vstab)
  call HistoryPut('vstab_low',  vstab_low)
  call HistoryPut('vstab_high', vstab_high)
  call HistoryPut('cp', cp)
  call HistoryPut('cldns', cldns)
  

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


!!!
!!!名前を入れる
!!!
  call netcdf_open(AtmosFile, ncid)
  
  name = "sname"
  dim_name1 = "size"
  dim_name2 = "spc"
  long_name = "name of chemical species"
  units = "1"
  call netcdf_redef(ncid, name, dim_name1, dim_name2, &
       & long_name, units, nf90_char)
  call netcdf_put(ncid, name, property%name)
  
  call netcdf_close(ncid)


contains

  subroutine output_gtool4_init
    use prm_common, only: planet, AtmosFile
    
    implicit none
    real(8)                 :: spc(smax)
    real(8)                 :: elm(emax)
    real(8)                 :: sizes(20)
    integer                 :: num  
    integer                 :: i 
    character(1000)         :: dataset
    
    !ファイルに書き出すべき初期値を変数に入れる
    write(dataset, *) planet, "; pres_ini=", pres(1), &
         & "; pres_end=", pres(step), &
         & "; temp_ini=", temp(1), &
         & "; step=", step
    num = len_trim(dataset)
    
    !大域属性, 座標軸の設定
    !"size" は化学種名を入れるためのもの
    call HistoryCreate(&
         & file=AtmosFile, &
         & title="Equilibrium Calculation", &
         & source=dataset(:num+1), &
         & institution="GFD_Dennou_Club oboro Project", &
         & dims=(/'pres', 'spc ', 'elm ', 'size'/), &
         & dimsizes=(/step, smax, emax, 20/), &
         & longnames=(/"Pressure        ", "chemical species", &
         &             "element         ", "size of name    "/), &
         & units=(/"Pa", "1 ", "1 ", "1 "/), &
         & 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 HistoryPut('pres', real(pres, 4))

    sizes = (/(i, i = 1, 20)/)
    call HistoryPut('size', real(sizes, 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="vtemp", &
         & dims=(/'pres'/), &
         & longname="Vertial Temperature", &
         & units='K', &
         & xtype="double" )
    
    !変数定義, モル
    call HistoryAddVariable( &
         & varname="mol", &           
         & dims=(/'pres', 'spc '/), &
         & longname="Chemical Species Abundance", &
         & units='mol', &
         & xtype="double" )

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

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

    !変数定義, 大気の分子量
    call HistoryAddVariable( &
         & varname="molwt", &           
         & dims=(/'pres'/), &
         & longname="Atmospheric Mol Weight", &
         & units='g/mol', &
         & xtype="double" )

    !変数定義, 乾燥大気の分子量
    call HistoryAddVariable( &
         & varname="molwt_dry", &           
         & dims=(/'pres'/), &
         & longname="Atmospheric Mol Weight", &
         & units='g/mol', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz_low", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz_high", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )

    !変数定義, 乾燥大気の温度減率
    call HistoryAddVariable( &
         & varname="dtdz_dry", &           
         & dims=(/'pres'/), &
         & longname="Dry Lapse Rate", &
         & units='K/km', &
         & xtype="double" )

    !変数定義, 仮温度減率
    call HistoryAddVariable( &
         & varname="vdtdz", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )

    !変数定義, 仮温度減率
    call HistoryAddVariable( &
         & varname="vdtdz_low", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )

    !変数定義, 仮温度減率
    call HistoryAddVariable( &
         & varname="vdtdz_high", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )

    !変数定義, 乾燥大気の仮温度減率
    call HistoryAddVariable( &
         & varname="vdtdz_dry", &           
         & dims=(/'pres'/), &
         & longname="Dry Lapse Rate", &
         & units='K/km', &
         & xtype="double" )

    !変数定義, 静的安定度
    call HistoryAddVariable( &
         & varname="stab", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )

    !変数定義, 静的安定度
    call HistoryAddVariable( &
         & varname="vstab", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )

    !変数定義, 静的安定度
    call HistoryAddVariable( &
         & varname="vstab_low", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )

    !変数定義, 静的安定度
    call HistoryAddVariable( &
         & varname="vstab_high", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )

    !変数定義, 比熱
    call HistoryAddVariable( &
         & varname="cp", &           
         & dims=(/'pres'/), &
         & longname="Specific Heat", &
         & units='J/ k mol', &
         & xtype="double" )

    !変数定義, 雲密度
    call HistoryAddVariable( &
         & varname="cldns", &           
         & dims=(/'pres', 'spc '/), &
         & longname="Cloud Density", &
         & units='g/l', &
         & xtype="double" )

  end subroutine output_gtool4_init
  
end program OboroAtmos
