@Michael Levy, @Keith Lindsay, here is a gx3v7 FESEDFLUX file:
/glade/work/mclong/cesm_inputdata/fesedflux_total_reduce_oxic_POP_gx3v7.c200618.nc
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
@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
no recollection. It might make sense for me to revist the workflow and check
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.
@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
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.
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?
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
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.
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.
You're finding bugs faster than I can fix them! :) So far
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 fieldsref_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 outStill to do (hopefully including everything mentioned above):
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.JRA025m
instead of rx1
IRON_FLUX
is accumulated, and make sure MOM isn't trying to apply any unit conversion since it's already in mksSome 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.
Updated list (I realized I had left one off the earlier version)
Done
_FillValue
for 2D fields with ref_depth
below ocean depthvolcello
and h
are written to 3D bgc streams on native grid, geolat
and geolon
are written to all bgc streamsTodo
(i,j)
loops in the driver [this should've been on the earlier list]IRON_FLUX
is accumulatedJRA025m
instead of rx1
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.
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.
I think I fixed both issues with IRON_FLUX
:
kg / m^2 / s
but MARBL wants input in nmol / cm^2 / s
so need to multiply by 1e8 / molw_Fe
, or roughly 1.8e6Looking 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?
Also, I'm starting to catch up to Keith's bug-reporting...
DONE:
bgc_daily_z
_FillValue
for 2D fields with ref_depth below ocean depthvolcello
and h
are written to 3D bgc streams on native grid, geolat
and geolon
are written to all bgc streams(i,j)
loops in the driverJRA025m
instead of rx1
IRON_FLUX
, and I think fixed the general computation as well (see above)TODO:
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.
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
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
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:
I think looking at Frank's fill tool is a good starting point, hopefully I can recreate the environment that runs this notebook...
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
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.
The notebook is still re-generating initial conditions, but I hope to get the two mom cases running later this afternoon or evening
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 :)
I think I've got FESEDFLUX
plumbed in correctly (thanks to @Keith Lindsay for all the help). There were two big things going wrong:
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 columnsFESEDFLUX
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:
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...
Thanks @Michael Levy! This sounds good to me.
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