Depending on the version of GDAL, there are a few different resample options available; see gdalwarp.
GDAL 1.10 or later using
average resampling, computes the weighted average of all non-NODATA contributing pixels
This isn't tested, but should look something like:
gdalwarp -t_srs EPSG:4326 -tr 0.5 0.66 -r average fine_one_sq_km.tif coarse_average.tif
Then to get the sum, multiply the average by the number of pixels of the fine resolution raster in one pixel of the coarse resolution raster, which hopefully is constant (you could assume it is).
GDAL 3.1 or later using
compute the weighted sum of all non-NODATA contributing pixels
This should look like this:
gdalwarp -t_srs EPSG:4326 -tr 0.5 0.66 -r sum fine_one_sq_km.tif coarse_sum.tif
sum for older versions) can be used to aggregate multidimensional sums. But this may rely on perfect matchings between grids.