!=======================================
! 2D cumulus model - kaminari
!  - subroutine ptemp_1stp
! 
! Author    : TAKAHASHI Koko
! Date      : 2003/11/05
! Note      : 1 $B%9%F%C%WL\$N290L(B
! 
!=======================================
subroutine 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)

  implicit none

  integer, intent(in) :: i, k
  integer, intent(in) :: im, km

  real(8), intent(in) :: nuh, nuv,cp, ml
  real(8), intent(in) :: dx, dz
  real(8), intent(in) :: dtb

  real(8), intent(in) :: prss_bs(0:im+1,0:km+1)
  real(8), intent(in) :: pi_bs(0:im+1,0:km+1)
  real(8), intent(in) :: dens_bs(0:im+1,0:km+1)
  real(8), intent(in) :: u(0:im+1,0:km+1)
  real(8), intent(in) :: omg(0:im+1,0:km+1)
  real(8), intent(in) :: qv(0:im+1,0:km+1)
  real(8), intent(in) :: qr(0:im+1,0:km+1)
  real(8), intent(in) :: qvs(0:im+1,0:km+1)

  real(8), intent(inout) :: ptemp(0:im+1,0:km+1)

  real(8)             :: ptemp_dif(0:im+1,0:km+1)
  real(8)             :: ptemp_adv(0:im+1,0:km+1)
  real(8)             :: ptemp_mic(0:im+1,0:km+1)
  real(8)             :: c(0:im+1,0:km+1)
  real(8)             :: er(0:im+1,0:km+1)
  real(8)             :: gamma(0:im+1,0:km+1)

  !--- $B290L$N?tCMG4@-9`(B
  ptemp_dif(i,k) = nuh*(                                    &
       &                ptemp(i+1,k) - 2.0D0*ptemp(i,k)     &
       &                + ptemp(i-1,k)                      &
       &               )/(dx**2.0D0)                        &
       &           + nuv*(                                  &
       &                  ptemp(i,k+1) - 2.0D0*ptemp(i,k)   &
       &                  + ptemp(i,k-1)                    &
       &                  )/(dz**2.0D0)

  !--- $B290L$N0\N.9`(B
  ptemp_adv(i,k) = (u(i+1,k) + u(i,k))/2.0D0                &
       &           *(ptemp(i+1,k) - ptemp(i-1,k))/(4.0D0*dx)&
       &          + (omg(i,k+1) + omg(i,k))/2.0D0           &
       &            *(ptemp(i,k+1) - ptemp(i,k-1))/(4.0D0*dz)&
       &          + ptemp_dif(i,k)

  !--- $B290L$N@8@.!&>CLG9`$N78?t72(B
  c(i,k) = 1.6D0 + 124.9D0*(dens_bs(i,k)*qr(i,k))**0.2046D0
  
  er(i,k) = (1 - qv(i,k)/qvs(i,k))*c(i,k)                  &
      &       *((dens_bs(i,k)*qr(i,k))**0.525D0)           &
      &       /(5.4e5 + 2.55e6/(prss_bs(i,k)*qvs(i,k)))

  gamma(i,k) = ml/(cp*pi_bs(i,k))

  !--- $B290L$N@8@.!&>CLG9`(B
  ptemp_mic(i,k) = - gamma(i,k)*er(i,k)

  !--- $B290L(B
  ptemp(i,k) = ptemp(i,k) + dtb*(- ptemp_adv(i,k) + ptemp_mic(i,k))

end subroutine ptemp_1stp
