!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Program Arare
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: arare.f90,v 1.1.1.1.2.15 2005/09/17 07:52:18 kitamo Exp $ 
!   * Tag Name: $Name: arare3m $
!   * Change History: 
!
!== Overview 
!
!$BHs@ENO3X%b%G%k(B deepconv/arare. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * $BJ}Dx<07O$O=`05=L7O(B.
!
!== Future Plans
!
                                                                 !=end
!-- (2005/04/21 $B>.9b@5;L(B) --
! $BA4BNE*$J%3%a%s%H(B
! - $B%W%m%0%i%`A4BN$NN.$l$O$3$l$G$h$$$H;W$o$l$k(B.
! - $BL?L>K!(B($B<jB3$-L>(B, $BJQ?t(B, $B2>0z?t(B)$B$O8!F$$NM>CO$"$j(B
! - $B%3%a%s%HJ8$r$I$NDxEY=q$-9~$`$+(B?
!   arare.f90 $B$+$i%3!<%I2r@b$N!V%a%$%s%W%m%0%i%`$NN.$l!W$,<+F0@8@.(B
!   $B$5$l$k$h$&$K$J$C$F$$$k$H$h$$$@$m$&(B. 
!  

program arare
                                                                 !=begin

  !== Dependency
    !-- (2005/04/21 $B>.9b@5;L(B) --
    ! use $B$5$l$kJQ?t$K$b%3%a%s%H$rIU$1$?$[$&$,$h$$$H;W$&$,(B, $B$&$k$5$$(B? 
    ! $B%3%a%s%H$O3F%5%V%k!<%A%s(B/$B%b%8%e!<%kFb$N5-=R$H0lCW$5$;$k(B.
    !
  use dc_trace, only: BeginSub, &      ! $B%G%P%C%0=PNO%5%V%k!<%A%s(B  
       &              EndSub,   &      ! $B%G%P%C%0=PNO%5%V%k!<%A%s(B  
       &              DbgMessage       ! $B%G%P%C%0=PNO%5%V%k!<%A%s(B  
  use dc_message                       ! only $BB0@-ITMW(B ?
  use fileset,  only: fileset_init,  & ! I/O $B%U%!%$%kL>=i4|2=%5%V%k!<%A%s(B  
       &              cfgfile,       & ! $B@_Dj%U%!%$%k(B
       &              ReStartFile,   & ! $BF~NOMQ%j%9%?!<%H%U%!%$%k(B
       &              NewReStartFile   ! $B=PNOMQ%j%9%?!<%H%U%!%$%k(B
  use debugset, only: debugset_init, & ! $B%G%P%C%0@_Dj=i4|2=%5%V%k!<%A%s(B
       &              DebugOn          ! $B%G%P%C%0=PNO%9%$%C%A(B
  use gridset,  only: gridset_init,  & ! $B3J;RE@>pJs=i4|2=%5%V%k!<%A%s(B
       &              DimXMin,       & ! x $BJ}8~$NG[Ns$N2<8B(B 
       &              DimXMax,       & ! x $BJ}8~$NG[Ns$N>e8B(B
       &              DimZMin,       & ! z $BJ}8~$NG[Ns$N2<8B(B
       &              DimZMax,       & ! z $BJ}8~$NG[Ns$N>e8B(B
       &              RegXMin,       & ! x $BJ}8~$NJ*M}NN0h$N2<8B(B
       &              RegXMax,       & ! x $BJ}8~$NJ*M}NN0h$N>e8B(B 
       &              RegZMin,       & ! z $BJ}8~$NJ*M}NN0h$N2<8B(B
       &              RegZMax,       & ! z $BJ}8~$NJ*M}NN0h$N>e8B(B 
       &              DelX,          & ! x $BJ}8~$N3J;R4V3V(B
       &              DelZ,          & ! z $BJ}8~$N3J;R4V3V(B
       &              s_Z              ! z $B:BI8(B
  use timeset,  only: timeset_init,  & ! $B;~4V@QJ,@_Dj=i4|2=%5%V%k!<%A%s(B
       &              NstepShort,    & ! $BD9$$%?%$%`%9%F%C%W$N%9%F%C%W?t(B
       &              NstepLong,     & ! $BC;$$%?%$%`%9%F%C%W$N%9%F%C%W?t(B
       &              DelTimeShort,  & ! $BD9$$%?%$%`%9%F%C%W(B
       &              DelTimeLong,   & ! $BC;$$%?%$%`%9%F%C%W(B
       &              TimeDisp         ! $B%U%!%$%k=PNO$9$k;~4V4V3V(B
  use nameset,  only: nameset_init     ! $B<B83L>=i4|2=%5%V%k!<%A%s(B
  use physset,  only: physset_init,  & ! $BJ*M}Dj?t$N=i4|2=%5%V%k!<%A%s(B
       &              PressSfc,      & ! $BCOI=LL05NO(B
       &              GasR,          & ! $B5$BNDj?t(B
       &              Grav             ! $B=ENO(B
  use linlib,   only: linlib_init      ! $B9TNs7W;;%i%$%V%i%j=i4|2=%5%V%k!<%A%s(B
  use bcset,    only: bcset_init       ! $B6-3&>r7o@_Dj%5%V%k!<%A%s(B
  use basicset, only: basicset_file_init, &  ! $B4pK\>l@_Dj%5%V%k!<%A%s(B
       &              ss_DensBasicZ,      &  ! $B4pK\>l$NL)EY(B
       &              ss_PotTempBasicZ,   &  ! $B4pK\>l$NL)EY(B
       &              ss_ExnerBasicZ,     &  ! $B4pK\>l$NL)EY(B
       &              ss_CpBasicZ,        &  ! $BDj05HfG.(B 
       &              ss_VelSoundBasicZ      ! $B2;B.(B
  use arareset, only: arareset_init    ! $B7W;;@_Dj$rFI$_9~$`%5%V%k!<%A%s(B
  use average,  only: ss_avr_fs,     & 
                        ! x $BJ}8~%U%i%C%/%93J;RE@$+$i%9%+%i!<3J;RE@$X$NJ?6QA`:n(B
                      ss_avr_sf    
                        ! z $BJ}8~%U%i%C%/%93J;RE@$+$i%9%+%i!<3J;RE@$X$NJ?6QA`:n(B
  use summation,only: sumxz            ! $BNN0hA4BN$NOB$r7W;;(B
  use gt4_history
                                                                 !=end
  !== $B0EL[$N7?@k8@6X;_(B
  implicit none

  !== $BJQ?tG[Ns(B
  real(8), allocatable  :: fs_VelX_as(:,:)    ! $B?eJ?B.EY(B      [t=(n-1)$B&$&S(B]
  real(8), allocatable  :: fs_VelX_ns(:,:)    ! $B?eJ?B.EY(B      [t=n$B&$&S(B]
  real(8), allocatable  :: sf_VelZ_as(:,:)    ! $B1tD>B.EY(B      [t=(n-1)$B&$&S(B]
  real(8), allocatable  :: sf_VelZ_ns(:,:)    ! $B1tD>B.EY(B      [t=n$B&$&S(B]
  real(8), allocatable  :: ss_Exner_ns(:,:)   ! $B%(%/%9%J!<4X?t(B[t=n$B&$&S(B]
  real(8), allocatable  :: ss_Exner_as(:,:)   ! $B%(%/%9%J!<4X?t(B[t=(n-1)$B&$&S(B]

  real(8), allocatable  :: ss_PotTemp_bl(:,:) ! $B290L(B [t=(n-1)$B&$(Bt]
  real(8), allocatable  :: ss_PotTemp_nl(:,:) ! $B290L(B [t=n$B&$(Bt    ]
  real(8), allocatable  :: ss_PotTemp_al(:,:) ! $B290L(B [t=(n+1)$B&$(Bt]
  real(8), allocatable  :: ss_Exner_bl(:,:)   ! $B%(%/%9%J!<4X?t(B [t=(n-1)$B&$(Bt]
  real(8), allocatable  :: ss_Exner_nl(:,:)   ! $B%(%/%9%J!<4X?t(B [t=n$B&$(Bt    ]
  real(8), allocatable  :: ss_Exner_al(:,:)   ! $B%(%/%9%J!<4X?t(B [t=(n+1)$B&$(Bt]
  real(8), allocatable  :: fs_VelX_bl(:,:)    ! $B?eJ?B.EY(B [t=(n-1)$B&$(Bt]
  real(8), allocatable  :: fs_VelX_nl(:,:)    ! $B?eJ?B.EY(B [t=n$B&$(Bt    ]
  real(8), allocatable  :: fs_VelX_al(:,:)    ! $B?eJ?B.EY(B [t=(n+1)$B&$(Bt]
  real(8), allocatable  :: sf_VelZ_bl(:,:)    ! $B1tD>B.EY(B [t=(n-1)$B&$(Bt]
  real(8), allocatable  :: sf_VelZ_nl(:,:)    ! $B1tD>B.EY(B [t=n$B&$(Bt    ]
  real(8), allocatable  :: sf_VelZ_al(:,:)    ! $B1tD>B.EY(B [t=(n+1)$B&$(Bt]

  real(8), allocatable  :: fs_AdvDiffX_nl(:,:)! $B?eJ?J}8~$N1?F0J}Dx<0$N0\N.9`$H3H;69`(B   [t=n$B&$(Bt]
  real(8), allocatable  :: sf_AdvDiffZ_nl(:,:)! $B1tD>J}8~$N1?F0J}Dx<0$N0\N.9`$H3H;69`(B   [t=n$B&$(Bt]
  
  real(8), allocatable  :: ss_Km_bl(:,:)      ! $B12G4@-78?t(B [t=(n-1)$B&$(Bt]
  real(8), allocatable  :: ss_Km_nl(:,:)      ! $B12G4@-78?t(B [t=n$B&$(Bt    ]
  real(8), allocatable  :: ss_Km_al(:,:)      ! $B12G4@-78?t(B [t=(n+1)$B&$(Bt]
  real(8), allocatable  :: ss_Kh_bl(:,:)      ! $B123H;678?t(B [t=(n-1)$B&$(Bt]
  real(8), allocatable  :: ss_Kh_nl(:,:)      ! $B123H;678?t(B [t=n$B&$(Bt    ]
  real(8), allocatable  :: ss_Kh_al(:,:)      ! $B123H;678?t(B [t=(n+1)$B&$(Bt]

  real(8), allocatable  :: sf_SfcTemp_al(:)      ! $BCOI=LL29EY(B [t=(n+1)$B&$(Bt    ]
  real(8), allocatable  :: sf_SfcHeatFlux(:)     ! $BCOI=LLG.%U%i%C%/%9(B
  real(8), allocatable  :: ff_SfcMomentumFlux(:) ! $BCOI=LL1?F0NL%U%i%C%/%9(B
  real(8), allocatable  :: ss_QRad(:,:)          ! $BJ|<M2CG.(B
  real(8), allocatable  :: ss_QDis(:,:)          ! $B;60o2CG.(B

  real(8), allocatable  :: ss_PotTempTendAdv_nl(:,:)  ! $B0\N.$K$h$k290LJQ2=(B
  real(8), allocatable  :: ss_PotTempTendSrc_nl(:,:)  ! $B@8@.$K$h$k290LJQ2=(B
  real(8), allocatable  :: ss_PotTempTendTurb_nl(:,:) ! $BMpN.3H;6$K$h$k290LJQ2=(B
  real(8), allocatable  :: s_PotTempTendAdvZ(:)  ! $B0\N.$K$h$k290LJQ2=(B($B?eJ?J?6Q(B)
  real(8), allocatable  :: s_PotTempTendSrcZ(:)  ! $B@8@.$K$h$k290LJQ2=(B($B?eJ?J?6Q(B)
  real(8), allocatable  :: s_PotTempTendTurbZ(:) ! $BMpN.3H;6$K$h$k290LJQ2=(B($B?eJ?J?6Q(B)
  real(8), allocatable  :: s_PotTempTendRadZ(:)  ! $BJ|<M2CG.$K$h$k290LJQ2=(B($B?eJ?J?6Q(B)
  real(8), allocatable  :: s_PotTempTendDisZ(:)  ! $B;60o2CG.$K$h$k290LJQ2=(B($B?eJ?J?6Q(B)
  real(8), allocatable  :: s_PotTempTendSfcZ(:)  ! $BCOI=LL2CG.$K$h$k290LJQ2=(B($B?eJ?J?6Q(B)

  real(8), allocatable  :: ss_Mass_nl(:,:)            ! $B<ANL(B
  real(8), allocatable  :: ss_KineticEnergy_nl(:,:)   ! $B1?F0%(%M%k%.!<(B
  real(8), allocatable  :: ss_PotentialEnergy_nl(:,:) ! $B%]%F%s%7%c%k%(%M%k%.!<(B
  real(8)               :: MassTotal             ! $BNN0hA4BN$N<ANL(B
  real(8)               :: KineticEnergyTotal    ! $BNN0hA4BN$N1?F0%(%M%k%.!<(B
  real(8)               :: PotentialEnergyTotal  ! $BNN0hA4BN$N%]%F%s%7%c%k%(%M%k%.!<(B

  real(8)               :: Time               ! $B7W;;7QB3;~4V(B
  real(8)               :: DelTime            ! $B;~4V%9%F%C%W(B
  real(8)               :: ReStartTime(2)     ! $B%j%9%?!<%H%U%!%$%k$N;~9o(B
  integer               :: t, tau             ! $B%k!<%W$r2s$9$?$a$NJQ?t(B
  integer               :: NStep_leapfrog     ! $B%k!<%W$r2s$92s?t(B

  integer               :: StartRunTime(8)    ! $B7W;;3+;O;~9o(B
  integer               :: EndRunTime(8)      ! $B7W;;=*N;;~9o(B
  integer               :: RunTime(8)         ! $B<B9T;~4V(B
  integer               :: i,j                ! $B%k!<%W$r2s$9$?$a$NJQ?t(B

  type(GT_HISTORY)   :: hst_pred, hst_diag    ! gt4_history $B$N%R%9%H%jL>(B

  !== $B=i4|2=(B

  !=== $B3+;O;~9o$N5-O?(B
  !fortran90 $B$NAH$_9~$_4X?t(B date_and_time $B$rMxMQ(B.
  call date_and_time(values = StartRunTime)

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

  !=== $B%G%P%0@_Dj$N=i4|2=(B
  !NAMELIST $B$+$i>pJs$rF@$F(B, $B%G%P%C%0=PNO%9%$%C%A$N@ZBX$($r9T$&(B.  
  call debugset_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)
  
  !=== $B;~4V@QJ,$N@_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)
  
  !=== $B<B83L>$N=i4|2=(B
  !netCDF $B$K=q$-9~$`<B83L>$r(B NAMELIST $B$+$iFI$_9~$`(B
  call nameset_init(cfgfile)
  
  !=== $BJ*M}Dj?t$N=i4|2=(B
  !NAMELIST $B$+$iJ*M}Dj?t$rFI$_9~$`(B
  call physset_init(cfgfile)
  
  !=== $B7W;;$N@_Dj$rFI$_9~$`(B
  !NameList $B%U%!%$%k$+$i(B arare $B<B9T$KI,MW$J%Q%i%a%?$r<hF@$9$k(B
  call arareset_init(cfgfile)

  !=== $B6-3&>r7o$r@_Dj(B
  !NAMELIST $B$+$i6-3&>r7o$N%?%$%W$r@_Dj$9$k(B. 
  call bcset_init(cfgfile)

  !=== $B=PNO%U%!%$%k$NDj5A(B
  call OpenDataFile(hst_pred)
  call OpenDiagDataFile(hst_diag)

!== $BJQ?t$N=i4|2=(B

  !=== $BG[Ns$N3d$jEv$F(B, $B=i4|2=(B
  call ArareAlloc

  !=== $B4pK\>l$NCM$r@_Dj(B
  !$B4pK\>l$NCM$r(B netCDF $B%U%!%$%k$+$i<hF@$9$k(B. 
  call basicset_file_init( ReStartFile )

  !=== $B>qMp>l$NCM$r@_Dj(B
  !$B>qMp>l$N=i4|CM$r(B netCDF $B%U%!%$%k$+$i<hF@$9$k(B. 
  call getdisturbvar( ReStartFile, ReStartTime,   &
       & fs_VelX_bl,     sf_VelZ_bl,  ss_Exner_bl,  &
       & ss_PotTemp_bl,  ss_Km_bl,    ss_Kh_bl,     &
       & fs_VelX_nl,     sf_VelZ_nl,  ss_Exner_nl,  &
       & ss_PotTemp_nl,  ss_Km_nl,    ss_Kh_nl      )

  !=== $B9TNs7W;;%i%$%V%i%j$N=i4|2=(B. 
  !$BMxMQ$9$k@~7A7W;;%Q%C%1!<%8$r(B NAMELIST $B$+$i@_Dj$7(B, 
  !$BG[Ns$NBg$-$5$r7hDj$9$k(B
  !  * $B0z?t$O%9%+%i!<$N(B z $BJ}8~$N7W;;NN0h$NBg$-$5(B
  call linlib_init(cfgfile, RegZMax - RegZMin)
  
  !=== $B;~9o$N=i4|2=(B
  Time = ReStartTime(2)          !$B%j%9%?!<%H3+;O;~9o(B
  DelTime = 2.0 * DelTimeLong    !leap-frog $B%9%-!<%`$G$N@QJ,;~4V(B

  !=== $B%k!<%W2s?t$N=i4|2=(B
  !$B;~9o(B 0 $BIC$+$i$N7W;;(B, $B$9$J$o$A(B ReStartTime(1) = 0 $B$G$O(B, 
  !$B;~9o(B 0 $BIC$+$i@QJ,;~4V(B TimeInt $B$^$G@QJ,$9$k(B. 
  !$B$3$N$H$-(B leap-frog $B%9%-!<%`$O(B NstepLong - 1 $B2s%k!<%W$5$;$k$3$H$H$J$k(B. 
  !
  !$B7QB37W;;(B, $B$9$J$o$A(B ReStartTime(1) /= 0 $B$G$O(B, 
  !RestartTime(2) $BIC$+$i@QJ,;~4V(B TimeInt $B$^$G@QJ,$9$k(B. 
  !$B$3$N$H$-(B leap-frog $B%9%-!<%`$O(B NstepLong1 $B2s%k!<%W$5$;$k$3$H$H$J$k(B. 
  !$B7QB37W;;$N>l9g$K$O%j%9%?!<%H%U%!%$%k$N;~9o$+$i7W;;$r9T$$$?$$$+$i$G$"$k(B.

  if ( ReStartTime(1) == 0.0 ) then 
     Nstep_leapfrog = NstepLong - 1
  else
     Nstep_leapfrog = NstepLong
  end if

  !=== $B;~9o(B 0 $B$+$i$N7W;;$G$O(B, $BI,MW$K1~$8$F=i4|CM$r%U%!%$%k=PNO(B
  if ( ReStartTime(1) == 0.0 ) then 
     call OutputDataFile( hst_pred, ReStartTime(1), &
          & fs_VelX_bl, sf_VelZ_bl, ss_Exner_bl, &
          & ss_PotTemp_bl, ss_Km_bl, ss_Kh_bl )

     if ( mod( ReStartTime(2), TimeDisp ) == 0 ) then 
        call OutputDataFile( hst_pred,  ReStartTime(2), &
             & fs_VelX_nl, sf_VelZ_nl, ss_Exner_nl, &
             & ss_PotTemp_nl, ss_Km_nl, ss_Kh_nl )
     end if
  end if


!== CFL $B>r7o$N%A%'%C%/(B

  call cflcheck


!== $BD9$$;~4V%9%F%C%W(B. leap-frog $B%9%-!<%`$rMxMQ(B
  do t = 1, Nstep_leapfrog
     Time = Time + DelTimeLong
     write(*,*) Time

     !=== $B12G4@-78?t(B, $B123H;678?t$N7W;;(B
     call EddyViscosity(                                     &
          & ss_Km_bl,                                        & 
          & DelTime,                                         & !$B;~4V4V3V(B
          & fs_VelX_nl, sf_VelZ_nl, ss_Km_nl,                & !$B0\N.(B
          & fs_VelX_bl, sf_VelZ_bl, ss_PotTemp_bl, ss_Km_bl, & !$B3H;6Ey(B
          & ss_Km_al, ss_Kh_al )

     !=== $BD9$$;~4V%9%F%C%W$G$N290L$N7W;;(B
     call Scalar(                                   & 
          & ss_PotTemp_bl,                          &
          & DelTime,                                & !$B;~4V4V3V(B
          & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl,  & !$B0\N.(B
          & ss_PotTemp_bl, ss_Kh_bl,                & !$BMpN.3H;6(B, $BEy(B
          & ss_PotTemp_al,                          &  
          & ss_PotTempTendAdv_nl,                  & !$B0\N.(B
          & ss_PotTempTendSrc_nl,                  & !$B@8@.(B
          & ss_PotTempTendTurb_nl)                    !$BMpN.3H;6(B

     !=== $BCOI=LL29EY$N7W;;(B    
     call SurfaceTemp(sf_SfcTemp_al)

     !=== $BCOLL%U%i%C%/%9$N7W;;(B    
     call SurfaceHeatFlux(ss_CpBasicZ,     &
          &               ss_DensBasicZ,   &
          &               ss_PotTempBasicZ, &
          &               ss_ExnerBasicZ,   &
          &               fs_VelX_bl,       &
          &               ss_PotTemp_bl,    &
          &               ss_Exner_bl,      &
          &               sf_SfcTemp_al,    &
          &               sf_SfcHeatFlux )

     !=== $BJ|<MNd5QN($N7W;;(B
     call Radiation(ss_Exner_bl, ss_ExnerBasicZ, ss_QRad)

     !=== $B;60o2CG.$N7W;;(B
     call DissipativeHeating(ss_Km_bl, ss_QDis)

     !=== $B%U%i%C%/%9(B, $BNd5Q(B, $B;60o2CG.$K$h$k290LJQ2=$r7W;;(B
     ss_PotTemp_al = ss_PotTemp_al + DelTime * ss_QRad

     ss_PotTemp_al = ss_PotTemp_al + DelTime * ss_QDis

     ss_PotTemp_al(:,RegZMin + 1) =   &
         &   ss_PotTemp_al(:,RegZMin + 1)        &
         & + DelTime * sf_SfcHeatFlux            & 
         &   / ( DelZ                            &
         &       * ss_DensBasicZ(:,RegZMin + 1)  & 
         &       * ss_CpBasicZ(:,RegZMin + 1)    &
         &      )

     !=== $B290LJQ2=$N4sM?$N?eJ?J?6QCM$r7W;;(B
     do i = RegXMin, RegXMax
        s_PotTempTendAdvZ(:) = s_PotTempTendAdvZ(:)  &
             & -  ss_PotTempTendAdv_nl(i,:) / (RegXMax - RegXMin + 1)

        s_PotTempTendSrcZ(:) = s_PotTempTendSrcZ(:)  &
             & -  ss_PotTempTendSrc_nl(i,:) / (RegXMax - RegXMin + 1)

        s_PotTempTendTurbZ(:) = s_PotTempTendTurbZ(:)  &
             & +  ss_PotTempTendTurb_nl(i,:) / (RegXMax - RegXMin + 1)

        s_PotTempTendRadZ(:) = s_PotTempTendRadZ(:)   &
             & + ss_QRad(i, :) / (RegXMax - RegXMin + 1)

        s_PotTempTendDisZ(:) = s_PotTempTendDisZ(:)   &
             & + ss_QDis(i, :) / (RegXMax - RegXMin + 1)

        s_PotTempTendSfcZ(RegZMin + 1) =   &
             & s_PotTempTendSfcZ(RegZMin + 1)   &
             & + sf_SfcHeatFlux(i)   &
             & / ( (RegXMax - RegXMin + 1)  &
             &     * DelZ   &
             &     * ss_DensBasicZ(i, RegZMin + 1)  &
             &     * ss_CpBasicZ(i, RegZMin + 1) )

     end do

     !=== $BD9$$;~4V%9%F%C%W$G$NB.EY(B u $B$N0\N.(B, $B3H;6(B, $BG4@-(B
     call AdvectDiffuseX(                     &
          & fs_VelX_nl, sf_VelZ_nl,           & !$B0\N.(B
          & fs_VelX_bl, sf_VelZ_bl, ss_Km_bl, & !$BMpN.3H;6(B, $B?tCMG4@-(B
          & fs_AdvDiffX_nl )
     
     !=== $BD9$$;~4V%9%F%C%W$G$NB.EY(B w $B$N0\N.(B, $B3H;6(B, $BG4@-(B
     call AdvectDiffuseZ(                          &
         & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, & !$B0\N.(B, $BIbNO(B
         & fs_VelX_bl, sf_VelZ_bl, ss_Km_bl,      & !$BMpN.3H;6(B, $B?tCMG4@-(B
         & sf_AdvDiffZ_nl )
     
    
     !=== $BC;$$%?%$%`%9%F%C%W$N=i4|CM$r:n$k(B
     fs_VelX_ns  = fs_VelX_bl
     sf_VelZ_ns  = sf_VelZ_bl
     ss_Exner_ns = ss_Exner_bl 
    
     !=== $BC;$$%?%$%`%9%F%C%W$G$N@QJ,(B
     do tau = 1, NstepShort 
        !==== $BB.EY(B u $B$N7W;;(B
        call VelocityX(                                             &
             & fs_VelX_ns,                                          &
             & DelTimeShort,                                        &
             & fs_AdvDiffX_nl, fs_VelX_ns, sf_VelZ_ns, ss_Exner_ns, &
             & fs_VelX_as )
        
        !==== $B%(%/%9%J!<4X?t$N7W;;(B
        call Exner(                                    &
             & sf_AdvDiffZ_nl, fs_VelX_ns, fs_VelX_as, &
             &   sf_VelZ_ns, ss_Exner_ns,              &
             & ss_Exner_as )

        !==== $BB.EY(B w $B$N7W;;(B
        call VelocityZ(                                &
             & sf_VelZ_ns,                             &
             & DelTimeShort,                           &
             & sf_AdvDiffZ_nl, fs_VelX_ns, sf_VelZ_ns, &
             &    ss_Exner_ns, ss_Exner_as,            &
             & sf_VelZ_as )

        !=== $BCOLL%U%i%C%/%9$N7W;;(B    
        call SurfaceMomentumFlux(ss_DensBasicZ,  &
                                 fs_VelX_ns,      &  ! ns ?
                                 ff_SfcMomentumFlux )

        !=== $B%U%i%C%/%9$K$h$kB.EYJQ2=$r7W;;(B
        fs_VelX_as(:, RegZMin + 1) = fs_VelX_as(:, RegZMin + 1)   &
             & + DelTimeShort * ff_SfcMomentumFlux                &
             &     / ( ss_DensBasicZ(:,RegZMin + 1) * DelZ )
        
        !==== $BC;$$;~4V%9%F%C%W$N%k!<%W$r2s$9$?$a$N=hCV(B
        ss_Exner_ns  = ss_Exner_as 
        fs_VelX_ns   = fs_VelX_as
        sf_VelZ_ns   = sf_VelZ_as
        
     end do
     
     !=== $B:G=*E*$JC;$$;~4V%9%F%C%W$G$NCM$rD9$$;~4V%9%F%C%W$G$NCM$H$_$J$9(B
     ss_Exner_al  = ss_Exner_as 
     fs_VelX_al   = fs_VelX_as
     sf_VelZ_al   = sf_VelZ_as     
     
     !=== $B%"%C%;%j%s$N;~4V%U%#%k%?(B
     call AsselinTimeFilter(ss_PotTemp_al, ss_PotTemp_nl, ss_PotTemp_bl )
     call AsselinTimeFilter(ss_Exner_al,   ss_Exner_nl,   ss_Exner_bl )
     call AsselinTimeFilter(fs_VelX_al,    fs_VelX_nl,    fs_VelX_bl )
     call AsselinTimeFilter(sf_VelZ_al,    sf_VelZ_nl,    sf_VelZ_bl )

     !=== $BD9$$;~4V%9%F%C%W$N%k!<%W$r2s$9$?$a$N=hCV(B
     ss_PotTemp_bl = ss_PotTemp_nl
     ss_Exner_bl   = ss_Exner_nl 
     fs_VelX_bl    = fs_VelX_nl
     sf_VelZ_bl    = sf_VelZ_nl
     ss_Km_bl      = ss_Km_nl
     ss_Kh_bl      = ss_Kh_nl
     ss_PotTemp_nl = ss_PotTemp_al
     ss_Exner_nl   = ss_Exner_al 
     fs_VelX_nl    = fs_VelX_al
     sf_VelZ_nl    = sf_VelZ_al
     ss_Km_nl      = ss_Km_al
     ss_Kh_nl      = ss_Kh_al

     !=== $BJ]B8NL$N7W;;(B

     ss_Mass_nl =  (  &
          &            ss_DensBasicZ   &
          &          - ss_DensBasicZ * ss_PotTemp_nl / ss_PotTempBasicZ  &
          &          + PressSfc * ss_CpBasicZ    &
          &            * ss_ExnerBasicZ**((ss_CpBasicZ - GasR)/ GasR) &
          &            * ss_Exner_nl  &
          &            / (ss_VelSoundBasicZ**2 * GasR)  &
          &        ) * DelX * DelZ
     
     MassTotal = sumxz(ss_Mass_nl)

     ss_KineticEnergy_nl =   &
          & ss_DensBasicZ  &
          & * (  &
          &       ss_avr_fs(fs_VelX_nl)**2  &
          &     + ss_avr_sf(sf_VelZ_nl)**2  &
          &   ) * DelX * DelZ

     KineticEnergyTotal = sumxz(ss_KineticEnergy_nl)

     do i = DimXMin, DimXMax
        ss_PotentialEnergy_nl(i,DimZMin:DimZMax) =  &
             & ss_DensBasicZ(i,DimZMin:DimZMax)  &
             & * (    &
             &       Grav * s_Z(DimZMin:DimZMax)   &
             &     + ss_CpBasicZ(i,DimZMin:DimZMax)   &
             &       * (ss_PotTemp_nl(i,DimZMin:DimZMax) + ss_PotTempBasicZ(i,DimZMin:DimZMax))  &
             &       * (ss_Exner_nl(i,DimZMin:DimZMax) + ss_ExnerBasicZ(i,DimZMin:DimZMax) )  &
             &   ) * DelX * DelZ
     end do

     PotentialEnergyTotal = sumxz(ss_PotentialEnergy_nl)

!     write(*,*) 'Mass', MassTotal
!     write(*,*) 'K.E.', KineticEnergyTotal
!     write(*,*) 'P.E.', PotentialEnergyTotal
!     write(*,*) 'T.E.', KineticEnergyTotal + PotentialEnergyTotal

     !=== $B%U%!%$%k=PNO(B
     if ( mod( Time, TimeDisp ) == 0 ) then 
        call DbgMessage(fmt="Time = %f", d=(/Time/))
        call OutputDataFile(hst_pred, Time, &
             &  fs_VelX_nl,    sf_VelZ_nl, ss_Exner_nl, &
             &  ss_PotTemp_nl, ss_Km_nl,   ss_Kh_nl )

!     write(*,*) Time

        !--- $B?GCGNL=PNO(B
        s_PotTempTendAdvZ = s_PotTempTendAdvZ * DelTimeLong / TimeDisp
        s_PotTempTendSrcZ = s_PotTempTendSrcZ * DelTimeLong / TimeDisp
        s_PotTempTendTurbZ = s_PotTempTendTurbZ * DelTimeLong / TimeDisp
        s_PotTempTendRadZ = s_PotTempTendRadZ * DelTimeLong / TimeDisp
        s_PotTempTendDisZ = s_PotTempTendDisZ * DelTimeLong / TimeDisp
        s_PotTempTendSfcZ = s_PotTempTendSfcZ * DelTimeLong / TimeDisp

        call OutputDiagDataFile(hst_diag, &
             & Time,   &
             & s_PotTempTendAdvZ,    &
             & s_PotTempTendSrcZ,    &
             & s_PotTempTendTurbZ,   &
             & s_PotTempTendRadZ,    &
             & s_PotTempTendDisZ,    &
             & s_PotTempTendSfcZ,    &
             & MassTotal,            &
             & KineticEnergyTotal,   &
             & PotentialEnergyTotal)

        s_PotTempTendAdvZ = 0.0d0
        s_PotTempTendSrcZ = 0.0d0
        s_PotTempTendTurbZ = 0.0d0
        s_PotTempTendRadZ = 0.0d0
        s_PotTempTendDisZ = 0.0d0
        s_PotTempTendSfcZ = 0.0d0
     end if
     
  end do
  
  !=== $B=PNO%U%!%$%k$N%/%m!<%:(B
  call CloseDataFile(hst_pred)
  call CloseDiagDataFile(hst_diag)

  !=== $B%j%9%?!<%H%U%!%$%k$N:n@.(B
  ReStartTime(1) = Time - DelTimeLong
  ReStartTime(2) = Time
  call MakeReStartFile( NewReStartFile, ReStartTime, &
       & fs_VelX_bl,    sf_VelZ_bl,  ss_Exner_bl,    &
       & ss_PotTemp_bl, ss_Km_bl,    ss_Kh_bl,       &
       & fs_VelX_nl,    sf_VelZ_nl,  ss_Exner_nl,    &
       & ss_PotTemp_nl, ss_Km_nl,    ss_Kh_nl        )


  !=== $B=*N;;~9o$N5-O?$H%W%m%0%i%`<B9T;~4V$N=q$-=P$7(B
  !fortran90 $B$NAH$_9~$_4X?t(B date_and_time $B$rMxMQ(B.
  call date_and_time(values = EndRunTime)

  RunTime = EndRunTime - StartRunTime
  do i = 6,7
     if ( RunTime(i) < 0 ) then
        RunTime(i - 1) = RunTime(i - 1) - 1
        RunTime(i)     = RunTime(i) + 60
     end if
  end do
  write(*,*) RunTime(5), 'hour', RunTime(6), 'min', RunTime(7), 'sec'

contains

!!!
!!!$BG[Ns$N3d$jEv$F(B, $B=i4|2=(B
!!!
  subroutine ArareAlloc

    call BeginSub("ArareAlloc", &
&                 fmt="%c",     &
&                 c1="Allocate and initialize allocatable variables.")

    allocate(ss_Exner_as(DimXMin:DimXMax, DimZMin:DimZMax),  &
         & fs_VelX_as(DimXMin:DimXMax, DimZMin:DimZMax)  , &
         & sf_VelZ_as(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Exner_ns(DimXMin:DimXMax, DimZMin:DimZMax), &
         & fs_VelX_ns(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_VelZ_ns(DimXMin:DimXMax, DimZMin:DimZMax), &
         & fs_VelX_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_VelZ_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Exner_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_PotTemp_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & fs_VelX_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_VelZ_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Exner_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_PotTemp_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & fs_VelX_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_VelZ_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Exner_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_PotTemp_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & fs_AdvDiffX_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_AdvDiffZ_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_SfcTemp_al(DimXMin:DimXMax),      &
         & sf_SfcHeatFlux(DimXMin:DimXMax),     &
         & ff_SfcMomentumFlux(DimXMin:DimXMax),   &
         & ss_QRad(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_QDis(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Km_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Km_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Km_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Kh_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Kh_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Kh_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_PotTempTendAdv_nl(DimXMin:DimXMax, DimZMin:DimZMax),   & 
         & ss_PotTempTendSrc_nl(DimXMin:DimXMax, DimZMin:DimZMax),   &
         & ss_PotTempTendTurb_nl(DimXMin:DimXMax, DimZMin:DimZMax),   &
         & s_PotTempTendAdvZ(DimZMin:DimZMax),   &
         & s_PotTempTendSrcZ(DimZMin:DimZMax),   &
         & s_PotTempTendTurbZ(DimZMin:DimZMax),   &
         & s_PotTempTendRadZ(DimZMin:DimZMax),  &
         & s_PotTempTendDisZ(DimZMin:DimZMax),  &
         & s_PotTempTendSfcZ(DimZMin:DimZMax),  &
         & ss_Mass_nl(DimXMin:DimXMax, DimZMin:DimZMax),  &
         & ss_KineticEnergy_nl(DimXmin:DimXMax, DimZMin:DimZMax),  &
         & ss_PotentialEnergy_nl(DimXMin:DimXMax, DimZMin:DimZMax) )

    fs_VelX_al = 0.0d0
    fs_VelX_nl = 0.0d0
    fs_VelX_bl = 0.0d0

    sf_VelZ_al = 0.0d0
    sf_VelZ_nl = 0.0d0
    sf_VelZ_bl = 0.0d0

    ss_Exner_al = 0.0d0
    ss_Exner_nl = 0.0d0
    ss_Exner_bl = 0.0d0
    
    ss_PotTemp_al = 0.0d0
    ss_PotTemp_nl = 0.0d0
    ss_PotTemp_bl = 0.0d0

    fs_VelX_as = 0.0d0
    fs_VelX_ns = 0.0d0

    sf_VelZ_as = 0.0d0
    sf_VelZ_ns = 0.0d0

    ss_Exner_as = 0.0d0
    ss_Exner_ns = 0.0d0

    fs_AdvDiffX_nl = 0.0d0
    sf_AdvDiffZ_nl = 0.0d0

    sf_SfcTemp_al      = 0.0d0
    sf_SfcHeatFlux     = 0.0d0 
    ff_SfcMomentumFlux = 0.0d0 
    ss_QRad = 0.0d0
    ss_QDis = 0.0d0

    ss_Km_bl = 0.0d0
    ss_Km_nl = 0.0d0
    ss_Km_al = 0.0d0
    ss_Kh_bl = 0.0d0
    ss_Kh_nl = 0.0d0
    ss_Kh_al = 0.0d0

    Time = 0.0d0

    !--- $B?GCGNL=PNO%U%!%$%kMQ$K;H$o$l$kJQ?tG[Ns(B

    ss_PotTempTendAdv_nl = 0.0d0
    ss_PotTempTendSrc_nl = 0.0d0
    ss_PotTempTendTurb_nl = 0.0d0
    s_PotTempTendAdvZ = 0.0d0
    s_PotTempTendSrcZ = 0.0d0
    s_PotTempTendTurbZ = 0.0d0
    s_PotTempTendRadZ = 0.0d0
    s_PotTempTendDisZ = 0.0d0
    s_PotTempTendSfcZ = 0.0d0

    ss_Mass_nl = 0.0d0
    ss_KineticEnergy_nl = 0.0d0
    ss_PotentialEnergy_nl = 0.0d0

    call EndSub("ArareAlloc")

  end subroutine ArareAlloc

end program arare

