!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Program Arare
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $ 
!   * Tag Name: $Name:  $
!   * 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

program arare
                                                                 !=begin
  !== Dependency
  use dc_trace, only: BeginSub, EndSub, DbgMessage
  use dc_message
  use fileset,  only: fileset_init, cfgfile, ReStartFile, NewReStartFile
  use debugset, only: debugset_init, DebugOn
  use gridset,  only: gridset_init, DimXMin, DimXMax, DimZMin, DimZMax, &
       &              RegZMin, RegZMax
  use timeset,  only: timeset_init, NstepShort, NstepLong, &
       &              DelTimeShort, DelTimeLong, TimeDisp
  use nameset,  only: nameset_init
  use physset,  only: physset_init
  use linlib,   only: linlib_init
  use bcset,    only: bcset_init
  use basicset, only: basicset_file_init
  use arareset, only: arareset_init 
                                                                 !=end
  implicit none

  real(8), allocatable  :: ss_Exner_as(:,:)
  real(8), allocatable  :: fs_VelX_as(:,:)
  real(8), allocatable  :: sf_VelZ_as(:,:)
  real(8), allocatable  :: ss_Exner_ns(:,:)
  real(8), allocatable  :: fs_VelX_ns(:,:)
  real(8), allocatable  :: sf_VelZ_ns(:,:)

  real(8), allocatable  :: fs_AdvDiffX_nl(:,:)
  real(8), allocatable  :: sf_AdvDiffZ_nl(:,:)
  
  real(8), allocatable  :: ss_PotTemp_bl(:,:)
  real(8), allocatable  :: ss_Exner_bl(:,:)
  real(8), allocatable  :: fs_VelX_bl(:,:)
  real(8), allocatable  :: sf_VelZ_bl(:,:)
  real(8), allocatable  :: ss_PotTemp_nl(:,:)
  real(8), allocatable  :: ss_Exner_nl(:,:)
  real(8), allocatable  :: fs_VelX_nl(:,:)
  real(8), allocatable  :: sf_VelZ_nl(:,:)
  real(8), allocatable  :: ss_PotTemp_al(:,:)
  real(8), allocatable  :: ss_Exner_al(:,:)
  real(8), allocatable  :: fs_VelX_al(:,:)
  real(8), allocatable  :: sf_VelZ_al(:,:)

  real(8), allocatable  :: ss_Km_bl(:,:)
  real(8), allocatable  :: ss_Km_nl(:,:)
  real(8), allocatable  :: ss_Km_al(:,:)
  real(8), allocatable  :: ss_Kh_bl(:,:)
  real(8), allocatable  :: ss_Kh_nl(:,:)
  real(8), allocatable  :: ss_Kh_al(:,:)

  real(8)               :: Time
  real(8)               :: DelTime
  real(8)               :: ReStartTime(2)
  integer               :: t, tau


!!!
!!!$B=i4|2=(B
!!!
  !--- I/O $B%U%!%$%kL>$N=i4|2=(B
  call fileset_init

  !--- $B%G%P%0@_Dj$N=i4|2=(B
  call debugset_init(cfgfile)

  !--- $B3J;RE@>pJs$N=i4|2=(B
  call gridset_init(cfgfile)
  
  !--- $B;~4V@QJ,$N@_Dj$N=i4|2=(B
  call timeset_init(cfgfile)
  
  !--- $B<B83L>$N=i4|2=(B
  call nameset_init(cfgfile)
  
  !--- $BJ*M}Dj?t$N=i4|2=(B
  call physset_init(cfgfile)
  
  !--- $B7W;;$N@_Dj$rFI$_9~$`(B
  call arareset_init(cfgfile)

  !--- $B6-3&>r7o$r@_Dj(B
  call bcset_init(cfgfile)

  !--- $B9TNs7W;;%i%$%V%i%j$N=i4|2=(B. 
  !--- $B0z?t$O%9%+%i!<$N(B z $BJ}8~$N7W;;NN0h$NBg$-$5(B
  call linlib_init(cfgfile, RegZMax - RegZMin)
  
  !--- $B=PNO%U%!%$%k$NDj5A(B
  call OpenDataFile

!!!
!!!$BJQ?t$N=i4|2=(B
!!!  
  !--- $BG[Ns$N3d$jEv$F(B, $B=i4|2=(B
  call ArareAlloc

  !--- $B4pK\>l$NCM$r@_Dj(B
  call basicset_file_init( ReStartFile )

  !--- $B>qMp>l$NCM$r@_Dj(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      )

  !--- $B;~9o$N=i4|2=(B
  Time = ReStartTime(2)
  DelTime = 2.0 * DelTimeLong


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


!!!
!!!$B?tCM@QJ,(B
!!!
  !--- $B=i4|>uBV$r%U%!%$%k=PNO(B
  if ( mod( ReStartTime(1), TimeDisp ) == 0 ) then 
     call OutputDataFile( ReStartTime(1), &
          & fs_VelX_bl, sf_VelZ_bl, ss_Exner_bl, &
          & ss_PotTemp_bl, ss_Km_bl, ss_Kh_bl )
  end if

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

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

     !--- $B12G4@-78?t(B, $B123H;678?t$N7W;;(B
     call EddyViscosity(                                     &
          & ss_Km_bl,                                        &
          & DelTime,                                         & 
          & 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,                                &
          & 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)
     
     !--- $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 )
        
        !--- $B%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*$JB.EY$rD9$$;~4V%9%F%C%W$G$N$b$N$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 )

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

     !--- $B%U%!%$%k=PNO(B
     if ( mod( Time, TimeDisp ) == 0 ) then 
        call DbgMessage(fmt="Time = %f", d=(/Time/))
        call OutputDataFile( Time, &
             &  fs_VelX_nl,    sf_VelZ_nl, ss_Exner_nl, &
             &  ss_PotTemp_nl, ss_Km_nl,   ss_Kh_nl )
     end if
     
  end do
  
  !--- $B=PNO%U%!%$%k$N%/%m!<%:(B
  call CloseDataFile

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

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), &
         & 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) )

    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

    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

    call EndSub("ArareAlloc")

  end subroutine ArareAlloc

end program arare

