!=====================================================================
! $BI9>2%b%G%k%a%$%s%k!<%A%s(B
!
! $BMzNr(B
!      2003/08/26: $B=$O@:n@.%3!<%I$+$i0\F0(B
!=====================================================================
program main 
  use prm,     only: IMAX, JMAX, KMAX,dt, dt_T             ! $BDj?t%Q%i%a!<%?(B
  use prmphys, only: s_a
  use com_var, only: get_frame, A_T               
  use Tbound, only: get_Ms_con, get_Ts_con, get_Ts_test1, &
       &            get_TS_test2, get_Ms_test2
  use com_step, only: bound_flg, &                ! $B6-3&CM%U%i%0(B
 &                    Ms_flg,  Ts_flg,  &         ! 
 &                    rs_flg,           &         ! $B%j%9%?!<%H%U%i%0(B 
 &                    monit2D_flg, monit3D_flg,  &! $B%b%K%?%U%i%0(B
 &                    exT_flg,                   &! $B29EY7W;;%U%i%0(B
 &                    com_step_ini,     &         ! $B%b%K%?=PNO@_Dj(B
 &                    com_step_run                ! $B%b%K%?(B on/off
  use file_switch, only : file_switch_on, file_switch_off
  use monit_out,   only : monit2D_out_tmp, monit3D_out_tmp, & 
       &                  monit2D_out,  monit3D_out,& 
 &                    rcNO_2D, rcNO_3D            ! $B%l%3!<%I?t(B

  use u_pack,  only: get_U_ana , get_F_ana , &    ! $BB.EY$K4X$9$k%b%8%e!<%k(B
 &                   get_U,                  &
 &                   get_W, get_Fgrd  
  use height, only : get_H, add_Ms, mk_thick
  use T_pack, only : get_T,   adjust_energy, energy2Temp, &
       &                     get_Qstrain  , homo_T



  
  implicit none  

  real(8)    :: Hsfc1(1:IMAX, 1:JMAX)            ! $BI9>29bEY(B($B;~9o(B: t)
  real(8)    :: Hthk1(1:IMAX, 1:JMAX)            ! $BI98|(B($B;~9o(B: t)
  real(8)    :: Temp1(1:IMAX, 1:JMAX, 1:KMAX)    ! $B29EYJ,I[(B($B;~9o(B: t)
  real(8)    :: energy(1:IMAX, 1:JMAX, 1:KMAX)    ! $B%(%M%k%.!<J,I[(B($B;~9o(B: t)
  real(8)    :: energy_fl(1:IMAX, 1:JMAX, 1:KMAX)    ! $B%(%M%k%.!<J,I[(B($B;~9o(B: t)
!  real(8)    :: energy_k(1:IMAX, 1:JMAX, 1:KMAX) ! $B1?F0%(%M%k%.!<(B($B;~9o(B: t)

  real(8)    :: Uvec1(0:IMAX, 1:JMAX, 1:KMAX)    ! $B?eJ?B.EY(B u($B?eJ?(B, $BH>@0?t(B)
!  real(8)    :: Uvec1_h(0:IMAX, 1:JMAX, 0:KMAX)    ! $B?eJ?B.EY(B u($B?eJ?(B, $BH>@0?t(B)
  real(8)    :: Usfc1(0:IMAX, 1:JMAX)            ! $BI=LLB.EY(B

  real(8)    :: Hsfc2(1:IMAX, 1:JMAX)            ! $BI9>29bEY(B($B;~9o(B: t+1)
  real(8)    :: Hthk2(1:IMAX, 1:JMAX)            ! $BI9>29bEY(B($B;~9o(B: t+1)
  real(8)    :: Temp2(1:IMAX, 1:JMAX, 1:KMAX)    ! $B29EYJ,I[(B($B;~9o(B: t+1)
  real(8)    :: TempH(1:IMAX, 1:JMAX, 1:KMAX)    ! $B29EYJ,I[(B($B;~9o(B: t+1)
  real(8)    :: Tpm(1:IMAX, 1:JMAX, 1:KMAX)    ! $B29EYJ,I[(B($B;~9o(B: t+1)

!  real(8)    :: Uvec2(0:IMAX, 1:JMAX, 1:KMAX) ! $B?eJ?B.EY(B u($B?eJ?(B, $BH>@0?t(B)
  real(8)    :: Wvec1(1:IMAX, 1:JMAX, 0:KMAX)   ! $B1tD>B.EY(B w($B1tD>(B, $BH>@0?t(B)


  real(8)    :: UFL(0:IMAX, 1:JMAX)              ! $B%U%i%C%/%9(B($BH>@0?t(B)$B$N:9(B
  real(8)    :: Uflgrd(0:IMAX, 1:JMAX, 1:KMAX)   ! $B%U%i%C%/%9(B($BH>@0?t(B)$B$N:9(B
  real(8)    :: Qstrain(1:IMAX, 1:JMAX, 1:KMAX)  ! strain heating
!  real(8)    :: sigmaXX(0:IMAX, 1:JMAX, 0:KMAX+1)  ! strain heating
!  real(8)    :: sigmaXZ(0:IMAX, 1:JMAX, 0:KMAX+1)  ! strain heating
 
  
  real(8)    :: Hbnd(1:IMAX, 1:JMAX)             ! $B4pHW9bEY(B($B;~9o(B: t)



  integer                 :: it                 ! do $BJQ?t(B
  integer                 :: loop_end           ! $B%k!<%W2s?t(B
  integer                 :: KT                 ! $B%j%9%?!<%H%U%!%$%k$K(B($B$+$i(B)
                                                ! $B=q$-(B($BFI$_(B)$B9~$`@QJ,2s?t(B

  integer    :: anl                             ! $B2r@O2r$+?tCM2r$+$N%U%i%0(B
  integer    :: gamma_flg                       ! $BCO3LG.N.NL$N(Bon/off
  integer    :: Apara_flg                       ! $B29EY0MB8%Q%i%a!<%?$N(B on/off
  integer                  :: k, i, j



  character(LEN = 100)    :: HEAD               ! $B%j%9%?!<%H%U%!%$%k$N%X%C%@!<(B
  integer                 :: IDATE              ! $B;~9o>pJs(B


  NAMELIST /NAM_TIME/ loop_end  ! loop_end *dt $B$,@QJ,4|4V(B
  NAMELIST /anl_flg / anl, gamma_flg
  NAMELIST /A_flg/ Apara_flg

  gamma_flg = 1
  Apara_flg = 1
  loop_end  = 20000
  anl = 0


  READ(5, NAM_TIME)
  WRITE(6, NAM_TIME)

  READ(5, anl_flg)
  WRITE(6,anl_flg)

  write(6, *)'dt=', dt/s_a, 'dt_T=', dt_T/s_a


!--------------------------------------------------------!
            ! 1. $B=i4|2=(B
!--------------------------------------------------------!
  ! $BF~NO(B, $B=PNO%j%9%?!<%H%U%!%$%k(B, $B%b%K%?%U%!%$%k$N(B open
  ! $B0z?t(B loop_end $B$O=PNO%j%9%?!<%H%U%!%$%kL>$r7h$a$k(B
  call file_switch_on(loop_end)   ! I 
  ! $B%j%9%?!<%H%U%!%$%k$NFI$_9~$_(B
  call rs_read(  &
       HEAD, IDATE, KT, Hsfc1, Hbnd, Temp1) ! O

  ! $B8|$5$rDj5A$9$k(B
  call mk_thick(Hsfc1, Hbnd, Hthk1)

  ! $B;~4V@_Dj(B($B@QJ,3+;O;~4V$H@QJ,2s?t(B)
  call TIME_ini(KT,        & ! I
       &  IDATE)             ! IO
  ! $B%b%K%?=PNO$N@_Dj(B($BJQ?t$NEPO?(B)
  write(6, *)'start monit_ini'
  call monit_ini
  write(6, *)'end monit_ini'
  ! $B%b%K%?(B, $B%j%9%?!<%H%U%!%$%k$N=PNO;~4|$K4X$9$k@_Dj(B
  ! $B6-3&CM$NFI$_9~$_4V3V$N@_Dj(B
  call com_step_ini
  ! $B:BI8$N@_Dj(B
  call get_frame


!--------------------------------------------------------!
            ! 2. $B;~4V@QJ,(B
!--------------------------------------------------------!
  do it = 1, loop_end
!     do j = 1, JMAX                        ! y $BJ}8~$N@QJ,(B($B$=$l$>$lFHN)(B)
     
        ! $B;~9oI=<((B
        call TIME_step(IDATE, it, & !I
             & KT)                  !O          
        ! $B=PNO$*$h$S6-3&CM$NFI$_9~$_$N(B on/off
        call com_step_run(it)       !I         

        !----------------------!
        ! 2.0 $B6-3&CM%G!<%?BeF~(B !
        !----------------------!
        if(Ms_flg == 1)then

!           call get_Ms_test2(it)
           call get_Ms_con
           
        end if
        if(Ts_flg == 1)then

!           call get_Ts_con
!           call get_Ts_test2(it, Hsfc1)
           
           call get_Ts_test1(Hsfc1)
        end if
        !----------------------!
        ! $B%"%$%=%9%?%7!<%b%G%k<B9T(B 
        !----------------------!
          if(bound_flg == 1)then
!            call bound(Hsfc1, Hbnd)! IO
          end if
          
        !----------------------!
        ! 2.1 H $B;;=P(B           !
        !----------------------!
          call homo_T(Hthk1, Temp1, Tpm, TempH)
          call A_T(Apara_flg, Temp1, Tpm)
     if(anl == 1)then
         call get_U_ana(Hsfc1, Hbnd, Hthk1, &     
         & Usfc1, Uvec1)               ! $B;~9o(B t $B$NB.EY(B u $B=PNO(B
         call get_F_ana(Hsfc1, Hbnd, Hthk1,  &  
          & Ufl)                        ! $B;~9o(B t $B$NB.EY(B u $B$NI98|@QJ,(B f
      else
         call get_U(Hsfc1, Hbnd, Hthk1, &     
             & Usfc1, Uvec1)          ! $B;~9o(B t $B$NB.EY(B u $B=PNO(B

         !         call get_U2(Hsfc1, Hbnd, &     
!             & Usfc1, Uvec1,Uvec1_h)   ! $B;~9o(B t $B$NB.EY(B u $B=PNO(B(half)
!         call get_F(Hsfc1, Hbnd, Usfc1, Uvec1, &  
!              & Ufl)                        ! $B;~9o(B t $B$NB.EY(B u $B$NI98|@QJ,(B f
         call get_Fgrd(Hsfc1, Hbnd, Hthk1, Usfc1, Uvec1, &  
              & Ufl, Uflgrd)                 ! $B;~9o(B t $B$NB.EY(B u $B$NI98|@QJ,(B f
         call get_W(Hsfc1, Hbnd, Hthk1, Uvec1, Usfc1, &
              & Wvec1)

      end if
!      call get_kinetic(Hsfc1, Hbnd, Hthk1, Uvec1, Usfc1, Wvec1, energy_k, sigmaXX, sigmaXZ)
      
      ! t+1 $B%9%F%C%WL\$N(B h $B$N7W;;(B($BN.$l$N4sM?(B)
      call get_H(Hthk1,  Ufl, &
           & Hthk2)                 
      !t+1 $B%9%F%C%WL\$N(B h $B$N7W;;(B($B$+$sM\$N4sM?(B)
      call add_Ms(Hthk2) 
      Hsfc2=Hthk2+Hbnd



      !----------------------!
      ! 2.3 T $B;;=P(B           !
      !----------------------!
     ! $B3J;R$ND4@a(B
      call adjust_energy(Temp1, Uflgrd,  Hthk1, Hthk2, &  ! i
           &  energy_fl, energy)                          ! o

      call energy2Temp(energy,  Hthk2,  &                 ! i 
              & Temp1)                                       ! o 

      if(exT_flg == 1)then
      ! $B29EY;;=P(B
         
      call get_Qstrain(Hsfc2, Hthk2, Qstrain)
!      write(6, *)'Qstrain(:, 2, :)', Qstrain(:, 2, :)

      call get_T(Hthk2, Temp1, Qstrain, Temp2, gamma_flg)

      do i =1, IMAX
         do j =1, JMAX
            do k =1, KMAX
               if(Temp2(i, j, k) < 0.0D0)then
                  write(6, *)'temperature is mistake', i, j, k, it
                  write(6, *)'Hsfc1(i, j, :)', Hsfc1(i, j)
                  write(6, *)'Hsfc2(i, j, :)', Hsfc2(i, j)
                  write(6, *)'temp1(i, j, :)', Temp1(i, j, :)
                  write(6, *)'energy(i, j, :)', energy(i, j, :)
                  stop 
               end if
            end do
         end do
      end do
   end if

        !----------------------!
      ! 2.4 $B%b%K%?!<=PNO(B     !
        !----------------------!
        if(monit2D_flg ==1)then
           write(6, *)'time =', KT
           call monit2D_out
!           call monit2D_out_tmp(Hsfc2)
!           call monit2D_out_tmp(UFL(1:IMAX, 1:JMAX))
!           call monit2D_out_tmp(Hbnd)
!           call monit2D_out_tmp(Ufl(1:IMAX, 1:JMAX))
           ! $B$H$j$"$($:(B,$B0z?t$rJQ?tL>$KM?$($F=V4VCM$r=PNO(B
           ! $B$3$NJ}<0$@$HJQ?tL>$@$1(B monit2D_out_tmp $B$r8F$S=P$5$J$$$H%@%a(B
           
!           call monit2D_out
           call monit_clear_2
        end if

        if(exT_flg == 1)then
           if(monit3D_flg ==1)then

           
!           call monit3D_out_tmp(Uvec1(1:IMAX, 1:JMAX, 1:KMAX)*s_a)
!           call monit3D_out_tmp(Wvec1(1:IMAX, 1:JMAX, 1:KMAX)*s_a)
!           call monit3D_out_tmp(Temp2(1:IMAX, 1:JMAX, 1:KMAX))
!           call monit3D_out_tmp(TempH(1:IMAX, 1:JMAX, 1:KMAX))
!           call monit3D_out_tmp(energy_fl(1:IMAX, 1:JMAX, 1:KMAX))
!           call monit3D_out_tmp(Qstrain(1:IMAX, 1:JMAX, 1:KMAX))
!           call monit3D_out_tmp(energy_k(1:IMAX, 1:JMAX, 1:KMAX))
!           call monit3D_out_tmp(sigmaXX(1:IMAX, 1:JMAX, 1:KMAX))
!           call monit3D_out_tmp(sigmaXZ(1:IMAX, 1:JMAX, 1:KMAX))

              call monit3D_out
              call monit_clear_3
        end if
     end if
        
        !----------------------------!
        !2.5 $B%j%9%?!<%H%U%!%$%k:n@.(B  !
        !----------------------------!
        if(rs_flg ==1 )then
           write(6, *)'rs_flg = 1'
           call rs_write(HEAD, IDATE, KT, Hsfc2, Hbnd, Temp2)
        end if

        !----------------------!
        ! $B8=:_$NJQ?t$X%3%T!<(B   !
        !----------------------!
        call reset2D(Hsfc2, Hsfc1)
        call reset2D(Hthk2, Hthk1)
        if(exT_flg == 1)then
           call reset3D(Temp2, Temp1)
        end if

     !   end do
     end do



!--------------------------------------------------------!
          ! $B=hM}$N=*N;(B
!--------------------------------------------------------!

!  call close_monit(30)
!  call close_monit(40)

end program main

