!=================================================
! 2D cumulus model - kaminari
!
! Author  : TAKAHASHI Koko
! Date    : 2003/11/28  $B:G=*99?7(B
!           2003/09/01  $B?75,:n@.(B
! Note    : - $BC10L$O(B MKS $B7O(B
!           - $BL5<!8505NO(B = $B%(%/%9%J!<4X?t(B
!           - 2$B<!85(B (x-z)
!           - $B6-3&MQ$NG[Ns(B i = -1,-2,im,im+1 k = -1,-2,km,km+1
!           - $BB&LL6-3&$O<~4|6-3&>r7o(B
!           - $B>qMp>l$N7W;;$K6-3&MQ$NG[Ns$rM?$($?$?$a(B,
!             $B4pK\>l$K$b6-3&MQ$NG[Ns$rM?$($?$,;HMQ$7$F$$$J$$(B
!             ($B4pK\>l$H>qMp>l$NG[Ns$N?t$r9g$o$;$k$?$a(B)
!           - $B=PNO%G!<%?7A<0(B: NetCDF $B7A<0(B
! 
!=================================================
program kaminari_2d

  !-------------------------------------------------
  ! $BJQ?t!&Dj?t$N@_Dj(B
  !-------------------------------------------------

  use gtool_history

  implicit none

  ! $B6u4V2rA|EY@_Dj(B
  integer(8), parameter :: im = 10      ! $B3J;RE@$N@_Dj(B (x)
  integer(8), parameter :: km = 10      ! $B3J;RE@$N@_Dj(B (z)

  ! $B:BI8JQ?t(B
  real(8), parameter :: xmin = 0.0d0    ! $B%b%G%kNN0h$N:8B&LL(B
  real(8), parameter :: xmax = 1.0d4    ! $B%b%G%kNN0h$N:8B&LL(B
  real(8), parameter :: dx = (xmax-xmin)/im ! x $B$N3J;R4V3V(B
  real(8), parameter :: zmin = 0.0d0    ! $B%b%G%kNN0h$NDl(B
  real(8), parameter :: zmax = 1.0d4    ! $B%b%G%kNN0h$N9b$5(B
  real(8), parameter :: dz = (zmax-zmin)/km ! z^{*} $B$N3J;R4V3V(B

  ! $B3J;R(B
  real(8) :: x(-2:im+1,-2:km+1)         ! $BCO7A$K1h$C$?:BI87O$G$N(B x $B:BI8(B
  real(8) :: z(-2:im+1,-2:km+1)         ! $BCO7A$K1h$C$?:BI87O$G$N(B z $B:BI8(B

  ! $B;~4VJQ?t(B
  integer(8), parameter :: nts = 100     ! $BC;$$;~4V%9%F%C%W?t(B
!  integer(8), parameter :: ntb = 5      ! $BD9$$;~4V%9%F%C%W?t(B
  real(8), parameter :: dts = 1.0d-2    ! $BC;$$;~4V4V3V(B d tau
!  real(8), parameter :: dtb = nts*dts   ! $BD9$$;~4V4V3V(B dt
  real(8), parameter :: ndisp = 1.0d0   ! $B=PNO4V3V(B

  integer(8) :: i, k                    ! do $BJQ?t(B/$B6u4V(B
  integer(8) :: its                     ! do $BJQ?t(B/$B;~4V(B
!  integer(8) :: itb                    ! do $BJQ?t(B/$B;~4V(B

  !--- $BJQ?t(B
  !--- $B4pK\>l$NJQ?t(B
  real(8) :: pi_bs(-2:im+1,-2:km+1)     ! $BL5<!8505NO(B
  real(8) :: prss_bs(-2:im+1,-2:km+1)   ! $B05NO(B
  real(8) :: dens_bs(-2:im+1,-2:km+1)   ! $BL)EY(B
  real(8) :: temp_bs(-2:im+1,-2:km+1)   ! $B29EY(B
  real(8) :: ptemp_bs(-2:im+1,-2:km+1)  ! $B290L(B
  real(8) :: vptemp_bs(-2:im+1,-2:km+1) ! $B2>290L(B
  real(8) :: qv_bs(-2:im+1,-2:km+1)     ! $B?e>x5$$N:.9gHf(B

  !--- $B>qMp$NJQ?t(B
  real(8) :: u(-2:im+1,-2:km+1)         ! x $BJ}8~B.EY(B
  real(8) :: w(-2:im+1,-2:km+1)         ! z $BJ}8~B.EY(B
  real(8) :: omg(-2:im+1,-2:km+1)       ! z^{*} $BJ}8~B.EY(B
  real(8) :: pi(-2:im+1,-2:km+1)        ! $BL5<!8505NO(B
  real(8) :: prss(-2:im+1,-2:km+1)      ! $B05NO(B
  real(8) :: temp(-2:im+1,-2:km+1)      ! $B29EY(B
  real(8) :: ptemp(-2:im+1,-2:km+1)     ! $B290L(B
!  real(8) :: qv(-2:im+1,-2:km+1)        ! $B?e>x5$$N:.9gHf(B
!  real(8) :: qc(-2:im+1,-2:km+1)        ! $B1@?e$N:.9gHf(B
!  real(8) :: qr(-2:im+1,-2:km+1)        ! $B1+?e$N:.9gHf(B
!  real(8) :: qvs(-2:im+1,-2:km+1)       ! $B?e>x5$$NK0OB:.9gHf(B
!  real(8) :: km_sub(-2:im+1,-2:km+1)    ! $B12:.9g78?t(B K_{m}
  !--- t - dt $B$G$NCM(B
  real(8) :: u_old(-2:im+1,-2:km+1)     ! x $BJ}8~B.EY(B
!  real(8) :: w_old(-2:im+1,-2:km+1)     ! z $BJ}8~B.EY(B
!  real(8) :: omg_old(-2:im+1,-2:km+1)     ! z $BJ}8~B.EY(B
  real(8) :: pi_old(-2:im+1,-2:km+1)    ! $BL5<!8505NO(B
!  real(8) :: ptemp_old(-2:im+1,-2:km+1)          ! $B290L(B
!  real(8) :: qv_old(-2:im+1,-2:km+1)    ! $B?e>x5$$N:.9gHf(B
!  real(8) :: qc_old(-2:im+1,-2:km+1)    ! $B1@?e$N:.9gHf(B
!!  real(8) :: qr_old(-2:im+1,-2:km+1)    ! $B1+?e$N:.9gHf(B
!  real(8) :: qvs_old(-2:im+1,-2:km+1)   ! $B?e>x5$$NK0OB:.9gHf(B
!  real(8) :: km_sub_old(-2:im+1,-2:km+1) ! $B%5%V%0%j%C%I1?F0%(%M%k%.!<(B
  !--- t - 2*dt $B$G$NCM(B
!  real(8) :: ptemp_old2(-2:im+1,-2:km+1) ! $B290L(B
!  real(8) :: qv_old2(-2:im+1,-2:km+1)   ! $B?e>x5$$N:.9gHf(B
!  real(8) :: qc_old2(-2:im+1,-2:km+1)   ! $B1@?e$N:.9gHf(B
!  real(8) :: qr_old2(-2:im+1,-2:km+1)   ! $B1+?e$N:.9gHf(B
!  real(8) :: qvs_old2(-2:im+1,-2:km+1)  ! $B?e>x5$$NK0OB:.9gHf(B
!  real(8) :: km_sub_old2(-2:im+1,-2:km+1) ! $B%5%V%0%j%C%I1?F0%(%M%k%.!<(B

  real(8) :: u_adv(-2:im+1,-2:km+1)     ! x $BJ}8~B.EY(B
  real(8) :: u_cori(-2:im+1,-2:km+1)    ! x $BJ}8~B.EY(B
  real(8) :: u_trb(-2:im+1,-2:km+1)     ! x $BJ}8~B.EY(B
  real(8) :: w_adv(-2:im+1,-2:km+1)     ! x $BJ}8~B.EY(B
  real(8) :: w_buoy(-2:im+1,-2:km+1)    ! x $BJ}8~B.EY(B
  real(8) :: w_trb(-2:im+1,-2:km+1)     ! x $BJ}8~B.EY(B

  !--- $B2;GH8:?j9`(B
  real(8) :: div(-2:im+1,-2:km+1)
  real(8) :: div_x_u(-2:im+1,-2:km+1)   ! x $BJ}8~HyJ,(B
  real(8) :: div_z_u(-2:im+1,-2:km+1)   ! z^{*} $BJ}8~HyJ,(B
  real(8) :: div_z_w(-2:im+1,-2:km+1)   ! z^{*} $BJ}8~HyJ,(B

  !--- $B$=$NB>$NJQ?t(B         
  real(8) :: g_sqrt(-2:im+1,-2:km+1)    ! $B%d%3%S%"%s$N(B1/2$B>h(B
  real(8) :: g_sqrt_x(-2:im+1,-2:km+1)  ! g_sqrt$B$r(Bx$BJ}8~$KJ?6Q2=(B
  real(8) :: g_sqrt_z(-2:im+1,-2:km+1)  ! g_sqrt$B$r(Bz^{*}$BJ}8~$KJ?6Q2=(B
  real(8) :: g_sqrt_xz(-2:im+1,-2:km+1) ! g_sqrt$B$r(Bz^{*}$BJ}8~$KJ?6Q2=(B
  real(8) :: g13(-2:im+1,-2:km+1)       ! $B%d%3%S%"%s$N(B13$B@.J,(B
  real(8) :: g13x(-2:im+1,-2:km+1)      ! g13$B$r(Bx$BJ}8~$KJ?6Q2=(B
  real(8) :: g13z(-2:im+1,-2:km+1)      ! g13$B$r(Bz^{*}$BJ}8~$KJ?6Q2=(B
  real(8) :: g13xz(-2:im+1,-2:km+1)     ! g13$B$r(Bxz^{*}$BJ}8~$KJ?6Q2=(B

  !--- $BJ*M}%Q%i%a!<%?(B
  ! $B%5%V%k!<%A%sFb$@$1$G;H$&Dj?t$b$3$3$GDj5A(B
  real(8), parameter :: prss_sfc = 1.0d5 ! $BCOI=LL$G$N05NO(B [Pa]
  real(8), parameter :: rd = 2.87d2     ! $B4%Ag6u5$$N5$BNDj?t(B
  real(8), parameter :: cp = 1.952d3    ! $BDj05HfG.(B
  real(8), parameter :: cv = 1.463d3    ! $BDj@QHfG.(B
  real(8), parameter :: grv = 9.8d0     ! $B=ENO(B
  real(8), parameter :: alp = 5.0d-4*dz/dts ! $BDj?t(B
!  real(8), parameter :: eps = 6.22d-1   ! $BDj?t(B
!  real(8), parameter :: lv = 5.8d2      ! $B>xH/$N@xG.(B $B4X?t$K$9$k(B?
!  real(8), parameter :: cm = 2d-1       ! $BHfNcDj?t(B
  real(8), parameter :: beta = 1.0d0     ! $B1tD>1"2rK!$N78?t(B
!  real(8), parameter :: nuh = 1.0d-5*dx*dx ! $B?eJ?J}8~$N?tCMG4@-9`$N78?t(B
!  real(8), parameter :: nuv = 1.0d-5*dz*dz ! $B1tD>J}8~$N?tCMG4@-9`$N78?t(B
  ! $B:.9g5wN%(B: parameter $B$@$H$Y$->h7W;;$9$k$H%(%i!<$K$J$k(B
!  real(8)            :: ml              ! $B:.9g5wN%(B
!  real(8),parameter :: ml = sqrt(dx*dz) ! $B:.9g5wN%(B
!  ml = sqrt(dx*dz)

  !-------------------------------------------------
  ! $BCO7A4XO"$N@_Dj(B
  !-------------------------------------------------

 do k = 0, km-1
     do i = 0, im-1
        !--- $B%d%3%S%"%s(B
        !--- 1 $B9TL\(B: $BF~NO(B, 2,3 $B9TL\(B: $B=PNO(B
        call jacobian(i,k,im,km,                            &
             &        g_sqrt,g_sqrt_x,g_sqrt_z,g_sqrt_xz,   &
             &        g13,g13x,g13z,g13xz)
     end do
  end do

  g_sqrt(:,:) = 1.0d0
  g_sqrt_x(:,:) = 1.0d0
  g_sqrt_z(:,:) = 1.0d0
  g_sqrt_xz(:,:) = 1.0d0

  !-------------------------------------------------
  ! x-z^{*} $B:BI8$N@_Dj(B
  !-------------------------------------------------
  !--- $B86E@(B (x(1,1),z(1,1))
  do k = -2, km+1
     do i = -2, im+1
     x(i,k) = dx*i
     z(i,k) = dz*k
     end do
  end do

  call output_gtool4_init

  !-------------------------------------------------
  ! $B4pK\>l(B
  !-------------------------------------------------
  !--- $BA4BN$N29EY$H?e>x5$$N:.9gHf$rM?$($k(B
  do k = 0, km-1
     do i = 0, im-1
     !--- $B29EY(B
     temp_bs(i,k) = 3.0d2! - 6.5d-3*z(i,k)

     !--- $B?e>x5$$N:.9gHf(B
     qv_bs(i,k) = 0.0d0
     end do
  end do

  !--- $BCOI=LL05NO$rM?$($k(B
  do i =0,im-1
     prss_bs(i,0) = prss_sfc
  end do

  !--- $BCOI=LL$G$N3FCM$N7W;;(B
  do i = 0, im-1
     !--- $BL5<!8505NO(B
     pi_bs(i,0) = (prss_bs(i,0)/prss_sfc)**(rd/cp)

     !--- $B290L(B
     ptemp_bs(i,0) = temp_bs(i,0)/pi_bs(i,0)

     !--- $B2>290L(B
     vptemp_bs(i,0) = ptemp_bs(i,0)*(1.0d0 + 6.1d-2*qv_bs(i,0))

     !--- $BL)EY(B
     dens_bs(i,0) = prss_sfc/rd*pi_bs(i,0)**(cv/rd)/vptemp_bs(i,0)
  end do

  !--- $BCOI=LL$h$j>e$G$N3FCM(B
  k = 1
  do i = 0, im-1
     !--- $B@E?e05J?9U$N<0(B -> $BL5<!8505NO(B
     pi_bs(i,k) = pi_bs(i,k-1)                           &
          &       - dz*g_sqrt(i,k-1)*grv/(cp*vptemp_bs(i,k-1))

     !--- $B05NO(B
     prss_bs(i,k) = prss_sfc*pi_bs(i,k)**(cp/rd)

     !--- $B290L(B
     ptemp_bs(i,k) = temp_bs(i,k)/pi_bs(i,k)

     !--- $B2>290L(B
     vptemp_bs(i,k) = ptemp_bs(i,k)*(1.0d0 + 6.1d-2*qv_bs(i,k))

     !--- $BL)EY(B
     dens_bs(i,k) = prss_sfc/rd*pi_bs(i,k)**(cv/cp)/vptemp_bs(i,k)
  end do

  !--- $BCOI=LL$h$j>e$G$N3FCM(B
  do k = 2, km-1
     do i = 0, im-1
        !--- $B@E?e05J?9U$N<0(B -> $BL5<!8505NO(B
        pi_bs(i,k) = pi_bs(i,k-2)                           &
             &       - 2.0d0*dz*g_sqrt(i,k-1)               &
             &         *grv/(cp*vptemp_bs(i,k-1))

        !--- $B05NO(B
        prss_bs(i,k) = prss_sfc*pi_bs(i,k)**(cp/rd)

        !--- $B290L(B
        ptemp_bs(i,k) = temp_bs(i,k)/pi_bs(i,k)

        !--- $B2>290L(B
        vptemp_bs(i,k) = ptemp_bs(i,k)*(1.0d0 + 6.1d-2*qv_bs(i,k))

        !--- $BL)EY(B
        dens_bs(i,k) = prss_sfc/rd*pi_bs(i,k)**(cv/rd)/vptemp_bs(i,k)
     end do
  end do

  !-------------------------------------------------
  ! $B4pK\>l$N6-3&>r7o(B
  !---  i = -1 : $BB&LL6-3&(B
  do k = -2,km+1
     qv_bs(-1,k) = 0.0d0
     prss_bs(-1,k) = prss_bs(im-2,k)
     pi_bs(-1,k) = pi_bs(im-2,k)
     temp_bs(-1,k) = temp_bs(im-2,k)
     ptemp_bs(-1,k) = ptemp_bs(im-2,k)
     vptemp_bs(-1,k) = vptemp_bs(km-2,k)
     dens_bs(-1,k) = dens_bs(im-2,k)
  end do

  !---  i = -2 : $BB&LL6-3&(B
  do k = -2,km+1
     qv_bs(-2,k) = 0.0d0
     prss_bs(-2,k) = prss_bs(im-3,k)
     pi_bs(-2,k) = pi_bs(im-3,k)
     temp_bs(-2,k) = temp_bs(im-3,k)
     ptemp_bs(-2,k) = ptemp_bs(im-3,k)
     vptemp_bs(-2,k) = vptemp_bs(im-3,k)
     dens_bs(-2,k) = dens_bs(im-3,k)
  end do

  !---  i = im-1 : $BB&LL6-3&(B
  do k = -2,km+1
     qv_bs(im-1,k) = 0.0d0
     prss_bs(im-1,k) = prss_bs(1,k)
     pi_bs(im-1,k) = pi_bs(1,k)
     temp_bs(im-1,k) = temp_bs(1,k)
     ptemp_bs(im-1,k) = ptemp_bs(1,k)
     vptemp_bs(im-1,k) = vptemp_bs(1,k)
     dens_bs(im-1,k) = dens_bs(1,k)
  end do

  !---  i = im : $BB&LL6-3&(B
  do k = -2,km+1
     qv_bs(im,k) = 0.0d0
     prss_bs(im,k) = prss_bs(1,k)
     pi_bs(im,k) = pi_bs(1,k)
     temp_bs(im,k) = temp_bs(1,k)
     ptemp_bs(im,k) = ptemp_bs(1,k)
     vptemp_bs(im,k) = vptemp_bs(1,k)
     dens_bs(im,k) = dens_bs(1,k)
  end do

  !---  i = im+1 : $BB&LL6-3&(B
  do k = -2,km+1
     qv_bs(im+1,k) = 0.0d0
     prss_bs(im+1,k) = prss_bs(2,k)
     pi_bs(im+1,k) = pi_bs(2,k)
     temp_bs(im+1,k) = temp_bs(2,k)
     ptemp_bs(im+1,k) = ptemp_bs(2,k)
     vptemp_bs(im+1,k) = vptemp_bs(2,k)
     dens_bs(im+1,k) = dens_bs(2,k)
  end do

  !---  k = -1 : $B2<It6-3&(B
  do i = 0,im-1
     qv_bs(i,-1) = 0.0d0
     prss_bs(i,-1) = prss_bs(i,0)
     pi_bs(i,-1) = pi_bs(i,0)
     temp_bs(i,-1) = temp_bs(i,0)
     ptemp_bs(i,-1) = ptemp_bs(i,0)
     vptemp_bs(i,-1) = vptemp_bs(i,0)
     dens_bs(i,-1) = dens_bs(i,0)
  end do

  !---  k = -2 : $B2<It6-3&(B
  do i = 0,im-1
     qv_bs(i,-2) = 0.0d0
     prss_bs(i,-2) = prss_bs(i,0)
     pi_bs(i,-2) = pi_bs(i,0)
     temp_bs(i,-2) = temp_bs(i,0)
     ptemp_bs(i,-2) = ptemp_bs(i,0)
     vptemp_bs(i,-2) = vptemp_bs(i,0)
     dens_bs(i,-2) = dens_bs(i,0)
  end do

  !---  k = km : $B>eIt6-3&(B
  do i = 0,im-1
     qv_bs(i,km) = 0.0d0
     prss_bs(i,km) = prss_bs(i,km-1)
     pi_bs(i,km) = pi_bs(i,km-1)
     temp_bs(i,km) = temp_bs(i,km-1)
     ptemp_bs(i,km) = ptemp_bs(i,km-1)
     vptemp_bs(i,km) = vptemp_bs(i,km-1)
     dens_bs(i,km) = dens_bs(i,km-1)
  end do

  !---  k = km+1 : $B>eIt6-3&(B
  do i = 0,im-1
     qv_bs(i,km+1) = 0.0d0
     prss_bs(i,km+1) = prss_bs(i,km-1)
     pi_bs(i,km+1) = pi_bs(i,km-1)
     temp_bs(i,km+1) = temp_bs(i,km-1)
     ptemp_bs(i,km+1) = ptemp_bs(i,km-1)
     vptemp_bs(i,km+1) = vptemp_bs(i,km-1)
     dens_bs(i,km+1) = dens_bs(i,km-1)
  end do

  !--- $B=PNO(B
  call output_base_gtool4

  !-------------------------------------------------
  ! $B>qMp>l$N7W;;(B
  !------------------------------------------------- 
  !-------------------------------------------------
  ! $B=i4|CM(B
  its = 0   ! itb = 0
  do k = 0,km-1
     do i = 0,im-1
        !--- $BB.EY(B, $B05NO(B, $B29EY$rM?$($k(B
        !--- 1$B9TL\F~NOJQ?t(B, 2$B9TL\=PNOJQ?t(B
        call init_value(i,k,im,km,prss_bs,temp_bs,   &
             &         u,w,prss,temp)
     end do
  end do

  !--- $BM?$($?=i4|CM$+$iB>$N=i4|CM$r7W;;(B
  do k = 0,km-1
     do i = 0,im-1
        !--- $BCO7A$K1h$C$?:BI87O$G$N1tD>B.EY(B omega
        !--- 1$B9TL\F~NOJQ?t(B, 2$B9TL\=PNOJQ?t(B
        call omgw(i,k,im,km,g_sqrt_z,g13x,u,w,              &
             &    omg)

        !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $B=PNO(B)
        call prss2pi(0,i,k,im,km,prss_sfc,rd,cp,prss,pi)

        !--- $B290L(B
        ptemp(i,k) = temp(i,k)/pi(i,k)
     end do
  end do

  !-------------------------------------------------
  ! $B6-3&>r7o(B
     !---  i = -1 : $BB&LL6-3&(B
     do k = -2,km+1
        u(-1,k) = u(im-2,k)
        w(-1,k) = w(im-2,k)
        omg(-1,k) = omg(im-2,k)
        pi(-1,k) = pi(im-2,k)
        ptemp(-1,k) = ptemp(im-2,k)
     end do
   
     !---  i = -2 : $BB&LL6-3&(B
     do k = -2,km+1
        u(-2,k) = u(im-3,k)
        w(-2,k) = w(im-3,k)
        pi(-2,k) = pi(im-3,k)
        ptemp(-2,k) = ptemp(im-3,k)
     end do
   
     !---  i = im-1 : $BB&LL6-3&(B
     do k = -2,km+1
        u(im-1,k) = u(0,k)
        w(im-1,k) = w(0,k)
        omg(im-1,k) = omg(0,k)
        pi(im-1,k) = pi(0,k)
        ptemp(im-1,k) = ptemp(0,k)
     end do
   
     !---  i = im : $BB&LL6-3&(B
     do k = -2,km+1
        u(im,k) = u(1,k)
        w(im,k) = w(1,k)
        omg(im,k) = omg(1,k)
        pi(im,k) = pi(1,k)
        ptemp(im,k) = ptemp(1,k)
     end do
   
     !---  i = im+1 : $BB&LL6-3&(B
     do k = -2,km+1
        u(im+1,k) = u(2,k)
        w(im+1,k) = w(2,k)
        pi(im+1,k) = pi(2,k)
        ptemp(im+1,k) = ptemp(2,k)
     end do

     !---  k = -1 : $B2<It6-3&(B
     do i = 0,im-1
        u(i,-1) = u(i,0)
        w(i,-1) = 0.0d0
        pi(i,-1) = pi(i,0)
        ptemp(i,-1) = ptemp(i,0)
     end do
   
     !---  k = -2 : $B2<It6-3&(B
     do i = 0,im-1
        u(i,-2) = u(i,0)
        w(i,-2) = 0.0d0
        pi(i,-2) = pi(i,0)
        ptemp(i,-2) = ptemp(i,0)
     end do
   
     !---  k = km : $B>eIt6-3&(B
     do i = 0,im-1
        u(i,km) = u(i,km-1)
        w(i,km) = 0.0d0
        pi(i,km) = pi(i,km-1)
        ptemp(i,km) = ptemp(i,km-1)
     end do
   
     !---  k = km+1 : $B>eIt6-3&(B
     do i = 0,im-1
        u(i,km+1) = u(i,km-1)
        w(i,km+1) = 0.0d0
        pi(i,km+1) = pi(i,km-1)
        ptemp(i,km+1) = ptemp(i,km-1)
     end do
   
  !--- $B=PNO(B
  call output_gtool4

  !---------------------------
  !--- $B1?F0J}Dx<0Cf$ND9$$;~4V4V3V$G7W;;$9$kItJ,(B
  do k = 0,km-2
     do i = 0,im-2
        u_adv(i,k) = 0.0d0
        u_cori(i,k) = 0.0d0
        u_trb(i,k) = 0.0d0

        w_adv(i,k) = 0.0d0
        w_buoy(i,k) = 0.0d0
        w_trb(i,k) = 0.0d0
     end do
  end do

  do k = -2,km
     do i = -2,im
        u_old(i,k) = u(i,k)
        pi_old(i,k) = pi(i,k)
     end do
  end do

  !-------------------------------------------------
  !--- $BC;$$;~4V4V3V(B($BD9$$;~4V4V3V$N(B1$B%9%F%C%WL\$^$G(B)
  do its = 0,nts-1

     !--- $B2;GH8:?j9`(B
     !--- 1$B9TL\F~NOJQ?t(B, 2$B9TL\=PNOJQ?t(B
     do k = -1,km
        do i = -1,im
           call div_diff(i,k,im,km,dx,dz,u,w,g_sqrt,g13xz, &
                &        div)
        end do
     end do

     do k = 0,km-1
        do i = 0,im-1
           call div_diff_avr(i,k,im,km,dx,dz,div,           &
                &            div_x_u,div_z_u,div_z_w)
        end do
     end do

     do k = 0,km-1
        do i = 0,im-2
           !--- $BCO7A$K1h$C$?:BI87O$G$N1tD>B.EY(B omega
           !--- 1$B9TL\F~NOJQ?t(B, 2$B9TL\=PNOJQ?t(B
           call omgw(i,k,im,km,g_sqrt_z,g13x,u,w,           &
                &    omg)

           !--- $BB.EY(B u
           !--- 1,2,3$B9TL\F~NOJQ?t(B, 4$B9TL\F~=PNOJQ?t(B
           call vel_u(i,k,im,km,cp,alp,dts,g13z,div_x_u,    &
                &     div_z_u,dx,dz,vptemp_bs,pi,u_adv,     &
                &     u_cori,u_trb,                         &
                &     u)
        end do
     end do

     do i = 0,im-2
        !--- $BL5<!8505NO(B pi (SSL2$B;HMQ(B)
        !--- 1,2,3,4$B9TL\F~NOJQ?t(B, $B!V(Bk$B!W$H(B5$B9TL\F~=PNOJQ?t(B
        call non_dim_pi(i,k,its,im,km,cp,cv,rd,alp,beta,dx, &
             &          dz,dts,pi_bs,vptemp_bs,dens_bs,u,w, &
             &          w_adv,w_buoy,w_trb,g_sqrt,g13xz,    &
             &          div_z_w,                            &
             &          pi)
     end do

     do k = -1,km-1
        do i = -1,im-1
           call div_diff(i,k,im,km,dx,dz,u_old,w,g_sqrt,g13xz, &
                &        div)
        end do
     end do

     !--- $B05NO(B ($BL5<!8505NO(B pi -> $B05NO(B prss $B=PNO(B)
     do k = 0,km-1
        do i = 0,im-2
           call prss2pi(1,i,km,im,km,prss_sfc,rd,cp,prss,pi)
        end do
     end do

     do k = 0,km-1
        do i = 0,im-2
           !--- $BB.EY(B w
          !--- 1,2,3$B9TL\F~NOJQ?t(B, 4$B9TL\F~=PNOJQ?t(B
           call vel_w(i,k,im,km,dts,cp,alp,beta,dz,vptemp_bs,pi, &
                &     pi_old,g_sqrt_z,div_z_w,w_adv,w_buoy,w_trb,&
                &     w)
        end do
     end do

     !-------------------------------------------------
     ! $B6-3&>r7o(B
     !---  i = -1 : $BB&LL6-3&(B
     do k = -2,km+1
        u(-1,k) = u(im-2,k)
        w(-1,k) = w(im-2,k)
        omg(-1,k) = omg(im-2,k)
        pi(-1,k) = pi(im-2,k)
        ptemp(-1,k) = ptemp(im-2,k)
     end do
   
     !---  i = -2 : $BB&LL6-3&(B
     do k = -2,km+1
        u(-2,k) = u(im-3,k)
        w(-2,k) = w(im-3,k)
        pi(-2,k) = pi(im-3,k)
        ptemp(-2,k) = ptemp(im-3,k)
     end do
   
     !---  i = im-1 : $BB&LL6-3&(B
     do k = -2,km+1
        u(im-1,k) = u(0,k)
        w(im-1,k) = w(0,k)
        omg(im-1,k) = omg(0,k)
        pi(im-1,k) = pi(0,k)
        ptemp(im-1,k) = ptemp(0,k)
     end do
   
     !---  i = im : $BB&LL6-3&(B
     do k = -2,km+1
        u(im,k) = u(1,k)
        w(im,k) = w(1,k)
        omg(im,k) = omg(1,k)
        pi(im,k) = pi(1,k)
        ptemp(im,k) = ptemp(1,k)
     end do
   
     !---  i = im+1 : $BB&LL6-3&(B
     do k = -2,km+1
        u(im+1,k) = u(2,k)
        w(im+1,k) = w(2,k)
        pi(im+1,k) = pi(2,k)
        ptemp(im+1,k) = ptemp(2,k)
     end do

     !---  k = -1 : $B2<It6-3&(B
     do i = 0,im-1
        u(i,-1) = u(i,0)
        w(i,-1) = 0.0d0
        pi(i,-1) = pi(i,0)
        ptemp(i,-1) = ptemp(i,0)
     end do
   
     !---  k = -2 : $B2<It6-3&(B
     do i = 0,im-1
        u(i,-2) = u(i,0)
        w(i,-2) = 0.0d0
        pi(i,-2) = pi(i,0)
        ptemp(i,-2) = ptemp(i,0)
     end do
   
     !---  k = km : $B>eIt6-3&(B
     do i = 0,im-1
        u(i,km) = u(i,km-1)
        w(i,km) = 0.0d0
        pi(i,km) = pi(i,km-1)
        ptemp(i,km) = ptemp(i,km-1)
     end do
   
     !---  k = km+1 : $B>eIt6-3&(B
     do i = 0,im-1
        u(i,km+1) = u(i,km-1)
        w(i,km+1) = 0.0d0
        pi(i,km+1) = pi(i,km-1)
        ptemp(i,km+1) = ptemp(i,km-1)
     end do
     !--- $B=PNO(B
     call output_gtool4


     !--- $B<!(B(t + dt)$B$N7W;;;~$K;HMQ(B: t $B$NCM(B

     do k = -2,km+1
        do i = -2,im+1
           u_old(i,k) = u(i,k)
           pi_old(i,k) = pi(i,k)
        end do
     end do

  end do

  call output_gtool4_close

contains

!-----------------------------------------------------------
!--- $B%5%V%k!<%A%s(B

  !--- gtool4 $B=PNO4XO"(B
  subroutine output_gtool4_init
    ! $B%R%9%H%j!<:n@.(B
    call HistoryCreate(                                     &
          & file='kaminari_2d.nc', title='2D cumulus model',&
          & source='kaminari - 2D cumulus model',           &
          & institution='GFD_Dennou Club deepconv project', &
          & dims=(/'x','z','t'/),                           &
          & dimsizes=(/int(im+4,4),int(km+4,4),0/),         &
          & longnames=(/'X-coordinate','Z-coordinate',      &
          &             '    Time    '/),                   &
          & units=(/'m','m','s'/), origin=0.0,              &
          & interval=real(ndisp))

    ! $BJQ?t=PNO(B
    call HistoryPut('x',x(-2:im+1,0))
    call HistoryPut('z',z(0,-2:km+1))

    ! $BJQ?tDj5A(B
    ! $B4pK\>l$N05NO(B
    call HistoryAddVariable(                            &
         & varname='prss_bs', dims=(/'x','z'/),             &
         & longname='pressure', units='Pa', xtype='double')
    ! $B4pK\>l$NL5<!8505NO(B
    call HistoryAddVariable(                            &
         & varname='pi_bs', dims=(/'x','z'/),             &
         & longname='nondimensional pressure', units='1',&
         & xtype='double')
    ! $B4pK\>l$NL)EY(B
    call HistoryAddVariable(                            &
         & varname='dens_bs', dims=(/'x','z'/),             &
         & longname='density', units='kg/m3',           &
         & xtype='double')
    ! $B4pK\>l$N?e>x5$:.9gHf(B
    call HistoryAddVariable(                            &
         & varname='qv_bs', dims=(/'x','z'/),             &
         & longname='mixing ratio of water vapor',      &
         &  units='g/g', xtype='double')
    ! $B4pK\>l$N29EY(B
    call HistoryAddVariable(                            &
         & varname='temp_bs', dims=(/'x','z'/),             &
         & longname='temperature', units='K', xtype='double')
    ! $B4pK\>l$N290L(B
    call HistoryAddVariable(                            &
         & varname='ptemp_bs', dims=(/'x','z'/),             &
         & longname='potential temperature (base)',         &
         & units='K', xtype='double')

    ! $BB.EY(B u
    call HistoryAddVariable(                                &
         & varname='u', dims=(/'x','z','t'/),           &
         & longname='u-velocity', units='m/s', xtype='double')
    ! $BB.EY(B w
    call HistoryAddVariable(                                &
         & varname='w', dims=(/'x','z','t'/),           &
         & longname='w-velocity', units='m/s', xtype='double')
    ! $B05NO(B
    call HistoryAddVariable(                                &
         & varname='prss', dims=(/'x','z','t'/),          &
         & longname='pressure', units='Pa', xtype='double')
    ! $BL5<!8505NO(B
    call HistoryAddVariable(                                &
         & varname='pi', dims=(/'x','z','t'/),          &
         & longname='nondimensional pressure', units='1', &
         & xtype='double')
    ! $B29EY(B
    call HistoryAddVariable(                                &
         & varname='temp', dims=(/'x','z','t'/),          &
         & longname='temperature', units='K', xtype='double')
    ! $B290L(B
    call HistoryAddVariable(                                &
         & varname='ptemp', dims=(/'x','z','t'/),          &
         & longname='potential temperature', units='K',     &
         & xtype='double')
    ! $B?e>x5$:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qv', dims=(/'x','z','t'/),          &
!         & longname='mixing ratio of water vapor',          &
!         &  units='g/g', xtype='double')
    ! $B1@?e:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qc', dims=(/'x','z','t'/),          &
!        & longname='mixing ratio of cloud drops',          &
!         & units='g/g', xtype='double')
    ! $B1+?e:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qr', dims=(/'x','z','t'/),          &
!         & longname='mixing ratio of rainwater',            &
!         & units='g/g', xtype='double')
    ! $BK0OB?e>x5$:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qvs', dims=(/'x','z','t'/),          &
!         & longname='saturation mixing ratio',          &
!         &  units='g/g', xtype='double')
    ! $B12:.9g78?t(B
!    call HistoryAddVariable(                                &
!         & varname='km_sub', dims=(/'x','z','t'/),         &
!         & longname='momentum eddy mixing coefficient',        &
!         & units='1', xtype='double')
    ! $B%5%V%0%j%C%I1?F0%(%M%k%.!<(B
!    call HistoryAddVariable(                                &
!         & varname='e_sub', dims=(/'x','z','t'/),         &
!         & longname='kinetic energy',        &
!         & units='1', xtype='double')
    ! $B%^%C%W%U%!%/%?!<(B
!    call HistoryAddVariable(                                &
!         & varname='g_sqrt', dims=(/'x','z','t'/),           &
!         & longname='G_SQRT', units='1', xtype='double')
!    call HistoryAddVariable(                                &
!         & varname='g_sqrt_z', dims=(/'x','z','t'/),           &
!         & longname='G_SQRT z', units='1', xtype='double')
!    call HistoryAddVariable(                                &
!         & varname='g13', dims=(/'x','z','t'/),           &
!         & longname='G13', units='1', xtype='double')
!    call HistoryAddVariable(                                &
!         & varname='g13x', dims=(/'x','z','t'/),           &
!         & longname='G13 x', units='1', xtype='double')
!    call HistoryAddVariable(                                &
!         & varname='g13z', dims=(/'x','z','t'/),           &
!         & longname='G13 z', units='1', xtype='double')
!    call HistoryAddVariable(                                &
 !         & varname='g13xz', dims=(/'x','z','t'/),           &
!         & longname='G13 xz', units='1', xtype='double')
  end subroutine output_gtool4_init

  subroutine output_base_gtool4
    call HistoryPut('prss_bs', prss_bs)
    call HistoryPut('pi_bs', pi_bs)
    call HistoryPut('temp_bs', temp_bs)
    call HistoryPut('ptemp_bs', ptemp_bs)
    call HistoryPut('qv_bs', qv_bs)
    call HistoryPut('dens_bs', dens_bs)
  end subroutine output_base_gtool4

  subroutine output_gtool4
    write(*,*) 'its =', its
    call HistoryPut('u', u)
    call HistoryPut('w', w)
    call HistoryPut('prss', prss)
    call HistoryPut('pi', pi)
    call HistoryPut('temp', temp)
    call HistoryPut('ptemp', ptemp)
!    call HistoryPut('qv', qv)
!    call HistoryPut('qc', qc)
!    call HistoryPut('qr', qr)
!    call HistoryPut('qvs', qvs)
!    call HistoryPut('km_sub', km_sub)
!    call HistoryPut('e_sub', e_sub)
!    call HistoryPut('g_sqrt', g_sqrt)
!    call HistoryPut('g_sqrt_z', g_sqrt_z)
!    call HistoryPut('g13', g13)
!    call HistoryPut('g13x', g13x)
!    call HistoryPut('g13z', g13z)
!    call HistoryPut('g13xz', g13xz)
  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program kaminari_2d
