FUNCTION GRID_COMBINE, A, B

;- Composite grid b with grid a (e.g. to combine days)

;- Check arguments
if (n_elements(a) eq 0) then message, 'Argument A is undefined'
if (n_elements(b) eq 0) then message, 'Argument B is undefined'

;- Check grid array dimensions
dims = size(a, /dimensions)
if (dims[0] ne 4L) or (dims[1] ne 814880L) then $
  message, 'Dimensions of A must be [4, 814880, N]'
dims = size(b, /dimensions)
if (dims[0] ne 4L) or (dims[1] ne 814880L) then $
  message, 'Dimensions of B must be [4, 814880, N]'
if (n_elements(a) ne n_elements(b)) then $
  message, 'A and B must have the same number of elements'

;- Save start time
t0 = systime(1)

;- Create the output array
output = a
  
;- Composite grid b with grid a
for dim = 0L, dims[2] - 1L do begin
  index = where(b[0, *, dim] gt 0.0, count)
  if (count gt 0) then begin
    output[0, index, dim] = a[0, index, dim] + b[0, index, dim]
    output[1, index, dim] = a[1, index, dim] + b[1, index, dim]
    output[2, index, dim] = a[2, index, dim] < b[2, index, dim]
    output[3, index, dim] = a[3, index, dim] > b[3, index, dim]
  endif
endfor

;- Print elapsed time
t1 = systime(1)
print, 'Elapsed time (sec): ', t1 - t0

;- Return the output array
return, output

END