Stream: MARBL

Topic: forcing


view this post on Zulip Matt Long (Feb 04 2021 at 22:29):

@Michael Levy, @Keith Lindsay, here is a gx3v7 FESEDFLUX file:
/glade/work/mclong/cesm_inputdata/fesedflux_total_reduce_oxic_POP_gx3v7.c200618.nc

view this post on Zulip Michael Levy (Feb 04 2021 at 22:48):

Cool, I'll do a quick run on POP and then (1) add the file to inputdata, (2) make a PR to update the default on master, and (3) let Jun know how to update the file in her runs

view this post on Zulip Michael Levy (Feb 04 2021 at 23:22):

@Matt Long -- Do you recall off the top of your head what you changed between the June 16th and June 18th versions of these files? I noticed that the default for tx0.1v3 is the June 16th file (although the June 18th one is in your work directory), and so that's what we are using in #0.1° JRA BGC Run

view this post on Zulip Matt Long (Feb 04 2021 at 23:23):

no recollection. It might make sense for me to revist the workflow and check

view this post on Zulip Michael Levy (Feb 04 2021 at 23:27):

It looks like we last talked about forcing in #0.1° JRA BGC Run > forcing on June 16th, but you mentioned

thinking more about this, I would like to have an opportunity to revisit the fesedflux field to ensure comparability with what we use for the 1° spin-up.

view this post on Zulip Michael Levy (Feb 05 2021 at 00:01):

@Matt Long My run crashed due to a NaN in Iron Sediment Flux. It looks like there is something wrong with the land mask in your forcing file -- I'm counting 3692 missing values in the surface layer, but there should only be 3689 missing values for the 11600 surface grid points. Log files for my run are in /glade/scratch/mlevy/c.e23b1.C1850ECO.T62_g37.new_fesedflux_file/run, and I just used xarray / numpy to compare your file to the topography file:

>>> import xarray as xr
>>> import numpy as np
>>> ds = xr.open_dataset('fesedflux_total_reduce_oxic_POP_gx3v7.c200618.nc')
>>> da = ds['FESEDFLUXIN'].isel(z_t=0)
>>> np.sum(np.isnan(da.data))
3692
>>> with open('../grid/topography_20100105.ieeei4', 'r') as f:
...     a = np.fromfile(f, dtype='>i4')
...
>>> sum(a == 0)
3689

view this post on Zulip Keith Lindsay (Feb 07 2021 at 17:11):

The FESEDFLUXIN variable in Matt's gx3v7 file has missing values for the Caspian Sea, while KMT is positive for the Caspian Sea for that grid.

view this post on Zulip Michael Levy (Mar 04 2021 at 23:01):

When Jun gave her update, I realized that we never really wrapped this up. @Matt Long, if it would be straightforward could you generate a gx3v7 forcing file that treats the Caspian as active ocean? Or should I update https://github.com/marbl-ecosys/marbl-forcing and do it myself?

view this post on Zulip Matt Long (Mar 04 2021 at 23:03):

I should be able to do it—but have been distracted with other things. We need to revisit the workflow too as I pulled out the FESEDFLUX into it's own repo: https://github.com/marbl-ecosys/forcing-Fe-sedflux

view this post on Zulip Keith Lindsay (Oct 06 2021 at 22:25):

In the latest MARBL MOM dev runs, I notice a steep increase in some BGC tracers in the Sea of Azov, a smaller sea on the north side of the Black Sea. For instance, DIC>8000, ALK>9000, SiO3>130 after 20 years.

It looks like we have large riverine BGC fluxes into this region. The large fluxes are in the dataset riv_nut.gnews_gnm.rx1_to_tx0.66v1_nnsm_e1000r300_190315.20210405.nc, in the directory /glade/work/mlevy/cesm_inputdata.

They don't appear to be present in riv_nut.gnews_gnm.JRA025m_to_tx0.66v1_nnsm_e333r100_190910.20210401.nc, in the same directory. I don't understand why the rx1_to_tx0.66v1 file has this feature that the JRA025m_to_tx0.66v1 doesn't have.

That said, the latter file seems preferable for these runs anyway, since we're using JRA forcing, which presumably has the runoff on the JRA025m grid.

view this post on Zulip Keith Lindsay (Oct 07 2021 at 02:55):

Looking at animations of the MARBL MOM oob run, it looks like primary productivity is declining over the course of the run, particularly in the Southern and Pacific Oceans. (I have annual means of some output streams in directory /glade/scratch/klindsay/mom_oob/ocn/proc/.) The Fe limitation factors are dropping, and the other limitation terms are not. IRON_FLUX is 8 orders of magnitude smaller in MOM than in POP. I haven't dug into the code yet to see if this is a diagnostic issue, a forcing issue, or both. I'm guessing both. This is an unusual diagnostic in POP, as its output units are MKS, as opposed to most other diagnostics whose units are CGS. The field is scaled before it is passed to POP's tavg accumulate subroutine. Perhaps something is off in the translation of this to MOM.

view this post on Zulip Michael Levy (Oct 07 2021 at 14:21):

You're finding bugs faster than I can fix them! :) So far

  1. I no longer create the bgc_daily_z stream (I had added {tracer}_SURF and {tracer}_zint_100m to the MOM output list but forgot to tell my scripts those are 2D fields
  2. MOM recognizes when 2D fields have a ref_depth below the bottom of the column, and replaces those values with 0. I think I found the missing value variable as well, so later this morning it should appropriately mask that out

Still to do (hopefully including everything mentioned above):

  1. Update which variables are written to which stream: we want geolat and geolon in every output file, and then I should pare down the additional variables that are added to streams with 3D variables (volcelo might be the only required field?) For 3D fields written on the native grid, I plan to include h as well to make it a little easier to reconstruct depth. Should I also include ocn_depth so we don't force the topmost interface to be 0? This would basically let us reconstruct SSH by summing h and subtracting the bathymetry.
  2. Switch river nutrients to the file mapped from JRA025m instead of rx1
  3. Look into how IRON_FLUX is accumulated, and make sure MOM isn't trying to apply any unit conversion since it's already in mks

view this post on Zulip Keith Lindsay (Oct 07 2021 at 15:55):

Some of the tracers, biomass and non-refractory DOM, in the MARBL MOM dev runs are negative. With MOM's numerics we should be able to avoid this. In this case, it looks like the negative values are in the initial condition file, which for these tracers looks like it was regridded from the spunup CESM2 IC. In POP, when we read tracer values from a non-restart file, we set negative values to 0. I suggest we also do that in the MOM MARBL driver, after these lines.

view this post on Zulip Michael Levy (Oct 07 2021 at 16:07):

Updated list (I realized I had left one off the earlier version)

Done

  1. No longer create bgc_daily_z
  2. Use _FillValue for 2D fields with ref_depth below ocean depth
  3. volcello and h are written to 3D bgc streams on native grid, geolat and geolon are written to all bgc streams

Todo

  1. Fix order of 2 (i,j) loops in the driver [this should've been on the earlier list]
  2. look into how IRON_FLUX is accumulated
  3. Set any negative initial tracer values to 0
  4. Switch river nutrients to the values mapped from JRA025m instead of rx1

view this post on Zulip Keith Lindsay (Oct 07 2021 at 16:20):

Looking at the MARBL MOM t061 IC, there are numerical artifacts close to topographic features in the tracers that were regridded from the POP grid to the MOM grid. For instance, there are isolated zeros around the northern coastline of Canada.

Did the regridding of these variables normalize the result if the source grid only partially covered the destination grid?

We might also need to apply a fill algorithm to get values set where the regridding doesn't get the job done.

I recall Frank updating a regrid/fill algorithm to take care of such concerns, and be fast, but I don't remember which repo it was in. It might have been pop_tools.

view this post on Zulip Keith Lindsay (Oct 07 2021 at 16:22):

Oh, and the metadata in that IC file has some issues. All of the fields that were regridded from the POP grid to the MOM grid have long_name and units from alkalinity.

view this post on Zulip Michael Levy (Oct 07 2021 at 20:35):

I think I fixed both issues with IRON_FLUX:

  1. As @Keith Lindsay pointed out, I was missing a few terms in the computation but I added them in
  2. The units were also wrong -- MOM is storing values from the coupler in kg / m^2 / s but MARBL wants input in nmol / cm^2 / s so need to multiply by 1e8 / molw_Fe, or roughly 1.8e6

Looking at the output from the model now, POP's IRON_FLUX from Jan 0001 has a range of (5e-14, 7e-8) while MOM is (1e-13, 3e-9). (MARBL does one last conversion, so this output has units of mmol / m^2 / s in the history files). I was a little surprised that MOM's max value was 1/20th the size of POP's max value, but does that seem like a reasonable change given differences in the grids, numerics, etc between the models?

view this post on Zulip Michael Levy (Oct 07 2021 at 20:44):

Also, I'm starting to catch up to Keith's bug-reporting...

DONE:

  1. No longer create bgc_daily_z
  2. Use _FillValue for 2D fields with ref_depth below ocean depth
  3. volcello and h are written to 3D bgc streams on native grid, geolat and geolon are written to all bgc streams
  4. Fixed order of 2 (i,j) loops in the driver
  5. Set any negative initial tracer values to 0 when not reading from restart file
  6. Switched river nutrients to the values mapped from JRA025m instead of rx1
  7. Got the units right for IRON_FLUX, and I think fixed the general computation as well (see above)

TODO:

  1. Investigate IC generation notebook, figure out how to get better values near topography for tracers mapped from previous POP run

view this post on Zulip Keith Lindsay (Oct 07 2021 at 20:45):

Peak values for atm dep usually occur along the coastline. Since these grids are different resolution, one grid might pick up an extreme coastal value that the other grid is missing. So looking at the range is tricky to interpret.

I suggest choosing an open ocean feature, .e.g, for IRON_FLUX, the low values along the equator in the Pacific, and use the mouse to determine typical values in the vicinity of the region on both grids. If those roughly agree then things are looking up.

view this post on Zulip Michael Levy (Oct 07 2021 at 21:10):

Thanks, that's good advice! I think I'm off by a factor of 10... setting the MOM range to [1e-12, 1e-8] looks a lot like POP's field with [1e-11, 1e-7]:

MOM-IRON_FLUX.png
POP-IRON_FLUX.png

view this post on Zulip Michael Levy (Oct 07 2021 at 22:09):

Found it - I was converting from mks -> cgs twice... I think this will improve both IRON_FLUX and DUST_FLUX, which might also have been off by a factor of 10

view this post on Zulip Michael Levy (Oct 08 2021 at 17:28):

Keith Lindsay said:

Looking at the MARBL MOM t061 IC, there are numerical artifacts close to topographic features in the tracers that were regridded from the POP grid to the MOM grid. For instance, there are isolated zeros around the northern coastline of Canada.

Did the regridding of these variables normalize the result if the source grid only partially covered the destination grid?

We might also need to apply a fill algorithm to get values set where the regridding doesn't get the job done.

I recall Frank updating a regrid/fill algorithm to take care of such concerns, and be fast, but I don't remember which repo it was in. It might have been pop_tools.

The notebook that generated these initial conditions is here... and looking closer, there is definitely some unexpected 0s in the fields mapped from POP. DOPr was the easiest for me to see, since most of the values are all pretty far from 0 in that field so the small values really stand out:

DOPr.png

I think looking at Frank's fill tool is a good starting point, hopefully I can recreate the environment that runs this notebook...

view this post on Zulip Michael Levy (Oct 08 2021 at 17:41):

Looking closer at the plots in the notebook, it seems like maybe we aren't constructing the POP values correctly -- instead of missing values, land has 0s? That would explain why the WOA data mapped over better. If that's the case, masking out the continents before doing the POP -> MOM map should improve the look of things

view this post on Zulip Michael Levy (Oct 08 2021 at 22:10):

I ended up using pop_tools.get_grid() to get POP's REGION_MASK and then using that to mask out the POP data... the DOPr field looks promising -- note that the color bar starts around 0.02 instead of right at 0. We already had calls to lateral and vertical fill routines, but I don't think they were doing anything because the source grid was happily mapping 0s from land rather than letting missing values propogate / get filled in later.

Updated-DOPr.png

The notebook is still re-generating initial conditions, but I hope to get the two mom cases running later this afternoon or evening

view this post on Zulip Michael Levy (Oct 09 2021 at 12:42):

I restarted mom_hybrid_z and created mom_oob.002 last night, both built with the latest code documented above and started with regenerated initial conditions. The hybrid case is on the big pe layout, so it should be finishing 4 years every 10 hours or so; the out-of-the-box run is on the old layout, so it’s running in 2-year chunks. I’m going to be out-of-touch until Monday night, but hopefully when I’m back in the office on Tuesday things look okay :)

view this post on Zulip Michael Levy (Oct 22 2021 at 14:46):

I think I've got FESEDFLUX plumbed in correctly (thanks to @Keith Lindsay for all the help). There were two big things going wrong:

  1. After reading in FESEDFLUX from disk, I was trying to accumulate all the flux from below the ocean floor into the layer containing bathyT... but it turned out that when bathyT was exactly equal to a layer interface depth on the forcing grid, it would accumulate into the layer below bathyT instead of above it. So we were losing all the sub-floor flux in many columns
  2. Keith pointed out that FESEDFLUX is the only diagnostic we have where we want to conserve the column sum when vertically interpolating. It turns out that there is a specific flag to use when registering a diagnostics (v_extensive = .true.) to conserve the column sum when remapping from the model vertical to whatever grid is being used when writing history files.

I also had a bug in my notebook that compares the column sums: I was writing the MOM output to a history file that used the POP vertical levels to make it easy to compare the input and output files. However, the MOM output has missing values below bathyT while the forcing file has full columns at every ocean cell (I believe the forcing file that was mapped from the POP horizontal to the MOM horizontal has 0s below KMT).

Lastly, I added a step to the generation of the forcing file. I now apply the POP land mask to the source file, and then use pop-tools to laterally fill the land before mapping to MOM (and the POP -> MOM mapping file does not know about the POP land mask). I had noticed some weird patterns of 0s in the MOM forcing file that I traced to areas of the globe that are land in POP but open ocean in MOM (e.g. between Australia and Papua New Guinea) and this seemed like a reasonable way to avoid those patterns.

Anyway, I think everything is resolved now and I've launched mom_oob.003 and mom_hybrid_z.003 with these updates. :fingers_crossed:

view this post on Zulip Michael Levy (Oct 22 2021 at 14:54):

Also, using 81 nodes per run (79 for POP, 1 for CICE, 1 for everyone else), I'm seeing ~10.5 SYPD for the oob run. That layout was good for ~9 with the hybrid grid, so both of these runs should be complete over the weekend...

view this post on Zulip Matt Long (Oct 22 2021 at 15:08):

Thanks @Michael Levy! This sounds good to me.

view this post on Zulip Michael Levy (Oct 22 2021 at 21:10):

FYI I just restarted both runs around 3:00 today... so they should both still finish over the weekend but now there will be a stream file with output using the POP vertical levels (-ish). Still expect them to finish tomorrow or Sunday, to be looked at Monday morning


Last updated: May 16 2025 at 17:14 UTC