Immediate 'bounding bracket for pH solution not found' error

[jcampbell2548]

Hello,

I'm attempting to run a CESM2.1.5 simulation with the following compset: 1850_CAM50_CLM40%CN_CICE_POP2%ECO_RTM_SGLC_SWAV in f19_g16 centered around the Mid Miocene (~15 million yrs ago). The simulation runs absolutely fine when MARBL is turned off. When I turn MARBL on, I inevitably run into this "bounding bracket for pH solution not found" error immediately after the ocean initialization and writing the pop.h.once.nc file. I have tried changing dt_count and my pe layout but that hasn't helped. The coordinates also don't seem to come from 1 specific region. My ecosys initial conditions file does use Miocene geography, as does my nitrogen deposition file. Here is an example of part of the cesm.log:

MARBL WARNING (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) x1,f = 0.1289921E-009 0.4220937E+072
MARBL WARNING (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) it = 2
MARBL WARNING (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) x2,f = 0.2044387E-006 0.1378627E+072
MARBL ERROR (marbl_co2calc_mod:drtsafe): bounding bracket for pH solution not found
MARBL ERROR (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) dic = 0.1975949E+004
MARBL WARNING (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) x1,f = 0.3358608E-008 0.2515583E+072
MARBL ERROR (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) ta = -0.5820493E+044
MARBL WARNING (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) x2,f = 0.2119139E-007 0.1923800E+072
MARBL ERROR (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) pt = 0.1767677E+078
Message from (lon, lat) ( 21.368, 25.022), which is global (i,j) (49, 232)
MARBL ERROR (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) sit = 0.1767677E+078
MARBL ERROR (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) temp = 0.8515730E+001
MARBL WARNING (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) it = 2
MARBL ERROR (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) salt = 0.3630198E+002
MARBL WARNING (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) x1,f = 0.3360479E-008 0.2511535E+072
MARBL WARNING (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) x2,f = 0.2120319E-007 0.1923635E+072
Message from (lon, lat) ( 57.856, 45.550), which is global (i,j) (83, 265)

Any advice is appreciated! I can send my case directory if it helps.

Thanks so much,
Julia

 

[mlevy]

Was this the issue you were talking about in the other thread? In this one, the model is reporting that total alkalinity is O(10^44) and total phosphate and silicate are both O(10^78):

MARBL ERROR (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) ta = -0.5820493E+044
MARBL ERROR (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) pt =  0.1767677E+078
MARBL ERROR (marbl_co2calc_mod:drtsafe): (marbl_co2calc_mod:drtsafe) sit =  0.1767677E+078

One issue with the cesm.log file is that lots of different processes write to the file simultaneously, so the output is somewhat jumbled. Is there a line like

Message from (lon, lat) (  21.368,  25.022), which is global (i,j) (49, 232)

Preceding the block you posted that might point to where on the grid these values are showing up? As mentioned in the other thread, my guess is that this is a mismatch in the grid in the initial condition file and the model. What computer are you doing this run on? If it's somewhere I have access, I'm happy to investigate more.

 

[jcampbell2548]

Yes, this is the issue I meant, and ah okay I will look for that! I'm running this on Derecho (down today) so will need to wait until tomorrow to check the logs. Oddly, one of the Miocene initial conditions files I've tried using (with actual values, not set to 0) worked on CESM1 with BGC on Cheyenne ~2 years ago, but now will not work on Derecho even when I try using CESM1 again. This CESM2 Miocene case is at /glade/work/juliacam/case/MMCO1xMARBL.21 if you want to take a look; I have other cases with CESM1-BGC and some Oligocene cases that have all resulted in the same error. The case "LOWP2xBGC" ran for 2 days which is the farthest I've gotten any of these but ends in the same error. Thank you for your help!

 

[mlevy]

I looked at /glade/work/juliacam/case/MMCO1xMARBL.21, and the initial condition file for that run is /glade/work/juliacam/MMCO/IC_ecosys_MIOv2.nc, which just has zeroes. As you've seen in your run, MARBL can't handle initializing everything to zero, so you need something more realistic. This is where it gets thorny -- the ocean BGC group doesn't have any expertise with the Miocene, so we don't know what a reasonable initialization would look like. One possibility would be to take level-by-level global / temporal means from the last year of a pre-industrial run, but that would be introducing the implicit assumption that the BGC inventory in the Miocene matches the recent preindustrial era. Another issue to be aware of is that you don't have any riverine inputs -- MARBL is set up to balance river fluxes with sediment burial, so without river fluxes you should set the burial coefficients to 0 to avoid long-term drift as sediment is buried but not replenished.

I did look at LOWP2xBGC, which has non-zero initial conditions in LOWP2xBGC.pop.r.0001-01-03-00000.nc, but there are some non-physical values:

* O2 approaches 10^12 in some cells
* Fe has missing values everywhere in the ocean

Lastly, I have some concerns about your Miocene grid -- it appears that there are several inland seas, and in the modern day grid we need to apply restoring / nudging to nitrate, silicate, phosphate, oxygen, and alkalinity to keep them reasonable in these Black Sea (as well as the Caspian Sea, when that was part of the CESM ocean grid).

So I'm happy to help you try to get this run going, but there are several science questions to consider which will play a role in how we configure the model that should be answered before going too far down any one path.

 

[jcampbell2548]

Thank you for taking a look!! I have tried running with non-zero initial conditions for the Miocene (see /glade/work/pacosta/PaleoBC/acostar/Mod_POP_Eocosys/Init.MIO_Ecosys.pop.h.climov4.nc as an IC file that worked ~2 years ago on Cheyenne) but it ends in the same error. Riverine inputs file is tricky to change to paleo with non-zero values so I will set burial coefficients to 0.

Do you know how I can apply restoring/nudging to those variables so they remain reasonable in the inland seas? I cannot change the Miocene geography at this point to erase those seas.

My purpose to running this simulation is to attempt to see phytoplankton biomass and primary productivity results in the Miocene at 280 ppm vs 560 ppm and conduct a model-data comparison with some collaborators. It's ok if several variables are unreasonable as long as variables like PIC/POC, DIC/DOC, and CaCO3 are somewhat realistic. I can also run the model for some time if needed, I was told by someone that if I'm able to run the model for a few decades then it's ok to start with somewhat unrealistic values. So it doesn't need to be perfect, I'm mostly looking globally at potential changes in marine carbon sequestration with doubled CO2 in this past warm interval.

Let me know what you think and thanks for looking into it!

 

[mlevy]

For restoring, look at the files /glade/campaign/cesm/cesmdata/inputdata/ocn/pop/gx1v6/forcing/ecosys_restore_inv_tau_POP_gx1v6_20170125.nc and /glade/campaign/cesm/cesmdata/inputdata/ocn/pop/gx1v6/forcing/ecosys_restore_POP_gx1v6_20220404.nc

The first file needs a field that is 1/[restoring time scale], which has units 1/s. It should be 0 in the open ocean (equivalent to "no restoring") and something like 1/(2*365*24*3600) = 1.59e-8 in the marginal seas [this sets a restoring time scale of 2 years and that seems to stabilize the model]. This particular file has two fields, but RESTORE_INV_TAU_MARGINAL_SEA_ONLY is the one to mimic. Notice that it has (i,j,k) indices, and you want the non-zero values to extend to the ocean floor.

The second file contains the fields you want to restore to NO3, SiO3, PO4, O2, ALK, and ALK_ALT_CO2

I went back to /glade/work/juliacam/case/MMCO1xMARBL.21, and actually wonder if maybe part of the problem is that it looks like tracer restoring is enabled but you are pointing to those two files that are set up on the modern grid... From CaseDocs/marbl_in, the last block of code shows that MARBL is going to request restoring for five of those six tracers (we don't restore O2):

! tracer restoring
tracer_restore_vars(1) = "PO4"
tracer_restore_vars(2) = "NO3"
tracer_restore_vars(3) = "SiO3"
tracer_restore_vars(4) = "ALK"
tracer_restore_vars(5) = "ALK_ALT_CO2"
tracer_restore_vars(6) = ""

CaseDocs/pop_in tells POP what to do when MARBL wants some tracers restored:

 restore_data_file_varnames = 'PO4','NO3','SiO3','O2','ALK','ALK'
 restore_data_filenames = '/glade/campaign/cesm/cesmdata/inputdata/ocn/pop/gx1v6/forcing/ecosys_restore_POP_
gx1v6_20170113.nc','/glade/campaign/cesm/cesmdata/inputdata/ocn/pop/gx1v6/forcing/ecosys_restore_POP_gx1v6_2
0170113.nc','/glade/campaign/cesm/cesmdata/inputdata/ocn/pop/gx1v6/forcing/ecosys_restore_POP_gx1v6_20170113
.nc','/glade/campaign/cesm/cesmdata/inputdata/ocn/pop/gx1v6/forcing/ecosys_restore_POP_gx1v6_20170113.nc','/
glade/campaign/cesm/cesmdata/inputdata/ocn/pop/gx1v6/forcing/ecosys_restore_POP_gx1v6_20170113.nc','/glade/c
ampaign/cesm/cesmdata/inputdata/ocn/pop/gx1v6/forcing/ecosys_restore_POP_gx1v6_20170113.nc'
 restore_inv_tau_input%file_fmt = 'nc'
 restore_inv_tau_input%file_varname = 'RESTORE_INV_TAU_MARGINAL_SEA_ONLY'
 restore_inv_tau_input%filename = '/glade/campaign/cesm/cesmdata/inputdata/ocn/pop/gx1v6/forcing/ecosys_restore_inv_tau_POP_gx1v6_20170125.nc'
 restore_inv_tau_input%scale_factor = 1.0

Before we go too far down the path of setting up the restoring, could you try a run with the non-zero initial conditions where you add the following to user_nl_marbl:

tracer_restore_vars(1) = ""
tracer_restore_vars(2) = ""
tracer_restore_vars(3) = ""
tracer_restore_vars(4) = ""
tracer_restore_vars(5) = ""

I was wondering why restoring would be necessary if we don't have river runoff flowing into these inland seas, and maybe the answer is "restoring isn't necessary, but we inadvertently turned it on and there might be locations restoring to missing values at a very fast time scale"

 

[jcampbell2548]

Wow.. Well, that worked! I got it to run for 1 year without problem after: setting the initial conditions file back to the one with values, adding the "tracer_restore_vars" lines of code in user_nl_marbl that you suggested, and setting burial coefficients to 0 as well. Thank you so much for your help, I really appreciate it!!!

 

[mlevy]

That's great to hear! Reach out if you run into more issues along the way :)