!=================================================
! 2D cumulus model - kaminari
!
! Author  : TAKAHASHI Koko
! Date    : 2003/12/06  $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 = 24      ! $B3J;RE@$N@_Dj(B (x)
  integer(8), parameter :: km = 20      ! $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 = 2.4d4    ! $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

  ! $B;~4VJQ?t(B
  integer(8), parameter :: nts = 100    ! $BC;$$;~4V%9%F%C%W?t(B
  integer(8), parameter :: ntb = 1      ! $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

  ! $B%Y%/%H%kMQ3J;R(B
  real(8) :: x_u(-2:im,-2:km+1)         ! $BCO7A$K1h$C$?:BI87O$G$N(B x $B:BI8(B
  real(8) :: z_u(-2:im+1,-2:km)         ! $BCO7A$K1h$C$?:BI87O$G$N(B z $B:BI8(B
  real(8) :: x_w(-2:im+1,-2:km)         ! $BCO7A$K1h$C$?:BI87O$G$N(B x $B:BI8(B
  real(8) :: z_w(-2:im,-2:km+1)         ! $BCO7A$K1h$C$?:BI87O$G$N(B z $B:BI8(B

  ! $B%9%+%i!<MQ3J;R(B
  real(8) :: x_s(-2:im,-2:km)           ! $BCO7A$K1h$C$?:BI87O$G$N(B x $B:BI8(B
  real(8) :: z_s(-2:im,-2:km)           ! $BCO7A$K1h$C$?:BI87O$G$N(B z $B:BI8(B

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

  !--- $B>qMp$NJQ?t(B
  real(8) :: u(-2:im,-2:km+1)           ! x $BJ}8~B.EY(B
  real(8) :: w(-2:im+1,-2:km)           ! z $BJ}8~B.EY(B
  real(8) :: omg(-2:im+1,-2:km)         ! z^{*} $BJ}8~B.EY(B

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

  !--- $B4pK\>l(B + $B>qMp>l$NJQ?t(B
  real(8) :: temp(-2:im,-2:km)      ! $B29EY(B
  real(8) :: prss(-2:im,-2:km)      ! $B05NO(B
  real(8) :: pi_all(-2:im,-2:km)    ! $BL5<!8505NO(B
  real(8) :: ptemp_all(-2:im,-2:km) ! $B290L(B
  real(8) :: dens(-2:im,-2:km)     ! $BL)EY(B

  !--- $B4pK\>l(B + $B>qMp>lJQ?t(B ($B7W;;NN0h$N$_(B)
  real(8) :: u_cal(0:im-2,0:km-1)         ! x $BJ}8~B.EY(B
  real(8) :: omg_cal(0:im-1,0:km-2)         ! z $BJ}8~B.EY(B

  real(8) :: prss_cal(0:im-2,0:km-2)      ! $B05NO(B
  real(8) :: pi_all_cal(0:im-2,0:km-2)    ! $BL5<!8505NO(B
  real(8) :: temp_cal(0:im-2,0:km-2)      ! $B29EY(B
  real(8) :: ptemp_all_cal(0:im-2,0:km-2) ! $B290L(B
  real(8) :: dens_cal(0:im-2,0:km-2)     ! $BL)EY(B

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

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

  !--- $B$=$NB>$NJQ?t(B         
  real(8) :: g_sqrt(-2:im,-2:km)    ! $B%d%3%S%"%s$N(B1/2$B>h(B
  real(8) :: g_sqrt_x(-2:im,-2:km)  ! g_sqrt$B$r(Bx$BJ}8~$KJ?6Q2=(B
  real(8) :: g_sqrt_z(-2:im,-2:km)  ! g_sqrt$B$r(Bz^{*}$BJ}8~$KJ?6Q2=(B
  real(8) :: g_sqrt_xz(-2:im,-2:km) ! g_sqrt$B$r(Bz^{*}$BJ}8~$KJ?6Q2=(B
  real(8) :: g13(-2:im,-2:km)       ! $B%d%3%S%"%s$N(B13$B@.J,(B
  real(8) :: g13x(-2:im,-2:km)      ! g13$B$r(Bx$BJ}8~$KJ?6Q2=(B
  real(8) :: g13z(-2:im,-2:km)      ! g13$B$r(Bz^{*}$BJ}8~$KJ?6Q2=(B
  real(8) :: g13xz(-2:im,-2:km)     ! 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 = 0.0d0     ! $B=ENO(B
  real(8), parameter :: alp = 5.0d-13*dx*dx/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-2*dx*dx ! $B?eJ?J}8~$N?tCMG4@-9`$N78?t(B
!  real(8), parameter :: nuv = 1.0d-2*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)

  !--- $B=PNO(B (gtool history) $BMQ(B
  type(GT_HISTORY) :: vector_u, vector_u_cal
  type(GT_HISTORY) :: vector_w, vector_w_cal
  type(GT_HISTORY) :: scalar, scalar_cal

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

  !-------------------------------------------------
  ! x-z^{*} $B:BI8$N@_Dj(B
  !-------------------------------------------------
  !--- $B%Y%/%H%kNL(B u
  do k = -2, km+1
     do i = -2, im
        x_u(i,k) = dx*(i + 5.0d-1)
        z_u(i,k) = dz*k
     end do
  end do

  !--- $B%Y%/%H%kNL(B u
  do k = -2, km
     do i = -2, im+1
        x_w(i,k) = dx*i
        z_w(i,k) = dz*(k + 5.0d-1)
     end do
  end do

  !--- $B%9%+%i!<NL(B
  do k = -2, km
     do i = -2, im
        x_s(i,k) = dx*(i + 5.0d-1)
        z_s(i,k) = dz*(k + 5.0d-1)
     end do
  end do

  call output_gtool4_init

  !-------------------------------------------------
  ! $B4pK\>l(B
  !-------------------------------------------------
  !--- 1 $B9TL\(B: $BF~NO(B, 2,3 $B9TL\(B: $B=PNO(B
  call base(im,km,dz,prss_sfc,rd,cp,cv,grv,g_sqrt,          &
       &    qv_bs,prss_bs,pi_bs,temp_bs,ptemp_bs,           &
       &    vptemp_bs,dens_bs)

  !-------------------------------------------------
  ! $B4pK\>l$N6-3&>r7o(B
  !--- 1 $B9TL\(B: $BF~NO(B, 2,3 $B9TL\(B: $B=PNO(B
  call boundary_bs(im,km,                                   &
       &           qv_bs,prss_bs,pi_bs,temp_bs,ptemp_bs,    &
       &           vptemp_bs,dens_bs)

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

  !-------------------------------------------------
  ! $B>qMp>l$N7W;;(B
  !------------------------------------------------- 

  !-------------------------------------------------
  ! $B=i4|CM(B
  !  its = 0,    itb = 0
  !--- $BB.EY(B, $B05NO(B, $B29EY$rM?$($k(B
  !--- 1$B9TL\F~NOJQ?t(B, 2$B9TL\=PNOJQ?t(B
  call init_value(im,km,dz,cp,rd,prss_sfc,pi_bs, &
     &                u,w,prss,pi,pi_all,temp,ptemp,vptemp)

  !--- $B6-3&>r7o$H$=$l$h$j30$NG[NsCM(B
  !--- 1$B9TL\F~NOJQ?t(B, 2$B9TL\=PNOJQ?t(B
  call boundary(im,km,                      &
       &        u,w,pi,ptemp)

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

        !--- $B4pK\>l(B+$B>qMp>l$NCM(B
        pi_all(i,k) = pi_bs(i,k) + pi(i,k)
        ptemp_all(i,k) = ptemp_bs(i,k) + ptemp(i,k)
     end do
  end do

  !--- $BCO7A$K1h$C$?:BI87O$G$N1tD>B.EY(B
  call omgw(im,km,g_sqrt_z,g13x,u,w,                        &
             &    omg)

  !--- $B7W;;NN0h$N$_=PNO$9$k$?$a(B
  call cal_value(im,km,u,omg,prss,pi_all,temp,ptemp_all,    &
       &         dens,                                      &
       &         u_cal,omg_cal,prss_cal,pi_all_cal,temp_cal,&
       &         ptemp_all_cal,dens_cal)

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

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

     call vel_w_tb(im,km,grv,ptemp_bs,ptemp,           &
          &        w_adv,w_buoy,w_trb)

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

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

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

        !--- $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(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,g_sqrt_z, &
             &          g13xz,div_z_w,                      &
             &          pi)

        !--- $BB.EY(B w
        !--- 1,2,3$B9TL\F~NOJQ?t(B, 4$B9TL\F~=PNOJQ?t(B
        call vel_w(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)

        !--- $B6-3&>r7o$H$=$l$h$j30$NG[NsCM(B
        call boundary(im,km,                                &
             &        u,w,pi,ptemp)

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

              !--- $B4pK\>l(B+$B>qMp>l$NCM(B
              pi_all(i,k) = pi_bs(i,k) + pi(i,k)
              ptemp_all(i,k) = ptemp_bs(i,k) + ptemp(i,k)
              
              vptemp(i,k) = ptemp_all(i,k)
              dens(i,k) = prss_sfc*pi_all(i,k)**(cv/rd)/(rd*vptemp(i,k))

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

        !--- $BCO7A$K1h$C$?:BI87O$G$N1tD>B.EY(B
        call omgw(im,km,g_sqrt_z,g13x,u,w,                  &
             &    omg)

        !--- $B7W;;NN0h$N$_=PNO$9$k$?$a(B
        call cal_value(im,km,u,omg,prss,pi_all,temp,        &
             &         ptemp_all,dens,                      &
             &         u_cal,omg_cal,prss_cal,pi_all_cal,   &
             &         temp_cal,ptemp_all_cal,dens_cal)

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

     end do

  end do

  call output_gtool4_close

contains

!-----------------------------------------------------------
!--- gtool4 $B=PNO4XO"(B

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

    ! $B%Y%/%H%kNL$NJQ?t=PNO(B
    call HistoryPut('x',x_u(-2:im,0),vector_u)
    call HistoryPut('z',z_u(0,-2:km+1),vector_u)

    ! $B%Y%/%H%kNL(B w
    call HistoryCreate(                                     &
          & file='kaminari_2d_w.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+3,4),0/),         &
          & longnames=(/'X-coordinate','Z-coordinate',      &
          &             '    Time    '/),                   &
          & units=(/'m','m','s'/), origin=0.0,              &
          & interval=real(ndisp), history=vector_w)

    ! $B%Y%/%H%kNL$NJQ?t=PNO(B
    call HistoryPut('x',x_w(-2:im+1,0),vector_w)
    call HistoryPut('z',z_w(0,-2:km),vector_w)

    ! $B%9%+%i!<NL(B
    call HistoryCreate(                                     &
          & file='kaminari_2d_scl.nc',                      &
          & title='2D cumulus model',                       &
          & source='kaminari - 2D cumulus model',           &
          & institution='GFD_Dennou Club deepconv project', &
          & dims=(/'x','z','t'/),                           &
          & dimsizes=(/int(im+3,4),int(km+3,4),0/),         &
          & longnames=(/'X-coordinate','Z-coordinate',      &
          &             '    Time    '/),                   &
          & units=(/'m','m','s'/), origin=0.0,              &
          & interval=real(ndisp),history=scalar)

    ! $B%9%+%i!<NL$NJQ?t=PNO(B
    call HistoryPut('x',x_s(-2:im,0),scalar)
    call HistoryPut('z',z_s(0,-2:km),scalar)

    !--- $B7W;;NN0h$N$_(B
    ! $B%Y%/%H%kNL(B u
    call HistoryCreate(                                     &
          & file='kaminari_2d_u_cal.nc',                    &
          & title='2D cumulus model',                       &
          & source='kaminari - 2D cumulus model',           &
          & institution='GFD_Dennou Club deepconv project', &
          & dims=(/'x','z','t'/),                           &
          & dimsizes=(/int(im-1,4),int(km,4),0/),           &
          & longnames=(/'X-coordinate','Z-coordinate',      &
          &             '    Time    '/),                   &
          & units=(/'m','m','s'/), origin=0.0,              &
          & interval=real(ndisp), history=vector_u_cal)

    ! $B%Y%/%H%kNL$NJQ?t=PNO(B
    call HistoryPut('x',x_u(0:im-2,0),vector_u_cal)
    call HistoryPut('z',z_u(0,0:km-1),vector_u_cal)

    ! $B%Y%/%H%kNL(B w
    call HistoryCreate(                                     &
          & file='kaminari_2d_w_cal.nc',                    &
          & title='2D cumulus model',                       &
          & source='kaminari - 2D cumulus model',           &
          & institution='GFD_Dennou Club deepconv project', &
          & dims=(/'x','z','t'/),                           &
          & dimsizes=(/int(im,4),int(km-1,4),0/),           &
          & longnames=(/'X-coordinate','Z-coordinate',      &
          &             '    Time    '/),                   &
          & units=(/'m','m','s'/), origin=0.0,              &
          & interval=real(ndisp), history=vector_w_cal)

    ! $B%Y%/%H%kNL$NJQ?t=PNO(B
    call HistoryPut('x',x_w(0:im-1,0),vector_w_cal)
    call HistoryPut('z',z_w(0,0:km-2),vector_w_cal)

    ! $B%9%+%i!<NL(B
    call HistoryCreate(                                     &
          & file='kaminari_2d_scl_cal.nc',                  &
          & title='2D cumulus model',                       &
          & source='kaminari - 2D cumulus model',           &
          & institution='GFD_Dennou Club deepconv project', &
          & dims=(/'x','z','t'/),                           &
          & dimsizes=(/int(im-1,4),int(km-1,4),0/),         &
          & longnames=(/'X-coordinate','Z-coordinate',      &
          &             '    Time    '/),                   &
          & units=(/'m','m','s'/), origin=0.0,              &
          & interval=real(ndisp),history=scalar_cal)

    ! $B%9%+%i!<NL$NJQ?t=PNO(B
    call HistoryPut('x',x_s(0:im-2,0),scalar_cal)
    call HistoryPut('z',z_s(0,0:km-2),scalar_cal)

    !----------------------
    ! $BJQ?tDj5A(B
    ! $B4pK\>l$N05NO(B
    call HistoryAddVariable(                                &
         & varname='prss_bs', dims=(/'x','z'/),             &
         & longname='pressure', units='Pa', xtype='double', &
         & history=scalar)
    ! $B4pK\>l$NL5<!8505NO(B
    call HistoryAddVariable(                                &
         & varname='pi_bs', dims=(/'x','z'/),               &
         & longname='nondimensional pressure', units='1',   &
         & xtype='double' ,history=scalar)
    ! $B4pK\>l$NL)EY(B
    call HistoryAddVariable(                                &
         & varname='dens_bs', dims=(/'x','z'/),             &
         & longname='density', units='kg/m3',               &
         & xtype='double', history=scalar)
    ! $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', history=scalar)
    ! $B4pK\>l$N29EY(B
    call HistoryAddVariable(                                &
         & varname='temp_bs', dims=(/'x','z'/),             &
         & longname='temperature', units='K',               &
         & xtype='double', history=scalar)
    ! $B4pK\>l$N290L(B
    call HistoryAddVariable(                                &
         & varname='ptemp_bs', dims=(/'x','z'/),            &
         & longname='potential temperature (base)',         &
         & units='K', xtype='double', history=scalar)

    ! $BB.EY(B u
    call HistoryAddVariable(                                &
         & varname='u', dims=(/'x','z','t'/),               &
         & longname='u-velocity', units='m/s',              &
         & xtype='double', history=vector_u)
    ! $BB.EY(B w
    call HistoryAddVariable(                                &
         & varname='omg', dims=(/'x','z','t'/),             &
         & longname='w-velocity', units='m/s',              &
         & xtype='double', history=vector_w)
    ! $B05NO(B
    call HistoryAddVariable(                                &
         & varname='prss', dims=(/'x','z','t'/),            &
         & longname='pressure', units='Pa', xtype='double', &
         & history=scalar)
    ! $BL5<!8505NO(B
    call HistoryAddVariable(                                &
         & varname='pi_all', dims=(/'x','z','t'/),          &
         & longname='nondimensional pressure', units='1',   &
         & xtype='double', history=scalar)
    ! $B29EY(B
    call HistoryAddVariable(                                &
         & varname='temp', dims=(/'x','z','t'/),            &
         & longname='temperature', units='K',               &
         & xtype='double', history=scalar)
    ! $B290L(B
    call HistoryAddVariable(                                &
         & varname='ptemp', dims=(/'x','z','t'/),           &
         & longname='potential temperature', units='K',     &
         & xtype='double', history=scalar)
    ! $BL)EY(B
    call HistoryAddVariable(                                &
         & varname='dens', dims=(/'x','z','t'/),            &
         & longname='density', units='kg/m^3',              &
         & xtype='double', history=scalar)
    ! $B?e>x5$:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qv', dims=(/'x','z','t'/),              &
!         & longname='mixing ratio of water vapor',          &
!         &  units='g/g', xtype='double', history=scalar)
    ! $B1@?e:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qc', dims=(/'x','z','t'/),              &
!        & longname='mixing ratio of cloud drops',           &
!         & units='g/g', xtype='double', history=scalar)
    ! $B1+?e:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qr', dims=(/'x','z','t'/),              &
!         & longname='mixing ratio of rainwater',            &
!         & units='g/g', xtype='double', history=scalar)
    ! $BK0OB?e>x5$:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qvs', dims=(/'x','z','t'/),             &
!         & longname='saturation mixing ratio',              &
!         &  units='g/g', xtype='double', history=scalar)
    ! $B12:.9g78?t(B
!    call HistoryAddVariable(                                &
!         & varname='km_sub', dims=(/'x','z','t'/),          &
!         & longname='momentum eddy mixing coefficient',     &
!         & units='1', xtype='double', history=scalar)
    ! $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', history=scalar)

    !--- $B7W;;NN0h$N$_(B
    ! $BB.EY(B u
    call HistoryAddVariable(                                &
         & varname='u', dims=(/'x','z','t'/),               &
         & longname='u-velocity', units='m/s',              &
         & xtype='double', history=vector_u_cal)
    ! $BB.EY(B w
    call HistoryAddVariable(                                &
         & varname='omg', dims=(/'x','z','t'/),             &
         & longname='w-velocity', units='m/s',              &
         & xtype='double', history=vector_w_cal)
    ! $B05NO(B
    call HistoryAddVariable(                                &
         & varname='prss', dims=(/'x','z','t'/),            &
         & longname='pressure', units='Pa', xtype='double', &
         & history=scalar_cal)
    ! $BL5<!8505NO(B
    call HistoryAddVariable(                                &
         & varname='pi_all', dims=(/'x','z','t'/),          &
         & longname='nondimensional pressure', units='1',   &
         & xtype='double', history=scalar_cal)
    ! $B29EY(B
    call HistoryAddVariable(                                &
         & varname='temp', dims=(/'x','z','t'/),            &
         & longname='temperature', units='K',               &
         & xtype='double', history=scalar_cal)
    ! $B290L(B
    call HistoryAddVariable(                                &
         & varname='ptemp', dims=(/'x','z','t'/),           &
         & longname='potential temperature', units='K',     &
         & xtype='double', history=scalar_cal)
    ! $BL)EY(B
    call HistoryAddVariable(                                &
         & varname='dens', dims=(/'x','z','t'/),            &
         & longname='density', units='kg/m^3',              &
         & xtype='double', history=scalar_cal)
    ! $B?e>x5$:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qv', dims=(/'x','z','t'/),              &
!         & longname='mixing ratio of water vapor',          &
!         &  units='g/g', xtype='double', history=scalar_cal)
    ! $B1@?e:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qc', dims=(/'x','z','t'/),              &
!        & longname='mixing ratio of cloud drops',           &
!         & units='g/g', xtype='double', history=scalar_cal)
    ! $B1+?e:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qr', dims=(/'x','z','t'/),              &
!         & longname='mixing ratio of rainwater',            &
!         & units='g/g', xtype='double', history=scalar_cal)
    ! $BK0OB?e>x5$:.9gHf(B
!    call HistoryAddVariable(                                &
!         & varname='qvs', dims=(/'x','z','t'/),             &
!         & longname='saturation mixing ratio',              &
!         &  units='g/g', xtype='double', history=scalar_cal)
    ! $B12:.9g78?t(B
!    call HistoryAddVariable(                                &
!         & varname='km_sub', dims=(/'x','z','t'/),          &
!         & longname='momentum eddy mixing coefficient',     &
!         & units='1', xtype='double', history=scalar_cal)
    ! $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', history=scalar_cal)
  end subroutine output_gtool4_init

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

  subroutine output_gtool4
    write(*,*) 'itb =', itb,'its =', its
    call HistoryPut('u', u, vector_u)
    call HistoryPut('omg', omg, vector_w)
    call HistoryPut('prss', prss, scalar)
    call HistoryPut('pi_all', pi_all, scalar)
    call HistoryPut('temp', temp, scalar)
    call HistoryPut('ptemp_all', ptemp_all, scalar)
    call HistoryPut('dens', dens, scalar)
!    call HistoryPut('qv', qv, scalar)
!    call HistoryPut('qc', qc, scalar)
!    call HistoryPut('qr', qr, scalar)
!    call HistoryPut('qvs', qvs, scalar)
!    call HistoryPut('km_sub', km_sub, scalar)
!    call HistoryPut('e_sub', e_sub, scalar)

    call HistoryPut('u', u_cal, vector_u_cal)
    call HistoryPut('omg', omg_cal, vector_w_cal)
    call HistoryPut('prss', prss_cal, scalar_cal)
    call HistoryPut('pi_all', pi_all_cal, scalar_cal)
    call HistoryPut('temp', temp_cal, scalar_cal)
    call HistoryPut('ptemp_all', ptemp_all_cal, scalar_cal)
    call HistoryPut('dens', dens_cal, scalar_cal)
!    call HistoryPut('qv', qv_cal, scalar_cal)
!    call HistoryPut('qc', qc_cal, scalar_cal)
!    call HistoryPut('qr', qr_cal, scalar_cal)
!    call HistoryPut('qvs', qvs_cal, scalar_cal)
!    call HistoryPut('km_sub', km_sub_cal, scalar_cal)
!    call HistoryPut('e_sub', e_sub_cal, scalar_cal)
  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose(vector_u)
    call HistoryClose(vector_w)
    call HistoryClose(scalar)
    call HistoryClose(vector_u_cal)
    call HistoryClose(vector_w_cal)
    call HistoryClose(scalar_cal)
  end subroutine output_gtool4_close

end program kaminari_2d
