Interpolation of bed stress for step_floats

[gli353]

Hi,
I'm currently modifying ROMS code for the floats time stepping. We want to allow particles to resuspend (rather than either rebound or be permanently stuck), according to a critical shear stress criterion. I know that I may use the interp_floats subroutine that comes with ROMS, but it may require modify the DRIFTER type a bit more, to add new fields etc., which for now I'm reluctant to do. So I wrote a quick-and-dirty interpolation in step_floats.F, as shown below:

But it seems like the bed stress computation is causing me troubles. I used ana_diag.h to print out the bed stress (and some other information not included in the standard float output) experienced by each particle in the model for each baroclinic time step (3s in my case), and found that:
(i) for many particles, the computed stress is often zero, though they are in the estuarine channel rather than at somewhere high and dry;
(ii) in some instances, I saw very large values of bstress , up to 10^5.

I know that the proper interpolation should involve a 2 more grid points than shown here (bilinear), yet I don't think this is responsible for the weird results I'm looking at. Is the coding correct? Am I using the correct mask (which I'm pretty sure...)? or, can here be something really wrong with the bed stress computed by ROMS, and has anyone seen these abnormal bstress values? Any suggestion will be appreciated!

The modified step_floats.F is now also attached. The modifications are in Line 1023 ~ Line 1118. But the only abnormality I have seen so far is related to the bed stress. In the file, theta_cr is set to 0 for debugging purpose, but this is not what we intend.

Kind Regards,
Gaoyang

[jcwarner]

i think 0 is at the lower left psi point. maybe Hernan is better to answer.
if so,then int(track(ixgrd)) would give the left side of the box at a u-point, and int(track(ixgrd))+1 would be the right side u-point.

but you have
i1 = MIN(INT (track(ixgrd,nfp1,1)+0.5_r8) ,Lm(ng)+1)
i2 = MAX(i1 + 1,0)
so if the track x-pos is 12.7, you would get
i1 =int(12.7+.5) = int(13.2)=13
and i2=14.
that seems off.

[gli353]

i think 0 is at the lower left psi point. maybe Hernan is better to answer.
if so,then int(track(ixgrd)) would give the left side of the box at a u-point, and int(track(ixgrd))+1 would be the right side u-point.

but you have
i1 = MIN(INT (track(ixgrd,nfp1,1)+0.5_r8) ,Lm(ng)+1)
i2 = MAX(i1 + 1,0)
so if the track x-pos is 12.7, you would get
i1 =int(12.7+.5) = int(13.2)=13
and i2=14.
that seems off.

Thanks for your reply. So the grid numbers associated with floats are for Psi points? I thought they are for the rho points, so it may work as shown in the diagram:

2.5 is the location of the u-point in between the rho points indexed as 2 and 3. So for a float at 2.7, MIN(INT (track(ixgrd,nfp1,1)+0.5_r8) ,Lm(ng)+1) gives 3, which is also the index of the u-point to its left (shown in the blue circle), and i2 = 4 is to its right, and I use these two values for interpolation.

I think I may need a bit more clarification on the indexing of floats?

[arango]

I don't have time to look at this now. However, as a guideline, the indexing in the I-direction for the U- and PSI-points are identical. Similarly, the indexing in the J-direction for the V- and PSI-points are identical. Thus, if you look at interpolating at PSI-points are a combination of both U and V grids methods.

You may look at the plot for the fractional coordinate system (XI, ETA), which is different in 4D-Var.

[gli353]

Here must be something wrong with my interpolation... But I've adopted another approach now; I just added two more fields to the track array in DRIFTER to store the bed stress, interpolated using interp_float subroutine that comes with the ROMS code. Preliminarily, the results seem to make sense to me.
If I have any new problems, I will ask again. Thanks for all the suggestions

[wilkin]

Another approach you might consider is offline particle tracking using ROMSPath https://gmd.copernicus.org/articles/15/4297/2022/

ROMSpath uses almost exactly the same floats algorithm internal to ROMS. Importantly, the velocity is expressed in ROMS native xi,eta,s coordinates using u,v,omega, so it does not incur the errors that other offline particle algorithms do with ROMS output by interpolating and averaging in geographic coordinates. The only real difference is a 4th-order Runge-Kutta time step instead of the 3rd-order in ROMS floats.

ROMSpath was conceived to make it easy to add behavior to floats that represent larvae, but in your case you could add behavior for sediment. An environmental cue we explored with ROMSpath simulations was the effect of turbulence on larval settlement, so many of the hooks you'd need are already there.

Plus, it supports floats traversing nested grids in both directions.

[gli353]

Thanks for pointing out in this direction, we are considering it as well. However, the major concern, for which I have not found an answer yet, is that tracking particles in 3D with turbulence would require me to output u,v,w,akt,salt,temp (the last two for computing density; if I tweak the code and output density itself, it would still be 5 variables) at all depth levels, and at a reasonably small time step like 5 minutes, that will mean 2.5~3 T of output in total, for one scenario. The storage we are granted by the HPC centre is a constraint at the moment... For now, my code seems to run well, but ROMSpath is something we are very curious about indeed, and may manage to run a simulation it and compare.

By the way, I would really appreciate it if you could share the link to the larval settlement research, if it has been published. Thanks!

Another approach you might consider is offline particle tracking using ROMSPath https://gmd.copernicus.org/articles/15/4297/2022/

ROMSpath uses almost exactly the same floats algorithm internal to ROMS. Importantly, the velocity is expressed in ROMS native xi,eta,s coordinates using u,v,omega, so it does not incur the errors that other offline particle algorithms do with ROMS output by interpolating and averaging in geographic coordinates. The only real difference is a 4th-order Runge-Kutta time step instead of the 3rd-order in ROMS floats.

ROMSpath was conceived to make it easy to add behavior to floats that represent larvae, but in your case you could add behavior for sediment. An environmental cue we explored with ROMSpath simulations was the effect of turbulence on larval settlement, so many of the hooks you'd need are already there.

Plus, it supports floats traversing nested grids in both directions.

[wilkin]

A paper describing larval transport results using ROMSpath in a nested grid configuration is here:
Garwood, J.C., Fuchs, H.L., Gerbi, G.P., Hunter, E.J., Chant, R.J. and Wilkin, J.L., 2022. Estuarine retention of larvae: Contrasting effects of behavioral responses to turbulence and waves. Limnology and Oceanography, 67(4), pp.992-1005. https://aslopubs.onlinelibrary.wiley.co ... /lno.12052

The ROMS outputs that the article used are here:
https://www.seanoe.org/data/00833/94520/

As a general message to the ROMS user community reading this forum, I highly recommend SEANOE https://www.seanoe.org as a repository for serving your model outputs with a citable DOI. They are very quick to respond to submissions, and they accept not just output file uploads but also the option to give a dataset description that points to web services for a full data set -- when the run output is too voluminous for it to be practical for SEANOE to duplicate. See e.g. https://www.seanoe.org/data/00785/89673/.

SEANOE allow you to later edit the data set description (without altering the DOI) to note new papers that use the data, or if a newer version exists to superceed the original data.

They also send notifications when the dataset DOI is cited in publications ... great at promotion time so you can prove folks you don't know actually used your data!

[gli353]

A paper describing larval transport results using ROMSpath in a nested grid configuration is here:
Garwood, J.C., Fuchs, H.L., Gerbi, G.P., Hunter, E.J., Chant, R.J. and Wilkin, J.L., 2022. Estuarine retention of larvae: Contrasting effects of behavioral responses to turbulence and waves. Limnology and Oceanography, 67(4), pp.992-1005. https://aslopubs.onlinelibrary.wiley.co ... /lno.12052

The ROMS outputs that the article used are here:
https://www.seanoe.org/data/00833/94520/

As a general message to the ROMS user community reading this forum, I highly recommend SEANOE https://www.seanoe.org as a repository for serving your model outputs with a citable DOI. They are very quick to respond to submissions, and they accept not just output file uploads but also the option to give a dataset description that points to web services for a full data set -- when the run output is too voluminous for it to be practical for SEANOE to duplicate. See e.g. https://www.seanoe.org/data/00785/89673/.

SEANOE allow you to later edit the data set description (without altering the DOI) to note new papers that use the data, or if a newer version exists to superceed the original data.

They also send notifications when the dataset DOI is cited in publications ... great at promotion time so you can prove folks you don't know actually used your data!

Thanks for pointing out the existence of SEANOE! I think we will upload some of our results and code to it. Now resuspension looks fine

Gaoyang