There are several 3d values that tend to appear in the lateral boundary file: horizontal winds (u and v), the moist perturbation potential temperature, mixing ratio scalars (vapor, cloud, snow, ice, etc), and geopotential.
Each of these 3d fields is "coupled" or "mass weighted" with the 3d value of d(Pd)/d(eta). These d(Pd)/d(eta) values range from 0 to about 10^5, and are effectively used as scale factors. They DRAMATICALLY change the interpretation of the value of the original field.
Pd(i,k,j) = C3(k) * ( Pds(i,j) - Pt ) + C4(k) + Pt.
The eta array is the full-level values of the vertical coordinate (can be specified in the namelist, or can be computed automatically).
Pt is the model lid (Pa), and defined in the namelist.
The 1d constants (C3 and C4) are computed internally.
The dry surface pressure is Pds.
Then, d(Pd(i,k,j))/d(eta(k)) = C1(k) * ( Pcb(i,j) + Pc(i,j) ) + C2(k).
Pcb is the base column pressure (reference), and Pc is the perturbation column pressure (in the WRF model, these are the 2d mub and mu fields). Care must be taken to use the correct horizontal stagger for mub and mu for the different variables for momentum cell-face staggering vs mass point staggering. Also the internally computed and readily available half-layer vs full-level values of the 1d arrays must match the expected field (potential temperature vs geopotential, for example).
So, the boundary values (the BXE, BYE, BXS, BYS arrays) for the moist fields would be qv(i,k,j)*C1(k) * ( mub(i,j) + mu(i,j) ) + C2(k).
For the suffix:
B = Boundary
X/Y = X or Y lateral boundary
S/E = Start or End
The fields that include the "T" in the suffix are the tendencies of the lateral boundaries. If the data from ungrib is 3-hourly (then metgrid is 3-hourly, then the lateral BC data from the real program is 3-hourly).
The tendency array for the 3d array alpha would be alpha_tend = ( alpha(at hour t + 3h) - alpha(at hour t) ) / 10800 s.