!=================================================
! 2D cumulus model - kaminari
!
! Author  : TAKAHASHI Koko
! Date    : 2003/11/11
!           2003/11/10
!           2003/11/05
!           2003/11/04
!           2003/11/01
!           2003/10/31
!           2003/09/01
! 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, parameter :: im = 5          ! $B3J;RE@$N@_Dj(B (x)
  integer, parameter :: km = 5          ! $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 = 10000.0D0 ! $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 = 10000.0D0 ! $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 :: i, k                       ! do $BJQ?t(B/$B6u4V(B
  integer :: its                        ! do $BJQ?t(B/$B;~4V(B
!  integer :: 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) :: w_old(0:im+1,0:km+1)      ! $B1tD>B.EY(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

  !--- SSL2 $B9TNs7W;;$G$NJQ?t(B
  real(8) :: sbd(0:im+1,0:km+1)
  real(8) :: d(0:im+1,0:km+1)
  real(8) :: spd(0:im+1,0:km+1)
  real(8) :: b(0:im+1,0:km+1)

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

  !--- $B$=$NB>$NJQ?t(B         
  real(8) :: g12(0:im+1,0:km+1)         ! $B%d%3%S%"%s$N(B1/2$B>h(B
  real(8) :: g12z(0:im+1,0:km+1)        ! g12$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

  real(8) :: c_bs(0:im+1,0:km+1)        ! $B4pK\>l$G$N2;B.(B
  real(8) :: ssl_fctr(0:im+1,0:km+1)    ! $B9TNs@.J,$N78?t(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 = 100000.0D0  ! $BCOI=LL$G$N05NO(B [Pa]
  real(8), parameter :: rd = 287.0D0    ! $B4%Ag6u5$$N5$BNDj?t(B
  real(8), parameter :: cp = 1952.0D0   ! $BDj05HfG.(B
  real(8), parameter :: cv = 1463.0D0   ! $BDj@QHfG.(B
  real(8), parameter :: grv = 9.8D0     ! $B=ENO(B
  real(8), parameter :: alp = 0.50D0*dz/dts ! $BDj?t(B
  real(8), parameter :: eps = 6.22e-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 :: nuh = 0.00001D0 ! $B?eJ?J}8~$N?tCMG4@-9`$N78?t(B
  real(8), parameter :: nuv = 0.00001D0 ! $B1tD>J}8~$N?tCMG4@-9`$N78?t(B
  real(8), parameter :: ml = (dx*dz)**(1/2) ! $B:.9g5wN%(B
  real(8), parameter :: beta = 0.5D0    ! Crank-Nicolson $BK!$N78?t(B

!  real(8) :: f(1:im)                   ! $B%3%j%*%j%Q%i%a!<%?(B
!  real(8), parameter :: d2r =  3.141592654e0/180.e0
                                        ! $B!VEY!W$+$i!V%i%8%"%s!WJQ49(B
!  real(8), parameter :: omega = 7.292e-5   ! $BCO5e$N3QB.EY(B
!  real(8) :: sinlat                     ! sin($B0^EY(B)
!  real(8) :: lat(0:im)                  ! $B0^EY(B

  ! SSL2 $BMQ%Q%i%a!<%?(B
  integer       :: is
  integer       :: ip
  integer       :: vw
  integer       :: icon

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

  do k = 1, km
     do i = 1, im
        !--- $B%d%3%S%"%s(B
        call jacobian(i,k,im,km,g12,g13,g12z,g13x,g13z,g13xz)
     end do
  end do
  
  !-------------------------------------------------
  ! 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) = 300.0D0 - 0.0065D0*z(i,k)

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

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

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

     !--- $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 + 0.61D0*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
  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*g12(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 + 0.61D0*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*g12(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 + 0.61D0*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
  end do

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

  do i =0.im+1
     call boundary()
  end do

  !--- $B3Q(B
   u(0,0) = 0.0D0
   w(0,0) = 0.0D0
   omg(0,0) = 0.0D0
   prss(0,0) = prss_bs(1,1)
   ptemp(0,0) = ptemp_bs(1,1)
   qv(0,0) = qv_bs(1,1)
   qc(0,0) = 0.0D0
   qr(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)
   !--- $BK0OB:.9gHf(B
   call mxrt_satu(0,0,im,km,prss_bs,pi_bs,ptemp,qvs)

   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)
   ptemp(0,km+1) = ptemp_bs(1,km)
   qv(0,km+1) = qv_bs(1,km)
   qc(0,km+1) = 0.0D0
   qr(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)
   !--- $BK0OB:.9gHf(B
   call mxrt_satu(0,km+1,im,km,prss_bs,pi_bs,ptemp,qvs)

   u(im+1,0) = 0.0D0
   w(im+1,0) = 0.0D0
   omg(im+1,0) = 0.0D0
   prss(im+1,0) = prss_bs(1,km)
   ptemp(im+1,0) = ptemp_bs(1,km)
   qv(im+1,0) = qv_bs(1,km)
   qc(im+1,0) = 0.0D0
   qr(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)
   !--- $BK0OB:.9gHf(B
   call mxrt_satu(im+1,0,im,km,prss_bs,pi_bs,ptemp,qvs)

   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)
   ptemp(im+1,km+1) = ptemp_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
   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)
   !--- $BK0OB:.9gHf(B
   call mxrt_satu(im+1,km+1,im,km,prss_bs,pi_bs,ptemp,qvs)

  !---  i = 0 : $BB&LL6-3&(B
  do k = 1,km
     u(0,k) = 0.0D0
     w(0,k) = 0.0D0
     omg(0,k) = 0.0D0
     prss(0,k) = prss_bs(1,k)
     ptemp(0,k) = ptemp_bs(1,k)
     qv(0,k) = qv_bs(1,k)
     qc(0,k) = 0.0D0
     qr(0,k) = 0.0D0
     e_sub(0,k) = 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,k,im,km,prss_sfc,rd,cp,prss,pi)
     !--- $B290L(B
     ptemp(0,k) = temp(0,k)/pi(0,k)
     !--- $BK0OB:.9gHf(B
     call mxrt_satu(0,k,im,km,prss_bs,pi_bs,ptemp,qvs)
  end do

  !---  i = im+1 : $BB&LL6-3&(B
  do k = 1,km
     u(im+1,k) = 0.0D0
     w(im+1,k) = 0.0D0
     omg(im+1,k) = 0.0D0
     prss(im+1,k) = prss_bs(im,k)
     ptemp(im+1,k) = ptemp_bs(im,k)
     qv(im+1,k) = qv_bs(im,k)
     qc(im+1,k) = 0.0D0
     qr(im+1,k) = 0.0D0
     e_sub(im+1,k) = 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,k,im,km,prss_sfc,rd,cp,prss,pi)
     !--- $B290L(B
     ptemp(im+1,k) = temp(im+1,k)/pi(im+1,k)
     !--- $BK0OB:.9gHf(B
     call mxrt_satu(im+1,k,im,km,prss_bs,pi_bs,ptemp,qvs)
  end do

  !---  k = 0 : $B2<It6-3&(B
  do i = 1,im
     u(i,0) = 0.0D0
     w(i,0) = 0.0D0
     omg(i,0) = 0.0D0
     prss(i,0) = prss_bs(i,1)
     ptemp(i,0) = ptemp_bs(i,1)
     qv(i,0) = qv_bs(i,1)
     qc(i,0) = 0.0D0
     qr(i,0) = 0.0D0
     e_sub(i,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,i,0,im,km,prss_sfc,rd,cp,prss,pi)
     !--- $B290L(B
     ptemp(i,0) = temp(i,0)/pi(i,0)
     !--- $BK0OB:.9gHf(B
     call mxrt_satu(i,0,im,km,prss_bs,pi_bs,ptemp,qvs)
  end do

  !---  k = km+1 : $B>eIt6-3&(B
  do i = 1,im
     u(i,km+1) = 0.0D0
     w(i,km+1) = 0.0D0
     omg(i,km+1) = 0.0D0
     prss(i,km+1) = prss_bs(i,km)
     ptemp(i,km+1) = ptemp_bs(i,km)
     qv(i,km+1) = qv_bs(i,km)
     qc(i,km+1) = 0.0D0
     qr(i,km+1) = 0.0D0
     e_sub(i,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,i,km+1,im,km,prss_sfc,rd,cp,prss,pi)
     !--- $B290L(B
     ptemp(i,km+1) = temp(i,km+1)/pi(i,km+1)
     !--- $BK0OB:.9gHf(B
     call mxrt_satu(i,km+1,im,km,prss_bs,pi_bs,ptemp,qvs)
  end do

  !-------------------------------------------------
  ! $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,g12z,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

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

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

  !-------------------------------------------------
  !--- $BC;$$;~4V4V3V(B($BD9$$;~4V4V3V$N(B1$B%9%F%C%WL\$^$G(B)

  do its = 1,nts-1

     do k = 1,km
        do i = 1,im
           call omgw(i,k,im,km,g12z,g13x,u,w,omg)
           call div3(i,k,im,km,dx,dz,g12,g13xz,u,w,div,     &
                &    divx,divz)

           !--- $BB.EY(B u
           call vel_u(i,k,im,km,cp,alp,dts,g13z,divx,divz,  &
                &    dx,dz,vptemp_bs,pi,u_adv,u_cori,u_trb,u)

        end do
     end do

     !--- SSL2 $B;HMQ(B: $BD>8rD>@~:BI87O(B($BCO7A$J$7(B)$B$G:n$C$F$"$k(B
     do k = 1,km
        do i = 1,im
           !--- $B2;B.(B
           c_bs(i,k) = cp*rd/cv*pi_bs(i,k)*vptemp_bs(i,k)
           !--- $B9TNs@.J,$N78?t(B
           ssl_fctr(i,k) = (c_bs(i,k)*dts/dens_bs(i,k)      &
                &           /ptemp_bs(i,k))**2/dz**2
        end do
     end do

     !--- $B05NOJ}Dx<0$r2r$/$?$a$N9TNs(B: $BBP3QItJ,(B
     k = 1
     do i = 1,im
        d(i,k) = 1                                       &
             &   + ssl_fctr(i,k)                         &
             &     *dens_bs(i,k)*ptemp_bs(i,k)**2
     end do

     do k = 2,km-1
        do i = 1,im
           d(i,k) = 1                                       &
                &   + ssl_fctr(i,k)                         &
                &     *(                                    &
                &       dens_bs(i,k)*ptemp_bs(i,k)**2       &
                &       + dens_bs(i,k-1)*ptemp_bs(i,k-1)**2 &
                &      )
        end do
     end do

     k = km
     do i = 1,im
        d(i,k) = 1                                       &
             &   + ssl_fctr(i,k)                         &
             &     *dens_bs(i,k-1)*ptemp_bs(i,k-1)**2
     end do

     !--- $B05NOJ}Dx<0$r2r$/$?$a$N9TNs(B: $B2<J!BP3QItJ,(B
     do k = 2,km
        do i = 1,im
           sbd(i,k) = - ssl_fctr(i,k)                      &
                &       *dens_bs(i,k-1)*ptemp_bs(i,k-1)**2
        end do
     end do

     !--- $B05NOJ}Dx<0$r2r$/$?$a$N9TNs(B: $B>eJ!BP3QItJ,(B
     do k = 1,km-1
        do i = 1,im
           spd(i,k) = - ssl_fctr(i,k)                       &
                &       *(dens_bs(i,k)*ptemp_bs(i,k)**2)
        end do
     end do

     do i = 1,im
        !--- $BDj?t%Y%/%H%k(B
        b(i,1) = pi(i,1)                                    &
             &   - c_bs(i,1)**2/(cp*ptemp_bs(i,1))          &
             &     *dts*(u(i+1,1) - u(i,1))/dx              &
             &   - c_bs(i,1)**2                             &
             &     /(cp*dens_bs(i,1)*ptemp_bs(i,1)**2)*dts  &
             &       *(                                     &
             &         w(i,2)                               &
             &         + (                                  &
             &            - w_adv(i,2)+ w_buoy(i,2)         &
             &            + w_trb(i,2)                      &
             &           )*dts                              &
             &         - dens_bs(i,1)*ptemp_bs(i,1)         &
             &           *(                                 &
             &             w(i,1)                           &
             &             + (                              &
             &                - w_adv(i,1)+ w_buoy(i,1)     &
             &                + w_trb(i,1)                  &
             &               )*dts                          &
             &            )                                 &
             &         )/dz                                 &
             &   - (c_bs(i,1)*dts)**2                       &
             &     /(dens_bs(i,1)*ptemp_bs(i,1)**2)/dz      &
             &     *dens_bs(i,1)*ptemp_bs(i,1)              &
             &      *(- w_adv(i,1)+ w_buoy(i,1)+ w_trb(i,1))&
             &      /cp
        do k = 2,km-1
           b(i,k) = pi(i,k)                                 &
                &   - dts*c_bs(i,k)**2/(cp*ptemp_bs(i,k))   &
                &     *(u(i+1,k) - u(i,k))/dx               &
                &   - dts*c_bs(i,k)**2                      &
                &     /(cp*dens_bs(i,k)*ptemp_bs(i,k)**2)   &
                &     *(                                    &
                &       dens_bs(i,k+1)*ptemp_bs(i,k+1)*dts  &
                &       *(                                  &
                &         w(i,k+1)                          &
                &         + (                               &
                &            - w_adv(i,k+1)+ w_buoy(i,k+1)  &
                &            + w_trb(i,k+1)                 &
                &           )*dts                           &
                &        )                                  &
                &       - dens_bs(i,k-1)*ptemp_bs(i,k-1)    &
                &         *(                                &
                &           w(i,k-1)                        &
                &           + (                             &
                &              - w_adv(i,k-1)+ w_buoy(i,k-1)&
                &              + w_trb(i,k-1)               &
                &             )*dts                         &
                &          )                                &
                &      )/(2.0D0*dz)
        end do
        b(i,km) = pi(i,km)                                  &
             &   - c_bs(i,km)**2/(cp*ptemp_bs(i,km))        &
             &     *dts*(u(i+1,km) - u(i,km))/dx            &
             &   - c_bs(i,km)**2                            &
             &     /(cp*dens_bs(i,km)*ptemp_bs(i,km)**2)*dts&
             &       *(                                     &
             &         w(i,2)                               &
             &         + (                                  &
             &            - w_adv(i,2)+ w_buoy(i,2)         &
             &            + w_trb(i,2)                      &
             &           )*dts                              &
             &         - dens_bs(i,km)*ptemp_bs(i,km)       &
             &           *(                                 &
             &             w(i,km)                          &
             &             + (                              &
             &                - w_adv(i,km)+ w_buoy(i,km)   &
             &                + w_trb(i,km)                 &
             &               )*dts                          &
             &            )                                 &
             &         )/dz                                 &
             &   + (c_bs(i,km)*dts)**2                      &
             &     /(dens_bs(i,km)*ptemp_bs(i,km)**2)/dz    &
             &     *dens_bs(i,km)*ptemp_bs(i,km)            &
             &     *(- w_adv(i,km)+ w_buoy(i,km)+ w_trb(i,km))&
             &      /cp

        !--- $BL5<!8505NO(B
        call dltx(sbd,d,spd,km,b,0.0,2,is,ip,vw,icon)
        write(*,*) is,ip,vw,icon
     end do

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

           !--- $BB.EY(B w
           w(i,k) = 0.0D0
           call vel_w(i,k,im,km,dts,cp,alp,beta,dz,ptemp_bs,&
                &     pi,pi_old,g12z,divz,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,g12z,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,g12,g12z,&
             &              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
!        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

!     do k = 1,km
!        do i = 1,im
!           call omgw(i,k,im,km,g12z,g13x,u,w,omg)
!           call div3(i,k,im,km,dx,dz,g12,g13xz,u,w,div,     &
!                &    divx,divz)

           !--- $BB.EY(B u
!           call vel_u(i,k,im,km,cp,alp,dts,g13z,divx,divz,  &
!                &    dx,dz,vptemp_bs,omg,pi,u)

           !--- SSL2 $B;HMQ(B
!           call 
           !--- $B05NO(B
!           pi(i,k) = 1.0D0
           !--- $BL5<!8505NO(B pi -> $B05NO(B prss $BJQ49(B
!           call prss2pi(1,i,k,im,km,prss_sfc,rd,cp,prss,pi)

           !--- $BB.EY(B w
!           w(i,k) = 0.0D0

           !--- $B2>290L(B
!           call ver_ptemp(i,k,im,km,ptemp,qv,vptemp)

           !--- $BL)EY(B
!           call density(i,k,im,km,prss_sfc,rd,cp,cv,pi,     &
!                &       vptemp,dens)

!        end do
!     end do

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

!  end do

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

!     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,g12,g12z,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 1$B%9%F%C%W0J9_(B)
!     do its = 1,nts-1

!        do k = 1,km
!           do i = 1,im
!              call omgw(i,k,im,km,g12z,g13x,u,w,omg)
!              call div3(i,k,im,km,dx,dz,g12,g13xz,u,w,div,     &
!                   &    divx,divz)

              !--- $BB.EY(B u
!              call vel_u(i,k,im,km,cp,alp,dts,g13z,divx,divz,   &
!                   &    dx,dz,vptemp_bs,omg,pi,u)

              !--- $BB.EY(B w
!              w(i,k) = 0.0D0

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

              !--- $B2>290L(B
!              call ver_ptemp(i,k,im,km,ptemp,qv,vptemp)

              !--- $BL)EY(B
!              call density(i,k,im,km,prss_sfc,rd,cp,cv,pi,     &
!                   &       vptemp,dens)

!           end do
!        end do

        !--- $B=PNO(B
!        call output_gtool4
!     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=(/im+2,km+2,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-vorticity', units='m/s', xtype='double')
    ! $BB.EY(B w
    call HistoryAddVariable(                                &
         & varname='w', dims=(/'x','z','t'/),           &
         & longname='w-vorticity', 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='g12', dims=(/'x','z','t'/),           &
!         & longname='G12', units='1', xtype='double')
!    call HistoryAddVariable(                                &
!         & varname='g12z', dims=(/'x','z','t'/),           &
!         & longname='G12 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
!    write(6,*) 'i =',i,', k =',k
    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
!    write(6,*) 'its =',its,', i =',i,', k =',k
    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('g12', g12)
!    call HistoryPut('g12z', g12z)
!    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
