!= Program ArareCAPE
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: arare-cape.f90,v 1.2 2010-03-05 03:54:28 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
! 非静力学モデル deepconv/arare 2 次データ作成プログラム
!
!
program ArareCAPE

  !----- モジュール読み込み ------

  !-----   型宣言, 文字列処理   ----
  use dc_types,       only : STRING
  use dc_string,      only : StoA

  !-----   メッセージ出力   -----
  use dc_message,     only: MessageNotify

  !  コマンドライン引数解釈
  use argset,        only : argset_init  

  !-----    管理モジュール   -----
  !  化学量計算モジュール
  use chemcalc, only: chemcalc_init, xz_SvapPress, xz_LatentHeat, ReactHeatNH4SH
  use chemdata, only: chemdata_init

  !  入出力ファイル名管理モジュール
  use fileset,       only : fileset_init

  !  デバッグ出力管理モジュール
  use debugset,      only : debugset_init

  !  時間管理モジュール
  use timeset,       only : timeset_init, TimeDisp, TImeInt

  !  格子点管理モジュール 
  use gridset,       only : gridset_init, DelZ, &
    &                       DimXMin, DimXMax, DimZMin, DimZMax, SpcNum, &
    &                       RegXMax, RegZMax, RegXMin, RegZMin, s_Z, MarginZ

  !  基本場設定モジュール
  use basicset,      only : basicset_init,  &
    &                       MolWtDry, MolWtWet, CvDry, CpDry, PressBasis, GasRDry,    &
    &                       SpcWetID, Tropopause, &
    &                       xza_MixRtBasicZ, xz_PotTempBasicZ, xz_ExnerBasicZ, xz_DensBasicZ, &
    &                       xz_TempBasicZ, xz_PressBasicZ, Grav

  !  湿潤ルーチン設定モジュール
  use moistset,      only: moistset_init,   &
    &                      GasNum, CloudNum, RainNum, CondNum, RactNum, &
    &                      IdxG,   IdxC,   IdxR,                        &
    &                      IdxCG,  IdxCC,  IdxCR,                       &
    &                      IdxNH3, IdxH2S, IdxNH4SHc, IdxNH4SHr         

  !-----    下請けモジュール   -----
  !微分平均演算
  use average

  !  断熱上昇気塊の温度減率計算用モジュール
  use ECCM,         only : eccm_wet2
  
  ! 境界条件
  use boundary,      only : BoundaryXCyc_xz, BoundaryZSym_xz

  !暗黙の型宣言禁止
  implicit none
  
  !変数定義
  real(8), allocatable :: xz_PotTemp(:,:)
  real(8), allocatable :: xz_PotTempAll(:,:)
  real(8), allocatable :: z_PotTempAll(:)
  real(8), allocatable :: z_VPotTempAll(:)
  real(8), allocatable :: z_DVPotTempAllDt(:)
  real(8), allocatable :: z_DPotTempAllDt(:)
  real(8), allocatable :: z_PotTempMoist(:,:)
  real(8), allocatable :: z_VPotTempMoist(:,:)
  real(8), allocatable :: z_VPotTempDown(:,:)

  real(8), allocatable :: xz_Temp(:,:)
  real(8), allocatable :: xz_TempAll(:,:)
  real(8), allocatable :: z_TempAll(:)
  real(8), allocatable :: z_TempBasicZ(:)
  real(8), allocatable :: z_TempMoist(:)

  real(8), allocatable :: z_VTempAll(:)
  real(8), allocatable :: z_VTempDown(:,:)
  real(8), allocatable :: z_VTempMoist(:,:)
  real(8), allocatable :: z_VTempMoist2(:,:)
  real(8), allocatable :: z_DelTempMoist(:,:)
  real(8), allocatable :: z_DelPotTempMoist(:,:)
  real(8), allocatable :: z_DelVTempMoist(:,:)

  real(8), allocatable :: xz_Exner(:,:)
  real(8), allocatable :: z_ExnerAll(:)
  real(8), allocatable :: z_ExnerMoist(:)
  real(8), allocatable :: z_DensMoist(:)
  real(8), allocatable :: z_DensBasicZ(:)
  real(8), allocatable :: z_CAPE(:,:)
  real(8), allocatable :: CAPE(:)
  real(8), allocatable :: z_CloudWorkFunc(:,:)
  real(8), allocatable :: z_CloudWorkFuncTemp(:,:)
  real(8), allocatable :: z_CloudWorkFuncMolWt(:,:)
  real(8), allocatable :: z_CloudWorkFuncMolWt2(:,:)
  real(8), allocatable :: z_CloudWorkFunc2(:,:)
  real(8), allocatable :: z_CloudWorkFunc3(:,:)
  real(8), allocatable :: CloudWorkFunc(:)
  real(8), allocatable :: CloudWorkFuncTemp(:)
  real(8), allocatable :: CloudWorkFuncMolWt(:)  
  real(8), allocatable :: CloudWorkFuncMolWt2(:)  
  real(8), allocatable :: CloudWorkFunc2(:)  
  real(8), allocatable :: CloudWorkFunc3(:)  

  real(8), allocatable :: xza_MixRt(:,:,:)
  real(8), allocatable :: xza_MixRtAll(:,:,:)
  real(8), allocatable :: za_MixRtAll(:,:)
  real(8), allocatable :: za_MixRtMax(:,:)
  real(8), allocatable :: za_MixRtMoist(:,:,:)
  real(8), allocatable :: za_MixRtMoist2(:,:,:)
  real(8), allocatable :: za_MolFrMoist(:,:,:)
  real(8), allocatable :: za_MolFrAll(:,:)

  real(8)              :: AnalTime
  real(8)              :: Hum
  real(8), allocatable :: Altitude(:)
  real(8), allocatable :: Altitude2(:)

  real(8), allocatable :: xz_PressAll(:,:)
  real(8), allocatable :: z_PressAll(:)
  real(8), allocatable :: z_PressMoist(:)

  real(8), allocatable :: z_EffMolWt(:)
  real(8), allocatable :: z_EffMolWtMoist(:)
  real(8), allocatable :: z_EffMolWtMoist2(:)
  real(8), allocatable :: za_MixRtDivMolWt(:,:)
  real(8), allocatable :: z_LatentHeat(:)

  real(8), allocatable :: z_MoistCond(:)
  real(8), allocatable :: z_MoistNH4SH(:)
  real(8), allocatable :: z_DryENG(:,:)
  real(8), allocatable :: z_MoistENG(:,:)

  real(8), allocatable :: z_MolWtAll(:)
  real(8), allocatable :: z_MolWtMoist(:)
  real(8), allocatable :: z_StabTempAll(:,:)
  real(8), allocatable :: z_StabTempMoist(:,:)
  real(8), allocatable :: z_StabMolWtAll(:,:)
  real(8), allocatable :: z_StabMolWtMoist(:,:)
  real(8), allocatable :: z_StabAll(:,:)
  real(8), allocatable :: z_StabMoist(:,:)

  integer              :: Idx
  integer              :: Idx2
  integer              :: Idx3
  integer              :: t, s, i, tt, step, j, l, m, k
  integer, parameter   :: jmax=40
!  integer, parameter   :: jmax=6
!  integer, parameter   :: jmax=1
  character(100)       :: cfgfile
  character(100)       :: info

  !コマンドライン引数の解釈
  !  NAMELIST ファイル名の読み込み
  call argset_init(cfgfile)

  !物質特性の初期化
  call chemdata_init()

  !時刻に関する設定の初期化
  !  NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. 
  call timeset_init(cfgfile)

  !格子点情報の初期化
  !  NAMELIST から情報を得て, 格子点を計算する
  call gridset_init(cfgfile)
  
  !化学計算ルーチンの初期化
  call chemcalc_init()

  !基本場の情報の初期化
  !  NAMELIST から情報を得て, 基本場を設定する.
  call basicset_init(cfgfile)
  
  !I/O ファイル名の初期化
  !  NAMELIST ファイル名を指定し, deepconv/arare の
  !  出力ファイル名を NAMELIST から得る
  call fileset_init(cfgfile)
  
  !基本場の取得
  call AnalFile_BasicZ_Get( )

  !湿潤ルーチンの共有変数の初期化
  call moistset_init()

  !内部変数の初期化. とりあえずゼロを入れて値を確定させておく. 
  call ArareAlloc

  ! 高度を入力
!  altitude(1:jmax) = (/0.0d3,50.0d3,60.0d3,70.0d3,75.0d3,80.0d3,85.0d3,90.0d3,95.0d3,100.0d3,105.0d3,110.0d3,115.0d3,120.0d3/)
!  altitude(1:jmax) = (/111.0d3,113.0d3,115.0d3,117.0d3,119.0d3,121.0d3/)
!  altitude(1:jmax) = (/0.0d3/)

  SearchIdx: do m = 1, DimZMax
     if (sum(xza_MixRtBasicZ(1,m,:),1) < sum(xza_MixRtBasicZ(1,1,:),1) * 0.95d0) then 
        Idx = m
        exit SearchIdx
     end if
  end do SearchIdx

  SearchIdx3: do m = 1, DimZMax
     if (xza_MixRtBasicZ(1,m,IdxNH3) < xza_MixRtBasicZ(1,1,IdxNH3) * 0.5d0) then 
        Idx3 = m
        exit SearchIdx3
     end if
  end do SearchIdx3

  do m = 1, jmax 
     Altitude(m) = s_Z(Idx + m - jmax / 2)
  end do
  write(*,*) "####### Alt ", Altitude, "########"


  !NetCDF ファイルの作成と基本場の取得  
  call AnalFile_MoistENG_Open(  )
    
  !----------------------------------------------------------------------
  ! 解析値の計算開始
  !   時刻を進めながら, ファイルの値を得る
  !----------------------------------------------------------------------
!  do t = 1, int( TimeInt / TimeDisp ) + 1
  step = int( TimeInt / TimeDisp )
  STEPS: do while (step >= 600)
   write(*,*) step
   if (mod(step, 2) == 1) then 
       exit STEPS
    end if
    step = step / 2 
  end do STEPS
  write(*,*) "**** STEPS ***** ", int( TimeInt / TimeDisp ), " ==> ", step

  TIMESTEP: do tt = 0, int( TimeInt / TimeDisp) -1
     if ( mod(tt, int(TimeInt / TimeDisp / step)) /= 0 ) then 
        cycle TIMESTEP
     end if
     t = tt + 1
     
    !----------------------------------------------------------------
    ! ヒストリファイルを開き値を得る. 
    !----------------------------------------------------------------
    call AnalFile_Get( t, AnalTime, xz_PotTemp, xz_Exner, xza_MixRt )
    write(*,*) t, AnalTime    
    
    !----------------------------------------------------------------
    ! 初期化
    !----------------------------------------------------------------
    xz_PressAll = PressBasis * ((xz_Exner + xz_ExnerBasicZ ) ** (CpDry / GasRDry))
    z_PressAll = z_MeanX_xz( xz_PressAll )

    xz_PotTempAll = xz_PotTemp + xz_PotTempBasicZ
    z_PotTempAll  = z_MeanX_xz( xz_PotTempAll )
    
    xz_TempAll = ( xz_Exner + xz_ExnerBasicZ ) * xz_PotTempAll   
    z_TempAll = z_MeanX_xz( xz_TempAll )
    z_TempBasicZ = z_MeanX_xz( xz_TempBasicZ )
    z_ExnerAll = z_MeanX_xz( xz_Exner + xz_ExnerBasicZ ) 
    z_DensBasicZ = z_MeanX_xz( xz_DensBasicZ )

    xza_MixRtAll = xza_MixRt + xza_MixRtBasicZ
    za_MixRtAll  = za_MeanX_xza( xza_MixRtAll )    

    !----------------------------------------------------------------
    ! 密度, 仮温度の計算
    !----------------------------------------------------------------
    za_MixRtDivMolWt = 0.0d0
    do s = 1, GasNum
       za_MixRtDivMolWt(:,IdxG(s)) = za_MixRtAll(:,IdxG(s)) / MolWtWet(IdxG(s))
    end do

    z_EffMolWt = &
      & (1.0d0 + sum(za_MixRtAll,2)) &
!      & (1.0d0 + sum(za_MixRtAll(:,1:GasNum),2)) &
      & / (MolWtDry * ((1.0d0 / MolWtDry) + sum(za_MixRtDivMolWt,2))) 
    
    z_VTempAll = z_TempAll / z_EffMolWt
    z_VPotTempAll = z_PotTempAll / z_EffMolWt

    z_DVPotTempAllDt = 0.0d0
    do k = DimZMin+1, DimZMax-1
       z_DVPotTempAllDt(k) = (z_VPotTempAll(k+1) - z_VPotTempAll(k-1)) / DelZ * 0.5d0
    end do

    z_DPotTempAllDt = 0.0d0
    do k = DimZMin+1, DimZMax-1
       z_DPotTempAllDt(k) = (z_PotTempAll(k+1) - z_PotTempAll(k-1)) / DelZ * 0.5d0 
    end do

    za_MixRtMax = 0.0d0
    do s = 1, SpcNum
       do k = DimZMin+1, DimZMax-1
          za_MixRtMax(k,s) = maxval( xza_MixRtAll(:,k,s))
       end do
    end do

    
    !---------------------------------------------------------------
    ! CAPE 
    !---------------------------------------------------------------    
    do j = 1, jmax
       ! 入力された高度に対応する座標添字を求める
       Idx  = minloc(s_Z, 1, s_Z > Altitude(j)) - 2 - MarginZ
       Idx2 = minloc(s_Z, 1, s_Z > Tropopause) - 2 - MarginZ  !対流圏
       Altitude2(j) = s_Z(Idx)
       
       write(*,*) "#######", Idx, Idx2, "###########"
       write(*,*) "#######", Altitude(j), "=>", Altitude2(j),  "###########"
       write(*,*) "#######", Tropopause,  "=>", s_Z(Idx2), "###########"

       ! 湿潤断熱線の計算
       Hum = 1.0d0
       z_TempMoist    = z_TempAll 
       z_PressMoist   = z_PressAll
       do s = 1, SpcNum
          za_MolFrAll(:,s) = z_MeanX_xz(xza_MixRtAll(:,:,s)) / MolWtWet(s) * MolWtDry 
          za_MolFrMoist(:,s,j) = za_MolFrAll(:,s)
       end do
       
       !湿潤断熱線を引く
       call ECCM_Wet2(Idx, Hum, z_TempMoist, z_PressMoist, za_MolFrMoist(:,:,j))
       
       !モル比を混合比に変換
       za_MixRtMoist(:,:,j) = 0.0d0
       do s = 1, SpcNum
          za_MixRtMoist(:,s,j) = za_MolFrMoist(:,s,j) * MolWtWet(s) / MolWtDry
       end do
       ! Moist2 では, 雲と雨の量は最大値を利用
       za_MixRtMoist2(:,:,j) = za_MixRtMoist(:,:,j)
       do s = GasNum+1, SpcNum
          za_MixRtMoist2(:,s,j) = za_MixRtMax(:,s)
       end do
       ! 雲の量は蒸気の減少分. 雨は取りあえずゼロ.
       do s = 1, GasNum
          do k = DimZMin+1, DimZMax-10
             za_MixRtMoist(k,IdxCC(s),j) = za_MixRtMoist(k-1,IdxCG(s),j) - za_MixRtMoist(k,IdxCG(s),j)
          end do
       end do
       
       !気塊の分子量効果
       za_MixRtDivMolWt = 0.0d0
       do s = 1, GasNum
          za_MixRtDivMolWt(:,IdxG(s)) = za_MixRtMoist(:,IdxG(s),j) / MolWtWet(IdxG(s))
       end do

       z_EffMolWtMoist = &
            & (1.0d0 + sum(za_MixRtMoist(:,:,j),2)) &
!            & (1.0d0 + sum(za_MixRtMoist(:,1:GasNum,j),2)) &
            & / (MolWtDry * ((1.0d0 / MolWtDry) + sum(za_MixRtDivMolWt,2))) 
       z_EffMolWtMoist2 = &
            & (1.0d0 + sum(za_MixRtMoist2(:,:,j),2)) &
!            & (1.0d0 + sum(za_MixRtMoist2(:,1:GasNum,j),2)) &
            & / (MolWtDry * ((1.0d0 / MolWtDry) + sum(za_MixRtDivMolWt,2))) 
       
       !気塊の仮温度
       z_VTempMoist(:,j)    = z_TempMoist / z_EffMolWtMoist
       z_VTempMoist2(:,j)   = z_TempMoist / z_EffMolWtMoist2
       z_DelTempMoist(:,j)  = z_TempMoist - z_TempAll
       z_DelVTempMoist(:,j) = z_VTempMoist(:,j) - z_VTempAll
       
       !気塊の温位. 圧力は周囲と同じ = 乾燥成分でおおよそ決まる
       z_PotTempMoist(:,j)    = z_TempMoist * (PressBasis / z_PressAll) ** (GasRDry / CpDry) 
       z_DelPotTempMoist(:,j) = z_PotTempMoist(:,j) - z_PotTempAll
       z_VPotTempMoist(:,j)   = z_PotTempMoist(:,j) / z_EffMolWtMoist
       
       !気塊のエクスナー関数
       z_ExnerMoist = z_TempMoist / z_PotTempMoist(:,j)
       
!       !気塊の密度
!       z_DensMoist = &
!            & PressBasis * (z_ExnerMoist ** (CvDry / GasRDry)) &
!            &  / (GasRDry * z_PotTempMoist(:,j) / z_EffMolWtMoist)
       
       ! H2O と NH3 の凝結の潜熱. 
       z_LatentHeat = z_MeanX_xz( xz_LatentHeat( SpcWetID(IdxCC(1)), xz_MeanX_z(z_TempMoist) ) )
       
       !CAPE
       z_CAPE(:,j)      = Grav * (z_VTempMoist(:,j) - z_VTempAll) / z_VTempAll 
       ! まともな計算
       z_CloudWorkFunc(:,j)  = z_DensBasicZ * Grav * (z_VTempMoist(:,j)  - z_VTempAll) / z_VTempAll 
       ! 計算結果から得られた雲の量を利用
       z_CloudWorkFunc2(:,j) = z_DensBasicZ * Grav * (z_VTempMoist2(:,j) - z_VTempAll) / z_VTempAll 
       ! 温度の寄与のみ
       z_CloudWorkFuncTemp(:,j) = (z_TempMoist - z_TempAll) / z_TempAll * z_DensBasicZ * Grav 
       ! 分子量の寄与のみ
       z_CloudWorkFuncMolWt(:,j) = - (sum(za_MixRtMoist(:,1:SpcNum,j), 2) - sum(za_MixRtAll(:,1:GasNum), 2)) * z_DensBasicZ * Grav        
       ! 分子量の寄与のみ Moist2 を利用
       z_CloudWorkFuncMolWt2(:,j)= - (sum(za_MixRtMoist2(:,1:SpcNum,j),2) - sum(za_MixRtAll(:,1:GasNum), 2)) * z_DensBasicZ * Grav        
       
       CAPE(j) = 0
       CloudWorkFunc(j) = 0       
       CloudWorkFunc2(j) = 0       
       CloudWorkFuncTemp(j) = 0       
       CloudWorkFuncMolWt(j) = 0       
       CloudWorkFuncMolWt2(j) = 0       
       do l = RegZMin+1, Idx2-2
          CAPE(j) = CAPE(j) + z_CAPE(l,j) * DelZ 
          CloudWorkFunc(j) = CloudWorkFunc(j) + (z_CloudWorkFunc(l,j) + z_CloudWorkFunc(l+1,j)) * DelZ * 5.0d-1
          CloudWorkFunc2(j)= CloudWorkFunc2(j)+ (z_CloudWorkFunc2(l,j)+ z_CloudWorkFunc2(l+1,j))* DelZ * 5.0d-1
          CloudWorkFuncTemp(j) = CloudWorkFuncTemp(j) + (z_CloudWorkFuncTemp(l,j)+z_CloudWorkFuncTemp(l+1,j)) * DelZ * 5.0d-1 
          CloudWorkFuncMolWt(j) = CloudWorkFuncMolWt(j) + (z_CloudWorkFuncMolWt(l,j) + z_CloudWorkFuncMolWt(l+1,j)) * DelZ * 5.0d-1
          CloudWorkFuncMolWt2(j)= CloudWorkFuncMolWt2(j)+ (z_CloudWorkFuncMolWt2(l,j)+ z_CloudWorkFuncMolWt2(l+1,j)) * DelZ * 5.0d-1
       end do
       
       !---------------------------------------------------------------
       ! 湿潤静的エネルギー
       !---------------------------------------------------------------
       ! 乾燥静的エネルギー CpT+gz
       z_DryEng(:,j) = CpDry * z_TempMoist + Grav * s_Z

       !NH4SH の生成反応. H2S は全て NH3 と反応して NH4SH になると考える. 
       ! H2S と反応した分だけ NH3 を減少させる. 
       z_MoistNH4SH = ReactHeatNH4SH * za_MixRtMoist(:,IdxH2S,j) * MolWtWet(IdxNH4SHc) / MolWtWet(IdxH2S)
       za_MixRtMoist(:,IdxNH3,j) = za_MixRtMoist(:,IdxNH3,j) -  za_MixRtMoist(:,IdxH2S,j) &  
            &          * MolWtWet(IdxNH3) / MolWtWet(IdxH2S)  
       ! H2O と NH3 の凝結の潜熱. 
       z_MoistCond = 0.0d0
       do s = 1, CondNum
          z_LatentHeat = z_MeanX_xz( xz_LatentHeat( SpcWetID(IdxCC(s)), xz_MeanX_z(z_TempMoist) ) )
          z_MoistCond = z_MoistCond + z_LatentHeat * za_MixRtMoist(:,IdxCG(s),j) 
       end do
       z_MoistEng(:,j) = z_DryEng(:,j) + z_MoistCond + z_MoistNH4SH

       !NH4SH の蒸発熱. 下部境界から持ち上げた気塊の NH4SH 凝結量を利用. 
!       z_MoistNH4SH = ReactHeatNH4SH * (za_MixRtMoist(:,IdxNH4SHc,1) + za_MixRtMoist(:,IdxNH4SHr,1))
       z_MoistNH4SH = ReactHeatNH4SH * (za_MixRtMax(:,IdxNH4SHc) + za_MixRtMax(:,IdxNH4SHr))
       ! H2O と NH3 の蒸発. 下部境界から持ち上げた気塊の凝結量を利用. 
       z_MoistCond = 0.0d0
       do s = 1, CondNum
          z_LatentHeat = z_MeanX_xz( xz_LatentHeat( SpcWetID(IdxCC(s)), xz_MeanX_z(z_TempMoist) ) )
!          z_MoistCond = z_MoistCond + z_LatentHeat * (za_MixRtMoist(:,IdxCC(s),1) + za_MixRtMoist(:,IdxCR(s),1))
          z_MoistCond = z_MoistCond + z_LatentHeat * (za_MixRtMax(:,IdxCC(s)) + za_MixRtMax(:,IdxCR(s)))
       end do

       !とある高度の温位で fix. さらに凝結物 (下部境界の気塊を上昇させた場合) を蒸発させる. 
       do k = DimZMin, DimZMax
!          z_VTempDown(k,j) = (z_PotTempAll(Idx) * z_ExnerAll(k) - (z_MoistNH4SH(k) + z_MoistCond(k)) / CpDry) / z_EffMolWt(k)
          z_VTempDown(k,j) = z_VTempAll(Idx)
       end do
!       z_CloudWorkFunc3(:,j) = z_DensBasicZ * Grav * (z_VTempDown(:,j)  - z_VTempAll) / z_VTempAll 
       z_CloudWorkFunc3(:,j) = z_DensBasicZ * Grav * (z_VTempAll(Idx)  - z_VTempAll) / z_VTempAll 
       !積分はとある高度から下で計算
       CloudWorkFunc3(j) = 0       
       do l = RegZMin, Idx
          CloudWorkFunc3(j)= CloudWorkFunc3(j)+ (z_CloudWorkFunc3(l,j)+ z_CloudWorkFunc3(l+1,j)) * DelZ * 5.0d-1
       end do
!       z_VPotTempDown(:,j) = z_VPotTempAll - (z_MoistNH4SH(Idx) - z_MoistCond(Idx)) / CpDry / z_ExnerAll(Idx) / z_EffMolWt(Idx)
!       z_VPotTempDown(:,j) = - (z_MoistNH4SH(Idx) - z_MoistCond(Idx)) / CpDry / z_ExnerAll(Idx) / z_EffMolWt(Idx)
       z_VPotTempDown(:,j) = - (z_MoistNH4SH(:) + z_MoistCond(:)) / CpDry / z_ExnerAll(:) / z_EffMolWt(:)
!       CloudWorkFunc3(j) = 0       
!       CloudWorkFunc3(j)= sum(z_VPotTempDown(RegZMin:Idx,j), 1)

       !----------------------------------------------------------------
       ! 安定度
       !----------------------------------------------------------------
       z_MolWtAll = 0.0d0
       z_MolWtMoist = 0.0d0
       do k = DimZMin, DimZMax
          z_MolWtMoist(k) = dot_product(MolWtWet(1:GasNum), za_MolFrMoist(k,1:GasNum,j))
          z_MolWtAll(k)   = dot_product(MolWtWet(1:GasNum), za_MolFrAll(k,1:GasNum))
       end do

       do k = DimZMin + 1, DimZMax - 1
          z_StabTempAll(k,j) =                                                &
               &  Grav / z_TempAll(k)                                        &
               &  * ( (z_TempAll(k+1) - z_TempAll(k-1)) * 5.0d-1 / DelZ &
!               &     + Grav * z_EffMolWt(k) / CpDry )
               &     + Grav / CpDry )
          z_StabTempMoist(k,j) =                                           &
               & Grav / z_TempMoist(k)                        &
               &  * ( (z_TempMoist(k+1) - z_TempMoist(k-1)) * 5.0d-1 / DelZ &
!               &     + Grav * z_EffMolWtMoist(k) / CpDry ) 
               &     + Grav / CpDry ) 
          z_StabMolWtAll(k,j) =                                          &
!               & - Grav * (z_MolWtAll(k+1) - z_MolWtAll(k-1)) * 5.0d-1 / DelZ / z_MolWtAll(k)
               & - Grav * (z_MolWtAll(k+1) - z_MolWtAll(k-1)) * 5.0d-1 / DelZ / MolWtDry
          z_StabMolWtMoist(k,j) =                                          &
!               & - Grav * (z_MolWtMoist(k+1) - z_MolWtMoist(k-1)) * 5.0d-1 / DelZ / z_MolWtMoist(k)
               & - Grav * (z_MolWtMoist(k+1) - z_MolWtMoist(k-1)) * 5.0d-1 / DelZ / MolWtDry
       end do
       z_StabTempAll(DimZMin,j) = z_StabTempAll(DimZMin+1,j) 
       z_StabTempMoist(DimZMax,j) = z_StabTempMoist(DimZMax-1,j)
       z_StabMolWtAll(DimZMin,j) = z_StabMolWtAll(DimZMin+1,j) 
       z_StabMolWtMoist(DimZMax,j) = z_StabMolWtMoist(DimZMax-1,j)
       
       z_StabMoist(:,j) = z_StabTempMoist(:,j) + z_StabMolWtMoist(:,j)
       z_StabAll(:,j)   = z_StabTempAll(:,j) + z_StabMolWtAll(:,j)
       
       where (z_StabMoist < 1.0d-7) 
          z_StabMoist = 1.0d-7
       end where
       where (z_StabAll < 1.0d-7) 
          z_StabAll = 1.0d-7
       end where
       
    end do
    !----------------------------------------------------------------    
    ! ファイル出力
    !----------------------------------------------------------------
    call AnalFile_MoistENG_OutPut( AnalTime )
    
 end do TIMESTEP
  
  !----------------------------------------------------------------    
  ! ファイルを閉じる
  !----------------------------------------------------------------    
  call AnalFile_Close

    
contains
  
  subroutine AnalFile_Close( )
    !
    !解析ファイルのクローズ
    !
    use gt4_history, only: HistoryClose
    use fileset, only: gt_hist

    call HistoryClose
    
  end subroutine AnalFile_Close

  
  subroutine AnalFile_MoistENG_Open( )
    !
    !解析ファイルの定義
    !
    
    use gt4_history, only: HistoryCreate, HistoryPut, HistoryAddVariable
    use fileset, only: gt_hist, exptitle, expsrc, expinst, HistoryFilePrefix    
    use gridset, only: FileNX, FileNZ, s_X, s_Z, &
      &                FileXMin, FileXMax, FileZMin, FileZMax
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    character(150)  :: File
    integer, allocatable:: spc(:)
    integer :: s

    !ファイル名の指定
    File = trim(HistoryFilePrefix)// '_CAPE.nc'
    write(*,*) "******* OUTPUT FILE: ", File

    !-----------------------------------------------------------
    ! ヒストリー作成
    !-----------------------------------------------------------
    call HistoryCreate(                              &
      & file = File,                                 &
      & title = exptitle,                            &
      & source = expsrc,                             &
      & institution = expinst,                       &
      & dims=(/'z  ','s  ','alt','t  '/),            &
      & dimsizes=(/FileNZ, SpcNum, jmax, 0/),                &
      & longnames=(/'Z-coordinate',                  &
      &             'Species     ',                  &
      &             'Altitude    ',                  &
      &             'Time        '/),                &
      & units=(/'m','1','m','s'/), origin=0.0,       &
      & interval=0.0 )
    
    !-----------------------------------------------------------  
    ! 軸の出力
    !-----------------------------------------------------------
    allocate( spc(SpcNum) )
    do s = 1, SpcNum
       spc(s) = s
    end do
    write(*,*) spc
    call HistoryPut('s  ', spc( 1:SpcNum ))
    call HistoryPut('alt', Altitude)
    call HistoryPut('z  ', s_Z( FileZMin: FileZMax ))

    !-----------------------------------------------------------  
    ! 解析用の変数の出力
    !-----------------------------------------------------------  
    call HistoryAddVariable(                                   &
      & varname='CAPEMolFr', dims=(/'z  ','s  ','alt','t  '/), &
      & longname='Mole Fraction (air parcel)',                 &
      & units='mol/mol', xtype='double' )

    call HistoryAddVariable(                                   &
      & varname='CAPEMixRt', dims=(/'z  ','s  ','alt','t  '/), &
      & longname='Mixing Ratio (air parcel)',           &
      & units='kg/kg', xtype='double' )
  
    call HistoryAddVariable(                                &
      & varname='CAPEPotTemp', dims=(/'z  ','alt','t  '/),  &
      & longname='Potential Temperature (air parcel)',      &
      & units='K', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='CAPEVPotTemp', dims=(/'z  ','alt','t  '/), &
      & longname='Virtual Potential Temperature (air parcel)', &
      & units='K', xtype='double' )

    call HistoryAddVariable(                                   &
      & varname='CAPEVTemp', dims=(/'z  ','alt','t  '/),       &
      & longname='Virtual Temperature (air parcel)',           &
      & units='K', xtype='double' )

    call HistoryAddVariable(                                   &
      & varname='CAPEVPotTempDown', dims=(/'z  ','alt','t  '/),       &
      & longname='Virtual Potential Temperature (air parcel)',           &
      & units='K', xtype='double' )

    call HistoryAddVariable(                                   &
      & varname='CAPEDelTemp', dims=(/'z  ','alt','t  '/),       &
      & longname='Del Potential Temperature (air parcel)',           &
      & units='K', xtype='double' )

    call HistoryAddVariable(                                   &
      & varname='CAPEDelPotTemp', dims=(/'z  ','alt','t  '/),       &
      & longname='Del Temperature (air parcel)',           &
      & units='K', xtype='double' )

    call HistoryAddVariable(                                   &
      & varname='CAPEDelVTemp', dims=(/'z  ','alt','t  '/),       &
      & longname='Del Virtual Temperature (air parcel)',           &
      & units='K', xtype='double' )

    call HistoryAddVariable(                                   &
      & varname='CAPEDVPotTempAllDt', dims=(/'z  ','t  '/),       &
      & longname='Gradient of Virtual Potential Temperature',           &
      & units='K/m', xtype='double' )

    call HistoryAddVariable(                                   &
      & varname='CAPEDPotTempAllDt', dims=(/'z  ','t  '/),       &
      & longname='Gradient of Potential Temperature',           &
      & units='K/m', xtype='double' )

    call HistoryAddVariable(                                         &
      & varname='CAPEMoistStaticEnergy', dims=(/'z  ','alt','t  '/), &
      & longname='Moist Static Energy (air parcel)',                 &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                       &
      & varname='CAPEDryStaticEnergy', dims=(/'z  ','alt','t  '/), &
      & longname='Dry Static Energy (air parcel)',                 &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='CAPE', dims=(/'z  ','alt','t  '/),        &
      & longname='Grav * DTv /Tv',                   &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='CAPEInt', dims=(/'alt', 't  '/),           &
      & longname='CAPE',                                    &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='CloudWorkFunc', dims=(/'z  ','alt','t  '/),        &
      & longname='Grav * rho * DTv /Tv',                   &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='CloudWorkFuncInt', dims=(/'alt', 't  '/),           &
      & longname='int Grav * rho * DTv /Tv',                                    &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='CloudWorkFuncTemp', dims=(/'z  ','alt','t  '/),        &
      & longname='Grav * rho * DTv /Tv (Temp)',                   &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='CloudWorkFuncTempInt', dims=(/'alt', 't  '/),           &
      & longname='int Grav * rho * DTv /Tv (Temp)',                                    &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='CloudWorkFuncMolWt', dims=(/'z  ','alt','t  '/),        &
      & longname='Grav * rho * DTv /Tv (MolWt)',                   &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='CloudWorkFuncMolWtInt', dims=(/'alt', 't  '/),           &
      & longname='int Grav * rho * DTv /Tv (MolWt)',                                    &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='CloudWorkFuncMolWt2', dims=(/'z  ','alt','t  '/),        &
      & longname='Grav * rho * DTv /Tv (MolWt)',                   &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='CloudWorkFuncMolWtInt2', dims=(/'alt', 't  '/),           &
      & longname='int Grav * rho * DTv /Tv (MolWt)',                                    &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='CloudWorkFunc2', dims=(/'z  ','alt','t  '/),        &
      & longname='Grav * rho * DTv /Tv (Temp + MolWt)',                   &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='CloudWorkFuncInt2', dims=(/'alt', 't  '/),           &
      & longname='int Grav * rho * DTv /Tv (Temp + MolWt)',                                    &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='CloudWorkFunc3', dims=(/'z  ','alt','t  '/),        &
      & longname='Grav * rho * DTv /Tv (Temp + MolWt)',                   &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='CloudWorkFuncInt3', dims=(/'alt', 't  '/),           &
      & longname='int Grav * rho * DTv /Tv (Temp + MolWt)',                                    &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='CAPEStabTempAll', dims=(/'z  ','alt','t  '/),        &
      & longname='stability from dT/dz',                   &
      & units='1/s/s', xtype='double' )
    call HistoryAddVariable(                               &
      & varname='CAPEStabTempMoist', dims=(/'z  ','alt','t  '/),        &
      & longname='stability from dT/dz (moist)',                   &
      & units='1/s/s', xtype='double' )
    call HistoryAddVariable(                               &
      & varname='CAPEStabMolWtAll', dims=(/'z  ','alt','t  '/),        &
      & longname='stability from dM/dz',                   &
      & units='1/s/s', xtype='double' )
    call HistoryAddVariable(                               &
      & varname='CAPEStabMolWtMoist', dims=(/'z  ','alt','t  '/),        &
      & longname='stability from dM/dz (moist)',                   &
      & units='1/s/s', xtype='double' )
    call HistoryAddVariable(                               &
      & varname='CAPEStabAll', dims=(/'z  ','alt','t  '/),        &
      & longname='stability',                   &
      & units='1/s/s', xtype='double' )
    call HistoryAddVariable(                               &
      & varname='CAPEStabMoist', dims=(/'z  ','alt','t  '/),        &
      & longname='stability (moist)',                   &
      & units='1/s/s', xtype='double' )
  end subroutine AnalFile_MoistENG_Open


  subroutine AnalFile_MoistENG_OutPut( Time )
    use gridset, only: FileZMin, FileZMax
    use gt4_history, only: HistoryPut

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in) :: Time    
 
    !----------------------------------------------------------------
    ! 値を出力
    !----------------------------------------------------------------
    call HistoryPut( 't', Time ) 
    
    call HistoryPut( 'CAPEMolFr', za_MolFrMoist(FileZMin:FileZMax,1:SpcNum,1:jmax) )

    call HistoryPut( 'CAPEPotTemp', z_PotTempMoist(FileZMin:FileZMax,1:jmax) )

    call HistoryPut( 'CAPEVPotTemp', z_VPotTempMoist(FileZMin:FileZMax,1:jmax) )

    call HistoryPut( 'CAPEVTemp', z_VTempMoist(FileZMin:FileZMax,1:jmax) )

    call HistoryPut( 'CAPEVPotTempDown', z_VPotTempDown(FileZMin:FileZMax,1:jmax) )

    call HistoryPut( 'CAPEDelTemp', z_DelTempMoist(FileZMin:FileZMax,1:jmax) )

    call HistoryPut( 'CAPEDelPotTemp', z_DelPotTempMoist(FileZMin:FileZMax,1:jmax) )

    call HistoryPut( 'CAPEDelVTemp', z_DelVTempMoist(FileZMin:FileZMax,1:jmax) )

    call HistoryPut( 'CAPEMixRt', za_MixRtMoist(FileZMin:FileZMax,1:SpcNum,1:jmax) )

    call HistoryPut( 'CAPEDVPotTempAllDt', z_DVPotTempAllDt(FileZMin:FileZMax) )

    call HistoryPut( 'CAPEDPotTempAllDt', z_DPotTempAllDt(FileZMin:FileZMax) )

    call HistoryPut( 'CAPEDryStaticEnergy', z_DryENG(FileZMin:FileZMax,1:jmax) )

    call HistoryPut( 'CAPEMoistStaticEnergy', z_MoistENG(FileZMin:FileZMax,1:jmax) )

    call HistoryPut( 'CAPE', z_CAPE(FileZMin:FileZMax,1:jmax) )

    call HistoryPut( 'CAPEInt', CAPE(1:jmax) )

    call HistoryPut( 'CloudWorkFunc', z_CloudWorkFunc(FileZMin:FileZMax,1:jmax) )

    call HistoryPut( 'CloudWorkFuncInt', CloudWorkFunc(1:jmax) )

    call HistoryPut( 'CloudWorkFuncTemp', z_CloudWorkFuncTemp(FileZMin:FileZMax,1:jmax) )    

    call HistoryPut( 'CloudWorkFuncTempInt', CloudWorkFuncTemp(1:jmax) )

    call HistoryPut( 'CloudWorkFuncMolWt', z_CloudWorkFuncMolWt(FileZMin:FileZMax,1:jmax) )    

    call HistoryPut( 'CloudWorkFuncMolWtInt', CloudWorkFuncMolWt(1:jmax) )

    call HistoryPut( 'CloudWorkFuncMolWt2', z_CloudWorkFuncMolWt2(FileZMin:FileZMax,1:jmax) )    

    call HistoryPut( 'CloudWorkFuncMolWtInt2', CloudWorkFuncMolWt2(1:jmax) )

    call HistoryPut( 'CloudWorkFunc3', z_CloudWorkFunc3(FileZMin:FileZMax,1:jmax) )    

    call HistoryPut( 'CloudWorkFuncInt3', CloudWorkFunc3(1:jmax) )

    call HistoryPut( 'CloudWorkFunc2', z_CloudWorkFunc2(FileZMin:FileZMax,1:jmax) )    

    call HistoryPut( 'CloudWorkFuncInt2', CloudWorkFunc2(1:jmax) )

    call HistoryPut( 'CAPEStabTempAll', z_StabTempAll(FileZMin:FileZMax,1:jmax) )    

    call HistoryPut( 'CAPEStabTempMoist', z_StabTempMoist(FileZMin:FileZMax,1:jmax) )    

    call HistoryPut( 'CAPEStabMolWtAll', z_StabMolWtAll(FileZMin:FileZMax,1:jmax) )    

    call HistoryPut( 'CAPEStabMolWtMoist', z_StabMolWtMoist(FileZMin:FileZMax,1:jmax) )    

    call HistoryPut( 'CAPEStabAll', z_StabAll(FileZMin:FileZMax,1:jmax) )    

    call HistoryPut( 'CAPEStabMoist', z_StabMoist(FileZMin:FileZMax,1:jmax) )    
   
  end subroutine AnalFile_MoistENG_OutPut


!!!------------------------------------------------------------------------------!!!
  subroutine AnalFile_Get( i, AnalTime, xz_PotTemp, xz_Exner, xza_MixRt )
 
    use dc_string
    use gt4_history,   only : HistoryGet
    use fileset,       only : HistoryFile
    use gridset,       only : DimXMin, DimXMax, DimZMin, DimZMax, SpcNum, &
      &                       FileXMin, FileXMax, FileZMin, FileZMax
    use basicset,      only : SpcWetSymbol
    use boundary,      only : BoundaryXCyc_xz, BoundaryZSym_xz,   &
      &                       BoundaryXCyc_xza, BoundaryZSym_xza, &
      &                       BoundaryXCyc_pz, BoundaryZSym_pz,   &
      &                       BoundaryXCyc_xr, BoundaryZAntiSym_xr

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer, intent(in)  :: i    
    real(8), intent(out) :: AnalTime
    real(8), intent(out) :: xz_Exner(DimXMin:DimXMax,DimZMin:DimZMax)
    real(8), intent(out) :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)
    real(8), intent(out) :: xza_MixRt(DimXMin:DimXMax,DimZMin:DimZMax,1:SpcNum)
    character(30)        :: name               !変数名
    character(10)        :: step
    
    step = 't=^' // adjustl(toChar(i))

    !-------------------------------------------------------------
    !Get a Value from netCDF File
    !-------------------------------------------------------------
    name = "t"
    call HistoryGet( HistoryFile(1), name, AnalTime, step )

    !-------------------------------------------------------------    
    ! Get a Value from netCDF File 
    !-------------------------------------------------------------
    name = "Exner"
    call HistoryGet( HistoryFile(1), name, xz_Exner(FileXMin:FileXMax, FileZMin:FileZMax), range=step )
    call BoundaryXCyc_xz( xz_Exner )
    call BoundaryZSym_xz( xz_Exner )
    
    name = "PotTemp"
    call HistoryGet( HistoryFile(2), name, xz_PotTemp(FileXMin:FileXMax, FileZMin:FileZMax), range=step )
    call BoundaryXCyc_xz( xz_PotTemp )
    call BoundaryZSym_xz( xz_PotTemp )
    
    do s = 1, SpcNum
      name = trim(SpcWetSymbol(s))
      call HistoryGet( HistoryFile(8+s), name, xza_MixRt(FileXMin:FileXMax, FileZMin:FileZMax, s), range=step )
    end do
    call BoundaryXCyc_xza( xza_MixRt )
    call BoundaryZSym_xza( xza_MixRt )

  end subroutine AnalFile_Get


  subroutine AnalFile_BasicZ_Get( )
    use dc_string
    use gt4_history,   only: HistoryGet
    use fileset,       only: HistoryFile
    use basicset,      only: BasicSetArray_Init, SpcWetSymbol
    use gridset,       only: DimXMin, DimXMax, DimZMin, DimZMax, SpcNum, &
      &                      FileXMin, FileXMax, FileZMin, FileZMax
    use boundary,      only: BoundaryXCyc_xz, BoundaryZSym_xz,     &
      &                      BoundaryXCyc_xza,BoundaryZSym_xza, &
      &                      BoundaryXCyc_pz, BoundaryZSym_pz,     &
      &                      BoundaryXCyc_xr, BoundaryZAntiSym_xr

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8) :: xz_DensBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xz_PotTempBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xz_ExnerBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xza_MixRtBZ(DimXMin:DimXMax, DimZMin:DimZMax,1:SpcNum)
    real(8) :: xz_EffMolWtBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xz_PressBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xz_TempBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xz_VelSoundBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    character(30)        :: name               !変数名
    
    !-------------------------------------------------------------
    ! 基本場の取得
    !-------------------------------------------------------------
    name = "TempBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_TempBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_TempBZ )
    call BoundaryZSym_xz( xz_TempBZ )

    name = "PressBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_PressBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_PressBZ )
    call BoundaryZSym_xz( xz_PressBZ )

    name = "ExnerBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_ExnerBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_ExnerBZ )
    call BoundaryZSym_xz( xz_ExnerBZ )

    name = "PotTempBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_PotTempBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_PotTempBZ )
    call BoundaryZSym_xz( xz_PotTempBZ )

    name = "DensBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_DensBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_DensBZ )
    call BoundaryZSym_xz( xz_DensBZ )

    name = "VelSoundBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_VelSoundBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_VelSoundBZ )
    call BoundaryZSym_xz( xz_VelSoundBZ )

    name = "EffMolWtBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_EffMolWtBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_EffMolWtBZ )
    call BoundaryZSym_xz( xz_EffMolWtBZ )

    do s = 1, SpcNum
      name = trim(SpcWetSymbol(s))//'BasicZ'
      call HistoryGet( HistoryFile(7), name, xza_MixRtBZ(FileXMin:FileXMax, FileZMin:FileZMax, s) )
    end do

    call BoundaryXCyc_xza( xza_MixRtBZ )
    call BoundaryZSym_xza( xza_MixRtBZ )

    !----------------------------------------------------------
    ! BasicSet モジュールに値を設定
    !----------------------------------------------------------
    call BasicSetArray_Init(                       &
      & xz_PressBZ,   xz_ExnerBZ,   xz_TempBZ,     &
      & xz_PotTempBZ, xz_DensBZ,    xz_VelSoundBZ, &
      & xza_MixRtBZ,  xz_EffMolWtBZ )

  end subroutine AnalFile_BasicZ_Get


  subroutine ArareAlloc

    use gridset,       only: DimXMin, DimXMax, DimZMin, DimZMax, SpcNum

    !基本場, 擾乱場の取得.
    allocate( &
      &       Altitude(1:jmax), &
      &       Altitude2(1:jmax), &
      &       xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax ),             &
      &       xz_PotTempAll(DimXMin:DimXMax, DimZMin:DimZMax ),          &
      &        z_PotTempAll(DimZMin:DimZMax),                          &
      &       z_PotTempMoist(DimZMin:DimZMax,1:jmax ),          &
      &       z_VPotTempMoist(DimZMin:DimZMax,1:jmax ),          &
      &       z_VPotTempAll(DimZMin:DimZMax),          &
      &       z_DVPotTempAllDt(DimZMin:DimZMax),          &
      &       z_DPotTempAllDt(DimZMin:DimZMax),          &
      &       xz_Temp(DimXMin:DimXMax, DimZMin:DimZMax ),             &
      &       xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax ),          &
      &       z_TempAll(DimZMin:DimZMax ),          &
      &       z_TempBasicZ(DimZMin:DimZMax ),          &
      &       z_TempMoist(DimZMin:DimZMax ),          &
      &       z_VTempAll(DimZMin:DimZMax ),          &
      &       z_VTempMoist(DimZMin:DimZMax, 1:jmax ),          &
      &       z_VTempDown(DimZMin:DimZMax, 1:jmax ),          &
      &       z_VPotTempDown(DimZMin:DimZMax, 1:jmax ),          &
      &       z_DelTempMoist(DimZMin:DimZMax, 1:jmax ),          &
      &       z_DelPotTempMoist(DimZMin:DimZMax, 1:jmax ),          &
      &       z_DelVTempMoist(DimZMin:DimZMax, 1:jmax ),          &
      &       z_VTempMoist2(DimZMin:DimZMax, 1:jmax ),          &
      &       xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax ),               &
      &       z_ExnerMoist(DimZMin:DimZMax ),               &
      &       z_ExnerAll(DimZMin:DimZMax ),               &
      &       z_DensMoist( DimZMin:DimZMax ),               &
      &       z_DensBasicZ( DimZMin:DimZMax ),               &
      &       z_CAPE( DimZMin:DimZMax, jmax),               &
      &       CAPE(jmax),               &
      &       z_CloudWorkFunc( DimZMin:DimZMax, jmax),               &
      &       CloudWorkFunc(jmax),               &
      &       z_CloudWorkFuncTemp( DimZMin:DimZMax, jmax),               &
      &       CloudWorkFuncTemp(jmax),               &
      &       z_CloudWorkFuncMolWt( DimZMin:DimZMax, jmax),               &
      &       CloudWorkFuncMolWt(jmax),               &
      &       z_CloudWorkFuncMolWt2( DimZMin:DimZMax, jmax),               &
      &       CloudWorkFuncMolWt2(jmax),               &
      &       z_CloudWorkFunc2( DimZMin:DimZMax, jmax),               &
      &       CloudWorkFunc2(jmax),               &
      &       z_CloudWorkFunc3( DimZMin:DimZMax, jmax),               &
      &       CloudWorkFunc3(jmax),               &
      &       xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ),    &
      &       xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ), &
      &       za_MixRtAll(DimZMin:DimZMax, 1:SpcNum ), &
      &       za_MixRtMax(DimZMin:DimZMax, 1:SpcNum ), &
      &       za_MixRtMoist(DimZMin:DimZMax, 1:SpcNum, 1:jmax ), &
      &       za_MixRtMoist2(DimZMin:DimZMax, 1:SpcNum, 1:jmax ), &
      &       za_MolFrMoist(DimZMin:DimZMax, 1:SpcNum, 1:jmax ), &
      &       xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      &       z_PressAll(DimZMin:DimZMax ),           &
      &       z_PressMoist(DimZMin:DimZMax ),           &
      &       z_EffMolWt(DimZMin:DimZMax ),           &
      &       z_EffMolWtMoist(DimZMin:DimZMax ),           &
      &       z_EffMolWtMoist2(DimZMin:DimZMax ),           &
      &       za_MixRtDivMolWt(DimZMin:DimZMax,SpcNum),&
      &       z_LatentHeat(DimZMin:DimZMax),   &
      !
      &       z_MoistCond(DimZMin:DimZMax),   &
      &       z_MoistNH4SH(DimZMin:DimZMax),   &
      &       z_MoistENG(DimZMin:DimZMax,jmax),   &
      &       z_DryENG(DimZMin:DimZMax,jmax),   &
      !
      & za_MolFrAll(DimXMin:DimXMax, DimZMin:DimZMax),   &
      & z_MolWtAll(DimZMin:DimZMax),   &
      & z_MolWtMoist(DimZMin:DimZMax),   &
      & z_StabTempAll(DimZMin:DimZMax,jmax),   &
      & z_StabTempMoist(DimZMin:DimZMax,jmax),   &
      & z_StabMolWtAll(DimZMin:DimZMax,jmax),   &
      & z_StabMolWtMoist(DimZMin:DimZMax,jmax),   &
      & z_StabAll(DimZMin:DimZMax,jmax),   &
      & z_StabMoist(DimZMin:DimZMax,jmax)   &
      &   )

  end subroutine ArareAlloc

  function z_MeanX_xz (xz_Var)
    real(8) :: xz_Var(DimXMin:DimXMax,DimZMin:DimZMax)
    real(8) :: z_MeanX_xz(DimZMin:DimZMax)

    z_MeanX_xz = sum(xz_Var(RegXMin+1:RegXMax,:), 1) / (RegXMax - RegXMin)

  end function z_MeanX_xz

  function xz_MeanX_z (z_Var)
    real(8) :: z_Var(DimZMin:DimZMax)
    real(8) :: xz_MeanX_z(DimXMin:DimXMax,DimZMin:DimZMax)
    integer :: i

    do i = DimXMin, DimXMax
       xz_MeanX_z(i,DimZMin:DimZMax) = z_Var(DimZMin:DimZMax)
    end do
    
  end function xz_MeanX_z

  function za_MeanX_xza (xza_Var)
    real(8) :: xza_Var(DimXMin:DimXMax,DimZMin:DimZMax, 1:SpcNum)
    real(8) :: za_MeanX_xza(DimZMin:DimZMax, 1:SpcNum)
    integer :: s

    do s = 1, SpcNum
       za_MeanX_xza(:,s) = sum(xza_Var(RegXMin+1:RegXMax,:,s), 1) / (RegXMax - RegXMin)
    end do

  end function za_MeanX_xza


end program ArareCAPE
