Re write obcs tides

[jm-c]

What changes does this PR introduce?

Implement suggestions from issue #617, renaming OB tidal forcing variables and improving
efficiency of tidal component addition. Also allows to add tides to OB tangential velocity.

What is the current behaviour?

1. not very efficient (see point 1, issue Improve OBC Tidal forcing implementation #617 )
2. not very flexible regarding number of tidal-components to set in binary input files (see point 2, issue Improve OBC Tidal forcing implementation #617 )
3. some confusing names, with "normal/tangential" velocity component used for Tides but "u/v" component used for the remaining OBCS pkg (see point 3, issue Improve OBC Tidal forcing implementation #617 )

What is the new behaviour

Address all 3 points above.
Also enable to specify the tangential tidal velocity at the OB (An other implementation of PR #73).

Does this PR introduce a breaking change?

Names of Tidal forcing parameters have changed. A check and stop has been added to obcs_readparams.F if old-names were found in data.obcs.

Other information:

1. loading and setting of OB tidal forcing has been moved from obcs_init_variables.F to obcs_init_fixed.F since the tidal components amplitude and phase are kept constant during one simulation.
2. point (2) above is now implemented (using updated code from PR Update mdsio slice I/O routines #747 which has now been merged in).
3. the new (more efficient) calculation of tidal flow introduces differences at machine-truncation level ;
   a temporary CPP option "OLD_OBCS_TIDES" was allowing to recover the original results (but got removed).
4. The updated example test seaice_obcs.tides from PR Yoshi Nakayama code for tangential component of tides #73 has been moved here with minor renaming and smaller binary files (since only 4 components are used).

Suggested addition to tag-index

o pkg/obcs (tides):

• more efficient implementation of OB barotropic tidal velocity components,
  with both Cos & Sin of phase times Amplitude stored in common block ;
• more flexible specification of OB tidal-components (only read-in components
  with non-zero tidal period) ;
• allows to specify tangential component of barotropic tidal velocity at OB ;
• rename all OBC tidal component input parameter (period, amplitude and
  phase of each OB barotropic tidal velocity component input files) ;
• upgrade OB tidal example (in seaice_obcs/input.tides) with fewer tidal
  components (only 4) but including few tangential components.

[jm-c]

Two comments:

1. I would like to remove completely the CPP option OLD_OBCS_TIDES that I added in this branch once we are happy with the new implementation. So we don't have to document this new CPP option.
2. the reference output for the seaice_obcs test with tides (output.tides.txt) has been imported from PR Yoshi Nakayama code for tangential component of tides #73, and it shows that this test passes. Might update this ref. output before merging this PR into master branch, so that it will be cleaner.

[jm-c]

@menemenlis Any chancee you will find a little bit of time to review this PR ? Thanks.

[antnguyen13]

A quick q: generally we either use amp & phase or split into sin and cos components , e.g., y = Acos(2pi/T t) + Bsin(2pi/T t) = sqrt(A^2+B^2) cos (2pi/T t - tan^{-1}(B/A)) . Am I correct you're switching from the latter (old, am, ph) to the former (cos and sin)?

Also, in the "old" way of am/ph, it was only for the "normal" component, e.g., if eastern boundary, we get u = normal. What you're redoing now, not only to use the sin and cos form (thus getting away from am, ph), but that there'll be an additional sin + cos for the tangential component? Thus the way you've coded now, with cos+sin for u and same for v, naturally one will be "normal" and one will be tangential, right? Thus bypassing the mod code of PR#73 where the "t" am and ph were added?

[antnguyen13]

I think the way the new variables are created, it's a bit confusing to understand this line: i'm used to understanding the form am*cos(2pi/T (t - phase)) , which is what the old code was. Now, it takes the same form, but with variables OBE_uTideCs (<-- to do with cosine projection onto u-component?) replacing "am" and OBE_uTideSn (<-- to do with sine projection onto u-component?) replace "ph" . Are the new variables identical in definition? or are they just place holders?

[jm-c]

The lines you are pointed to are within "#ifdef OLD_OBCS_TIDES" (lines 63 - 124) and I don't disagree with you. This is the reason I added this comment in Aug 22, 2023:

Two comments:
I would like to remove completely the CPP option OLD_OBCS_TIDES that I added in this branch once we are happy with the new implementation. So we don't have to document this new CPP option.

[mjlosch]

If you only read this part of the code, then these lines appear to make no sense. But the fields are in "dual-use" and if OLD_OBCS_TIDES is defined, they have a different meaning: then OB?_uTideCs = OB?am and OBN_uTideSn = OB?ph and this corresponds to what one would expect.

I think the only way to avoid the admittedly confusing code is to keep variables with the old names just for the OLD_OBCS_TIDES code. Doing so only makes sense (too me), if we plan on keeping the OLD_OBCS_TIDES code for bitwise reproducibility. But as I understand, this is not the plan.

The way I understand this PR, one can reproduce the old results by making sure that OB?_uTideCs = OB?_vTideCs =OB?am and OBN_uTideSn = OBN_vTideSn = OB?ph by specifying the appropriate input fields (subject to roundoff error), so the OLD_OBCS_TIDES code is not necessary and can be removed. This would also remove the confusion.

Neither the old code nor the new code is yet properly documented anywhere. I'd vote for removing OLD_OBCS_TIDES code and document the new code properly (and maybe how to recover old results with the new code).

[jm-c]

@antnguyen13 Thanks for reviewing this and regarding you questions:

A quick q: generally we either use amp & phase or split into sin and cos components , e.g., y = Acos(2pi/T t) + Bsin(2pi/T t) = sqrt(A^2+B^2) cos (2pi/T t - tan^{-1}(B/A)) . Am I correct you're switching from the latter (old, am, ph) to the former (cos and sin)?

Also, in the "old" way of am/ph, it was only for the "normal" component, e.g., if eastern boundary, we get u = normal. What you're redoing now, not only to use the sin and cos form (thus getting away from am, ph), but that there'll be an additional sin + cos for the tangential component? Thus the way you've coded now, with cos+sin for u and same for v, naturally one will be "normal" and one will be tangential, right? Thus bypassing the mod code of PR#73 where the "t" am and ph were added?

You are right, and this match this PR description:

Implement suggestions from issue #617, renaming OB tidal forcing variables and improving
efficiency of tidal component addition. Also allows to add tides to OB tangential velocity.

And to reiterate, we continue to specify (in input files) the amplitude and phase (no changes here) but we store (in common block) the pair product: Amplitude x cos,sin [2piPhase/Period]

[mjlosch]

In general, the changes look good to me. I have these comments:

• At first sight, it is complicated to relate the new code to the old one, but that's maybe not a big problem.
• Related to the previous comment: With this PR the results of seaice_obcs.tides are very different. This may be intentional, but would it be helpful to create results that would reproduce results of the old code? I tried by using the old OB?am/ph.seaice_obcs and adjusting data.obc like this:

 OBCS_tidalPeriod  = 44714.16,43200.,45569.88,43081.92,86164.2,92949.48,86637.24,96726.24,1180295.64,2380716.,
 OBS_vTidamFile='OBSam.seaice_obcs',
 OBS_vTidphFile='OBSph.seaice_obcs',
 OBN_vTidamFile='OBNam.seaice_obcs',
 OBN_vTidphFile='OBNph.seaice_obcs',
 OBE_uTidamFile='OBEam.seaice_obcs',
 OBE_uTidphFile='OBEph.seaice_obcs',
 OBW_uTidamFile='OBWam.seaice_obcs',
 OBW_uTidphFile='OBWph.seaice_obcs',
 useOBCStides = .TRUE.,

and got nearly identical results (13 matching digits on my MacBook, when it was 14 with upstream/master).

• If we can illustrate or describe somewhere how one can reproduce the old results with the new code without the cpp-flag OLD_OBCS_TIDES, then we can certainly get rid off this cpp-flag and the corresponding code (and should).

[mjlosch]

Is there a good reason for the truncated naming? Why not call them OBN_uTideAmpltFile and OBN_uTidePhaseFile or OBN_uTidalAmpltFile and OBN_uTidalPhaseFileetc.

[jm-c]

All the other reference to field name within OB file-name parameter are short or very short (e.g., "u" instead of "uVel", "S" instead of "Salt"), so I would prefer to keep these ones short too (as were the previous names, with "am" and "ph" for Amplitude and Phase).

[jm-c]

@mjlosch Thanks for the review and the coments. Will think about file-names and other things but regarding this:

Related to the previous comment: With this PR the results of seaice_obcs.tides are very different. This may be intentional, but would it be helpful to create results that would reproduce results of the old code?

The "new" output.tides.txt is in fact not that new since it is the same as the updated one in PR #73 (and from the discussion there it seems that the "tangential-component" needed to be part of an updated test). So we could:

1. compared the 2 updated "data.obcs" from the 2 PR to see how one would transpose the old setting to new setting (it's mostly renaming file name setting); or
2. convert the original "data.obcs" to reproduce old results with the caveat that (a) it does not test the "tangential-component" setting and (b) there are many random number in the 10 tidal-component original input files, so may be updating this set-up is not a bad thing.

[mjlosch]

I think it is enough to document somewhere (maybe in seaice_obcs/README, which could also be renamed to README.md to make it a markdown file, or even in seaice_obcs/input.tides/data.obcs itself) that with this PR everything (namelist names, variable names) changed, but that it is possible to recover the old results obtained with old code and

 useOBCStides = .TRUE.,
#components   = M2       S2     N2       K2       K1      O1       P1       Q1       Mf         Mm
#periods (hr) = 12.4206  12     12.6583  11.9672  23.9345 25.8193  24.0659  26.8684  327.8599   661.31
 tidalPeriod  = 44714.16,43200.,45569.88,43081.92,86164.2,92949.48,86637.24,96726.24,1180295.64,2380716.,
 OBSamFile='OBSam.seaice_obcs',
 OBNamFile='OBNam.seaice_obcs',
 OBEamFile='OBEam.seaice_obcs',
 OBWamFile='OBWam.seaice_obcs',
 OBSphFile='OBSph.seaice_obcs',
 OBNphFile='OBNph.seaice_obcs',
 OBEphFile='OBEph.seaice_obcs',
 OBWphFile='OBWph.seaice_obcs',

by new code and

 useOBCStides = .TRUE.,
#components   = M2       S2     N2       K2       K1      O1       P1       Q1       Mf         Mm
#periods (hr) = 12.4206  12     12.6583  11.9672  23.9345 25.8193  24.0659  26.8684  327.8599   661.31
 OBCS_tidalPeriod  = 44714.16,43200.,45569.88,43081.92,86164.2,92949.48,86637.24,96726.24,1180295.64,2380716.,
 OBS_vTidamFile='OBSam.seaice_obcs',
 OBS_vTidphFile='OBSph.seaice_obcs',
 OBN_vTidamFile='OBNam.seaice_obcs',
 OBN_vTidphFile='OBNph.seaice_obcs',
 OBE_uTidamFile='OBEam.seaice_obcs',
 OBE_uTidphFile='OBEph.seaice_obcs',
 OBW_uTidamFile='OBWam.seaice_obcs',
 OBW_uTidphFile='OBWph.seaice_obcs',

This would provide some sort of recipe for recovering old results in different contexts, if required.
And we do not need to keep the *.seaice_obcs files. If necessary, they can always be recovered from an old checkpoint.

[jm-c]

@mjlosch I like your suggestion, I will do one or the other (REAME or data.obcs). Also will remove completely OLD_OBCS_TIDES code.

[jm-c]

@mjlosch I removed all the OLD_OBCS_TIDES code and add comment in README with reference to a small file for updating data.obcs content. You can check if it's ok.
@antnguyen13 I think I addressed your comments, and if not, please let us know.

And @menemenlis we are still waiting for you approval (you seems to agreee with issue #617, and this PR is mostly implementing that).

[mjlosch]

I'll push small modifications to the README.me. I can also add some links files in the read (e.g. input.tides/update_TideFileName.sed), so one can click on this on the web-page, if desired.

[jm-c]

@mjlosch the README.md now looks even better ! thanks.

[menemenlis]

It looks good to me, ready for merge, but a few minor comments for your consideration:

1. perhaps add a little more detail in MITgcm/verification/seaice_obcs/README.md, for example:

Note: naming of tidal input files (OB*File) have been changed and augmented to accomodate tangential flows in PR #752
(see input.tides/update_TideFileName.sed on how to update data.obcs to map the old normal flow components to the new names)

2. perhaps add a comment in OBCS_FIELDS.h and obcs_init_fixed.F to make clear that OB[N,S,E,W]_[u,v]Tide[Cs,Sn] arrays contain amplitude and phase during initialization before being converted to Sine and Cosine, for example:

C Note that during initialization, the Cs arrays contain Amplitude
C and the Sn arrays contain Phase

[jm-c]

@menemenlis Thanks for the review. I made little adjustment to the README.md file in seaice_obcs following your suggestion.
Regarding your second point, since the only place where

the Cs arrays contain Amplitude and the Sn arrays contain Phase

is locally in obcs_init_fixed.F, it might add more confusion if this comment is found in OBCS_FIELDS.h. And in obcs_init_fixed.F we already have:

454 IF ( useOBCStides ) THEN
455 C-- Read from files Barotropic Tidal Amplitude and Phase:
...
536 C-- Set Tidal coeff (= Amplit x COS & SIN of Phase ) from Amplitude and Phase:

so it's not obvious to me that we need to repeat/re-phrase it there. So at the end, would prefer to leave it as it is now.

[jm-c]

I've jsut updated this PR description (including suggested addition to tag-index) so, from my side, this PR is ready to be merged (as I understood, @menemenlis approval extends to @antnguyen13 approval) and will merge it in the coming days unless someone wants to add something.