Self-consistently modeling baryonic effects in survey-scale cosmological simulations has become increasingly important as the diversity, precision, and statistical reach of modern observations continue to improve. The advent of exascale computing now enables a new generation of simulations that couple these physical processes across full-sky volumes with excellent statistical sampling of large-scale structure tracers such as galaxies, groups, and clusters. To support these efforts, we extend the CRK-HACC framework, a GPU-accelerated cosmological hydrodynamics code, with a suite of astrophysical subgrid models that simulate radiative cooling, star formation, stellar evolution, and AGN feedback within a numerically robust formulation optimized for scalability on modern exascale architectures. The models were selected and calibrated to reproduce observed galaxy stellar mass functions over the redshift range $0 < z < 2$ and cluster populations probed by cosmological surveys, capturing the large-scale baryonic evolution relevant for multi-wavelength, cross-correlated analyses. We describe the implementation and calibration of these models and demonstrate their consistency with observed galaxy population statistics and modern hydrodynamic simulations, establishing the baseline for exascale efforts that extend this framework to survey-scale volumes.