ROMS different src version leading to unalike results

[Timh37]

Dear all,

I'm running a ROMS configuration for the northwestern European shelf region. The configuration has 1/4° by 1/4° resolution, and some smoothing of the bathymetry is applied to avoid horizontal pressure gradient errors. It's forced by ERA_interim atmospheric forcing and GLORYS lateral ocean boundary conditions. We're focusing on sea level output. I've received the input files from another researcher and encountered some problems setting this configuration up on the HPC environment of my institution. The model is running OK, but I'm getting dissimilar results in terms of time-mean sea level.

I have triple-checked that all my input settings agree with the files that I received. I'm using the latest version of ROMS at the time of downloading, i.e. version 3.7 SVN revision 966M. When running for 10 years, I get the following time-mean sea level (left panel), which is different from the 'original' results obtained by the colleague that tuned the model using version 3.7 revision 783 (middle panel). I included an observational product, albeit over a longer time period, for comparison (right panel):



In all cases the time-mean sea level displayed are anomalies with respect to the regional mean.

I proceeded downloading the revision 783 src code on my HPC environment, and ran the experiment again. The result was exactly the same as the middle panel. So, it seems like the solution of the original model is strongly affected by updates in the ROMS version. It seems odd that a different version can have so much effect. We would like to find out if there are particular reasons that could cause the differences that we're seeing? Any help would be greatly appreciated, and I'm happy to provide more information is needed.

Kind regards,
Tim

[kate]

You can hunt this down by using the equivalent of "git bisection", where you check out the version in the middle and see which camp it falls into. Pick a simulation of only a day or so to compare to. You can iteratively try intermediate versions until you find out which commit changed your results, then just look at that code change to see what was done to change answers.

[arango]

The ROMS basic numerical kernel has not changed much in a few years. However, configuration, CPP options, and parameters in the standard input have been modified and enhanced. You need to be sure that your header file (application.h) and standard input script (roms.in) is equivalent between the two versions. Usually, users ignore the updates to roms.in. Sometimes, CPP options are replaced by standard input parameters. All of these changes are explained in detail in trac for each code update. Except for bug corrections, ROMS is backward compatible, and solutions should be similar within round off with the same version of the Fortran compiler.

Because of coupling, the time clocks in ROMS are now very strict and expanded for each forcing file. One thing that I usually do is to compare the standard output between different solutions. I use KDiff3 to compare files. It is a wonderful free tool and easy to install/ You also can compare the *.f90 in the Build subdirectory files between the two solutions. It is a faster way to analyze the solutions.

[Timh37]

Thanks Kate & Arrango, both great suggestions. I have one additional question regarding the .h and .in files - I have tried to copy the .h and .in files for the older revision code directly to the environment with the newer revision code, and ran the experiment again. I thought this way I would make sure that these files are identical between the two different versions. This still gives different results between the older and newer versions. Can we conclude from this that the different results indeed come from changes in the code or bug fixes? Or is the test I'm doing not sensible?

[wilkin]

There are a lot of changes since version 783 (3 or more years ago?) so it is going to be tedious to track down.

As Kate suggests, you can do intermediate updates to successively newer versions and see when changes in your solution arise. This will be a bit messy because you may need to thoughtfully update your ocean.in (now roms.in) on each step.

For guidance on what changes occurred with each revision you can read the corresponding log message in trac at https://www.myroms.org/projects/src/log/ and jump past revisions that are extremely unlikely to have had any impact. Note that the ticket number is not the revision number. @rev is the first column ... go to the right most column "src:ticket ..." to read the revision notes.

You can also go to the trac browser page
https://www.myroms.org/projects/src/browser
and click on "View changes ..."

You'll need to fill in the fields so that "from" is something like
/trunk/ROMS/Nonlinear
and of course the revision numbers 966 and 783.

This will throw up a huge number of differences, so again you might want to use small ranges of change to make it easier to follow.

To answer your last question ... Having identical application .h files will not guarantee the same solution. Where bugs have been fixed we can't retain full backward compatibility.

Also, though you have checked that your inputs are the same, a useful check on yourself is whether ROMS itself reports that the inputs were the same. The log file and global attributes in output netcdf files list the CPP options actually used as opposed to what you think you used. Looking for changes there might give you some guidance on tracking this down.

[Timh37]

Thanks John for the response. After testing with different versions I think I've pinpointed the change in solution (based on values of zeta) to version 839 to 840. This corresponds to ticket #723: https://www.myroms.org/projects/src/ticket/723.

Ticket 723 concerns the following bug correction:

Several miscellaneous corrections are applied:

Using potential/neutral (pden) density in the tracer diffusion isopycnic tensor in routines t3dmix2_iso.h and t3dmix4_iso.h. It was using in situ density (rho). Oops, I missed this one. I just discovered by looking at those routines. I thought that we were using pden.
Correct the extraction of free-surface in shallow-water nesting applications. We need instead in nesting.F routine get_refine:

Added wet/dry mask during penetration of solar radiation into the water column in pre_step3d.F

So apparently my solution of zeta is affected by this bug correction? I've proceeded to replace the corrected .F files in my src version one by one. It turns out that the change in t3dmix2_iso.h is the file that affects the solution. Is it logical that this changes the mean dynamic topography posted earlier in this thread by such a large amount? And should I interpret this as all solutions obtained with versions before 840 somehow being wrong? What is the right philosophy of using earlier versions than the most recent one given these bug fixes?

[Timh37]

An addition to the above:

I've rerun the experiment with rvn 966 (more recent version), now with mixing along geopotential lines instead of along isopycnals. The results I obtain are now very similar to the results obtained with version rvn 783, unlike when using rvn 966 with isopycnal mixing (first post).

Could this hint to some issues with the bug fixed version of t3dmix2_iso.h, since I wouldn't expect such big differences between the two different mixing schemes?

[wilkin]

The rotation of the mixing tensor to "isopycnal" coordinates should really be done with respect to the slope of potential density surfaces referenced to the local depth, not the sea surface. I checked on this yesterday with Trevor McDougall, just to be sure. So neither option (not in-situ density or potential density relative to sea surface) is really correct.

We never implemented rotation in local potential density because this requires the cost of another equation of state calculation.

Below 1000 m the in-situ density begins to align quite strongly with geopotentials, so this is consistent with your experience:

I've rerun the experiment with rvn 966 (more recent version), now with mixing along geopotential lines instead of along isopycnals. The results I obtain are now very similar to the results obtained with version rvn 783, unlike when using rvn 966 with isopycnal mixing (first post).

I don't know of many users actively running with the isopycnal rotated tensor. Geopotential mixing is the more widely used approach. If you have relatively high resolution then much of the lateral mixing is actually be accomplished by resolved eddy stirring, and the model dynamics will effectively do this in density coordinates because it's resolved advection (not parameterized diffusion).

Also, make sure your lateral diffusion coefficient is not too large.