Regridding two different curvilinear grids onto a common rectilinear grid with R

Question

I have numerical model outputs (in this example, ocean surface temperature) from two different models that cover different spatial domains. Both sets of outputs are on their own curvilinear grid (which has been a learning curve!). The larger domain model extent incorporates that of the smaller domain model. I am trying to create a raster (rectilinear) grid where the outputs (data) from the two models are seamlessly combined, so that the data from the smaller domain model (which also has a higher native resolution) replaces/over-writes the data from the larger domain (and lower resolution) model. I think I have managed to do this with a convoluted approach, using the stars R package and the final step accomplished with gdalwarp in a bash terminal (see screenshot of the final raster loaded in QGIS). However, my questions are:

1. Can the two curvilinear grids be warped/transformed onto the rectilinear grid within R (without having to warp them to rectilinear grids and saving those to implement the final step with gdalwarp)?

2. Is there a better way to approach this problem? I might have gone down the wrong rabbit hole.

I have provided a GitLab repository with the below R code and the files needed to run the code, which demonstrates my approach. I don't have the skills to build minimal example curvilinear grids and apologies for the lengthy script.

I should acknowledge code and (limited) understanding from pages that include:
https://r-spatial.org/book/07-Introsf.html
https://r-spatial.github.io/stars/articles/stars4.html#curvilinear-grids
https://www.mattreusswig.com/post/use-stars-to-visualize-curvilinear-netcdf-rasters/
https://gis.stackexchange.com/questions/390148/handle-curvilinear-rotated-grid-netcdf-file-in-r

#####################################################################################
##
## Purpose of script:
## develop a reproducible example of trying to add data from a larger-domain woesII
## (1/12th degree resolution) curvilinear grid, as well as a smaller-domain woesIII
## (1/36th degree resolution) curvilinear grid to a common (and higher resolution)
## rectilinear grid, so that the woesIII data replace the woesII data where the former
## exist (and the data from the two original grids are seemless in the output).
##
#####################################################################################
### packages & functions

require(stars)

#####################################################################################
### settings

### filenames & directories (note files are read from and saved to an 'outputs' folder, you may need to replace throughout if you did not create the same directory)

grd_fl <- "outputs/receiving_grid_2km.tif" 

wIII_fl <- "outputs/temp_woesIII_singleday.nc"

wII_fl <- "outputs/temp_woesII_singleday.nc"

varname <- "temp"

#####################################################################################
### load data

### load the new grid file
grd <- read_stars(grd_fl)
grd

### load the two (curvilinear) model output grids:
## woesII data
wII <- stars::read_ncdf(wII_fl, var = varname)
wII # not recognised as a curvilinear grid (unlike when I opened the original model files), so incorporate the below to tell it:  
coords_wII  <- stars::read_ncdf(wII_fl, var = c("lon_rho", "lat_rho")) #
## Convert the lat lon coordinates into matrices with the same dimensions as the raster data.
x_wII = matrix(coords_wII[[1]], nrow(coords_wII[[1]]), ncol(coords_wII[[1]]))
y_wII = matrix(coords_wII[[2]], nrow(coords_wII[[2]]), ncol(coords_wII[[2]]))
## tell stars this is a curvilinear grid with coords from x & y
wII <- st_as_stars(wII, curvilinear = list(longitude = x_wII, latitude = y_wII))
wII # now it recognises it as a curvilinear grid

## woesIII data
wIII <- stars::read_ncdf(wIII_fl, var = varname)
wIII # not recognised as a curvilinear grid (unlike when I opened the original model files), so incorporate the below to tell it:
coords_wIII  <- stars::read_ncdf(wIII_fl, var = c("lon_rho", "lat_rho")) #
## Convert the lat lon coordinates into matrices with the same dimensions as the raster data.
x_wIII = matrix(coords_wIII[[1]], nrow(coords_wIII[[1]]), ncol(coords_wIII[[1]]))
y_wIII = matrix(coords_wIII[[2]], nrow(coords_wIII[[2]]), ncol(coords_wIII[[2]]))
## tell stars this is a curvilinear grid with coords from x & y
wIII <- st_as_stars(wIII, curvilinear = list(longitude = x_wIII, latitude = y_wIII))
wIII # now it recognises it as a curvilinear grid

## in order to recognise whether the below methods are working, I'm temporarily going to alter the data of wIII (so that clear discontinuities will be visible where the two datasets meet in the final product, avoiding any uncertainty of whether the output is a combination of the two datasets).
wIII <- wIII*1.5
wIII # 50% higher temperatures

#####################################################################################
### try adding them to the new grd

newgrd <- st_warp(wII, grd)
image(newgrd) # seems to have worked to add the larger extent woesII grid
newgrd <- st_warp(wIII, newgrd) # also tried with , options="overwrite=F"
image(newgrd) # inserts the wIII grid but doesn't keep the surrounding data from wII, replacing those with NA
## is there a way to insert the wIII data without losing the wII data (in areas outside the extent of the wIII data)?

### trying with gdal_utils; note this was guessing; I'm not familiar with gdal_utils()
aa <- gdal_utils(util = "warp", source = c(wII, wIII), destination = grd) # I tried various 'along' arguments, but didn't manage to get it right

#####################################################################################
### I did seemingly manage to warp the two grids onto one higher resolution grid with 'gdalwarp' (implemented from a bash terminal as shown below); but to do that, the two curvilinear grids need to be saved as rectilinear raster grids first, as far as I can make out:
### so, warping to rectilinear grids and saving those outputs as follows:

## create rectilinear grids with same number of gridcells as the curvilinear grids
grdII <- st_bbox(wII) |>
  st_as_sfc() |>
  st_transform('OGC:CRS84') |>
  st_bbox() |>
  st_as_stars(nx = dim(wII)["longitude"], ny = dim(wII)["latitude"])

grdIII <- st_bbox(wIII) |>
  st_as_sfc() |>
  st_transform('OGC:CRS84') |>
  st_bbox() |>
  st_as_stars(nx = dim(wIII)["longitude"], ny = dim(wIII)["latitude"])

## and warp the data from the curvilinear grids onto the rectilinear grids
wII_rectilinear <- st_warp(wII, grdII)
wIII_rectilinear <- st_warp(wIII, grdIII)

## save the files:

write_stars(wII_rectilinear, "outputs/temp_woesII_singleday_rectilinear.tif")

write_stars(wIII_rectilinear, "outputs/temp_woesIIIx1.5_singleday_rectilinear.tif")

and then running the below two lines in a bash terminal (will throw errors in R) after navigating to the outputs folder

cp receiving_grid_2km.tif woesII_2km_grid.tif # copy the receiving grid as it gets over-written in the next line  
gdalwarp temp_woesII_singleday_rectilinear.tif temp_woesIIIx1.5_singleday_rectilinear.tif woesII_2km_grid.tif