I guess know what is going on here. Please look the codes below:
#if 1
! This is the hydrostatic equation used in the model after the small timesteps. In
! the model, grid%al (inverse density) is computed from the geopotential.
IF (grid%hypsometric_opt == 1) THEN
DO kk = 2,kte
k = kk - 1
grid%ph_2(i,kk,j) = grid%ph_2(i,kk-1,j) - &
grid%dnw(kk-1) * ( ((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_2(i,j)))*grid%al(i,kk-1,j) &
+ (grid%c1h(k)*grid%mu_2(i,j))*grid%alb(i,kk-1,j) )
grid%ph0(i,kk,j) = grid%ph_2(i,kk,j) + grid%phb(i,kk,j)
END DO
ELSE IF (grid%hypsometric_opt == 2) THEN
! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
! Note that al*p approximates Rd*T and dLOG(p) does z.
! Here T varies mostly linear with z, the first-order integration produces better result.
grid%ph_2(i,1,j) = grid%phb(i,1,j)
DO k = 2,kte
pfu = grid%c3f(k )*grid%MU0(i,j)+grid%c4f(k ) + grid%p_top
pfd = grid%c3f(k-1)*grid%MU0(i,j)+grid%c4f(k-1) + grid%p_top
phm = grid%c3h(k-1)*grid%MU0(i,j)+grid%c4h(k-1) + grid%p_top
grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
END DO
DO k = 1,kte
grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
END DO
END IF
#else
! Get the perturbation geopotential from the 3d height array from WPS.
DO k = 2,kte
grid%ph_2(i,k,j) = grid%ph0(i,k,j)*g - grid%phb(i,k,j)
END DO
#endif
and change the "#if 1" to "#if 0", then recompile WRF. Please let em know whether this works.