All my calculations involve quantities like < a'b' > where a' is the fluctuation from some temporal mean < a > and the < > symbol is the time-averaging operator. To compute such statistics, I first interpolate all ROMS fields to a new vertical grid that is fixed in time to avoid averaging over different depths at the same (lon, lat) because ROMS has a time-varying vertical grid. For some of these statistics, I often need to compute horizontal gradients, i.e., d/dx or d/dy. What is the best way to do this?
Method 1: Compute d/dx and d/dy on the original ROMS grid using pm and pn. Then interpolate the gradients onto the new vertical grid to enable spatial/time averaging at the same depth. Presumably, the difference in vertical locations for the same k index between two neighboring horizontal locations is small enough to justify this approach.
Method 2: Interpolate all ROMS fields to the time-invariant vertical grid first and then compute the horizontal gradients. In this approach, can I use the same pm and pn with the new vertical grid? Or do I need to recompute pm, pn after the vertical interpolation?