!=======================================
! 2D cumulus model - kaminari
!  - subroutine ptn_temp
! 
! Author    : TAKAHASHI Koko
! Date      : 2003/11/25  $B:G=*99?7(B
!             2003/11/25  $B?75,:n@.(B
! Note      : 2 $B%9%F%C%WL\0J9_$N290L(B
! 
!=======================================
!--- 1,2,3$B9TL\F~NOJQ?t(B, 4$B9TL\F~=PNOJQ?t(B
subroutine ptn_temp(i,k,im,km,nuh,nuv,cp,ml,dx,dz,dtb,u,    &
     &              g_sqrt,g13x,g13z,g13xz,                 &
     &              prss_bs,pi_bs,dens_bs,omg,qv,qr,qvs,    &
     &              ptemp_old2,km_sub_old2,                 &
     &              ptemp)

  implicit none

  integer(8), intent(in) :: i, k
  integer(8), 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) :: g_sqrt(-2:im+1,-2:km+1)
  real(8), intent(in) :: g13x(-2:im+1,-2:km+1)
  real(8), intent(in) :: g13z(-2:im+1,-2:km+1)
  real(8), intent(in) :: g13xz(-2:im+1,-2:km+1)

  real(8), intent(in) :: prss_bs(-2:im+1,-2:km+1)
  real(8), intent(in) :: pi_bs(-2:im+1,-2:km+1)
  real(8), intent(in) :: dens_bs(-2:im+1,-2:km+1)

  real(8), intent(in) :: u(-2:im+1,-2:km+1)
  real(8), intent(in) :: omg(-2:im+1,-2:km+1)
  real(8), intent(in) :: qv(-2:im+1,-2:km+1)
  real(8), intent(in) :: qr(-2:im+1,-2:km+1)
  real(8), intent(in) :: qvs(-2:im+1,-2:km+1)

  real(8), intent(in) :: ptemp_old2(-2:im+1,-2:km+1)
  real(8), intent(in) :: km_sub_old2(-2:im+1,-2:km+1)

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

  real(8)             :: ptemp_dif(-2:im+1,-2:km+1)
  real(8)             :: ptemp_adv(-2:im+1,-2:km+1)
  real(8)             :: ptemp_mic(-2:im+1,-2:km+1)
  real(8)             :: ptemp_trb(-2:im+1,-2:km+1)
  real(8)             :: c(-2:im+1,-2:km+1)
  real(8)             :: er(-2:im+1,-2:km+1)
  real(8)             :: gamma(-2:im+1,-2: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))/(2.0d0*dx)&
       &           + (omg(i,k+1) + omg(i,k))/2.0d0          &
       &             *(ptemp(i,k+1) - ptemp(i,k-1))         &
       &             /(2.0d0*dz)&
       &           + ptemp_dif(i,k)

  !--- $B290L$N@8@.!&>CLG9`$N78?t72(B
  c(i,k) = 1.6d0 + 1.249d2*(dens_bs(i,k)*qr(i,k))**2.046d-1
  
  er(i,k) = (1.0d0 - qv(i,k)/qvs(i,k))*c(i,k)               &
      &      *((dens_bs(i,k)*qr(i,k))**5.25d-1)             &
      &      /dens_bs(i,k)                                  &
      &      /(5.4d5 + 2.55d6/(prss_bs(i,k)*qvs(i,k)))

  !--- $B290L$N0\N.9`(B
  ptemp_trb(i,k) = 3.0d0                                    &
       &           *(                                       &
       &             (km_sub_old2(i+1,k) + km_sub_old2(i,k))&
       &             /2.0d0                                 &
       &             *(                                     &
       &               (ptemp_old2(i+1,k) - ptemp_old2(i,k))&
       &               /dx                                  &
       &               + g13z(i+1,k)                        &
       &                 *(                                 &
       &                   ptemp_old2(i+1,k+1)              &
       &                   + ptemp_old2(i,k+1)              &
       &                   - ptemp_old2(i+1,k-1)            &
       &                   - ptemp_old2(i,k-1)              &
       &                  )/(4.0d0*dz)                      &
       &               )                                    &
       &             - (                                    &
       &                km_sub_old2(i,k)                    &
       &                + km_sub_old2(i-1,k)                &
       &               )/2.0d0                              &
       &               *(                                   &
       &                 (                                  &
       &                  ptemp_old2(i,k)                   &
       &                  - ptemp_old2(i-1,k)               &
       &                 )/dx                               &
       &                 + g13z(i,k)                        &
       &                   *(                               &
       &                     ptemp_old2(i,k+1)              &
       &                     + ptemp_old2(i-1,k+1)          &
       &                     - ptemp_old2(i,k-1)            &
       &                     - ptemp_old2(i-1,k-1)          &
       &                    )/(4.0d0*dz)                    &
       &                 )                                  &
       &            )/dx                                    &
       &            + 3.0d0*g13xz(i,k)                      &
       &             *(                                     &
       &               (                                    &
       &                km_sub_old2(i,k+1)                  &
       &                + km_sub_old2(i,k)                  &
       &                )/2.0d0                             &
       &                *(                                  &
       &                  (                                 &
       &                   ptemp_old2(i+1,k+1)              &
       &                   + ptemp_old2(i+1,k)              &
       &                   - ptemp_old2(i-1,k+1)            &
       &                   - ptemp_old2(i-1,k)              &
       &                  )/(4.0d0*dx)                      &
       &                 + g13x(i,k+1)                      &
       &                   *(                               &
       &                     ptemp_old2(i,k+1)              &
       &                     - ptemp_old2(i,k)              &
       &                    )/dz                            &
       &                 )                                  &
       &               - (                                  &
       &                  km_sub_old2(i,k)                  &
       &                  + km_sub_old2(i,k-1)              &
       &                  )/2.0d0                           &
       &                  *(                                &
       &                    (                               &
       &                     ptemp_old2(i+1,k)              &
       &                     + ptemp_old2(i+1,k-1)          &
       &                     - ptemp_old2(i-1,k)            &
       &                     - ptemp_old2(i-1,k-1)          &
       &                    )/(4.0d0*dx)                    &
       &                   + g13x(i,k)                      &
       &                     *(                             &
       &                       ptemp_old2(i,k)              &
       &                       - ptemp_old2(i,k-1)          &
       &                      )/dz                          &
       &                   )                                &       
       &              )                                     &
       &            + 3.0d0/g_sqrt(i,k)                     &
       &              *(                                    &
       &                1.0d0/g_sqrt(i,k)*km_sub_old2(i,k)  &
       &                *(                                  &
       &                  ptemp_old2(i,k+1)                 &
       &                  - ptemp_old2(i,k)                 &
       &                 )/dz                               &
       &                - 1.0d0/g_sqrt(i,k-1)               &
       &                  *km_sub_old2(i,k-1)               &
       &                  *(                                &
       &                    ptemp_old2(i,k)                 &
       &                    - ptemp_old2(i,k-1)             &
       &                   )/dz                             &
       &               )/dz

  !--- $B%,%s%^(B
  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_old2(i,k)                              &
       &       + 2.0d0*dtb                                  &
       &         *(                                         &
       &           - ptemp_adv(i,k) + ptemp_mic(i,k)        &
       &           + ptemp_trb(i,k)                         &
       &          )

end subroutine ptn_temp
