!=======================================
! 2D cumulus model - kaminari
!  - subroutine adv_s
! 
! Author    : TAKAHASHI Koko
! Date      : 2003/12/19 $B?75,:n@.(B
! Note      : $B%9%+%i!<0\N.(B
! 
!=======================================
program adv

  use gtool_history

  implicit none

  integer(8) :: i,k
  integer(8) :: its
!  integer(8), parameter :: nts = 10   ! $BC;$$;~4V%9%F%C%W?t(B
  integer(8), parameter :: ntb = 1500   ! $BC;$$;~4V%9%F%C%W?t(B

  integer(8), parameter :: im = 23      ! $B3J;RE@$N@_Dj(B (x)
  integer(8), parameter :: km = 19      ! $B3J;RE@$N@_Dj(B (z)
!  real(8), parameter :: dts = 1.0d0     ! $BC;$$;~4V4V3V(B d tau
  real(8), parameter :: dtb = 5.0d0   ! $BD9$$;~4V4V3V(B dt

  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+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
  real(8), parameter :: nuh = 1.0d-3*dx*dx/dtb ! $B?eJ?J}8~$N?tCMG4@-9`$N78?t(B
  real(8), parameter :: nuv = 5.0d-3*dz*dz/dtb ! $B1tD>J}8~$N?tCMG4@-9`$N78?t(B
!  real(8), parameter :: nuh = 0.0d0
!  real(8), parameter :: nuv = 0.0d0

  real(8) :: f(-1:im+1,-1:km+1)
  real(8) :: f_old(-1:im+1,-1:km+1)
  real(8) :: f_old2(-1:im+1,-1:km+1)
  real(8) :: f_cal(0:im,0:km)

  real(8) :: u(-1:im+1,-1:km+1)
  real(8) :: omg(-1:im+1,-1:km+1)

  real(8) :: f_adv(-1:im+1,-1:km+1)
  real(8) :: f_dif(-1:im+1,-1:km+1)
  real(8) :: x(-1:im+1,-1:km+1)
  real(8) :: z(-1:im+1,-1:km+1)
  real(8) :: x_cal(-1:im+1,-1:km+1)
  real(8) :: z_cal(-1:im+1,-1:km+1)

  !--- $B:BI8<4(B
  do k = -1,km+1
     do i = -1,im+1
        x(i,k) = dx*i
        z(i,k) = dz*k
     end do
  end do

  do k = 0,km
     do i = 0,im
        x_cal(i,k) = dx*i
        z_cal(i,k) = dz*k
     end do
  end do

        call HistoryCreate(                                    &
             & file='iryu.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(dtb))

        ! $B%Y%/%H%kNL$NJQ?t=PNO(B
        call HistoryPut('x',x_cal(0:im,0))
        call HistoryPut('z',z_cal(0,0:km))

        call HistoryAddVariable(                                &
             & varname='f', dims=(/'x','z','t'/),               &
             & longname='scalar', units='m/s',              &
             & xtype='double')

 do k = -1,km+1
    do i = -1,im+1
        u(i,k) = 1.0d1
        omg(i,k) = 0.0d0
     end do
  end do

  do k = -1,km+1
     do i = -1,im+1
        f(i,k) = exp(-(x(i,k) - 1.0d4/2.0d0)**2.0d0/5.0d6 &
             &       -(z(i,k) - 1.0d4/2.0d0)**2.0d0/5.0d6)
        f_old(i,k) = f(i,k)
!        f_old2(i,k) = f_old(i,k)
     end do
  end do

 do k = 0,km
    do i = 0,im
       f_cal(i,k) = f(i,k)
     end do
  end do

  write(*,*) 'its =', its
  call HistoryPut('f', f_cal)

do its = 1,ntb
 
 do k = 0,km
    do i = 0,im
        !--- $B?tCMG4@-9`(B
        f_dif(i,k) = nuh*(f(i+1,k) - 2.0d0*f(i,k) + f(i-1,k))     &
             &       /dx**2.0d0                                   &
             &       + nuv*(f(i,k+1) - 2.0d0*f(i,k) + f(i,k-1))   &
             &         /dz**2.0d0

        !--- $B0\N.9`(B
        f_adv(i,k) = (u(i+1,k) + u(i,k))/2.0d0              &
             &           *(f(i+1,k) - f(i-1,k))/(2.0d0*dx)  &
             &           + (omg(i,k+1) + omg(i,k))/2.0d0    &
             &             *(f(i,k+1) - f(i,k-1))/(2.0d0*dz)&
             &           + f_dif(i,k)

        f(i,k) = f_old2(i,k) - 2.0d0*dtb*f_adv(i,k)
        f_old2(i,k) = f_old(i,k)
        f_old(i,k) = f(i,k)
     end do
  end do

 do k = 0,km
    do i = 0,im
       f_cal(i,k) = f(i,k)
     end do
  end do

  write(*,*) 'its =', its
  call HistoryPut('f', f_cal)

end do

call HistoryClose

end program adv
