!=================================================
! 2D cumulus model - kaminari
!
! Author  : TAKAHASHI Koko
! Date    : 2003/11/20  $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
!             - **(0:im+1,0), **(0:im+1,km+1)
!               **(0,0:km+1), **(im+1,0:km+1)
!           - $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)
!           - $BIA2h$9$k$H$-$O6-3&ItJ,$r:o=|$9$k$3$H(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-1) ! 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-1)  ! z^{*} $B$N3J;R4V3V(B

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

  ! $B;~4VJQ?t(B
  integer, parameter :: nts = 5           ! $BC;$$;~4V%9%F%C%W?t(B
!  integer, parameter :: ntb = 5         ! $BD9$$;~4V%9%F%C%W?t(B
  real(8), parameter :: dts = 1.0d0     ! $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(0:im+1,0:km+1)       ! $BL5<!8505NO(B
  real(8) :: prss_bs(0:im+1,0:km+1)     ! $B05NO(B
  real(8) :: dens_bs(0:im+1,0:km+1)     ! $BL)EY(B
  real(8) :: temp_bs(0:im+1,0:km+1)     ! $B29EY(B
  real(8) :: ptemp_bs(0:im+1,0:km+1)    ! $B290L(B
  real(8) :: vptemp_bs(0:im+1,0:km+1)    ! $B2>290L(B
  real(8) :: qv_bs(0:im+1,0:km+1)       ! $B?e>x5$$N:.9gHf(B

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

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

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

  !--- $B$=$NB>$NJQ?t(B         
  real(8) :: g_sqrt(0:im+1,0:km+1)         ! $B%d%3%S%"%s$N(B1/2$B>h(B
  real(8) :: g_sqrt_z(0:im+1,0:km+1)        ! g_sqrt$B$r(Bz^{*}$BJ}8~$KJ?6Q2=(B
  real(8) :: g13(0:im+1,0:km+1)         ! $B%d%3%S%"%s$N(B13$B@.J,(B
  real(8) :: g13x(0:im+1,0:km+1)        ! g13$B$r(Bx$BJ}8~$KJ?6Q2=(B
  real(8) :: g13z(0:im+1,0:km+1)        ! g13$B$r(Bz^{*}$BJ}8~$KJ?6Q2=(B
  real(8) :: g13xz(0:im+1,0: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-6*dz/dts ! $BDj?t(B
  real(8), parameter :: eps = 6.22d-1  ! $BDj?t(B
  real(8), parameter :: lv = 580      ! $B>xH/$N@xG.(B $B4X?t$K$9$k(B?
  real(8), parameter :: cm = 0.2d0     ! $BHfNcDj?t(B
  real(8), parameter :: ml = (dx*dz)**(1/2) ! $B:.9g5wN%(B
  real(8), parameter :: nuh = 1.0d-5 ! $B?eJ?J}8~$N?tCMG4@-9`$N78?t(B
  real(8), parameter :: nuv = 1.0d-5 ! $B1tD>J}8~$N?tCMG4@-9`$N78?t(B
  real(8), parameter :: beta = 1.0d0    ! $B1tD>1"2rK!$N78?t(B

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

 do k = 1, km
     do i = 1, im
        !--- $B%d%3%S%"%s(B
        call jacobian(i,k,im,km,g_sqrt,g13,g_sqrt_z,g13x,g13z,g13xz)
     end do
  end do
  g_sqrt(:,:) = 1.0d0
  g_sqrt_z(:,:) = 1.0d0

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

  call output_gtool4_init

  !-------------------------------------------------
  ! $B4pK\>l(B
  !-------------------------------------------------
  !--- $BA4BN$N29EY$H?e>x5$$N:.9gHf$rM?$($k(B
  do k = 1, km
     do i = 1, im
     !--- $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
  prss_bs(1:im,1) = prss_sfc

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

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

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

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

  !--- $BCOI=LL$h$j>e$G$N3FCM(B
  k = 2
  do i = 1, im
     !--- $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 = 3, km
     do i = 1, im
        !--- $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

  ! $B6-3&>r7o(B
  !--- $B2<It6-3&(B  : *(0:im+1,0), $B>eIt6-3&(B  : *(0:im+1,km+1)
  !--- $B:8B&LL6-3&(B: *(0,1:km),   $B1&B&LL6-3&(B: *(im+1,1:km)
  !--- $B$"$H$G4JC12=$G$-$?$i$d$k(B

  !--- $B3Q(B
  prss_bs(0,0) = prss_bs(1,1)
  temp_bs(0,0) = temp_bs(1,1)
  qv_bs(0,0) = qv_bs(1,1)
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  pi_bs(0,0) = (prss_bs(0,0)/prss_sfc)**(rd/cp)
  !--- $B290L(B
  ptemp_bs(0,0) = temp_bs(0,0)/pi_bs(0,0)
  !--- $B2>290L(B
  vptemp_bs(0,0) = ptemp_bs(0,0)*(1.0d0 + 6.1d-2*qv_bs(0,0))
  !--- $BL)EY(B
  dens_bs(0,0) = prss_sfc/rd*pi_bs(0,0)**(cv/rd)            &
       &         /vptemp_bs(0,0)

  prss_bs(im+1,0) = prss_bs(im,1)
  temp_bs(im+1,0) = temp_bs(im,1)
  qv_bs(im+1,0) = qv_bs(im,1)
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  pi_bs(im+1,0) = (prss_bs(im+1,0)/prss_sfc)**(rd/cp)
  !--- $B290L(B
  ptemp_bs(im+1,0) = temp_bs(im+1,0)/pi_bs(im+1,0)
  !--- $B2>290L(B
  vptemp_bs(im+1,0) = ptemp_bs(im+1,0)                      &
       &              *(1.0d0 + 6.1d-2*qv_bs(im+1,0))
  !--- $BL)EY(B
  dens_bs(im+1,0) = prss_sfc/rd*pi_bs(im+1,0)**(cv/rd)      &
       &            /vptemp_bs(im+1,0)

  prss_bs(0,km+1) = prss_bs(1,km)
  temp_bs(0,km+1) = temp_bs(1,km)
  qv_bs(0,km+1) = qv_bs(1,km)
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  pi_bs(0,km+1) = (prss_bs(0,km+1)/prss_sfc)**(rd/cp)
  !--- $B290L(B
  ptemp_bs(0,km+1) = temp_bs(0,km+1)/pi_bs(0,km+1)
  !--- $B2>290L(B
  vptemp_bs(0,km+1) = ptemp_bs(0,km+1)                      &
       &              *(1.0d0 + 6.1d-2*qv_bs(0,km+1))
  !--- $BL)EY(B
  dens_bs(0,km+1) = prss_sfc/rd*pi_bs(0,km+1)**(cv/rd)      &
       &            /vptemp_bs(0,km+1)

  prss_bs(im+1,km+1) = prss_bs(im,km)
  temp_bs(im+1,km+1) = temp_bs(im,km)
  qv_bs(im+1,km+1) = qv_bs(im,km)
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  pi_bs(im+1,km+1) = (prss_bs(im+1,km+1)/prss_sfc)**(rd/cp)
  !--- $B290L(B
  ptemp_bs(im+1,km+1) = temp_bs(im+1,km+1)/pi_bs(im+1,km+1)
  !--- $B2>290L(B
  vptemp_bs(im+1,km+1) = ptemp_bs(im+1,km+1)                &
       &                 *(1.0d0 + 6.1d-2*qv_bs(im+1,km+1))
  !--- $BL)EY(B
  dens_bs(im+1,km+1) = prss_sfc/rd*pi_bs(im+1,km+1)**(cv/rd)&
       &               /vptemp_bs(im+1,km+1)

  !---  k = 0 : $B2<It6-3&(B
  prss_bs(1:im,0) = prss_bs(1:im,1)
  temp_bs(1:im,0) = temp_bs(1:im,1)
  qv_bs(1:im,0) = qv_bs(1:im,1)
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  pi_bs(1:im,0) = (prss_bs(1:im,0)/prss_sfc)**(rd/cp)
  !--- $B290L(B
  ptemp_bs(1:im,0) = temp_bs(1:im,0)/pi_bs(1:im,0)
  !--- $B2>290L(B
  vptemp_bs(1:im,0) = ptemp_bs(1:im,0)                         &
       &           *(1.0d0 + 6.1d-2*qv_bs(1:im,0))
  !--- $BL)EY(B
  dens_bs(1:im,0) = prss_sfc/rd*pi_bs(1:im,0)**(cv/rd)         &
       &         /vptemp_bs(1:im,0)

  !---  k = km+1 : $B>eIt6-3&(B
  prss_bs(1:im,km+1) = prss_bs(1:im,km)
  temp_bs(1:im,km+1) = temp_bs(1:im,km)
  qv_bs(1:im,km+1) = qv_bs(1:im,km)
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  pi_bs(1:im,km+1) = (prss_bs(1:im,km+1)/prss_sfc)**(rd/cp)
  !--- $B290L(B
  ptemp_bs(1:im,km+1) = temp_bs(1:im,km+1)/pi_bs(1:im,km+1)
  !--- $B2>290L(B
  vptemp_bs(1:im,km+1) = ptemp_bs(1:im,km+1)                   &
       &              *(1.0d0 + 6.1d-2*qv_bs(1:im,km+1))
  !--- $BL)EY(B
  dens_bs(1:im,km+1) = prss_sfc/rd*pi_bs(1:im,km+1)**(cv/rd)   &
       &            /vptemp_bs(1:im,km+1)

  !---  i = 0 : $BB&LL6-3&(B
  prss_bs(0,1:km) = prss_bs(1,1:km)
  temp_bs(0,1:km) = temp_bs(1,1:km)
  qv_bs(0,1:km) = qv_bs(1,1:km)

  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  pi_bs(0,1:km) = (prss_bs(0,1:km)/prss_sfc)**(rd/cp)
  !--- $B290L(B
  ptemp_bs(0,1:km) = temp_bs(0,1:km)/pi_bs(0,1:km)
  !--- $B2>290L(B
  vptemp_bs(0,1:km) = ptemp_bs(0,1:km)                         &
       &           *(1.0d0 + 6.1d-2*qv_bs(0,1:km))
  !--- $BL)EY(B
  dens_bs(0,1:km) = prss_sfc/rd*pi_bs(0,1:km)**(cv/rd)         &
       &         /vptemp_bs(0,1:km)

  !---  i = im+1 : $BB&LL6-3&(B
  prss_bs(im+1,1:km) = prss_bs(im,1:km)
  temp_bs(im+1,1:km) = temp_bs(im,1:km)
  qv_bs(im+1,1:km) = qv_bs(im,1:km)
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  pi_bs(im+1,1:km) = (prss_bs(im+1,1:km)/prss_sfc)**(rd/cp)
  !--- $B290L(B
  ptemp_bs(im+1,1:km) = temp_bs(im+1,1:km)/pi_bs(im+1,1:km)
  !--- $B2>290L(B
  vptemp_bs(im+1,1:km) = ptemp_bs(im+1,1:km)                   &
       &              *(1.0d0 + 6.1d-2*qv_bs(im+1,1:km))
  !--- $BL)EY(B
  dens_bs(im+1,1:km) = prss_sfc/rd*pi_bs(im+1,1:km)**(cv/rd)   &
       &            /vptemp_bs(im+1,1:km)

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

  !-------------------------------------------------
  ! $B>qMp>l$N7W;;(B
  !------------------------------------------------- 
  !-------------------------------------------------
  ! $B6-3&>r7o(B
  !--- $B2<It6-3&(B  : *(0:im+1,0), $B>eIt6-3&(B  : *(0:im+1,km+1)
  !--- $B:8B&LL6-3&(B: *(0,1:km),   $B1&B&LL6-3&(B: *(im+1,1:km)

  !--- $B3Q(B
  u(0,0) = 0.0d0
  w(0,0) = 0.0d0
  omg(0,0) = 0.0d0
  prss(0,0) = prss_bs(1,1)
  temp(0,0) = temp_bs(1,1)
  qv(0,0) = qv_bs(1,1)
  qc(0,0) = 0.0d0
  qr(0,0) = 0.0d0
  qvs(0,0) = 0.0d0
  e_sub(0,0) = 0.0d0
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  call prss2pi(0,0,0,im,km,prss_sfc,rd,cp,prss,pi)
  !--- $B290L(B
  ptemp(0,0) = temp(0,0)/pi(0,0)

  u(im+1,0) = 0.0d0
  w(im+1,0) = 0.0d0
  omg(im+1,0) = 0.0d0
  prss(im+1,0) = prss_bs(im,1)
  temp(im+1,0) = temp_bs(im,1)
  qv(im+1,0) = qv_bs(im,1)
  qc(im+1,0) = 0.0d0
  qr(im+1,0) = 0.0d0
  qvs(im+1,0) = 0.0d0
  e_sub(im+1,0) = 0.0d0
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  call prss2pi(0,im+1,0,im,km,prss_sfc,rd,cp,prss,pi)
  !--- $B290L(B
  ptemp(im+1,0) = temp(im+1,0)/pi(im+1,0)

  u(0,km+1) = 0.0d0
  w(0,km+1) = 0.0d0
  omg(0,km+1) = 0.0d0
  prss(0,km+1) = prss_bs(1,km)
  temp(0,km+1) = temp_bs(1,km)
  qv(0,km+1) = qv_bs(1,km)
  qc(0,km+1) = 0.0d0
  qr(0,km+1) = 0.0d0
  qvs(0,km+1) = 0.0d0
  e_sub(0,km+1) = 0.0d0
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  call prss2pi(0,0,km+1,im,km,prss_sfc,rd,cp,prss,pi)
  !--- $B290L(B
  ptemp(0,km+1) = temp(0,km+1)/pi(0,km+1)

  u(im+1,km+1) = 0.0d0
  w(im+1,km+1) = 0.0d0
  omg(im+1,km+1) = 0.0d0
  prss(im+1,km+1) = prss_bs(im,km)
  temp(im+1,km+1) = temp_bs(im,km)
  qv(im+1,km+1) = qv_bs(im,km)
  qc(im+1,km+1) = 0.0d0
  qr(im+1,km+1) = 0.0d0
  qvs(im+1,km+1) = 0.0d0
  e_sub(im+1,km+1) = 0.0d0
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
  call prss2pi(0,im+1,km+1,im,km,prss_sfc,rd,cp,prss,pi)
  !--- $B290L(B
  ptemp(im+1,km+1) = temp(im+1,km+1)/pi(im+1,km+1)

  !---  k = 0 : $B2<It6-3&(B
  u(1:im,0) = 0.0d0
  w(1:im,0) = 0.0d0
  omg(1:im,0) = 0.0d0
  prss(1:im,0) = prss_bs(1:im,1)
  temp(1:im,0) = temp_bs(1:im,1)
  qv(1:im,0) = qv_bs(1:im,1)
  qc(1:im,0) = 0.0d0
  qr(1:im,0) = 0.0d0
  qvs(1:im,0) = 0.0d0
  e_sub(1:im,0) = 0.0d0
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  do i = 1,im
     !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
     call prss2pi(0,i,0,im,km,prss_sfc,rd,cp,prss,pi)
  end do
  !--- $B290L(B
  ptemp(1:im,0) = temp(1:im,0)/pi(1:im,0)

  !---  k = km+1 : $B>eIt6-3&(B
  u(1:im,km+1) = 0.0d0
  w(1:im,km+1) = 0.0d0
  omg(1:im,km+1) = 0.0d0
  prss(1:im,km+1) = prss_bs(1:im,km)
  temp(1:im,km+1) = temp_bs(1:im,km)
  qv(1:im,km+1) = qv_bs(1:im,km)
  qc(1:im,km+1) = 0.0d0
  qr(1:im,km+1) = 0.0d0
  qvs(1:im,km+1) = 0.0d0
  e_sub(1:im,km+1) = 0.0d0
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  do i = 1,im
     !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
     call prss2pi(0,i,km+1,im,km,prss_sfc,rd,cp,prss,pi)
  end do
  !--- $B290L(B
  ptemp(1:im,km+1) = temp(1:im,km+1)/pi(1:im,km+1)

  !---  i = 0 : $BB&LL6-3&(B
  u(0,1:km) = 0.0d0
  w(0,1:km) = 0.0d0
  omg(0,1:km) = 0.0d0
  prss(0,1:km) = prss_bs(1,1:km)
  temp(0,1:km) = temp_bs(1,1:km)
  qv(0,1:km) = qv_bs(1,1:km)
  qc(0,1:km) = 0.0d0
  qr(0,1:km) = 0.0d0
  e_sub(0,1:km) = 0.0d0
  !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  do k = 1,km
     !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
     call prss2pi(0,0,k,im,km,prss_sfc,rd,cp,prss,pi)
  end do
  !--- $B290L(B
  ptemp(0,1:km) = temp(0,1:km)/pi(0,1:km)

  !---  i = im+1 : $BB&LL6-3&(B
  u(im+1,1:km) = 0.0d0
  w(im+1,1:km) = 0.0d0
  omg(im+1,1:km) = 0.0d0
  prss(im+1,1:km) = prss_bs(im,1:km)
  temp(im+1,1:km) = temp_bs(im,1:km)
  qv(im+1,1:km) = qv_bs(im,1:km)
  qc(im+1,1:km) = 0.0d0
  qr(im+1,1:km) = 0.0d0
  qvs(im+1,1:km) = 0.0d0
  e_sub(im+1,1:km) = 0.0d0
    !--- $BM?$($?6-3&>r7o$+$iB>$N6-3&$G$NCM$r7W;;(B
  do k = 1,km
     !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(B)
     call prss2pi(0,im+1,k,im,km,prss_sfc,rd,cp,prss,pi)
  end do
  !--- $B290L(B
  ptemp(im+1,1:km) = temp(im+1,1:km)/pi(im+1,1:km)

  !-------------------------------------------------
  ! $B=i4|CM(B
  its = 0!, itb = 0
  do k = 1,km
     do i = 1,im
        !--- $B05NO(B, $B29EY(B, $B?e>x5$:.9gHf(B, $B1@?e:.9gHf(B, $B1+?e:.9gHf(B
        !--- $B%5%V%0%j%C%I1?F0%(%M%k%.!<$rM?$($k(B
        call init_value(i,k,im,km,u,w,prss_bs,temp_bs,  &
             &          prss,temp,qv,qc,qr,e_sub)
     end do
  end do

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

        !--- $BL5<!8505NO(B ($B05NO(B prss -> $BL5<!8505NO(B pi $BJQ49(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)

        !--- $BK0OB:.9gHf(B
        call mxrt_satu(i,k,im,km,prss_bs,pi_bs,ptemp,qvs)
     end do
  end do

  !--- $B1?F0J}Dx<0Cf$ND9$$;~4V4V3V$G7W;;$9$kItJ,(B
  do k = 1,km
     do i = 1,im
        call vel_u_tb(i,k,im,km,nuh,nuv,dx,dz,u,omg,        &
             &        u_adv,u_cori,u_trb)

        call vel_w_tb(i,k,im,km,nuh,nuv,grv,dx,dz,ptemp_bs, &
             &        qv_bs,u,w,omg,ptemp,qv,qc,qr,w_adv,   &
             &        w_buoy,w_trb)
     end do
  end do

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

  !--- t - dt $B$NCM$K$9$k(B
  do k = 1,km
     do i = 1,im
        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
     call div_diff(i,k,im,km,dx,dz,u,w,g_sqrt,g13xz,        &
          &        div_x_u,div_z_u,div_z_w)

     do k = 1,km
        do i = 1,im
           call omgw(i,k,im,km,g_sqrt_z,g13x,u,w,omg)

           !--- $BB.EY(B u
           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 = 1,im
        !--- $BL5<!8505NO(B pi (SSL2$B;HMQ(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)

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

     do k = 1,km
        do i = 1,im
           !--- $BB.EY(B w
           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

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

     !--- $B<!(B(t + dt)$B$N7W;;;~$K;HMQ(B: t $B$NCM(B
     do k = 1,km
        do i = 1,im
           pi_old(i,k) = pi(i,k)
        end do
     end do

  end do

  !-------------------------------------------------
  !--- $BD9$$;~4V4V3V$N(B 1$B%9%F%C%WL\(B
!  itb = 1
  do k = 1,km
     do i = 1,im       
        call omgw(i,k,im,km,g_sqrt_z,g13x,u,w,omg)

        !--- $B290L(B
        call ptemp_1stp(i,k,im,km,nuh,nuv,cp,ml,dx,dz,dtb,  &
             &          u,prss_bs,pi_bs,dens_bs,omg,qv,qr,  &
             &          qvs,ptemp)

        !--- $B%5%V%0%j%C%I1?F0%(%M%k%.!<(B
        call kener_sub_1stp(i,k,im,km,nuh,nuv,eps,lv,rd,cp, &
             &              cm,ml,grv,dtb,dx,dz,ptemp_bs,u, &
             &              w,omg,temp,ptemp,qv,qc,g_sqrt,g_sqrt_z,&
             &              g13x,g13z,g13xz,e_sub,km_sub)

        !--- $B12:.9g78?t(B
        km_sub(i,k) = cm*ml*e_sub(i,k)**(1/2)

        !--- $B?e>x5$$N:.9gHf(B
        qv(i,k) = 0.00d0

        !--- $B1@?e$N:.9gHf(B
        qc(i,k) = 0.00d0

        !--- $B1+?e$N:.9gHf(B
        qr(i,k) = 0.00d0

     end do
  end do

!  do k = 1,km
!     do i = 1,im
!        ptemp_old2(i,k) = ptemp_old(i,k)
!        ptemp_old(i,k) = ptemp(i,k)
!        e_sub_old2(i,k) = e_sub_old(i,k)
!        e_sub_old(i,k) = e_sub(i,k)
!     end do
!  end do

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

     !--- $B2;GH8:?j9`(B
     call div_diff(i,k,im,km,dx,dz,u,w,g_sqrt,g13xz,        &
          &        div_x_u,div_z_u,div_z_w)

     do k = 1,km
        do i = 1,im
           call omgw(i,k,im,km,g_sqrt_z,g13x,u,w,omg)

           !--- $BB.EY(B u
           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 = 1,im
        !--- $BL5<!8505NO(B pi (SSL2$B;HMQ(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)

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

     do k = 1,km
        do i = 1,im
           !--- $BB.EY(B w
           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

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

     !--- $B<!(B(t + dt)$B$N7W;;;~$K;HMQ(B: t $B$NCM(B
     do k = 1,km
        do i = 1,im
           pi_old(i,k) = pi(i,k)
        end do
     end do

  end do

  !-------------------------------------------------
  !--- $BD9$$;~4V4V3V$N(B 2$B%9%F%C%W0J9_(B
!  do itb = 2,ntb

!     do k = 1,km
!        do i = 1,im

        !--- $B290L(B: $B=q$-49$(I,MW(B
!        call ptemp(i,k,im,km,nuh,nuv,cp,ml,dx,dz,dtb,       &
!             &          u,prss_bs,pi_bs,dens_bs,omg,qv,qr,  &
!             &          qvs,ptemp)

        !--- $B%5%V%0%j%C%I1?F0%(%M%k%.!<(B: $B=q$-49$(I,MW(B
!        call kener_sub(i,k,im,km,nuh,nuv,eps,lv,rd,cp,cm,  &
!             &        ml,grv,dtb,dx,dz,ptemp_bs,u,w,omg,   &
!             &        temp,ptemp,qv,qc,g_sqrt,g_sqrt_z,g13x,g13z, &
!             &        g13xz,e_sub,km_sub)

        !--- $B?e>x5$$N:.9gHf(B
!        qv(i,k) = 0.00d0

        !--- $B1@?e$N:.9gHf(B
!        qc(i,k) = 0.00d0

        !--- $B1+?e$N:.9gHf(B
!        qr(i,k) = 0.00d0

!        end do
!     end do

     !------------------------
     !--- $BC;$$;~4V4V3V(B($BD9$$;~4V4V3V$N(B 2$B%9%F%C%W0J9_(B)
!     do its = 0,nts-1

!     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+2,4),int(km+2,4),0/),         &
          & longnames=(/'X-coordinate','Z-coordinate',      &
          &             '    Time    '/),                   &
          & units=(/'m','m','s'/), origin=0.0,              &
          & interval=real(ndisp*dts))

    ! $BJQ?t=PNO(B
    call HistoryPut('x',x(0:im+1,1))
    call HistoryPut('z',z(1,0: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')
    ! $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')
    ! $B12:.9g78?t(B
    call HistoryAddVariable(                                &
         & varname='km_sub', dims=(/'x','z','t'/),         &
         & longname='momentum eddy mixing coefficient',        &
         & 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('dens_bs', dens_bs)
    call HistoryPut('qv_bs', qv_bs)
    call HistoryPut('temp_bs', temp_bs)
    call HistoryPut('ptemp_bs', ptemp_bs)
  end subroutine output_base_gtool4

  subroutine output_gtool4
    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('e_sub', e_sub)
    call HistoryPut('km_sub', km_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
