!= Program ArareAnal
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: arare-analysis.f90,v 1.4 2011-02-28 11:39:32 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
! $BHs@ENO3X%b%G%k(B deepconv/arare 2 $B<!%G!<%?:n@.%W%m%0%i%`(B
!
!
program ArareAnal

  !----- $B%b%8%e!<%kFI$_9~$_(B ------

  !-----   $B7?@k8@(B, $BJ8;zNs=hM}(B   ----
  use dc_types,       only : STRING
  use dc_string,      only : StoA

  !-----   $B%a%C%;!<%8=PNO(B   -----
  use dc_message,     only: MessageNotify

  !  $B%3%^%s%I%i%$%s0z?t2r<a(B
  use argset,        only : argset_init  

  !-----    $B4IM}%b%8%e!<%k(B   -----
  !  $B2=3XNL7W;;%b%8%e!<%k(B
  use chemcalc, only: chemcalc_init, xz_SvapPress, xz_LatentHeat, ReactHeatNH4SH
  use chemdata, only: chemdata_init

  !  $BF~=PNO%U%!%$%kL>4IM}%b%8%e!<%k(B
  use fileset,       only : fileset_init

  !  $B%G%P%C%0=PNO4IM}%b%8%e!<%k(B
  use debugset,      only : debugset_init

  !  $B;~4V4IM}%b%8%e!<%k(B
  use timeset,       only : timeset_init, TimeDisp, TImeInt

  !  $B3J;RE@4IM}%b%8%e!<%k(B 
  use gridset,       only : gridset_init, DelX, DelZ, &
    &                       DimXMin, DimXMax, DimZMin, DimZMax, SpcNum, &
    &                       RegXMax, RegZMax, RegXMin, RegZMin

  !  $B4pK\>l@_Dj%b%8%e!<%k(B
  use basicset,      only : basicset_init, SpcWetSymbol, &
    &                       MolWtDry, MolWtWet, CpDry, PressBasis, GasRDry,    &
    &                       SpcWetID, SpcWetMolfr, xz_DensBasicZ,    &
    &                       xza_MixRtBasicZ, xz_PotTempBasicZ, xz_ExnerBasicZ, &
    &                       xz_TempBasicZ, xz_PressBasicZ, Grav, xz_EffMolWtBasicZ

  !  $B<>=a%k!<%A%s@_Dj%b%8%e!<%k(B
  use moistset,      only: moistset_init,   &
    &                      GasNum, CloudNum, RainNum, CondNum, RactNum, &
    &                      IdxG,   IdxC,   IdxR,                        &
    &                      IdxCG,  IdxCC,  IdxCR,                       &
    &                      IdxNH3, IdxH2S, IdxNH4SHc, IdxNH4SHr         

  use ChemCalc, only: xz_LatentHeat,    &!$B@xG.(B
    &                 ReactHeatNH4SH     !NH4SH $B$NH?1~G.(B

  !  $B@Q;;CM4IM}%b%8%e!<%k(B
  use StorePotTemp,      only : StorePotTemp_init
  use StoreMixRt,        only : StoreMixRt_init
  use StoreBuoy,         only : StoreBuoy_init
  use StoreStab,         only : StoreStab_init

  !-----    $B2<@A$1%b%8%e!<%k(B   -----
  !$BHyJ,J?6Q1i;;(B
  use average

  !  $BCGG.>e>:5$2t$N29EY8:N(7W;;MQ%b%8%e!<%k(B
  use ECCM,         only : eccm_molfr, eccm_stab, eccm_wet2
  
  ! $B6-3&>r7o(B
  use boundary,      only : BoundaryXCyc_xz, BoundaryZSym_xz

  !  $B<>=a5$2t$NIbNO7W;;MQ%b%8%e!<%k(B
  use MoistBuoyancy,only : MoistBuoy_Init, xr_BuoyMolWt, xr_BuoyDrag
  use DynFunc,       only : xr_Buoy
  
  use average, only:  xz_avr_xr
  use differentiate_center2, only: xr_dz_xz

  !$B0EL[$N7?@k8@6X;_(B
  implicit none
  
  !$BJQ?tDj5A(B
  real(8), allocatable :: xz_TempMoist(:,:)
  real(8), allocatable :: xz_PressMoist(:,:)
  real(8), allocatable :: xza_MolFrMoist(:,:,:)

  real(8), allocatable :: xz_PotTemp(:,:)
  real(8), allocatable :: xz_PotTempAll(:,:)
  real(8), allocatable :: xz_VTempAll(:,:)
  real(8), allocatable :: xz_Exner(:,:)
  real(8), allocatable :: pz_VelX(:,:)
  real(8), allocatable :: xr_VelZ(:,:)
  real(8), allocatable :: xza_MixRt(:,:,:)
  real(8), allocatable :: xza_MixRtAll(:,:,:)
  real(8), allocatable :: xza_MixRtAll2(:,:,:)
  real(8), allocatable :: xza_MixRtSat(:,:,:)
  real(8), allocatable :: xza_MixRtEq(:,:,:)
  real(8), allocatable :: xza_MixRtMoist(:,:,:)
  real(8), allocatable :: xza_Humidity(:,:,:)
  real(8), allocatable :: za_MolFrEq(:,:)
  real(8), allocatable :: xz_VPotTempAll(:,:)
  real(8), allocatable :: xz_EqPotTempAll(:,:)
  real(8), allocatable :: xz_Dens(:,:)
  real(8)              :: AnalTime
  real(8)              :: Hum

  real(8), allocatable :: xz_TempAll(:,:)
  real(8), allocatable :: xz_PressAll(:,:)
  real(8), allocatable :: xz_EffMolWt(:,:)
  real(8), allocatable :: xz_EffMolWtMoist(:,:)
  real(8), allocatable :: xza_MixRtDivMolWt(:,:,:)
  real(8), allocatable :: xza_LatentHeat(:,:,:)

  real(8), allocatable :: xz_EqConst(:,:)
  real(8), allocatable :: xz_EqConstSat(:,:)

  real(8), allocatable :: xza_MolFrAll(:,:,:)
  real(8), allocatable :: xz_MolWtWet(:,:)
  real(8), allocatable :: xz_MolWtMoist(:,:)
  real(8), allocatable :: xz_Stab(:,:)
  real(8), allocatable :: xz_StabTemp(:,:)
  real(8), allocatable :: xz_StabMolWt(:,:)
  real(8), allocatable :: xz_StabMoist(:,:)
  real(8), allocatable :: xz_StabTempMoist(:,:)
  real(8), allocatable :: xz_StabMolWtMoist(:,:)

  real(8), allocatable :: xz_StrmFunc(:,:)
  real(8), allocatable :: pq_StrmFunc(:,:)
  real(8), allocatable :: xr_Dens(:,:)
  real(8), allocatable :: pz_Dens(:,:)
  real(8), allocatable :: MolFrIni(:)

!  real(8), allocatable :: xza_CloudDens(:,:,:)
!  real(8), allocatable :: xz_CloudDens(:,:)
  real(8), allocatable :: xz_CloudMixRtEq(:,:)
  real(8), allocatable :: xza_CloudMixRtEq(:,:,:)
  real(8), allocatable :: xz_CloudMixRtMoist(:,:)
  real(8), allocatable :: xza_CloudMixRtMoist(:,:,:)

  real(8)              :: Humidity
  integer              :: t, s, i, k, tt, step
  integer              :: Idx
  character(100)       :: cfgfile


  !$B%3%^%s%I%i%$%s0z?t$N2r<a(B
  !  NAMELIST $B%U%!%$%kL>$NFI$_9~$_(B
  call argset_init(cfgfile)

  !$BJ*<AFC@-$N=i4|2=(B
  call chemdata_init()

  !$B;~9o$K4X$9$k@_Dj$N=i4|2=(B
  !  NAMELIST $B$+$iI,MW$J>pJs$rFI$_<h$j(B, $B;~4V4XO"$NJQ?t$N@_Dj$r9T$&(B. 
  call timeset_init(cfgfile)

  !$B3J;RE@>pJs$N=i4|2=(B
  !  NAMELIST $B$+$i>pJs$rF@$F(B, $B3J;RE@$r7W;;$9$k(B
  call gridset_init(cfgfile)

  !$B2=3X7W;;%k!<%A%s$N=i4|2=(B
  call chemcalc_init()

  !$B4pK\>l$N>pJs$N=i4|2=(B
  !  NAMELIST $B$+$i>pJs$rF@$F(B, $B4pK\>l$r@_Dj$9$k(B.
  call basicset_init(cfgfile)

  !I/O $B%U%!%$%kL>$N=i4|2=(B
  !  NAMELIST $B%U%!%$%kL>$r;XDj$7(B, deepconv/arare $B$N(B
  !  $B=PNO%U%!%$%kL>$r(B NAMELIST $B$+$iF@$k(B
  call fileset_init(cfgfile)

  !$B4pK\>l$N<hF@(B
  call AnalFile_BasicZ_Get( )

  !$B<>=a%k!<%A%s$N6&M-JQ?t$N=i4|2=(B
  call moistset_init()

  !$B@Q;;CM$rJ]4I$9$k$?$a$N%b%8%e!<%k$N=i4|2=(B
  !  NAMELIST $B$+$i>pJs$rF@$F(B, $B4pK\>l$r@_Dj$9$k(B.
  call StorePotTemp_init( )
  call StoreMixRt_init( )
  call StoreBuoy_init( )  
  call StoreStab_init( )  

  !$BFbItJQ?t$N=i4|2=(B. $B$H$j$"$($:%<%m$rF~$l$FCM$r3NDj$5$;$F$*$/(B. 
  call ArareAlloc

  !$BJ,;RNL$KBP$9$kIbNO7W;;%k!<%A%s$N=i4|2=(B
  call MoistBuoy_Init()             

  write(*,*) "OK1"
  
  !NetCDF $B%U%!%$%k$N:n@.$H4pK\>l$N<hF@(B  
  call AnalFile_Open( )
    

  write(*,*) "OK0"

  !----------------------------------------------------------------------
  ! $B4pK\>l$N29EY!&05NO$KBP$7$F<>EY(B 100 % $B$N;~$N:.9gHf$r5a$a$k(B
  ! $B4pK\>l$NJ*<ANL$O;~4VJQ2=$7$J$$$N$G(B, $B$3$N0LCV$G7W;;(B. 
  ! $BA4$F$N;~9o$KBP$7$FF1$8CM$,F~$k$3$H$K$J$k(B. 
  !----------------------------------------------------------------------
  Hum = 1.0d0
  do i = DimXMin, DimXMax
     call ECCM_MolFr( SpcWetMolFr(1:SpcNum), Hum, xz_TempBasicZ(i,:), &
          &              xz_PressBasicZ(i,:), za_MolFrEq )
  end do
  do i = DimXMin, DimXMax
     do s = 1, SpcNum
        xza_MixRtEq(i,:,s) = za_MolFrEq(:,s) * MolWtWet(s) / MolWtDry
     end do
  end do

  write(*,*) "*** s, MolWtWet(s) ***"
  do s = 1, GasNum
     write(*,*) s, MolWtWet(s)
  end do
  
  xz_CloudMixRtEq = 0.0d0
  do k = RegZMin, RegZMax
     do s = 1, GasNum
        xza_CloudMixRtEq(:,k,s) = xza_MixRtEq(:,k-1,s)  - xza_MixRtEq(:,k,s)
        xz_CloudMixRtEq(:,k) = xz_CloudMixRtEq(:,k) + xza_CloudMixRtEq(:,k,s) 
     end do
  end do
! WL1973 $BE*$J7W;;J}K!(B. 
!  do k = RegZMin, RegZMax
!     do s = 1, GasNum
!        xza_CloudDens(:,k,s) = (xza_MixRtEq(:,k-1,s)  - xza_MixRtEq(:,k,s)) * xz_PressBasicZ(:,k) / ( DelZ * Grav)
!        xz_CloudDens(:,k)    = xz_CloudDens(:,k) + xza_CloudDens(:,k,s)
!     end do
!       xz_CloudMixRt(:,k) = xz_CloudDens(:,k) / xz_DensBasicZ(:,k) ! $BHf<>$@$1$I(B
!  end do

  ! $B<>=aCGG.@~$N7W;;(B
  Hum = 1.0d0
  Idx = RegZMin
  ! $B=i4|CM$NMQ0U(B
  xz_TempMoist   = xz_TempBasicZ
  xz_PressMoist  = xz_PressBasicZ
  do i = DimXMin, DimXMax
     xza_MolFrMoist(i,:,:) = za_MolFrEq(:,:)
  end do
  do i = DimXMin, DimXMax
     call ECCM_Wet2(Idx, Hum, xz_TempMoist(i,:), xz_PressMoist(i,:), xza_MolFrMoist(i,:,:))
  end do
  do i = DimXMin, DimXMax
     do k = DimZMin, DimZMax
        do s = 1, SpcNum
           xza_MixRtMoist(i,k,s) = xza_MolFrMoist(i,k,s) * MolWtWet(s) / MolWtDry
        end do
     end do
  end do
  do k = DimZMin, DimZMax
     do i = DimXMin, DimXMax
        xz_MolWtMoist(i,k) = dot_product( MolWtWet(1:GasNum), xza_MolFrMoist(i,k,1:GasNum) )
     end do
  end do
  
  xz_StabTempMoist =                                           &
       &         Grav / xz_TempMoist                           &
       &           * (   xz_avr_xr( xr_dz_xz( xz_TempMoist) )  &
       &               + Grav / CpDry ) 
  xz_StabMolWtMoist =                                          &
       &       - Grav * xz_avr_xr( xr_dz_xz( xz_MolWtMoist ) ) &
       &         /  MolWtDry
  xz_StabMoist = xz_StabTempMoist + xz_StabMolWtMoist
  
  xz_CloudMixRtMoist = 0.0d0
  do k = RegZMin, RegZMax
     do s = 1, GasNum
        xza_CloudMixRtMoist(:,k,s) = xza_MixRtMoist(:,k-1,s)  - xza_MixRtMoist(:,k,s)
        xz_CloudMixRtMoist(:,k) = xz_CloudMixRtMoist(:,k) + xza_CloudMixRtMoist(:,k,s) 
     end do
  end do  

  !----------------------------------------------------------------------
  ! $B2r@OCM$N7W;;3+;O(B
  !   $B;~9o$r?J$a$J$,$i(B, $B%U%!%$%k$NCM$rF@$k(B
  !   $BE,59%G!<%?$r4V0z$/(B ($B:G=i$N%G!<%?$O=PNO$7$J$$(B, $B:GBg$G(B 500 $B%9%F%C%WJ,(B). 
  !----------------------------------------------------------------------
!  do t = 1, int( TimeInt / TimeDisp ) + 1
  step = int( TimeInt / TimeDisp )
!!  STEPS: do while (step >= 500)
!!   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
    
    !----------------------------------------------------------------
    ! $B%R%9%H%j%U%!%$%k$r3+$-CM$rF@$k(B. 
    !----------------------------------------------------------------
    call AnalFile_Get(                                           &
      &       t,  AnalTime,                                      &
      &       xz_PotTemp, xz_Exner, pz_VelX, xr_VelZ, xza_MixRt )
    write(*,*) t, AnalTime    
  
    !----------------------------------------------------------------
    ! $B=i4|2=(B
    !----------------------------------------------------------------
    xza_MixRtAll = xza_MixRt + xza_MixRtBasicZ

    !----------------------------------------------------------------
    ! $B29EY!&05NO$N7W;;(B
    !----------------------------------------------------------------
    xz_PressAll = &
      & PressBasis * ((xz_Exner + xz_ExnerBasicZ ) ** (CpDry / GasRDry))
    
    xz_PotTempAll  = xz_PotTemp + xz_PotTempBasicZ
    xz_TempAll = ( xz_Exner + xz_ExnerBasicZ ) * xz_PotTempAll

    !----------------------------------------------------------------
    ! $BL)EY(B, $B2>290L$N7W;;(B
    !----------------------------------------------------------------
    xza_MixRtDivMolWt = 0.0d0
    do s = 1, GasNum
      xza_MixRtDivMolWt(:,:,IdxG(s)) = xza_MixRtAll(:,:,IdxG(s)) / MolWtWet(IdxG(s))
    end do

    xz_EffMolWt = &
      & (1.0d0 + sum(xza_MixRtAll,3)) &
      & / (MolWtDry * ((1.0d0 / MolWtDry) + sum(xza_MixRtDivMolWt,3))) 

    xz_Dens = xz_PressAll * xz_EffMolWt / (GasRDry * xz_TempAll)
    xz_VPotTempAll = xz_PotTempAll / xz_EffMolWt
    xz_VTempAll = xz_TempAll / xz_EffMolWt

    !----------------------------------------------------------------
    ! $BN.@~4X?t(B
    !----------------------------------------------------------------
    xr_Dens = xr_avr_xz( xz_Dens )
    pz_Dens = pz_avr_xz( xz_Dens )
    pq_StrmFunc = 0.0d0
    xz_StrmFunc = 0.0d0

    do i = RegXMax, DimXMin+1, -1    
      pq_StrmFunc(i-1,:) = pq_StrmFunc(i,:) - DelX * xr_Dens(i,:) * xr_VelZ(i,:)
    end do
    do k = RegZMax, DimZMin+1, -1
      pq_StrmFunc(:,k-1) = pq_StrmFunc(:,k) + DelZ * pz_Dens(:,k) * pz_VelX(:,k)
    end do

    xz_StrmFunc = pq_StrmFunc 
    xz_StrmFunc(DimXMin+1 : DimXMax, DimZMin+1 : DimZMax) =           &
      &  (                                                            &
      &     pq_StrmFunc(DimXMin   : DimXMax-1, DimZMin   : DimZMax-1) &
      &   + pq_StrmFunc(DimXMin+1 : DimXMax,   DimZMin   : DimZMax-1) &
      &   + pq_StrmFunc(DimXMin   : DimXMax-1, DimZMin+1 : DimZMax)   &
      &   + pq_StrmFunc(DimXMin+1 : DimXMax,   DimZMin+1 : DimZMax)   &
      &   ) * 2.5d-1

    call BoundaryXCyc_xz( xz_StrmFunc )
    call BoundaryZSym_xz( xz_StrmFunc )

    where (xz_StrmFunc < 1.0d-20 .AND. xz_StrmFunc > - 1.0d-20) 
      xz_StrmFunc = 1.0d-20
    end where

    !----------------------------------------------------------------
    ! $BK0OB>x5$05$HJ?9UDj?t(B
    !----------------------------------------------------------------
    do s = 1, CondNum
      !$BK0OB>x5$05(B
      xza_MixRtSat(:,:,IdxCG(s)) =                        &
        & xz_SvapPress(SpcWetID(IdxCC(s)), xz_TempAll)    &
        &  * MolWtWet(IdxCC(s))                           &
        &  / (MolWtDry * xz_PressAll)                     

      !$BAjBP<>EY(B
      xza_Humidity(:,:,IdxCG(s)) =                        &
        & xza_MixRtAll(:,:,IdxCG(s)) / xza_MixRtSat(:,:,IdxCG(s)) * 100.0d0
    end do
    
    !$BJ?9UDj?t(B
    if ( IdxNH3 /= 0 ) then     
      xz_EqConstSat = 61.781d0 - 10834.0d0 / xz_TempAll - dlog(1.0d2)
      
      xz_EqConst =                                        &
        & dlog(                                           &
        &   ( xz_PressAll ** 2.0d0 )                      &
        &   * max( xza_MixRtAll(:,:,IdxNH3), 1.0d-20 )    &
        &   * max( xza_MixRtAll(:,:,IdxH2S), 1.0d-20 )    &
        &   * ( MolWtDry ** 2.0d0)                        &
        &   / ( MolWtWet(IdxNH3) )                        &
        &   / ( MolWtWet(IdxH2S) )                        &
        & )
    end if


    !----------------------------------------------------------------
    ! $BAjEv290L(B
    !----------------------------------------------------------------
    ! $B6u5$2t$rCGG.>e>:$5$;$?;~$K<B8=$9$k29EY!&J*<AJ,I[$r5a$a$k(B. 
!    MolFrIni = SpcWetMolFr(1:SpcNum) 
!    Humidity = 1.0
!    call ECCM_Wet( MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr )
    
!    xz_EqPotTempAll = xz_PotTempAll * dexp( xz_LatentHeat * xz_MixRt / (CpDry * xz_TempAll))    
!    xa_SatEqPotTempAll = xz_PotTempAll * dexp( xz_LatentHeat * xz_MixRt / (CpDry * xz_TempAll))

    !----------------------------------------------------------------
    ! $B0BDjEY(B
    !----------------------------------------------------------------
    do s = 1, SpcNum
      xza_MolFrAll(:,:,s) = xza_MixRtAll(:,:,s) * MolWtDry / MolWtWet(s) 
    end do
    
    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        xz_MolWtWet(i,k) = dot_product( MolWtWet(1:GasNum), xza_MolFrAll(i,k,1:GasNum) )
      end do
    end do
    
    xz_StabTemp =                                           &
      &         Grav / xz_TempAll                           &
      &           * (   xz_avr_xr( xr_dz_xz( xz_TempAll ) ) &
      &               + Grav / CpDry ) 
    xz_StabMolWt =                                          &
      &       - Grav * xz_avr_xr( xr_dz_xz( xz_MolWtWet ) ) &
      &         /  MolWtDry    
    xz_Stab = xz_StabTemp + xz_StabMolWt

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

    !----------------------------------------------------------------    
    ! $B%U%!%$%k=PNO(B
    !----------------------------------------------------------------
    call AnalFile_OutPut( AnalTime )    

 end do TIMESTEP
 
    
  !----------------------------------------------------------------    
  ! $B%U%!%$%k$rJD$8$k(B
  !----------------------------------------------------------------    
  call AnalFile_Close
 
contains
  
  subroutine AnalFile_Open( )
    !
    !$B2r@O%U%!%$%k$NDj5A(B
    !
    use gtool_history, only: HistoryCreate, HistoryPut, HistoryAddVariable
    use fileset, only: exptitle, expsrc, expinst, HistoryFilePrefix    
    use gridset, only: FileNX, FileNZ, s_X, s_Z, &
      &                FileXMin, FileXMax, FileZMin, FileZMax, SpcNum
    use basicset, only: SpcWetSymbol
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none

    !$BJQ?tDj5A(B
    character(150)  :: File
    integer        :: s
    
    File = trim(HistoryFilePrefix)// '_Analysis.nc'

    write(*,*) "A" 

    !-----------------------------------------------------------
    ! $B%R%9%H%j!<:n@.(B
    !-----------------------------------------------------------
    call HistoryCreate(                              &
      & file = trim(File),                           &
      & title = exptitle,                            &
      & source = expsrc,                             &
      & institution = expinst,                       &
      & dims=(/'x','z','t'/),                        &
      & dimsizes=(/FileNX, FileNZ, 0/),              &
      & longnames=(/'X-coordinate',                  &
      &             'Z-coordinate',                  &
      &             'Time        '/),                &
      & units=(/'m','m','s'/), origin=0.0,           &
      & interval=0.0)

    write(*,*) "B" 

    !-----------------------------------------------------------  
    ! $B<4$N=PNO(B
    !-----------------------------------------------------------
    call HistoryPut('x', s_X( FileXMin: FileXMax ))

    write(*,*) "C" 

    call HistoryPut('z', s_Z( FileZMin: FileZMax ))

    write(*,*) "D" 

    !-----------------------------------------------------------  
    ! $B2r@OMQ$NJQ?t$N=PNO(B
    !-----------------------------------------------------------  
    !$B29EY(B
    call HistoryAddVariable(                          &
      & varname='TempAll', dims=(/'x','z','t'/),      &
      & longname='temperature',                       &
      & units='K', xtype='double' )

    !$B2>290L(B
    call HistoryAddVariable(                         &
      & varname='VTempAll', dims=(/'x','z','t'/),    &
      & longname='virtual temperature',              &
      & units='K', xtype='double')

    !$B290L(B
    call HistoryAddVariable(                         &
      & varname='PotTempAll', dims=(/'x','z','t'/),  &
      & longname='potential temperature',            &
      & units='K', xtype='double' )

    !$B2>290L(B
    call HistoryAddVariable(                         &
      & varname='VPotTempAll', dims=(/'x','z','t'/), &
      & longname='virtual potential temperature',    &
      & units='K', xtype='double')

    !$B05NO(B
    call HistoryAddVariable(                          &
      & varname='PressAll', dims=(/'x','z','t'/),     &
      & longname='Pressure',                          &
      & units='Pa', xtype='double' )
    
    !$BL)EY(B
    call HistoryAddVariable(                          &
      & varname='Dens', dims=(/'x','z','t'/),         &
      & longname='Density',                           &
      & units='kg/m/m"', xtype='double' )

    !$B@EE*0BDjEY(B
    call HistoryAddVariable(                          &
      & varname='Stab', dims=(/'x','z','t'/),         &
      & longname='Stability',                         &
      & units='1/s/s"', xtype='double' )

    write(*,*) "J" 

    !$B@EE*0BDjEY(B
    call HistoryAddVariable(                          &
      & varname='StabTemp', dims=(/'x','z','t'/),     &
      & longname='Stability (Temp)',                  &
      & units='1/s/s"', xtype='double' )

    !$B@EE*0BDjEY(B
    call HistoryAddVariable(                          &
      & varname='StabMolWt', dims=(/'x','z','t'/),    &
      & longname='Stability (MolWt)',                 &
      & units='1/s/s"', xtype='double' )

    !$B@EE*0BDjEY(B
    call HistoryAddVariable(                          &
      & varname='StabMoist', dims=(/'x','z','t'/),         &
      & longname='Stability',                         &
      & units='1/s/s"', xtype='double' )

    !$B@EE*0BDjEY(B
    call HistoryAddVariable(                          &
      & varname='StabTempMoist', dims=(/'x','z','t'/),     &
      & longname='Stability (Temp)',                  &
      & units='1/s/s"', xtype='double' )

    write(*,*) "I" 

    !$B@EE*0BDjEY(B
    call HistoryAddVariable(                          &
      & varname='StabMolWtMoist', dims=(/'x','z','t'/),    &
      & longname='Stability (MolWt)',                 &
      & units='1/s/s"', xtype='double' )

    ! NH4SH $B$N05J?9UDj?t(B ($BK0OB(B)
    call HistoryAddVariable(                          &
      & varname='EqConstSat',                         &
      & dims=(/'x','z','t'/),                         &
      & longname='Saturated Equilibrium Constant',    &
      & units='1', xtype='double' )

    ! NH4SH $B$N05J?9UDj?t(B 
    call HistoryAddVariable(                          &
      & varname='EqConst',                            &
      & dims=(/'x','z','t'/),                         &
      & longname='Equibrium Constant',      &
      & units='1', xtype='double' )

    write(*,*) "H" 

    !$BN.@~4X?t(B 
    call HistoryAddVariable(                              &
      & varname='StrmFunc', dims=(/'x','z','t'/),         &
      & longname='Stream Function',                       &
      & units='kg/m/m', xtype='double')

    !$B1@L)EY(B
    call HistoryAddVariable(                              &
      & varname='CloudMixRtMoist', dims=(/'x','z','t'/),         &
      & longname='Cloud Density',                       &
      & units='kg/m/m', xtype='double')

    call HistoryAddVariable(                              &
      & varname='CloudMixRtEq', dims=(/'x','z','t'/),         &
      & longname='Cloud Mixing ratio',                       &
      & units='kg/m/m', xtype='double')

    do s = 1, SpcNum
      !$BG.J?9U>uBV$G$N:.9gHf(B
      call HistoryAddVariable(                               &
        & varname='D'//trim(SpcWetSymbol(s))//'DtAll',             &
        & dims=(/'x','z', 't'/),                             &
        & longname= trim(SpcWetSymbol(s))//' Mixing Ratio',  &
        & units='kg/kg', xtype='double')
    end do
    
    do s = 1, GasNum
      !$BG.J?9U>uBV$G$N:.9gHf(B
      call HistoryAddVariable(                                    &
        & varname='D'//trim(SpcWetSymbol(IdxG(s)))//'DtEq',             &
        & dims=(/'x','z', 't'/),                                       &
        & longname='Equilibrium '//trim(SpcWetSymbol(IdxG(s)))//' Mixing Ratio',&
        & units='kg/kg', xtype='double' )
    end do

    write(*,*) "G" 

    do s = 1, CondNum
      !$BK0OB:.9gHf(B
      call HistoryAddVariable(                          &
        & varname='D'//trim(SpcWetSymbol(IdxCG(s))) // '_Sat',      &
        & dims=(/'x','z','t'/),                         &
        & longname='Saturated vapor pressure',          &
        & units='kg/kg', xtype='double')
      
      !$BAjBP<>EY(B
      call HistoryAddVariable(                          &
        & varname='D'//trim(SpcWetSymbol(IdxCG(s))) // '_Hum',      &
        & dims=(/'x','z','t'/),                         &
        & longname='Humidity',                          &
        & units='%', xtype='double')
   end do

    write(*,*) "F" 

   do s = 1, GasNum
      !$B1@NL(B
      call HistoryAddVariable(                          &
        & varname='D'//trim(SpcWetSymbol(s)) // '_MixRtMoist',      &
        & dims=(/'x','z','t'/),                         &
        & longname='Cloud mixing ratio',                &
        & units='kg/kg', xtype='double')
    end do

    write(*,*) "E" 
    
  end subroutine AnalFile_Open


  subroutine AnalFile_Close( )
    !
    !$B2r@O%U%!%$%k$N%/%m!<%:(B
    !
    use gtool_history, only: HistoryClose

    call HistoryClose
    
  end subroutine AnalFile_Close
  
  
  subroutine AnalFile_OutPut( Time )

    use gtool_history, only: HistoryPut
    use gridset, only: FileXMin, FileXMax, FileZMin, FileZMax

    !$B0EL[$N7?@k8@6X;_(B
    implicit none
    
    !$BJQ?tDj5A(B
    real(8), intent(in) :: Time    
 
    !----------------------------------------------------------------
    ! $BCM$r=PNO(B
    !----------------------------------------------------------------
    call HistoryPut( 't', Time )

    call HistoryPut( 'TempAll',       &
      &            xz_TempAll(FileXMin:FileXMax, FileZMin:FileZMax) )

    call HistoryPut( 'VTempAll',       &
      &            xz_VTempAll(FileXMin:FileXMax, FileZMin:FileZMax) )
    
    call HistoryPut( 'PotTempAll',    &
      &            xz_PotTempAll(FileXMin:FileXMax, FileZMin:FileZMax))
   
    call HistoryPut( 'VPotTempAll',   &
      &            xz_VPotTempAll(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'PressAll',      &
      &            xz_PressAll(FileXMin:FileXMax, FileZMin:FileZMax))
    
    call HistoryPut( 'Stab',          &
      &            xz_Stab(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'StabTemp',          &
      &            xz_StabTemp(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'StabMolWt',          &
      &            xz_StabMolWt(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'StabMoist',          &
      &            xz_StabMoist(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'StabTempMoist',          &
      &            xz_StabTempMoist(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'StabMolWtMoist',          &
      &            xz_StabMolWtMoist(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'Dens',          &
      &            xz_Dens(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut(  'EqConstSat',   &
      & xz_EqConstSat(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'EqConst',       &
      & xz_EqConst(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'StrmFunc',      &
      &            xz_StrmFunc(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'CloudMixRtEq',      &
      &            xz_CloudMixRtEq(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'CloudMixRtMoist',      &
      &            xz_CloudMixRtMoist(FileXMin:FileXMax, FileZMin:FileZMax))

    do s = 1, SpcNum
      call HistoryPut( trim(SpcWetSymbol(s))//'DtAll',  &
        &           xza_MixRtAll(FileXMin:FileXMax, FileZMin:FileZMax, s))
    end do

    do s = 1, GasNum
      call HistoryPut( trim(SpcWetSymbol(IdxG(s)))//'DtEq',  &
        &           xza_MixRtEq(FileXMin:FileXMax, FileZMin:FileZMax, IdxG(s)))
    end do

    do s = 1, CondNum
      call HistoryPut(     &
        & trim(SpcWetSymbol(IdxCG(s))) // '_Sat', &
        & xza_MixRtSat(FileXMin:FileXMax, FileZMin:FileZMax, IdxCG(s)))

      call HistoryPut(     &
        & trim(SpcWetSymbol(IdxCG(s))) // '_Hum', &
        & xza_Humidity(FileXMin:FileXMax, FileZMin:FileZMax, IdxCG(s)))
   end do

   do s = 1, GasNum
      call HistoryPut(     &
        & trim(SpcWetSymbol(s)) // '_MixRtMoist', &
        & xza_CloudMixRtMoist(FileXMin:FileXMax, FileZMin:FileZMax, s))
   end do

  end subroutine AnalFile_OutPut
  


!!!------------------------------------------------------------------------------!!!
  subroutine AnalFile_Get(                                      &
    &       i, AnalTime,                                        &
    &       xz_PotTemp, xz_Exner, pz_VelX, xr_VelZ, xza_MixRt )
 
    use dc_string
    use gtool_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

    !$B0EL[$N7?@k8@6X;_(B
    implicit none

    !$BJQ?tDj5A(B
    integer, intent(in)  :: i    
    real(8), intent(out) :: AnalTime
    real(8), intent(out) :: pz_VelX(DimXMin:DimXMax,DimZMin:DimZMax)
    real(8), intent(out) :: xr_VelZ(DimXMin:DimXMax,DimZMin:DimZMax)
    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               !$BJQ?tL>(B
    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 )

    name = "VelX"
    call HistoryGet( HistoryFile(3), name, pz_VelX(FileXMin:FileXMax, FileZMin:FileZMax), range=step )    
    call BoundaryXCyc_pz( pz_VelX )
    call BoundaryZSym_pz( pz_VelX )
    
    name = "VelZ"
    call HistoryGet( HistoryFile(4), name, xr_VelZ(FileXMin:FileXMax, FileZMin:FileZMax), range=step )
    call BoundaryXCyc_xr( xr_VelZ )
    call BoundaryZAntiSym_xr( xr_VelZ )

!    name = "Km"
!    call HistoryGet( HistoryFile(5), name, xz_Km(FileXMin:FileXMax, FileZMin:FileZMax), range=step )
!    xz_Km = xz_BoundaryXCyc_xz( xz_Km )
!    xz_Km = xz_BoundaryZSym_xz( xz_Km )

    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 gtool_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

    !$B0EL[$N7?@k8@6X;_(B
    implicit none

    !$BJQ?tDj5A(B
    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               !$BJQ?tL>(B
    

    !-------------------------------------------------------------
    ! $B4pK\>l$N<hF@(B
    !-------------------------------------------------------------
    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 $B%b%8%e!<%k$KCM$r@_Dj(B
    !----------------------------------------------------------
    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

    !$B4pK\>l(B, $B>qMp>l$N<hF@(B.
    allocate( &
      &       xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax ),             &
      &       xz_PotTempAll(DimXMin:DimXMax, DimZMin:DimZMax ),          &
      &       xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax ),               &
      &       pz_VelX(DimXMin:DimXMax, DimZMin:DimZMax ),                &
      &       xr_VelZ(DimXMin:DimXMax, DimZMin:DimZMax ),                &
      &       xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ),    &
      &       xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ), &
      &       xza_MixRtAll2(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ), &
      &       xza_MixRtSat(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ), &
      &       xza_MixRtEq(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ),  &
      &       xza_MixRtMoist(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ),  &
      &       xza_Humidity(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ), &
      &       za_MolFrEq(DimZMin:DimZMax, 1:SpcNum ),                    &
      &       xz_VTempAll(DimXMin:DimXMax, DimZMin:DimZMax ),            &
      &       xz_VPotTempAll(DimXMin:DimXMax, DimZMin:DimZMax ),         &
      &       xz_EqPotTempAll(DimXMin:DimXMax, DimZMin:DimZMax ),        &
      &       xz_Dens(DimXMin:DimXMax, DimZMin:DimZMax),                 &
      &       xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax ),            &
      &       xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      &       xz_EqConstSat(DimXMin:DimXMax, DimZMin:DimZMax ),         &
      &       xz_EqConst(DimXMin:DimXMax, DimZMin:DimZMax ),            &
      !
      &       xz_EffMolWt(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      &       xz_EffMolWtMoist(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      &       xza_MixRtDivMolWt(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum),&
      &       xza_LatentHeat(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum),   &
      !
      &       xza_MolFrAll(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum), &
      &       xz_MolWtWet(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_MolWtMoist(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_Stab(DimXMin:DimXMax,DimZMin:DimZMax),       &
      &       xz_StabTemp(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_StabMolWt(DimXMin:DimXMax,DimZMin:DimZMax),  &
      &       xz_StabMoist(DimXMin:DimXMax,DimZMin:DimZMax),       &
      &       xz_StabTempMoist(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_StabMolWtMoist(DimXMin:DimXMax,DimZMin:DimZMax),  &
!      &       xza_CloudDens(DimXMin:DimXMax,DimZMin:DimZMax,1:SpcNum),  &
!      &       xz_CloudDens(DimXMin:DimXMax,DimZMin:DimZMax),  &
      &       xz_CloudMixRtMoist(DimXMin:DimXMax,DimZMin:DimZMax),  &
      &       xza_CloudMixRtMoist(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum),  &
      &       xz_CloudMixRtEq(DimXMin:DimXMax,DimZMin:DimZMax),  &
      &       xza_CloudMixRtEq(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum),  &
      !
      &       xz_TempMoist(DimXMin:DimXMax,DimZMin:DimZMax),  &
      &       xz_PressMoist(DimXMin:DimXMax,DimZMin:DimZMax),  &
      &       xza_MolFrMoist(DimXMin:DimXMax,DimZMin:DimZMax,1:SpcNum),  &
      !
      &       MolFrIni(SpcNum), &
      &       xz_StrmFunc(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      &       pq_StrmFunc(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      &       xr_Dens(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      &       pz_Dens(DimXMin:DimXMax, DimZMin:DimZMax )       )

    xz_StrmFunc = 0.0d0
    xza_MixRtSat = 0.0d0
    xza_MixRtDivMolWt = 0.0d0

  end subroutine ArareAlloc

end program ArareAnal
