\section{\src{dyn_vi_rhow_solver}}

\subsection{Description}

Kernel \src{dyn_vi_rhow_solver} is taken from the original subroutine
\src{vi_rhow_solver} in \NICAM.
%
This subrouine is originally defined in \src{mod_vi}.
%
Subroutine \src{vi_rhow_solver} is to solve the tridiagonal matrix
equations related to the vertical implicit scheme. See
\autoref{s:vert_coord} for the detail of this calculation.


\subsection{Discretization and code}

Argument lists and local variables definition part of this subroutine is
as follows.

\begin{LstF90}[name=vi_rhow_solver]
!-----------------------------------------------------------------------------
!> Tridiagonal matrix solver
subroutine vi_rhow_solver( &
     rhogw,  rhogw_pl,  &
     rhogw0, rhogw0_pl, &
     preg0,  preg0_pl,  &
     rhog0,  rhog0_pl,  &
     Srho,   Srho_pl,   &
     Sw,     Sw_pl,     &
     Spre,   Spre_pl,   &
     dt                 )
  !$ use omp_lib
  implicit none

  real(RP), intent(inout) :: rhogw    (ADM_gall   ,ADM_kall,ADM_lall   ) ! rho*w          ( G^1/2 x gam2 ), n+1
  real(RP), intent(inout) :: rhogw_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)

  real(RP), intent(in)    :: rhogw0   (ADM_gall   ,ADM_kall,ADM_lall   ) ! rho*w          ( G^1/2 x gam2 )
  real(RP), intent(in)    :: rhogw0_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), intent(in)    :: preg0    (ADM_gall   ,ADM_kall,ADM_lall   ) ! pressure prime ( G^1/2 x gam2 )
  real(RP), intent(in)    :: preg0_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), intent(in)    :: rhog0    (ADM_gall   ,ADM_kall,ADM_lall   ) ! rho            ( G^1/2 x gam2 )
  real(RP), intent(in)    :: rhog0_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), intent(in)    :: Srho     (ADM_gall   ,ADM_kall,ADM_lall   ) ! source term for rho  at the full level
  real(RP), intent(in)    :: Srho_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), intent(in)    :: Sw       (ADM_gall   ,ADM_kall,ADM_lall   ) ! source term for rhow at the half level
  real(RP), intent(in)    :: Sw_pl    (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), intent(in)    :: Spre     (ADM_gall   ,ADM_kall,ADM_lall   ) ! source term for pres at the full level
  real(RP), intent(in)    :: Spre_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), intent(in)    :: dt

  real(RP) :: Sall    (ADM_gall,   ADM_kall)
  real(RP) :: Sall_pl (ADM_gall_pl,ADM_kall)
  real(RP) :: beta    (ADM_gall   )
  real(RP) :: beta_pl (ADM_gall_pl)
  real(RP) :: gamma   (ADM_gall,   ADM_kall)
  real(RP) :: gamma_pl(ADM_gall_pl,ADM_kall)

  integer  :: gall, kmin, kmax, lall
  real(RP) :: grav
  real(RP) :: CVovRt2 ! Cv / R / dt**2
  real(RP) :: alpha

  integer  :: g, k, l
  integer  :: gstr, gend
  !$ integer  :: n_per_thread
  !$ integer  :: n_thread
  !---------------------------------------------------------------------------

\end{LstF90}

Here \src{rhogw} is $\rho \times w$ with metric $G^{1/2} \gamma^2$
multiplied at new time step $n+1$.
%
\src{rhogw0}, \src{preg0}, \src{rhog0}, are 
$\rho \times w$, pressure, $\rho$
with metric multiplied at time step $n$, respectively.
%
\src{Srho}, \src{Sw}, \src{Spre} are
source term for $\rho$ at the full level,
source term for $\rho \times w$ at the half level,
source term for pressure at the full level, respectively.
%
Other arguments with suffix \src{_pl} are for the pole region.
%
\src{dt} is a time step for fast-mode.

Among local variables,
\src{alpha} is the flag for non-hydrostatic/hydrostatic.
In this kernel, set to $1$ in \src{problem_size.inc}.




Main part of this subroutine is as follows.

\begin{LstF90}[name=vi_rhow_solver,firstnumber=last]
  call DEBUG_rapstart('____vi_rhow_solver')

  gall = ADM_gall
  kmin = ADM_kmin
  kmax = ADM_kmax
  lall = ADM_lall

  grav    = CONST_GRAV
  CVovRt2 = CONST_CVdry / CONST_Rdry / (dt*dt)
  alpha   = real(NON_HYDRO_ALPHA,kind=RP)

  !$omp parallel default(none),private(g,k,l), &
  !$omp private(gstr,gend,n_thread,n_per_thread) &
  !$omp shared(gall,kmin,kmax,lall,rhogw,rhogw0,preg0,rhog0,Srho,Sw,Spre,dt,Sall,beta,gamma,Mu,Mc,Ml, &
  !$omp        GRD_afact,GRD_bfact,GRD_rdgzh,VMTR_GSGAM2H,VMTR_RGAM,VMTR_RGAMH,VMTR_RGSGAM2,VMTR_RGSGAM2H,grav,alpha,CVovRt2)
  gstr = 1
  gend = gall
  !$ n_thread     = omp_get_num_threads()
  !$ n_per_thread = gall / n_thread + int( 0.5_RP + sign(0.5_RP,mod(gall,n_thread)-0.5_RP) )
  !$ gstr         = n_per_thread * omp_get_thread_num() + 1
  !$ gend         = min( gstr+n_per_thread-1, gall )

  do l = 1, lall
     ! calc Sall
     do k = kmin+1, kmax
     do g = gstr, gend
        Sall(g,k) = (   ( rhogw0(g,k,  l)*alpha + dt * Sw  (g,k,  l) ) * VMTR_RGAMH  (g,k,  l)**2             &
                    - ( ( preg0 (g,k,  l)       + dt * Spre(g,k,  l) ) * VMTR_RGSGAM2(g,k,  l)                &
                      - ( preg0 (g,k-1,l)       + dt * Spre(g,k-1,l) ) * VMTR_RGSGAM2(g,k-1,l)                &
                      ) * dt * GRD_rdgzh(k)                                                                   &
                    - ( ( rhog0 (g,k,  l)       + dt * Srho(g,k,  l) ) * VMTR_RGAM(g,k,  l)**2 * GRD_afact(k) &
                      + ( rhog0 (g,k-1,l)       + dt * Srho(g,k-1,l) ) * VMTR_RGAM(g,k-1,l)**2 * GRD_bfact(k) &
                      ) * dt * grav                                                                           &
                    ) * CVovRt2
     enddo
     enddo

     ! boundary conditions
     do g = gstr, gend
        rhogw(g,kmin,  l) = rhogw(g,kmin,  l) * VMTR_RGSGAM2H(g,kmin,  l)
        rhogw(g,kmax+1,l) = rhogw(g,kmax+1,l) * VMTR_RGSGAM2H(g,kmax+1,l)
        Sall (g,kmin+1)   = Sall (g,kmin+1) - Ml(g,kmin+1,l) * rhogw(g,kmin,  l)
        Sall (g,kmax  )   = Sall (g,kmax  ) - Mu(g,kmax,  l) * rhogw(g,kmax+1,l)
     enddo

     !---< solve tri-daigonal matrix >

     ! condition at kmin+1
     k = kmin+1
     do g = gstr, gend
        beta (g)     = Mc(g,k,l)
        rhogw(g,k,l) = Sall(g,k) / beta(g)
     enddo

     ! forward
     do k = kmin+2, kmax
     do g = gstr, gend
        gamma(g,k)   = Mu(g,k-1,l) / beta(g)
        beta (g)     = Mc(g,k,l) - Ml(g,k,l) * gamma(g,k) ! update beta
        rhogw(g,k,l) = ( Sall(g,k) - Ml(g,k,l) * rhogw(g,k-1,l) ) / beta(g)
     enddo
     enddo

     ! backward
     do k = kmax-1, kmin+1, -1
     do g = gstr, gend
        rhogw(g,k  ,l) = rhogw(g,k  ,l) - gamma(g,k+1) * rhogw(g,k+1,l)
        rhogw(g,k+1,l) = rhogw(g,k+1,l) * VMTR_GSGAM2H(g,k+1,l) ! return value ( G^1/2 x gam2 )
     enddo
     enddo

     ! boundary treatment
     do g = gstr, gend
        rhogw(g,kmin  ,l) = rhogw(g,kmin  ,l) * VMTR_GSGAM2H(g,kmin  ,l)
        rhogw(g,kmin+1,l) = rhogw(g,kmin+1,l) * VMTR_GSGAM2H(g,kmin+1,l)
        rhogw(g,kmax+1,l) = rhogw(g,kmax+1,l) * VMTR_GSGAM2H(g,kmax+1,l)
     enddo
  enddo
  !$omp end parallel
\end{LstF90}

Inside of the long $l$-loop(l.72), the first $k$-loop(l.73) calcurates total source
term \src{Sall}.
%
The section after setting the boundary condition at \src{kmin} and \src{kmax}
solves tri-diagonal matrix.
%
Its first part(l.97 to l.111) is for forward elimination and the second
part(l.114 to l.119) is for backward substitution.

Note that elements of the tridiagonal matrix \src{Mc}, \src{Mu} and
\src{Ml} are calculated in advance by other subroutine in the original
module, and they are read from input data file in this kernel program.


The last part of this subroutine is for the pole region, doing almost
the same process as the normal region.

\begin{LstF90}[name=vi_rhow_solver,firstnumber=last,breaklines=false]
  if ( ADM_have_pl ) then
     do l = 1, ADM_lall_pl
        do k  = ADM_kmin+1, ADM_kmax
        do g = 1, ADM_gall_pl
           Sall_pl(g,k) = (   ( rhogw0_pl(g,k,  l)*alpha + dt * Sw_pl  (g,k,  l) ) * VMTR_RGAMH_pl  (g,k,  l)**2             &
                          - ( ( preg0_pl (g,k,  l)       + dt * Spre_pl(g,k,  l) ) * VMTR_RGSGAM2_pl(g,k,  l)                &
                            - ( preg0_pl (g,k-1,l)       + dt * Spre_pl(g,k-1,l) ) * VMTR_RGSGAM2_pl(g,k-1,l)                &
                            ) * dt * GRD_rdgzh(k)                                                                            &
                          - ( ( rhog0_pl (g,k,  l)       + dt * Srho_pl(g,k,  l) ) * VMTR_RGAM_pl(g,k,  l)**2 * GRD_afact(k) &
                            + ( rhog0_pl (g,k-1,l)       + dt * Srho_pl(g,k-1,l) ) * VMTR_RGAM_pl(g,k-1,l)**2 * GRD_bfact(k) &
                            ) * dt * grav                                                                                    &
                          ) * CVovRt2
        enddo
        enddo

        do g = 1, ADM_gall_pl
           rhogw_pl(g,ADM_kmin,  l) = rhogw_pl(g,ADM_kmin,  l) * VMTR_RGSGAM2H_pl(g,ADM_kmin,  l)
           rhogw_pl(g,ADM_kmax+1,l) = rhogw_pl(g,ADM_kmax+1,l) * VMTR_RGSGAM2H_pl(g,ADM_kmax+1,l)
           Sall_pl (g,ADM_kmin+1)   = Sall_pl (g,ADM_kmin+1) - Ml_pl(g,ADM_kmin+1,l) * rhogw_pl(g,ADM_kmin,  l)
           Sall_pl (g,ADM_kmax  )   = Sall_pl (g,ADM_kmax  ) - Mu_pl(g,ADM_kmax,  l) * rhogw_pl(g,ADM_kmax+1,l)
        enddo

        k = ADM_kmin+1
        do g = 1, ADM_gall_pl
           beta_pl (g)     = Mc_pl(g,k,l)
           rhogw_pl(g,k,l) = Sall_pl(g,k) / beta_pl(g)
        enddo

        do k = ADM_kmin+2, ADM_kmax
        do g = 1, ADM_gall_pl
           gamma_pl(g,k)   = Mu_pl(g,k-1,l) / beta_pl(g)
           beta_pl (g)     = Mc_pl(g,k,l) - Ml_pl(g,k,l) * gamma_pl(g,k) ! update beta
           rhogw_pl(g,k,l) = ( Sall_pl(g,k) - Ml_pl(g,k,l) * rhogw_pl(g,k-1,l) ) / beta_pl(g)
        enddo
        enddo

        ! backward
        do k = ADM_kmax-1, ADM_kmin+1, -1
        do g = 1, ADM_gall_pl
           rhogw_pl(g,k  ,l) = rhogw_pl(g,k  ,l) - gamma_pl(g,k+1) * rhogw_pl(g,k+1,l)
           rhogw_pl(g,k+1,l) = rhogw_pl(g,k+1,l) * VMTR_GSGAM2H_pl(g,k+1,l) ! return value ( G^1/2 x gam2 )
        enddo
        enddo

        ! boundary treatment
        do g = 1, ADM_gall_pl
           rhogw_pl(g,ADM_kmin  ,l) = rhogw_pl(g,ADM_kmin  ,l) * VMTR_GSGAM2H_pl(g,ADM_kmin  ,l)
           rhogw_pl(g,ADM_kmin+1,l) = rhogw_pl(g,ADM_kmin+1,l) * VMTR_GSGAM2H_pl(g,ADM_kmin+1,l)
           rhogw_pl(g,ADM_kmax+1,l) = rhogw_pl(g,ADM_kmax+1,l) * VMTR_GSGAM2H_pl(g,ADM_kmax+1,l)
        enddo
     enddo
  endif

  call DEBUG_rapend('____vi_rhow_solver')

  return
end subroutine vi_rhow_solver
\end{LstF90}

\subsection{Input data and result}

Max/min/sum of input/output data of the kernel subroutine are output as
a log.
%
Below is an example of \src{$IAB_SYS=Ubuntu-gnu-ompi} case.


\begin{LstLog}
 ### Input ###
 +check[rhogw_prev      ] max=  1.0733675425174699E-14,min= -2.5212277839542070E-18,sum=  2.1443525022460023E-14
 +check[rhogw_prev_pl   ] max=  1.0733675425174699E-14,min= -3.0929371635479732E-15,sum=  7.6415689931329447E-15
 +check[check_rhogw     ] max=  6.3321830580406713E-01,min= -5.5759247708875415E-01,sum=  3.9071354372140144E+02
 +check[check_rhogw_pl  ] max=  9.9023353298473032E-02,min= -7.4852014128754477E-03,sum=  2.1725659627743767E+00
 +check[rhogw0          ] max=  1.2675708942695349E-03,min= -4.2678569215437437E-03,sum= -2.6525715444058587E-02
 +check[rhogw0_pl       ] max=  9.9481836269239257E-04,min= -4.2678569215437437E-03,sum= -8.3801034534239580E-02
 +check[preg0           ] max=  1.4047194007847271E+01,min= -2.6115802085359952E+01,sum= -1.0687083035942905E+05
 +check[preg0_pl        ] max=  1.3572400586711277E+01,min= -2.6115802085359952E+01,sum=  3.4640462790019097E+02
 +check[Srho            ] max=  3.3317357550311205E-04,min= -3.9691441465327759E-04,sum= -2.7468694177252501E-01
 +check[Srho_pl         ] max=  4.3441154740222974E-05,min= -2.2770470374092324E-05,sum= -3.8685936849254728E-05
 +check[Sw              ] max=  2.4651449651180712E-04,min= -3.2516247664737818E-04,sum= -1.6092196259156019E-01
 +check[Sw_pl           ] max=  2.8812604579719903E-04,min= -5.5442162236740700E-04,sum= -1.1046915056782920E-02
 +check[Spre            ] max=  3.6234075910381584E+01,min= -3.7498614557267203E+01,sum= -2.4763893276893483E+04
 +check[Spre_pl         ] max=  4.4863109100507081E+00,min= -2.0915571645010926E+00,sum=  1.0496606720811164E+00
 +check[Mc              ] max=  1.7037935104756485E+00,min=  0.0000000000000000E+00,sum=  8.2325595392291399E+05
 +check[Mc_pl           ] max=  1.4569727548844720E+00,min=  2.3177911966477871-310,sum=  6.0165678214473996E+02
 +check[Ml              ] max=  0.0000000000000000E+00,min= -8.3019108752152282E-01,sum= -3.9659607793632336E+05
 +check[Ml_pl           ] max=  3.2282782818256430E+02,min= -6.9953788050738130E-01,sum=  3.4604237063825449E+03
 +check[Mu              ] max=  0.0000000000000000E+00,min= -8.6888851023079583E-01,sum= -4.2616616191494197E+05
 +check[Mu_pl           ] max=  9.7233873292775215E+01,min= -7.5264063505825385E-01,sum=  8.5169724995286754E+02
 ### Output ###
 +check[rhogw           ] max=  6.3321830580406713E-01,min= -5.5759247708875415E-01,sum=  3.9071354372140144E+02
 +check[rhogw_pl        ] max=  9.9023353298473032E-02,min= -7.4852014128754477E-03,sum=  2.1725659627743767E+00
 ### Validation : point-by-point diff ###
 +check[check_rhogw     ] max=  0.0000000000000000E+00,min=  0.0000000000000000E+00,sum=  0.0000000000000000E+00
 +check[check_rhogw_pl  ] max=  0.0000000000000000E+00,min=  0.0000000000000000E+00,sum=  0.0000000000000000E+00
 *** Finish kernel
\end{LstLog}

Check the lines below \src{``Validation : point-by-point diff''} line,
that shows difference between calculated output array and
pre-calculated reference array.
These should be zero or enough small to be acceptable.

There are sample output log files in \file{reference/} 
in each kernel program directory, for reference purpose.



\subsection{Sample of perfomance result}

Here's an example of the performance result part of the log output.
Below is an example executed with the machine environment described in \autoref{s:measuring_env}.
%
Note that in this program kernel part is iterated one time.

\begin{LstLog}
 *** Computational Time Report
 *** ID=001 : MAIN_dyn_vi_rhow_solver          T=     0.016 N=      1
 *** ID=002 : ____vi_rhow_solver               T=     0.016 N=      1
\end{LstLog}