!c Description:
!c   大気科学の諸々の物理量の計算
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

program OboroAtmos
  use netcdf
  use gt4_history
  use prm_phys, only: gasr
  use prm_common, only: obrinit, property, AdbtFile, AtmosFile, &
       & grav, smax, emax
  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

  implicit none
  real(8), allocatable        :: pres(:)       !圧力
  real(8), allocatable        :: entr(:)       !エントロピー
  real(8), allocatable        :: temp(:)       !温度
  real(8), allocatable        :: mol(:,:)      !化学種のモル数
  real(8), allocatable        :: emol(:,:)     !元素のモル数
  real(8), allocatable        :: molfr(:,:)    !化学種のモル比 
  real(8), allocatable        :: molwt(:)      !分子量
  real(8), allocatable        :: dTdz(:)       !温度減率
  real(8), allocatable        :: dtdz_dry(:)   !温度減率
  real(8), allocatable        :: dMdz(:)       
  real(8), allocatable        :: dzdp(:)       
  real(8), allocatable        :: dT(:)       
  real(8), allocatable        :: dM(:)       
  real(8), allocatable        :: dX(:,:)       
  real(8), allocatable        :: dp(:)       
  real(8), allocatable        :: stab(:)       !静的安定度
  real(8), allocatable        :: cp(:)         !比熱
  real(8), allocatable        :: cldns(:,:)    !雲密度
  real(8), allocatable        :: cp_tmp(:,:)   !比熱(作業配列)
  real(8), allocatable        :: mol_tmp(:,:)
  real(8), allocatable        :: molfr_tmp(:,:)
  real(8), allocatable        :: emol_tmp(:,:)
  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         !名前

  integer                     :: w
  integer                     :: l
  real(8), allocatable        :: height(:)
  real(8), allocatable        :: lambda(:)
  real(8), allocatable        :: Mv(:)
  real(8), allocatable        :: Md(:)
  real(8), allocatable        :: dPw(:)
  real(8), allocatable        :: stab_low(:)
  real(8), allocatable        :: stab_high(:)
  integer                     :: csw
  real(8)                     :: IntStabDz
  real(8)                     :: NH_height(1000)
  real(8)                     :: NH_Stab(1000)
  
!!!
!!!データ呼び出し
!!!
  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), &
       & dMdz(step), &
       & dTdz(step), &
       & dTdz_dry(step), &
       & dzdp(step), &
       & stab(step), &
       & cldns(step, smax), &
       & cp(step), &
       & cp_tmp(step, smax), &
       & height(step), &
       & dT(step), &
       & dM(step), &
       & dX(step, smax), &
       & dp(step) )
  
  molwt = 0.0d0
  dMdz  = 0.0d0
  dTdz  = 0.0d0; dTdz_dry = 0.0d0
  dzdp  = 0.0d0 
  stab  = 0.0d0
  cldns = 0.0d0
  cp    = 0.0d0; cp_tmp = 0.0d0
  dT    = 0.0d0
  dM    = 0.0d0
  dX    = 0.0d0
  dp    = 0.0d0
  height = 0.0d0
  
!!!
!!!気塊の分子量の計算. 
!!!
  do num = 1, step 
     call molweight(molfr(num, :), property, molwt(num))
  end do
  
  molwt     = molwt * 1.0d-3   !単位を変換. kg/mol に
  
  
!!!
!!!気塊の比熱の計算. 
!!!
  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,:))

!     write(*,*) num, "temp", temp(num), "cp", cp_tmp(num, 3), cp_tmp(num, 1)
     
  end do
  
  
!!!
!!!乾燥断熱温度減率の計算
!!!
  dtdz_dry = - ( molwt * grav ) / cp


!!!
!!!鉛直微分
!!!
  !=== 圧力差の計算
  do num = 2, step - 1
     dp(num) = ( pres(num + 1) - pres(num - 1) ) / 2.0d0
  end do

  !=== 温度差の計算
  do num = 2, step - 1
     dT(num) = ( Temp(num + 1) - Temp(num - 1) ) / 2.0d0
  end do

  !=== 分子量差の計算
  do num = 2, step - 1
     dM(num) = ( MolWt(num + 1) - MolWt(num - 1) ) / 2.0d0
  end do

  !=== モル比の変化
  do num = 2, step - 1
     dX(num,:) = ( molfr(num+1,:) - molfr(num-1,:) ) / 2.0D0
  end do
  
  !=== 断熱減率の計算
  dTdz = - dT * pres * grav * molwt / ( dp * gasr * temp )
  
  !=== 分子量減率の計算
  dMdz = - dM * pres * grav * molwt / ( dp * gasr * temp )

  !=== 高度変化
  dzdp = - ( gasr * temp ) / ( molwt * pres * grav)

  !=== 高度積分
  do num = 2, step - 1
     height(num) = height(num-1) + dzdp(num) * dp(num)
  end do



!!!
!!!静的安定度
!!!
  stab = grav * ( dtdz - dtdz_dry ) / temp  - grav *  dMdz / molwt 

  where ( stab < 1.0d-7 )
     stab = 1.0d-7
  end where


  !(N dz) の積分
  IntStabDz = 0.0d0
  NH_Stab   = 0.0d0
  NH_height = 0.0d0

  Int: do num = 2, step - 2  ! step-1, step での Z の値は決まっていない. 

     if (pres(num) < 1.0d5)   exit Int
     if (Stab(num) <= 1.0d-7) cycle Int
        
     IntStabDz =   intStabDz &
          &      + sqrt( Stab(num) ) * ( height(num+1) - height(num-1) ) / 2.0d0
     !値の保管
     NH_Stab(num)   = Stab(num)
     NH_height(num) = height(num)

     write(*,*) num, pres(num), Stab(num),height(num), intStabDz

  end do Int

  !値の出力
  write(*,*) "Maximum value of N2: ", maxval(NH_Stab, 1)
  write(*,*) "Cloud Layer depth: ",   maxval(NH_height, 1) - minval(NH_height, 1, NH_height > 0.0d0)
  write(*,*) "Integral of Stability over height (NH): ", &
       & IntStabDz
  write(*,*) "Wave Speed: ", &
       & IntStabDz / 106.2384731188433d0 * 130.0d0, & 
       & "( ", IntStabDz / 106.2384731188433d0, ")"
  write(*,*) "Wave Speed: ", &
       & IntStabDz / 117.9405045091689d0 * 130.0d0, & 
       & "( ", IntStabDz / 117.9405045091689d0, ")"



!!!
!!!静的安定度の解析解(上限, 下限) 
!!!水についてのみ考えていることに注意.
!!!
  allocate( lambda(step), Mv(step), Md(step), dPw(step), &
       & stab_low(step), stab_high(step) )

  csw = 0
  w = 3;  l = 7
  Mv = 18.0d-3; Md = 2.2d-3
  dPw       = 0.0d0
  lambda    = 1.0d-40
  stab_low  = 0.0d0
  stab_high = 0.0d0

  do num = 2, step - 1
     dPw(num) = &
          & ( pres(num+1) * molfr(num+1,w) - pres(num-1) * molfr(num-1,w) ) &
          & / 2.0d0
  end do


  do num = 2, step - 1
     if ( mol(num, l) > 1.0d-40  ) then 
        
        lambda(num) = &
             & ( dPw(num) * gasr * ( temp(num)**2.0d0 ) ) &
             &   / ( pres(num) * molfr(num,w) * dT(num) )

!        write(*,*) "lambda", lambda(num)
        
        stab_low(num) = &
             & molwt(num) * ( grav ** 2.0d0 ) &
             & * ( lambda(num) / ( gasr * temp(num) ) - Cp(num) / gasr ) &
             & * ( lambda(num) / ( Cp(num) * temp(num) ) &
             &     + ( Mv(num) - Md(num) ) / molwt(num) ) &
             & * molfr(num,w) &
             & / ( Cp(num) * Temp(num) )
        
        stab_high(num) = &
             & ( molwt(num) * ( grav ** 2.0d0 ) &
             &   * ( 1.0d0 - ( ( Cp(num) * Temp(num) ) / lambda(num)  ) ) &
             &   / ( Cp(num) * Temp(num) ) ) & 
             & - grav * dMdz(num) / molwt(num) 
     end if
  end do

  where ( stab_low < 1.0d-7 ) 
     stab_low = 1.0d-7
  end where
  where ( stab_high < 1.0d-7 ) 
     stab_high = 1.0d-7
  end where


  
!!!
!!!雲量の計算
!!!
  do snum = 1, smax_gas  
     cldns(:,snum) = abs((property(snum)%molwt * dX(:, snum) * pres**2.0D0) &
          & / (gasr * temp * dp))      ! * 1.0D-3
  end do

  
!!!
!!! dtdz の単位を K/km へ変換
!!!
  dtdz      = dtdz     * ( - 1.0d3 )
  dtdz_dry  = dtdz_dry * ( - 1.0d3 )

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


!!!
!!!物理量の書き出し
!!!
  
  call HistoryPut('entr', entr(2:step-1) )
  call HistoryPut('temp', temp(2:step-1) )
  call HistoryPut('mol', mol(2:step-1, :) )
  call HistoryPut('emol', emol(2:step-1, :) )
  call HistoryPut('molfr', molfr(2:step-1, :) )
  call HistoryPut('molwt', molwt(2:step-1) )
  call HistoryPut('dtdz', dtdz(2:step-1) )
  call HistoryPut('dtdz_dry', dtdz_dry(2:step-1) )
  call HistoryPut('cp', cp(2:step-1) )
  call HistoryPut('cldns', cldns(2:step-1, :) )
  call HistoryPut('stab', stab(2:step-1) )
  call HistoryPut('stab_low', stab_low(2:step-1) )
  call HistoryPut('stab_high', stab_high(2:step-1) )
  call HistoryPut('lambda', lambda(2:step-1) )
  call HistoryPut('N',      sqrt( stab(2:step-1) ) )
  

!!!
!!!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)
    i   = step - 2
    
    !大域属性, 座標軸の設定
    !"size" は化学種名を入れるためのもの
    call HistoryCreate(&
         & file=AtmosFile, &
         & title="Equilibrium Calculation", &
         & source=dataset(:num+1), &
         & institution="GFD_Dennou_Club oboro Project", &
         & dims=(/'pres', 'spc ', 'elm ', 'size', 'z   '/), &
         & dimsizes=(/i, smax, emax, 20, i/), &
         & longnames=(/"Pressure        ", "chemical species", &
         &             "element         ", "size of name    ", &
         &             "height          " /), &
         & units=(/"Pa", "1 ", "1 ", "1 ", "km"/), &
         & 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(2:step-1), 4))
    call HistoryAddAttr( 'pres', 'positive', 'down' )

    sizes = (/(i, i = 1, 20)/)
    call HistoryPut('size', real(sizes, 4))    

    call HistoryPut('z', real(height(2:step-1) * 1.0d-3, 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=(/'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="dtdz", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz_dry", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & 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" )

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


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


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

    !変数定義, 潜熱
    call HistoryAddVariable( &
         & varname="lambda", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )

    !変数定義, ブラントバイサラ振動数
    call HistoryAddVariable( &
         & varname="N", &           
         & dims=(/'z'/), &
         & longname="Brunt-Vaisala Frequency", &
         & units='s^-1', &
         & xtype="double" )

  end subroutine output_gtool4_init
  
end program OboroAtmos
