!= Program ArareAnal
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: arareanal.f90,v 1.26 2009-02-28 08:19:11 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@O(B   -----
  use dc_args,        only : ARGS, Open, Debug, Help, Strict, Close, Option
  

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

  !  $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

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

  !  $B@Q;;CM4IM}%b%8%e!<%k(B
  use storeset,      only : storeset_init
  use storeset2,     only : store2set_init


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

  !-----       $BNO3X2aDx(B        -----
  !  $BNO3X2aDx7W;;MQ4X?t%b%8%e!<%k(B
  use DynFunc,       only :  xr_Buoy

  !-----       $BJ*M}2aDx(B        -----
!  !  $B<>=aK0OBD4@aK!7W;;MQ%b%8%e!<%k(B
!  use MoistAdjust,  only : MoistAdjust_Init, MoistAdjustSvapPress, MoistAdjustNH4SH

  !  $B1@J*M}%Q%i%a%?%j%<!<%7%g%s(B
  use WarmRainPrm,  only : WarmRainPrm_Init, WarmRainPrm_prm

  !  $B<>=a5$2t$NIbNO7W;;MQ%b%8%e!<%k(B
  use MoistBuoyancy,only : MoistBuoy_Init, xr_BuoyMolWt, xr_BuoyDrag
  
  !  $BCGG.>e>:5$2t$N29EY8:N(7W;;MQ%b%8%e!<%k(B
  use ECCM,         only : eccm_init, eccm_stab, eccm_molfr
  

  !$B0EL[$N7?@k8@6X;_(B
  implicit none
  
  !$BJQ?tDj5A(B
  real(8), allocatable :: xz_PotTemp(:,:)
  real(8), allocatable :: xz_PotTempAll(:,:)
  real(8), allocatable :: xz_Exner(:,:)
  real(8), allocatable :: pz_VelX(:,:)
  real(8), allocatable :: xr_VelZ(:,:)
!  real(8), allocatable :: xz_Km(:,:)
  real(8), allocatable :: xza_MixRt(:,:,:)
  real(8), allocatable :: xza_MixRtAll(:,:,:)
  real(8), allocatable :: xza_MixRtSat(:,:,:)
  real(8), allocatable :: xza_MixRtEq(:,:,:)
  real(8), allocatable :: za_MolFrEq(:,:)
  real(8), allocatable :: xz_Temp(:,:)
  real(8), allocatable :: xz_VPotTemp(:,:)
  real(8), allocatable :: xz_VPotTempAll(:,:)
  real(8), allocatable :: xz_Dens(:,:)
  real(8), allocatable :: xz_BuoyAll(:,:)
  real(8), allocatable :: xz_BuoyTemp(:,:)
  real(8), allocatable :: xz_BuoyMolWt(:,:)
  real(8), allocatable :: xz_BuoyDrag(:,:)
  real(8), allocatable :: xz_Stab(:,:)
  real(8), allocatable :: xz_StabTemp(:,:)
  real(8), allocatable :: xz_StabMolWt(:,:)
  real(8)              :: AnalTime
  real(8)              :: Hum

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

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

  real(8), allocatable :: xz_StrmFunc(:,:)
  real(8), allocatable :: xz_Cloud(:,:)
  real(8), allocatable :: xz_Rain(:,:)

  integer      :: LoopNum, LoopNum2
  integer      :: GasNum(10)
  integer      :: CloudNum(10)
  integer      :: RainNum(10)
  integer      :: NH3Num  
  integer      :: H2SNum  
  integer      :: NH4SHNum

  integer             :: t, s, i, k
  character(50)       :: cfgfile

  !-----   $B%3%^%s%I%i%$%s9T$-?t2r@OMQJQ?t(B    -----
  type(ARGS)           :: arg             ! $B%3%^%s%I%i%$%s0z?t>pJs(B
  logical              :: OPT_namelist    ! $B%3%^%s%I%i%$%s0z?tMQO@M}JQ?t(B
  character(STRING)    :: VAL_namelist    ! $B%3%^%s%I%i%$%s0z?t$NCM(B
  

  !----------------------------------------------------------------------
  ! $BJQ?t;2>H7?%b%8%e!<%k$N=i4|2=(B
  !----------------------------------------------------------------------
  !NAMELIST $B%U%!%$%k$N<hF@(B
  !  gtool5 $B%i%$%V%i%j$N(B dc_args $B%b%8%e!<%k$rMxMQ(B.
  !  $B;XDj2DG=$J%*%W%7%g%s$O0J2<$NDL$j(B.
  !
  !  -N=VAL, --namelist=VAL
  !    specify Namelist file (default is 'arare.conf'). 
  !
  !  -D=VAL, --debug=VAL
  !    call dc_trace#SetDebug (display a lot of messages for debug).
  !    VAL is unit number (default is standard output)
  !
  !  -h=VAL, -H=VAL, --help=VAL
  !    display this help and exit. VAL is unit number (default is
  !    standard output)
  !
  call Open(arg)
  call Option(arg, StoA('-N', '--namelist'), &
    &         OPT_namelist, VAL_namelist, &
    &         help="specify Namelist file (default is 'arare.conf')." )
                    ! "-N/--namelist" $B%*%W%7%g%s$N@_Dj(B
  call Debug(arg)   ! $B%G%P%C%0%*%W%7%g%s$N<+F0@_Dj(B
  call Help(arg)    ! $B%X%k%W%*%W%7%g%s$N<+F0@_Dj(B
  call Strict(arg)  ! $BL58z$J%*%W%7%g%s;XDj;~$K7Y9p$rI=<((B

  !"-N/-namelist" $B%*%W%7%g%s$N2r<a(B
  !  $BM?$($i$l$F$$$J$$>l9g$O%G%U%)%k%HCM(B (arare.conf) $B$r(B 
  !  NAMLIST $B%U%!%$%kL>$H$9$k(B.
  !
  if (OPT_namelist) then
     call MessageNotify( "M", "main", &
       &                 "Namelist file is '%c'", c1=trim(VAL_namelist) )
     cfgfile=trim(VAL_namelist)
  else
     call MessageNotify( "W", "main", &
       &                 "Namelist file is not specified." )
     call MessageNotify( "M", "main", &
       &                 "Use default Namelist file (arare.conf)." )
     cfgfile="arare.conf"
  end if

  call Close(arg)

  !$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)
  
  !$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)
  
  !  NAMELIST $B$+$i>pJs$rF@$F(B, $B4pK\>l$r@_Dj$9$k(B.
  call storeset_init( )
  call store2set_init( )
  
  !$BFbItJQ?t$N=i4|2=(B. $B$H$j$"$($:%<%m$rF~$l$FCM$r3NDj$5$;$F$*$/(B. 
  call ArareAlloc
  
  !NetCDF $B%U%!%$%k$N:n@.$H4pK\>l$N<hF@(B  
  call AnalFile_Open( )
  call AnalFile_OpenDens(  )

  !$B4pK\>l$N<hF@(B
  call AnalFile_BasicZ_Get( )
    
  !----------------------------------------------------------------------
  ! $B%Q%C%1!<%87?%b%8%e!<%k$N=i4|2=(B
  !   $B%G%U%)%k%H$NCM$+$iJQ99$9$kI,MW$N$"$k%k!<%A%s$N$_=i4|2=(B
  !----------------------------------------------------------------------
  call ECCM_Init
!  call Damping_Init( cfgfile )      !$BGH$N8:?j78?t$N=i4|2=(B
!  call NumDiffusion_Init()          !$B?tCM3H;69`$N=i4|2=(B
!  call Turbulence_Init()            !$BMpN.7W;;$N=i4|2=(B
!  call MoistAdjust_Init()           !$B<>=aK0OBD4@0K!$N=i4|2=(B
  call WarmRainPrm_Init( cfgfile )  !$BCH$+$$1+$N%Q%i%a%?%j%<!<%7%g%s$N=i4|2=(B
!  call Radiation_Init( cfgfile )    !$BJ|<M6/@)$N=i4|2=(B
  call MoistBuoy_Init()             !$BJ,;RNL$KBP$9$kIbNO7W;;%k!<%A%s$N=i4|2=(B

  call WarmRainPrm_prm( LoopNum, LoopNum2, RainNum, CloudNum, GasNum, NH3Num, H2SNum, NH4SHNum )

  !----------------------------------------------------------------------
  ! $B2r@OCM$N7W;;3+;O(B
  !   $B;~9o$r?J$a$J$,$i(B, $B%U%!%$%k$NCM$rF@$k(B
  !----------------------------------------------------------------------

  do t = 1, int( TimeInt / TimeDisp ) + 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
    !----------------------------------------------------------------
    xz_Cloud = 0.0d0
    xz_Rain  = 0.0d0
    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
    xz_Temp = xz_TempAll - xz_ExnerBasicZ * xz_PotTempBasicZ 
    
    !----------------------------------------------------------------
    ! $BL)EY(B, $B2>290L$N7W;;(B
    !----------------------------------------------------------------
    xza_MixRtDivMolWt = 0.0d0
    do s = 1, 3
      xza_MixRtDivMolWt(:,:,s) = xza_MixRtAll(:,:,s) / MolWtWet(s)
    end do

    xz_EffMolWt = &
      & (1.0d0 / MolWtDry) / ((1.0d0 / MolWtDry) + sum(xza_MixRtDivMolWt,3)) &
      & * (1.0d0 + sum(xza_MixRtAll,3))
    xz_Dens = xz_PressAll * xz_EffMolWt / (GasRDry * xz_TempAll)
    
    xz_VPotTemp    = xz_PotTemp / xz_EffMolWt
    xz_VPotTempAll = xz_PotTempAll / xz_EffMolWt
    
    !----------------------------------------------------------------
    ! $BN.@~4X?t(B
    !----------------------------------------------------------------
    do k = 1, DimZMax-1
      xz_StrmFunc(:,k+1) = xz_StrmFunc(:,k) - pz_VelX(:,k) * DelZ
    end do
    
    do i = 1, DimXMax-1
      xz_StrmFunc(i+1,:) = xz_StrmFunc(i,:) + xr_VelZ(i,:) * DelX
    end do
    
    !----------------------------------------------------------------
    ! $BIbNO(B
    !----------------------------------------------------------------
    xz_BuoyTemp  =  xz_avr_xr( xr_Buoy( xz_PotTemp ) )
    xz_BuoyMolWt =  xz_avr_xr( xr_BuoyMolWt( xza_MixRt ) )
    xz_BuoyDrag  =  xz_avr_xr( xr_BuoyDrag( xza_MixRt ) )
    xz_BuoyAll   =  xz_BuoyTemp + xz_BuoyMolWt + xz_BuoyDrag
    
    !----------------------------------------------------------------
    ! $B1@:.9gHf(B, $B1+:.9gHf$r$^$H$a$k(B
    !----------------------------------------------------------------
    do s = 1, LoopNum2
      xz_Cloud = xz_Cloud + xza_MixRt(:,:,CloudNum(s))
      xz_Rain  = xz_Rain  + xza_MixRt(:,:,RainNum(s))
    end do
    
    !----------------------------------------------------------------
    ! $BK0OB>x5$05$HJ?9UDj?t(B
    !----------------------------------------------------------------
    !$BK0OB>x5$05(B
    if (LoopNum /= 0 ) then 
      do s = 1, LoopNum
        xza_MixRtSat(:,:,GasNum(s)) =                       &
          & xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll) &
          &  * MolWtWet(CloudNum(s))                        &
          &  / (MolWtDry * xz_PressAll)                     
      end do
    end if
    
    !$BJ?9UDj?t(B
    if ( NH4SHNum /= 0 ) then     
      xz_EqConstSat = 61.781d0 - 10834.0d0 / xz_TempAll - dlog(1.0d2)
      
      xz_EqConst =                                        &
        & dlog(                                           &
        &   ( xz_PressAll ** 2.0d0 )                      &
        &   * max( xza_MixRtAll(:,:,NH3Num), 1.0d-20 )    &
        &   * max( xza_MixRtAll(:,:,H2SNum), 1.0d-20 )    &
        &   * ( MolWtDry ** 2.0d0)                        &
        &   / ( MolWtWet(NH3Num) )                        &
        &   / ( MolWtWet(H2SNum) )                        &
        & )
    end if

    !$B<>EY(B 100 % $B$N;~$N:.9gHf$r5a$a$k(B
    Hum = 1.0d0
    call ECCM_MolFr( SpcWetMolFr(1:SpcNum), Hum, xz_TempBasicZ(1,:), &
      &              xz_PressBasicZ(i,:), za_MolFrEq )
    
    !$BJ,;RNL$+$i:.9gHf$X$NJQ49(B. $B?eJ?J}8~$K$O0lMM(B
    do s = 1, SpcNum
      do i = DimXMin, DimXMax
        xza_MixRtEq(i,:,s) = za_MolFrEq(:,s) * MolWtWet(s) / MolWtDry
      end do
    end do

    
    !----------------------------------------------------------------
    ! $B0BDjEY$N7W;;(B
    !----------------------------------------------------------------
    call ECCM_Stab( xz_PotTemp, xz_Exner, xza_MixRt, &
      &             xz_Stab, xz_StabTemp, xz_StabMolWt)

    !----------------------------------------------------------------    
    ! $B%U%!%$%k=PNO(B
    !----------------------------------------------------------------

    if ( mod(t, 10) == 1 ) call AnalFile_OutPut( AnalTime )
    call AnalFile_OutPutDens( AnalTime )
    
  end do
  
  !----------------------------------------------------------------    
  ! $B%U%!%$%k$rJD$8$k(B
  !----------------------------------------------------------------    
  call AnalFile_Close

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

    call HistoryClose( gt_hist(1) )
    call HistoryClose( gt_hist(2) )
    
  end subroutine AnalFile_Close

  
  subroutine AnalFile_Open( )
    !
    !$B2r@O%U%!%$%k$NDj5A(B
    !
    use gtool_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
    use basicset, only: SpcWetSymbol
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none

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

    !-----------------------------------------------------------
    ! $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,                                &
      & history=gt_hist(1) )

    !-----------------------------------------------------------  
    ! $B<4$N=PNO(B
    !-----------------------------------------------------------
    call HistoryPut('x', s_X( FileXMin: FileXMax ), &
      & history=gt_hist(1) )
    call HistoryPut('z', s_Z( FileZMin: FileZMax ), &
      & history=gt_hist(1) )

    !-----------------------------------------------------------  
    ! $B2r@OMQ$NJQ?t$N=PNO(B
    !-----------------------------------------------------------  
    !$B2>290L$N>qMp(B
    call HistoryAddVariable(                                     &
      & varname='VPotTemp', dims=(/'x','z','t'/),                &
      & longname='disturbunce of virtual potential temperature', &
      & units='K', xtype='double' , &
      & history=gt_hist(1) )

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

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

    !$B29EY>qMp(B
    call HistoryAddVariable(                         &
      & varname='Temp', dims=(/'x','z','t'/),        &
      & longname='disturbunce of temperature',       &
      & units='K', xtype='double' , &
      & history=gt_hist(1) )

    !$B29EY(B
    call HistoryAddVariable(                          &
      & varname='TempAll', dims=(/'x','z','t'/),      &
      & longname='temperature',                       &
      & units='K', xtype='double' , &
      & history=gt_hist(1) )

    !$B05NO(B
    call HistoryAddVariable(                          &
      & varname='PressAll', dims=(/'x','z','t'/),     &
      & longname='Pressure',                          &
      & units='Pa', xtype='double' , &
      & history=gt_hist(1) )

    !$BL)EY(B
    call HistoryAddVariable(                          &
      & varname='Dens', dims=(/'x','z','t'/),         &
      & longname='Density',                           &
      & units='kg m|-2"', xtype='double' , &
      & history=gt_hist(1) )

    !-----------------------------------------------------------  
    ! $BIbNO4XO"(B    
    !-----------------------------------------------------------  
    !$BA4IbNO(B
    call HistoryAddVariable(                                        &
      & varname='BuoyAll', dims=(/'x','z','t'/),                    &
      & longname='Buoyancy',                                        &
      & units=' ', xtype='double' , &
      & history=gt_hist(1) )

    !$BIbNO(B($B29EY$N4sM?(B)
    call HistoryAddVariable(                                        &
      & varname='BuoyTemp', dims=(/'x','z','t'/),                   &
      & longname='Buoyancy (Temp)',                                 &
      & units=' ', xtype='double' , &
      & history=gt_hist(1) )

    !$BIbNO(B($BJ,;RNL$N4sM?(B)
    call HistoryAddVariable(                                        &
      & varname='BuoyMolWt', dims=(/'x','z','t'/),                  &
      & longname='Buoyancy (MolWt)',                                &
      & units=' ', xtype='double' , &
      & history=gt_hist(1) )

    !$BIbNO(B($B0z$-$:$j$N4sM?(B)
    call HistoryAddVariable(                                        &
      & varname='BuoyDrag', dims=(/'x','z','t'/),                   &
      & longname='Buoyancy (Drag)',                                 &
      & units=' ', xtype='double' , &
      & history=gt_hist(1) )

    !-----------------------------------------------------------  
    ! $B0BDjEY4XO"(B
    !-----------------------------------------------------------  
    !$B0BDjEY(B
    call HistoryAddVariable(                                        &
      & varname='Stab', dims=(/'x','z','t'/),                       &
      & longname='Static Stability',                                &
      & units='s|-2" ', xtype='double' , &
      & history=gt_hist(1) )

    !$B0BDjEY(B($B29EY$N4sM?(B)
    call HistoryAddVariable(                                        &
      & varname='StabTemp', dims=(/'x','z','t'/),                   &
      & longname='Static Stability (Temp)',                         &
      & units='s|-2" ', xtype='double' , &
      & history=gt_hist(1) )

    !$B0BDjEY(B($BJ,;RNL$N4sM?(B)
    call HistoryAddVariable(                                        &
      & varname='StabMolWt', dims=(/'x','z','t'/),                  &
      & longname='Static Stability (MolWt)',                        &
      & units='s|-2" ', xtype='double' , &
      & history=gt_hist(1) )

    !-----------------------------------------------------------  
    ! $BK0OB>x5$05$HJ?9UDj?t(B
    !-----------------------------------------------------------  
    do s = 1, 2
      !$BK0OB:.9gHf(B
      call HistoryAddVariable(                                          &
        & varname='D'//trim(SpcWetSymbol(s))//'Sat',                         &
        & dims=(/'x','z','t'/),                                         &
        & longname='Saturated '//trim(SpcWetSymbol(s))//' Mixing Ratio',&
        & units='kg kg|-1"', xtype='double' , &
        & history=gt_hist(1) )
    end do

    do s = 1, 3
      !$BG.J?9U>uBV$G$N:.9gHf(B
      call HistoryAddVariable(                                          &
        & varname='D'//trim(SpcWetSymbol(s))//'Eq',                          &
        & dims=(/'x','z','t'/),                                         &
        & longname='Equilibrium '//trim(SpcWetSymbol(s))//' Mixing Ratio',&
        & units='kg kg|-1"', xtype='double' , &
        & history=gt_hist(1) )
    end do
    
    call HistoryAddVariable(                                         &
      & varname='EqConstSat', dims=(/'x','z','t'/),               &
      & longname='Equilibrium Constant of NH4SH',                     &
      & units='kg kg|-1"', xtype='double' , &
      & history=gt_hist(1) )
      
    call HistoryAddVariable(                                         &
      & varname='EqConst', dims=(/'x','z','t'/),                  &
      & longname='log (P_H2O P_NH3) ',                               &
      & units='kg kg|-1"', xtype='double' , &
      & history=gt_hist(1) )

    !----------------------------------------------------------------
    ! $B1@$H1+(B
    !----------------------------------------------------------------
    call HistoryAddVariable(                                       &
      & varname='Cloud', dims=(/'x','z','t'/),                     &
      & longname='All Clouds Mixing Ratio',                        &
      & units='kg kg|-1"', xtype='double' , &
      & history=gt_hist(1) )

    call HistoryAddVariable(                                       &
      & varname='Rain', dims=(/'x','z','t'/),                      &
      & longname='All Rains Mixing Ratio',                         &
      & units='kg kg|-1"', xtype='double' , &
      & history=gt_hist(1) )

    !----------------------------------------------------------------
    ! $BN.@~4X?t(B
    !----------------------------------------------------------------
    call HistoryAddVariable(                          &
      & varname='StrmFunc', dims=(/'x','z','t'/),       &
      & longname='Stream Function',                   &
      & units='m|2" s', xtype='double' , &
      & history=gt_hist(1) )

  end subroutine AnalFile_Open



  subroutine AnalFile_OutPut( Time )

    use gtool_history, only: HistoryPut
    use fileset, only: gt_hist    
    use gridset, only: FileXMin, FileXMax, FileZMin, FileZMax
    use basicset, only: SpcWetSymbol

    !$B0EL[$N7?@k8@6X;_(B
    implicit none
    
    !$BJQ?tDj5A(B
    real(8), intent(in) :: Time    
 
    !----------------------------------------------------------------
    ! $BCM$r=PNO(B
    !----------------------------------------------------------------
    call HistoryPut( 't', Time, history=gt_hist(1) )
    
    call HistoryPut( 'Temp',          &
      &            xz_Temp(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    call HistoryPut( 'TempAll',       &
      &            xz_TempAll(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )
    
    call HistoryPut( 'VPotTemp',      &
      &            xz_VPotTemp(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    call HistoryPut( 'VPotTempAll',   &
      &            xz_VPotTempAll(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    call HistoryPut( 'PotTempAll',    &
      &            xz_PotTempAll(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )
    
    call HistoryPut( 'PressAll',      &
      &            xz_PressAll(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )
    
    call HistoryPut( 'Dens',          &
      &            xz_Dens(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )
    
    call HistoryPut( 'BuoyAll',       &
      &            xz_BuoyAll(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    call HistoryPut( 'BuoyTemp',      &
      &            xz_BuoyTemp(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )
    
    call HistoryPut( 'BuoyMolWt',     &
      &            xz_BuoyMolWt(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )
    
    call HistoryPut( 'BuoyDrag',      &
      &            xz_BuoyDrag(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    do s = 1, 2
      call HistoryPut( trim(SpcWetSymbol(s))//'Sat',  &
        &          xza_MixRtSat(FileXMin:FileXMax, FileZMin:FileZMax, s), &
        &            history=gt_hist(1) )
    end do

    do s = 1, 3
      call HistoryPut( trim(SpcWetSymbol(s))//'Eq',  &
        &          xza_MixRtEq(FileXMin:FileXMax, FileZMin:FileZMax, s), &
        &            history=gt_hist(1) )
    end do

    call HistoryPut( 'EqConstSat', &
      &            xz_EqConstSat(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    call HistoryPut( 'EqConst',    &
      &            xz_EqConst(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    call HistoryPut( 'Cloud',         &
      &            xz_Cloud(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    call HistoryPut( 'Rain',          &
      &            xz_Rain(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    call HistoryPut( 'Stab',          &
      &            xz_Stab(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )
    
    call HistoryPut( 'StabTemp',      &
      &            xz_StabTemp(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    call HistoryPut( 'StabMolWt',     &
      &            xz_StabMolWt(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    call HistoryPut( 'StrmFunc',      &
      &            xz_StrmFunc(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )

    call HistoryPut( 'Dens',          &
      &            xz_Dens(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(1) )
    
  end subroutine AnalFile_OutPut



  subroutine AnalFile_OpenDens( )
    !
    !$B2r@O%U%!%$%k$NDj5A(B
    !
    
    use gtool_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
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none

    !$BJQ?tDj5A(B
    integer        :: s
    character(50)  :: File

    !$B%U%!%$%kL>$N;XDj(B
    File = trim(HistoryFilePrefix)// '_Dens.nc'

    !-----------------------------------------------------------
    ! $B%R%9%H%j!<:n@.(B
    !-----------------------------------------------------------
    call HistoryCreate(                              &
      & file = 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,                                &
      & history=gt_hist(2) )
    
    !-----------------------------------------------------------  
    ! $B<4$N=PNO(B
    !-----------------------------------------------------------
    call HistoryPut('x', s_X( FileXMin: FileXMax ), history=gt_hist(2) )
    call HistoryPut('z', s_Z( FileZMin: FileZMax ), history=gt_hist(2) )

    !-----------------------------------------------------------  
    ! $B2r@OMQ$NJQ?t$N=PNO(B
    !-----------------------------------------------------------  
    !$BL)EY(B
    call HistoryAddVariable(                          &
      & varname='Dens', dims=(/'x','z','t'/),         &
      & longname='Density',                           &
      & units='kg m|-2"', xtype='double', history=gt_hist(2) )

  end subroutine AnalFile_OpenDens


  subroutine AnalFile_OutPutDens( Time )

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


    !$B0EL[$N7?@k8@6X;_(B
    implicit none
    
    !$BJQ?tDj5A(B
    real(8), intent(in) :: Time    
 
    !----------------------------------------------------------------
    ! $BCM$r=PNO(B
    !----------------------------------------------------------------
    call HistoryPut( 't', Time, history=gt_hist(2) )
    
    call HistoryPut( 'Dens',          &
      &            xz_Dens(FileXMin:FileXMax, FileZMin:FileZMax), &
      &            history=gt_hist(2) )
    
  end subroutine AnalFile_OutPutDens



!!!------------------------------------------------------------------------------!!!
  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 : xz_BoundaryXCyc_xz, xz_BoundaryZSym_xz,     &
      &                       xza_BoundaryXCyc_xza, xza_BoundaryZSym_xza, &
      &                       pz_BoundaryXCyc_pz, pz_BoundaryZSym_pz,     &
      &                       xr_BoundaryXCyc_xr, 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 = '^' // 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 )
    xz_Exner = xz_BoundaryXCyc_xz( xz_Exner )
    xz_Exner = xz_BoundaryZSym_xz( xz_Exner )
    
    name = "PotTemp"
    call HistoryGet( HistoryFile(2), name, xz_PotTemp(FileXMin:FileXMax, FileZMin:FileZMax), range=step )
    xz_PotTemp = xz_BoundaryXCyc_xz( xz_PotTemp )
    xz_PotTemp = xz_BoundaryZSym_xz( xz_PotTemp )

    name = "VelX"
    call HistoryGet( HistoryFile(3), name, pz_VelX(FileXMin:FileXMax, FileZMin:FileZMax), range=step )    
    pz_VelX = pz_BoundaryXCyc_pz( pz_VelX )
    pz_VelX = pz_BoundaryZSym_pz( pz_VelX )
    
    name = "VelZ"
    call HistoryGet( HistoryFile(4), name, xr_VelZ(FileXMin:FileXMax, FileZMin:FileZMax), range=step )
    xr_VelZ = xr_BoundaryXCyc_xr( xr_VelZ )
    xr_VelZ = xr_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

    xza_MixRt = xza_BoundaryXCyc_xza( xza_MixRt )
    xza_MixRt = xza_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: xz_BoundaryXCyc_xz, xz_BoundaryZSym_xz,     &
      &                      xza_BoundaryXCyc_xza, xza_BoundaryZSym_xza, &
      &                      pz_BoundaryXCyc_pz, pz_BoundaryZSym_pz,     &
      &                      xr_BoundaryXCyc_xr, 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) )
    xz_TempBZ  = xz_BoundaryXCyc_xz( xz_TempBZ )
    xz_TempBZ  = xz_BoundaryZSym_xz( xz_TempBZ )

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

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

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

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

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

    name = "EffMolWtBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_EffMolWtBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    xz_EffMolWtBZ = xz_BoundaryXCyc_xz( xz_EffMolWtBZ )
    xz_EffMolWtBZ = xz_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

    xza_MixRtBZ = xza_BoundaryXCyc_xza( xza_MixRtBZ )
    xza_MixRtBZ = xza_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_MixRtSat(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ), &
      &       xza_MixRtEq(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ),  &
      &       za_MolFrEq(DimZMin:DimZMax, 1:SpcNum ),                    &
      &       xz_Temp(DimXMin:DimXMax, DimZMin:DimZMax ),                &
      &       xz_VPotTemp(DimXMin:DimXMax, DimZMin:DimZMax ),            &
      &       xz_VPotTempAll(DimXMin:DimXMax, DimZMin:DimZMax ),         &
      &       xz_Dens(DimXMin:DimXMax, DimZMin:DimZMax),                 &
      !
!      &       xz_Km(DimXMin:DimXMax, DimZMin:DimZMax ),                  &
      !
      &       xz_Stab(DimXMin:DimXMax, DimZMin:DimZMax ),               &
      &       xz_StabTemp(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      &       xz_StabMolWt(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 ),           &
      &       xza_MixRtDivMolWt(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum),&
      !
      &       xz_StrmFunc(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      &       xz_Cloud(DimXMin:DimXMax, DimZMin:DimZMax ),              &
      &       xz_Rain(DimXMin:DimXMax, DimZMin:DimZMax ),               &
      !
      &       xz_BuoyAll(DimXMin:DimXMax, DimZMin:DimZMax ),            &
      &       xz_BuoyTemp(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      &       xz_BuoyMolWt(DimXMin:DimXMax, DimZMin:DimZMax ),          &
      &       xz_BuoyDrag(DimXMin:DimXMax, DimZMin:DimZMax )          )

    LoopNum  = 0 
    LoopNum2 = 0
    RainNum = 0
    CloudNum = 0
    GasNum   = 0
    NH3Num   = 0
    H2SNum   = 0
    NH4SHNum = 0
    xz_StrmFunc = 0.0d0
    xza_MixRtSat = 0.0d0
    xza_MixRtDivMolWt = 0.0d0

  end subroutine ArareAlloc

end program ArareAnal
