WIP: add a task for plotting the sea surface height vs. obs

[xylar]

A task for plotting the SSH vs. observations from AVISO.

This is a work in progress (WIP)

[xylar]

@vanroekel, I don't want to compete with your work but I wanted to show you what I would do as an SSH comparison (without variability, so far).

I figure no matter what we do, this could presumably act as a starting point.

[xylar]

Here's what I'm seeing with the QU240 test case that I use a lot on my laptop. The fact that the SSH looks bad doesn't really tell us anything, but at least the values are of the right order of magnitude...

[xylar]

Okay, so now that I see that what's expected from us by Dec. 15th is the variability of SSH, that makes things a bit more challenging. That isn't a metric that I think we want to be computing offline if there's a way to get it added to timeSeriesStatsMonthly (e.g. standard deviation?)

[milenaveneziani]

@xylar: we just need to output SSH^2 in the timeSeriesStats stream (as well as SSH). I know that @mark-petersen has added an AM to compute squared values and things like uT, etc to compute eddy fluxes, so all we need to do is to turn on that AM and then have the values output in timeSeriesStatsMonthly.

[vanroekel]

@milenaveneziani is that right? To get the variability don’t we also need the deviation from a mean as well? In my implementation I am computing the variance of ssh^2 essentially, so I guess if we had that field, computing the variance is simpler, although having ssh^2 in the time series stats would make implementation of this task much more straightforward in the mpas-analysis framework.

[xylar]

@vanroekel, if we had <ssh> and <ssh^2>, the deviation is sqrt(<ssh^2> - <ssh>^2), right? I don't think the variance of ssh^2 would get at the variability of ssh, unless I'm missing something more complex.

[xylar]

One issue is whether we have time to add something new (e.g. ssh^2) to timeSeriesStatsMonthly and then to re-run a high-res run and analyze the results before Dec. 15th. Or do we really need to make due with what we have already (e.g. high-frequency output)?

[vanroekel]

I strongly would suggest for the 15th we make sure with the high frequency output for now. It doesn't make sense to me to run the high res again for one variable.

Regarding ssh^2. Computing the variance of this quantity came from conversation with mat. I can't remember the reason for using the variance of ssh^2. I may have also misconstrued our discussions. I'll follow up with him again

[milenaveneziani]

yes, you can just do the rms of the variance <(SSH-<SSH>)^2> = <SSH^2> - <SSH>^2

Mark has added an AM (need to check how it's called, something about eddy stats) that computes squares of fields and cross-products like u*v, u*T. So we just need to turn that on and have the <SSH^2> saved in the monthly files. This will give us monthly variances, which we can use to compute annual climatologies.

[milenaveneziani]

oops, I didn't realize the conversation had newer comments..
We can't change the output that we already have (and no, we can't re-run a high-res run for just one variable), but we can change now as soon as possible. Sorry that I didn't think about this earlier, I assumed we had that eddy stats on by default in the high res.

I think it will be OK to rely on the high-frequency output for the first part of the run.

[vanroekel]

@milenaveneziani thanks for the clarification. It wasn't clear to me and I must have misunderstood my conversations with Mat about what was computed previously in v0. I can modify my python notebook quickly to make the required plots for the 15th.

[vanroekel]

@xylar and @milenaveneziani here is the output from my python script for ssh variability for MPAS vs. the AVISO1993-2013 dataset. Certainly room for formatting improvements, but I think the base need for the 15th can be fulfilled with this plot.

[milenaveneziani]

Agreed, this looks great already. My only suggestion would be to use the newer data set, for which I processed years from 1993 to 2015. I will transfer the new obs file to NERSC and anvil.

[milenaveneziani]

@vanroekel: I do not have permission to write into the obs directory on anvil (/lcrc/group/acme/mpas_analysis/observations). Could you please change that? Thanks.

[vanroekel]

@milenaveneziani there is full rwx permissions for the climate/collab groups on anvil. Are you not a part of either of those? Looking on anvil I don't see that you have a scratch space in /lcrc/group/acme. I'm a little hesitant to give write permissions to all.

[mark-petersen]

@vanroekel sorry I didn't make this connection before about the eddy product AM. I also forgot that the eddy product AM is not on by default in ACME. It should be turned on for eddying resolution. Here is the location of the ACME default:
components/mpas-o/bld/namelist_files/namelist_defaults_mpas-o.xml:779:<config_AM_eddyProductVariables_enable>.false.</config_AM_eddyProductVariables_enable>

@vanroekel I'm guessing this is untested in the 18to6, but it doesn't need any input files like the moc does, so shouldn't cause trouble.

The AM is hard-wired to compute SSH*2, u^2, v^2, u*T, v*T. The namelist and streams file would then look something like:

&AM_eddyProductVariables
    config_AM_eddyProductVariables_enable = .true.
    config_AM_eddyProductVariables_compute_interval = '0000-00-00_01:00:00'
    config_AM_eddyProductVariables_output_stream = 'eddyProductVariablesOutput'
    config_AM_eddyProductVariables_compute_on_startup = .true.
    config_AM_eddyProductVariables_write_on_startup = .false.
/

<stream name="eddyProductVariablesOutput"
        type="output"
        filename_template="mpaso.hist.am.eddyProductVariablesOutput.$Y-$M-$D.nc"
        filename_interval="01-00-00_00:00:00"
        clobber_mode="truncate"
        reference_time="01-01-01_00:00:00"
        output_interval="00-01-00_00:00:00"
        packages="eddyProductVariablesAMPKG" >

        <var name="xtime"/>
        <var name="SSHSquared"/>
        <var name="velocityZonalSquared"/>
        <var name="velocityMeridionalSquared"/>
        <var name="velocityZonalTimesTemperature"/>
        <var name="velocityMeridionalTimesTemperature"/>
</stream>

and those variables would also need to be put into the monthly time averaged file.

@vanroekel would you be willing to run a one-day test with the above changes in the RRS18to6 and check the timers? If it is reasonable, I'll make a PR to turn it default on in ACME for mid and high resolutions. Then I will point to that PR, as this issue is on the analysis repo.

[milenaveneziani]

ok, I transferred SSHstd_AVISO_1993-2015.nc to both nersc and anvil's observations directory.

[xylar]

It sounds like this PR should be put on hold for now while @vanroekel and @mark-petersen sort out the issues with the associated AMs and we reach agreement on what the observations are that we will use for SSH and SSH^2.

My preference would be to support seasonal variability in the observations as an option even if we don't use it very often. We have infrastructure in palace for easily choosing from a number of seasons so it seems unfortunate if the observations aren't there to make meaningful comparisons.

[milenaveneziani]

We have daily SSH Aviso data from 1993 to 2015, so definitely we can compute any season we want from the observations.

[xylar]

@milenaveneziani, definitely. We should just make sure we have monthly climatologies of SSH and rms(SSH) or SSH^2 that we can use in the analysis so we don't have to recompute them from daily data. It doesn't seem like we have that yet but I'm sure it wouldn't be that hard to compute. You're more experienced at that kind of thing than I am so maybe you could help me with it or at least give me tips on what you would do?

[milenaveneziani]

The AVISO data is Mean Sea Level Anomaly (MSLA), so from that we can only get SSH variability. For the mean SSH/dynamic height, we would probably want to use the Maximenko et al. data set that you mentioned in #278: would you agree @maltrud?
In terms of priority for this PR, which I think is the one we want to have in by ~Dec 15, we should focus on the annual SSH variability computed from the eddy stats AM, if possible (and otherwise use the high-frequency AM). That is the Tier 1a metric we have in the Ocean Metrics table.

After that, we can certainly add a seasonal cycle in SSH variability, and the mean SSH metric (for that we will need to download the observations and process them).

[xylar]

@milenaveneziani, I don't think there's any need for a PR by Dec. 15. We have committed to doing the analysis by then but not necessarily here in MPAS-Analysis. So the priority is figuring out what the analysis will be by Dec. 15th and not adding an analysis task. This is why I think this PR should be put on hold for now.

[milenaveneziani]

@xylar: do you have a newer branch for this PR? (I guess I meant 'for this task')

[xylar]

No, not yet but I think it makes sense to start over with var(ssh) instead of ssh.

…

On Jan 16, 2018 7:43 PM, "Milena Veneziani" ***@***.***> wrote: @xylar <https://github.com/xylar>: do you have a newer branch for this PR? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#279 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AD_EeJ-MWKyQfPrpEO1V57BLwNI7zqHmks5tLO3ngaJpZM4Qq2Af> .