!=======================================
! 2D cumulus model - kaminari
!  - subroutine non_dim_pi
! 
! Author    : TAKAHASHI Koko
! Date      : 2003/12/06  $B:G=*99?7(B
!             2003/11/19  $B?75,:n@.(B
! Note      : SSL2 $B;HMQ(B
!             $B5-9f$N0UL#(B/$BI=5-(B    : $B%N!<%H(B => $B$3$3(B
!             $BBP3Q@.J,(B           :    A   => d(k)
!             $B>eI{BP3Q@.J,(B       :    B   => spd(k)
!             $B2<I{BP3Q@.J,(B       :    C   => sbd(k)
!             $BDj?t%Y%/%H%k(B($B9TNs(B) :    D   => b(k)
! 
!=======================================
!--- 1,2,3$B9TL\F~NOJQ?t(B, $B!V(Bk$B!W$H(B4$B9TL\F~=PNOJQ?t(B
subroutine 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)

  implicit none

  integer(8), intent(in) :: its
  integer(8), intent(in) :: im, km

  real(8), intent(in) :: cp,cv,rd,alp,beta
  real(8), intent(in) :: dx, dz, dts

  real(8), intent(in) :: pi_bs(-2:im+1,-2:km+1)
  real(8), intent(in) :: vptemp_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) :: w(-2:im+1,-2:km+1)
  real(8), intent(in) :: w_adv(-2:im+1,-2:km+1)
  real(8), intent(in) :: w_buoy(-2:im+1,-2:km+1)
  real(8), intent(in) :: w_trb(-2:im+1,-2:km+1)

  real(8), intent(in) :: g_sqrt(-2:im+1,-2:km+1)
  real(8), intent(in) :: g_sqrt_z(-2:im+1,-2:km+1)
  real(8), intent(in) :: g13xz(-2:im+1,-2:km+1)
  real(8), intent(in) :: div_z_w(-2:im+1,-2:km+1)

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

  integer(8)          :: i, k

  real(8)              :: c_bs2(0:km-2)        ! $B4pK\>l$G$N2;B.(B
  real(8)              :: ssl_fctr(0:km-2)     ! $B9TNs@.J,$N78?t(B
  real(8)              :: ssl_fctr2(-1:km-2)    ! $B9TNs@.J,$N78?t(B

  !--- SSL2 $B9TNs7W;;MQ$NJQ?t(B
  real(8)              :: d(0:km-2)            ! $BBP3Q@.J,(B
  real(8)              :: sbd(0:km-3)        ! $B2<I{BP3Q@.J,(B
  real(8)              :: spd(0:km-3)        ! $B>eI{BP3Q@.J,(B
  real(8)              :: b(0:km-2)            ! $BDj?t%Y%/%H%k(B($B9TNs(B)

  ! SSL2 $B%i%$%V%i%jMQ%Q%i%a!<%?(B
  real(8) :: vw(0:km-1)
  integer(8) :: isw, is, ip(0:km-1), icon

  !--- $BO"N)(B1$B<!J}Dx<0$r2r$/$?$a$N%Q%i%a!<%?(B
  isw = 1
  
  do i = 0,im-2

     do k = 0,km-2
        !--- $B2;B.$N(B2$B>h(B
        c_bs2(k) = cp*rd/cv*pi_bs(i,k)*vptemp_bs(i,k)
        !--- $B9TNs@.J,$N78?t(B $B$=$N(B1
        ssl_fctr(k) = c_bs2(k)/g_sqrt(i,k)                     &
             &          *(beta*dts/(vptemp_bs(i,k)*dz))**2.0d0 &
             &          /dens_bs(i,k)
     end do

     !--- $B9TNs@.J,$N78?t(B $B$=$N(B2
     do k = -1,km-2
        ssl_fctr2(k) = dens_bs(i,k)                            &
             &         *vptemp_bs(i,k)**2.0d0/g_sqrt(i,k)
     end do

     !--- $B05NOJ}Dx<0$r2r$/$?$a$N9TNs(B: $BBP3QItJ,(B (A)
     ! k = 0
     d(0) = 1.0d0 + ssl_fctr(0)*ssl_fctr2(0)

     do k = 1,km-3
        d(k) = 1.0d0                                           &
             & + ssl_fctr(k)*(ssl_fctr2(k) + ssl_fctr2(k-1))
     end do

     ! k = km-2
     d(km-2) = 1.0d0 + ssl_fctr(km-2)*ssl_fctr2(km-3)

     do k = 0,km-3
        !--- $B05NOJ}Dx<0$r2r$/$?$a$N9TNs(B: $B>eI{BP3QItJ,(B (B)
        spd(k) = - ssl_fctr(k)*ssl_fctr2(k)

        !--- $B05NOJ}Dx<0$r2r$/$?$a$N9TNs(B: $B2<I{BP3QItJ,(B (C)
        sbd(k) = - ssl_fctr(k)*ssl_fctr2(k-1)
     end do

     !--- $BDj?t%Y%/%H%k(B ($B9TNs(B) (D)
     b(0) = pi(i,0)                                            &
          & - (1.0d0 - beta)/g_sqrt(i,0)*c_bs2(0)*dts          &
          &   /(cp*dens_bs(i,0)*vptemp_bs(i,0)**2.0d0)         &
          &   *(                                               &
          &     (dens_bs(i,1) + dens_bs(i,0))/2.0d0            &
          &     *(vptemp_bs(i,1) + vptemp_bs(i,0))/2.0d0       &
          &     *w(i,1)                                        &
          &     - (dens_bs(i,0) + dens_bs(i,-1))/2.0d0         &
          &       *(vptemp_bs(i,0) + vptemp_bs(i,-1))/2.0d0    &
          &      *w(i,0)                                       &
          &    )/dz                                            &
          &  - c_bs2(0)*dts/(cp*vptemp_bs(i,0))                &
          &    *(                                              &
          &      (u(i+1,0) - u(i,0))/dx                        &
          &      + g13xz(i,0)                                  &
          &        *(u(i+1,1) + u(i,1) - u(i+1,-1) - u(i,-1))  &
          &        /(4.0d0*dz)                                 &
          &     )                                              &
          & - beta/g_sqrt(i,0)*c_bs2(0)*dts                    &
          &   /(cp*dens_bs(i,0)*vptemp_bs(i,0)**2.0d0)         &
          &   *(                                               &
          &     (dens_bs(i,1) + dens_bs(i,0))/2.0d0            &
          &     *(vptemp_bs(i,1) + vptemp_bs(i,0))/2.0d0       &
          &     *(                                             &
          &       w(i,1)                                       &
          &       - cp*dts/g_sqrt_z(i,1)                       &
          &         *(vptemp_bs(i,1) + vptemp_bs(i,0))/2.0d0   &
          &         *(                                         &
          &           - alp*div_z_w(i,1)                       &
          &           + (1.0d0 - beta)*(pi(i,1) - pi(i,0))/dz  &
          &          )                                         &
          &       - (- w_adv(i,1) + w_buoy(i,1) + w_trb(i,1))  &
          &         *dts                                       &
          &       )                                            &
          &     - (dens_bs(i,0) + dens_bs(i,-1))/2.0d0         &
          &       *(vptemp_bs(i,0) + vptemp_bs(i,-1))/2.0d0    &
          &       *(                                           &
          &         w(i,0)                                     &
          &         - cp*dts/g_sqrt_z(i,0)                     &
          &           *(vptemp_bs(i,0) + vptemp_bs(i,-1))/2.0d0&
          &           *(                                       &
          &             - alp*div_z_w(i,0)                     &
          &             + (1.0d0 - beta)                       &
          &               *(pi(i,0) - pi(i,-1))/dz             &
          &            )                                       &
          &         - (                                        &
          &            - w_adv(i,0) + w_buoy(i,0) + w_trb(i,0) &
          &           )*dts                                    &
          &         )                                          &
          &     )/dz                                           &
          & - beta/g_sqrt(i,0)*c_bs2(0)*dts**2.0d0             &
          &   /(dens_bs(i,0)*vptemp_bs(i,0)**2.0d0*dz)         &
          &   *cp                                              &
          &   *(vptemp_bs(i,0) + vptemp_bs(i,-1)/2.0d0)        &
          &    **2.0d0/g_sqrt_z(i,0)                           &
          &   *(                                               &
          &     alp*div_z_w(i,0)                               &
          &     - (1.0d0 - beta)*(pi(i,0) - pi(i,-1))/dz       &
          &     + g_sqrt_z(i,0)                                &
          &       /cp                                          &
          &       /(vptemp_bs(i,0) + vptemp_bs(i,-1)/2.0d0)    &
          &       *(- w_adv(i,0) + w_buoy(i,0)+ w_trb(i,0))    &
          &     )

     do k = 1,km-3
        b(k) = pi(i,k)                                            &
             &   - (1.0d0 - beta)                                 &
             &     /g_sqrt(i,k)*c_bs2(k)*dts                      &
             &     /(cp*dens_bs(i,k)*vptemp_bs(i,k)**2.0d0)       &
             &     *(                                             &
             &       (dens_bs(i,k+1) + dens_bs(i,k))/2.0d0        &
             &       *(vptemp_bs(i,k+1) + vptemp_bs(i,k))/2.0d0   &
             &       *w(i,k+1)                                    &
             &       - (dens_bs(i,k) + dens_bs(i,k-1))/2.0d0      &
             &         *(vptemp_bs(i,k) + vptemp_bs(i,k-1))/2.0d0 &
             &         *w(i,k)                                    &
             &      )/dz                                          &
             &    - c_bs2(k)*dts/(cp*vptemp_bs(i,k))              &
             &      *(                                            &
             &        (u(i+1,k) - u(i,k))/dx                      &
             &        + g13xz(i,k)                                &
             &          *(                                        &
             &            u(i+1,k+1) + u(i,k+1)                   &
             &            - u(i+1,k-1) - u(i,k-1)                 &
             &           )/(4.0d0*dz)                             &
             &       )                                            &
             &   - beta/g_sqrt(i,k)*c_bs2(k)*dts                  &
             &     /(cp*dens_bs(i,k)*vptemp_bs(i,k)**2.0d0)       &
             &     *(                                             &
             &       (dens_bs(i,k+1) + dens_bs(i,k))/2.0d0        &
             &       *(vptemp_bs(i,k+1) + vptemp_bs(i,k))/2.0d0   &
             &       *(                                           &
             &         w(i,k+1)                                   &
             &         - cp*dts/g_sqrt_z(i,k+1)                   &
             &           *(vptemp_bs(i,k+1) + vptemp_bs(i,k))     &
             &           /2.0d0                                   &
             &           *(                                       &
             &             - alp*div_z_w(i,k+1)                   &
             &             + (1.0d0 - beta)                       &
             &               *(pi(i,k+1) - pi(i,k))/dz            &
             &            )                                       &
             &         - (                                        &
             &            - w_adv(i,k+1) + w_buoy(i,k+1)          &
             &            + w_trb(i,k+1)                          &
             &           )*dts                                    &
             &        )                                           &
             &       - (dens_bs(i,k) + dens_bs(i,k-1))/2.0d0      &
             &         *(vptemp_bs(i,k) + vptemp_bs(i,k-1))/2.0d0 &
             &         *(                                         &
             &           w(i,k)                                   &
             &           - cp*dts/g_sqrt_z(i,k)                   &
             &             *(vptemp_bs(i,k) + vptemp_bs(i,k-1))   &
             &             /2.0d0                                 &
             &             *(                                     &
             &               - alp*div_z_w(i,k-1)                 &
             &               + (1.0d0 - beta)                     &
             &                 *(pi(i,k) - pi(i,k-1))/dz          &
             &              )                                     &
             &           - (                                      &
             &              - w_adv(i,k) + w_buoy(i,k)            &
             &              + w_trb(i,k)                          &
             &             )*dts                                  &
             &           )                                        &
             &       )/dz
     end do

     b(km-2) = pi(i,km-2)                                      &
          &   - (1.0d0 - beta)/g_sqrt(i,km-2)*c_bs2(km-2)*dts  &
          &     /(cp*dens_bs(i,km-2)*vptemp_bs(i,km-2)**2.0d0) &
          &     *(                                             &
          &       (dens_bs(i,km-1) + dens_bs(i,km-2))/2.0d0    &
          &       *(vptemp_bs(i,km-1) + vptemp_bs(i,km-2))     &
          &       /2.0d0                                       &
          &       *w(i,km-1)                                   &
          &       - (dens_bs(i,km-2) + dens_bs(i,km-3))/2.0d0  &
          &         *(vptemp_bs(i,km-2) + vptemp_bs(i,km-3))   &
          &         /2.0d0                                     &
          &         *w(i,km-2)                                 &
          &      )/dz                                          &
          &    - c_bs2(km-2)*dts/(cp*vptemp_bs(i,km-2))        &
          &      *(                                            &
          &        (u(i+1,km-2) - u(i,km-2))/dx                &
          &        + g13xz(i,km-2)                             &
          &          *(                                        &
          &            u(i+1,km-1) + u(i,km-1)                 &
          &            - u(i+1,km-3) - u(i,km-3)               &
          &           )/(4.0d0*dz)                             &
          &       )                                            &
          &   - beta/g_sqrt(i,km-2)*c_bs2(km-2)*dts            &
          &     /(cp*dens_bs(i,km-2)*vptemp_bs(i,km-2)**2.0d0) &
          &     *(                                             &
          &       (dens_bs(i,km-1) + dens_bs(i,km-2))/2.0d0    &
          &       *(vptemp_bs(i,km-1) + vptemp_bs(i,km-2))     &
          &       /2.0d0                                       &
          &       *(                                           &
          &         w(i,km-1)                                  &
          &         - cp*dts/g_sqrt_z(i,km-1)                  &
          &           *(vptemp_bs(i,km-1) + vptemp_bs(i,km-2)) &
          &           /2.0d0                                   &
          &           *(                                       &
          &             - alp*div_z_w(i,km-1)                  &
          &             + (1.0d0 - beta)                       &
          &               *(pi(i,km-1) - pi(i,km-2))/dz        &
          &            )                                       &
          &         - (                                        &
          &            - w_adv(i,km-1) + w_buoy(i,km-1)        &
          &            + w_trb(i,km-1)                         &
          &           )*dts                                    &
          &         )                                          &
          &       - (dens_bs(i,km-2) + dens_bs(i,km-3))/2.0d0  &
          &         *(vptemp_bs(i,km-2) + vptemp_bs(i,km-3))   &
          &         /2.0d0                                     &
          &         *(                                         &
          &           w(i,km-2)                                &
          &           - cp*dts/g_sqrt_z(i,km-2)                &
          &             *(vptemp_bs(i,km-2) + vptemp_bs(i,km-3))&
          &             /2.0d0                                 &
          &             *(                                     &
          &               - alp*div_z_w(i,km-2)                &
          &               + (1.0d0 - beta)                     &
          &                 *(pi(i,km-2) - pi(i,km-3))/dz      &
          &              )                                     &
          &           - (                                      &
          &              - w_adv(i,km-2) + w_buoy(i,km-2)      &
          &              + w_trb(i,km-2)                       &
          &             )*dts                                  &
          &          )                                         &
          &         )/dz                                       &
          &   + beta/g_sqrt(i,km-2)*c_bs2(km-2)*dts**2.0d0     &
          &     /(dens_bs(i,km-2)*vptemp_bs(i,km-2)**2.0d0*dz) &
          &     *cp                                            &
          &     *(vptemp_bs(i,km-1)+vptemp_bs(i,k-2)/2.0d0)    &
          &      **2.0d0/g_sqrt_z(i,km-1)                      &
          &     *(                                             &
          &       alp*div_z_w(i,km-1)                          &
          &       - (1.0d0 - beta)*(pi(i,km-1) - pi(i,km-2))/dz&
          &       + g_sqrt_z(i,km-1)                           &
          &         /cp                                        &
          &         /(vptemp_bs(i,km-1)+vptemp_bs(i,k-2)/2.0d0)&
          &         *(                                         &
          &           - w_adv(i,km-1) + w_buoy(i,km-1)         &
          &           + w_trb(i,km-1)                          &
          &          )                                         &
          &       )

     !--- $BL5<!8505NO(B
     call dltx(sbd,d,spd,km-1,b,0.0d0,isw,is,ip,vw,icon)

     do k = 0,km-2
        pi(i,k) = b(k)
        !--- $B%(%i!<%A%'%C%/(B
        if (icon /= 0) then
           write(*,*) 'its =', its,', i =',i,', k =',k,        &
                &     ', Condition code:', icon
           stop
           !     else
           !        write(*,*) 'its =', its,', i =',i,', k =',k,        &
           !             &     ', Condition code:', icon
        end if
     end do

  end do

end subroutine non_dim_pi

