How can I set up the initial conditions’ tracer values for the MARBL CISO model?

[shiannder]

Hi,

I have a simulation run (jsl.009) with CESM 2.2 compset G1850ECOIAF out to 8866 yrs with carbon isotopes on because initialization of the carbon isotopes (DO13C, DI13C) is initialized too depleted and radiocarbon (DO14C, DI14C) is initialized too young.

After I simulated isotope variables out to steady state conditions, I extracted out 14 of the CISO variables below with units “mmol/m^3” and put them into the original initial condition netcdf file (ecosys_jan_IC_gx3v7_20180308.nc) and renamed it (updated_ecosys_jan_IC_gx3v7_20180308_new_ciso.nc). I made sure to keep same variable names and remove the time dimension to just have latitude, longitude, and depth dimensions.

'DI13C',
'DO13Ctot',
'DI14C',
'DO14Ctot',
'zootot13C',
'zootot14C',
'sp13C',
'sp14C',
'spCa13CO3',
'spCa14CO3',
'diat13C',
'diat14C',
'diaz13C',
'diaz14C'


In my user_nl_pop_ciso file, I put (init_ecosys_init_file = '/mnt/lustre/letscher/jsl1063/initcond/ciso8866y/updated_ecosys_jan_IC_gx3v7_20180308_new_ciso.nc') to call the right initial condition file.

I also start the run with (cp user_nl_pop_ciso ${SCRIPT_DIR_BASE}/${CASE}/user_nl_pop) to call the right user_nl_pop_ciso file.

In my new model run, the ocn log file for some reason initializes the CISO variables with non CISO variables (see below and attached). So the model grabs DIC and not DI13C, DOCtot and not DO13Ctot, etc.

How can I set up/add the initial conditions’ tracer values for the MARBL CISO model properly? Thank you for the support. I can also send my complete ocn log file and initial condition netcdf file if it helps.

 

[mlevy]

The way POP initializes tracers is a bit confusing; there is a namelist variable ciso_tracer_init_ext that is a derived type and instructs POP on how to initialize the CISO tracers. In the default case,

 ciso_tracer_init_ext(1)%file_varname = 'DIC'
 ciso_tracer_init_ext(1)%mod_varname = 'DI13C'
 ciso_tracer_init_ext(1)%scale_factor = 1.025
 ciso_tracer_init_ext(2)%file_varname = 'DOCtot'
 ciso_tracer_init_ext(2)%mod_varname = 'DO13Ctot'
 ciso_tracer_init_ext(2)%scale_factor = 1.0
 ciso_tracer_init_ext(3)%file_varname = 'zooC'
 ciso_tracer_init_ext(3)%mod_varname = 'zootot13C'
 ciso_tracer_init_ext(3)%scale_factor = 1.0
 ciso_tracer_init_ext(4)%file_varname = 'DIC'
 ciso_tracer_init_ext(4)%mod_varname = 'DI14C'
 ciso_tracer_init_ext(4)%scale_factor = 0.9225
 ciso_tracer_init_ext(5)%file_varname = 'DOCtot'
 ciso_tracer_init_ext(5)%mod_varname = 'DO14Ctot'
 ciso_tracer_init_ext(5)%scale_factor = 1.0
 ciso_tracer_init_ext(6)%file_varname = 'zooC'
 ciso_tracer_init_ext(6)%mod_varname = 'zootot14C'
 ciso_tracer_init_ext(6)%scale_factor = 1.0
 ciso_tracer_init_ext(7)%file_varname = 'spC'
 ciso_tracer_init_ext(7)%mod_varname = 'sp13C'
 ciso_tracer_init_ext(7)%scale_factor = 1.0
 ciso_tracer_init_ext(8)%file_varname = 'spC'
 ciso_tracer_init_ext(8)%mod_varname = 'sp14C'
 ciso_tracer_init_ext(8)%scale_factor = 1.0
 ciso_tracer_init_ext(9)%file_varname = 'spCaCO3'
 ciso_tracer_init_ext(9)%mod_varname = 'spCa13CO3'
 ciso_tracer_init_ext(9)%scale_factor = 1.0
 ciso_tracer_init_ext(10)%file_varname = 'spCaCO3'
 ciso_tracer_init_ext(10)%mod_varname = 'spCa14CO3'
 ciso_tracer_init_ext(10)%scale_factor = 1.0
 ciso_tracer_init_ext(11)%file_varname = 'diatC'
 ciso_tracer_init_ext(11)%mod_varname = 'diat13C'
 ciso_tracer_init_ext(11)%scale_factor = 1.0
 ciso_tracer_init_ext(12)%file_varname = 'diatC'
 ciso_tracer_init_ext(12)%mod_varname = 'diat14C'
 ciso_tracer_init_ext(12)%scale_factor = 1.0
 ciso_tracer_init_ext(13)%file_varname = 'diazC'
 ciso_tracer_init_ext(13)%mod_varname = 'diaz13C'
 ciso_tracer_init_ext(13)%scale_factor = 1.0
 ciso_tracer_init_ext(14)%file_varname = 'diazC'
 ciso_tracer_init_ext(14)%mod_varname = 'diaz14C'
 ciso_tracer_init_ext(14)%scale_factor = 1.0

I think the best way to read in the data from the initial condition file you created is to set

 ciso_tracer_init_ext(N)%file_varname = 'unknown'
 ciso_tracer_init_ext(N)%mod_varname = 'unknown'
 ciso_tracer_init_ext(N)%scale_factor = 1.0

in user_nl_pop for N in 1, 2, ..., 14. If that doesn't work, I can show you how to comment out the part of the Fortran code responsible for making these changes but that's a riskier fix (if you create a new case that reverts back to the default initial condition file, it won't be able to initialize the ciso tracers)

 

[shiannder]

Thanks Mike for the support,

Could you clarify what and where is the file that contains the nameslist variable "ciso_tracer_init_ext" that I need to change to match your code for
"ciso_tracer_init_ext(N)%file_varname = 'unknown' "?

Is it suppose to be written/pasted within the user_nl_pop file with the for loop?

In my user_nl_pop file, I currently have

• init_ecosys_init_file = '/mnt/lustre/letscher/jsl1063/initcond/ciso8866y/updated_ecosys_jan_IC_gx3v7_20180308_new_ciso.nc'

to reference the new initialization file. Am I starting with the new ciso_tracer_init_ext (N)%... here?


Also, does the for loop you mentioned in user_nl_pop, suppose to be a do loop for Fortran? I am also confused on how to write the for loop too.


For N in 1,2,3,4,5,6,7,8,9,10,11,12,13,14:
"what goes here...?"

 

[mlevy]

The lines should go in user_nl_pop; it sounds like you want

init_ecosys_init_file = '/mnt/lustre/letscher/jsl1063/initcond/ciso8866y/updated_ecosys_jan_IC_gx3v7_20180308_new_ciso.nc'
ciso_tracer_init_ext(1)%file_varname = 'unknown'
ciso_tracer_init_ext(1)%mod_varname = 'unknown'
ciso_tracer_init_ext(1)%scale_factor = 1.0
ciso_tracer_init_ext(2)%file_varname = 'unknown'
ciso_tracer_init_ext(2)%mod_varname = 'unknown'
ciso_tracer_init_ext(2)%scale_factor = 1.0

and so on through

ciso_tracer_init_ext(14)%file_varname = 'unknown'
ciso_tracer_init_ext(14)%mod_varname = 'unknown'
ciso_tracer_init_ext(14)%scale_factor = 1.0

(Don't actually put a "for" loop in any files, it was just a shorthand to say "you'll need to add 42 lines, but it's 14 copies of the same 3 lines with one small change from set to set)

 

[shiannder]

Hi Mike,

Thank you for the clarification.

I ran a new model (jsl.041) with the suggested code.

I got DependencyNeverSatisfied, so the model never ran.

I looked into the cesm.log.18471.231214-144151 file and found some errors attached. Do you think my netcdf file has values that are not compatible with the code?

 

[mlevy]

I think you might have a problem with overflows trying to read data that is masked out in your data set. Notice that one of your error messages is

(Task 4, block 1) Message from (lon, lat) ( 172.100, -77.981), which is global (i,j) (59, 2). Level: 1
(Task 4, block 1) MARBL ERROR (marbl_ciso_diagnostics_mod:store_diagnostics_ciso_interior): abs(CISO_Jint_13Ctot)= 0.461E+019 exceeds CISO_Jint_13Ctot_thres= 0.317E-011
(Task 4, block 1) MARBL ERROR (marbl_ciso_interior_tendency_mod:marbl_ciso_interior_tendency_compute): Error reported from store_diagnostics_ciso_interior
(Task 4, block 1) MARBL ERROR (marbl_interior_tendency_mod:marbl_interior_tendency_compute): Error reported from marbl_ciso_interior_tendency_compute()
(Task 4, block 1) MARBL ERROR (marbl_interface:interior_tendency_compute): Error reported from marbl_interior_tendency_compute()
(Task 4, block 1) MARBL ERROR (ecosys_driver:ecosys_driver_set_interior): Error reported from marbl_instances(1)%set_interior_forcing()
ERROR reported from MARBL library

If I compare the DIC field at this location to to the DI13C field, I can see that there is one more level of DIC:

>>> import numpy as np
>>> import xarray as xr
>>> ds = xr.open_dataset('updated_ecosys_jan_IC_gx3v7_20180308_new_ciso.nc')
>>> ds.isel(nlon=58,nlat=1).DIC.data
array([2172.043 , 2174.3196, 2182.2224, 2196.0862, 2214.3633, 2228.3574,
       2238.0679, 2247.7786, 2251.0413, 2254.3037, 2256.7415, 2258.3542,
       2259.9668, 2262.203 , 2264.439 , 2266.7336, 2269.1106, 2271.5542,
       2274.0986, 2276.7744, 2278.7842, 2280.7173, 2282.7952, 2284.9524,
       2286.3242, 2287.8264, 2289.2239, 2289.9646, 2290.7893, 2291.7134,
       2292.8425, 2294.371 , 2296.11  , 2298.3977, 2298.3977,       nan,
             nan,       nan,       nan,       nan,       nan,       nan,
             nan,       nan,       nan,       nan,       nan,       nan,
             nan,       nan,       nan,       nan,       nan,       nan,
             nan,       nan,       nan,       nan,       nan,       nan],
      dtype=float32)
>>> np.sum(np.isfinite(ds.isel(nlon=58,nlat=1).DIC.data))
35
>>> ds.isel(nlon=58,nlat=1).DI13C.data
array([2325.9336, 2327.2986, 2330.7856, 2335.0076, 2338.969 , 2343.0322,
       2345.9392, 2347.8699, 2349.096 , 2349.8875, 2350.3806, 2350.686 ,
       2350.8762, 2351.0125, 2351.118 , 2351.1997, 2351.2676, 2351.33  ,
       2351.3848, 2351.44  , 2351.4944, 2351.5464, 2351.599 , 2351.6619,
       2351.7332, 2351.813 , 2351.915 , 2352.0557, 2352.1802, 2352.3447,
       2352.5503, 2352.8428, 2353.2192, 2353.6404,       nan,       nan,
             nan,       nan,       nan,       nan,       nan,       nan,
             nan,       nan,       nan,       nan,       nan,       nan,
             nan,       nan,       nan,       nan,       nan,       nan,
             nan,       nan,       nan,       nan,       nan,       nan],
      dtype=float32)
>>> np.sum(np.isfinite(ds.isel(nlon=58,nlat=1).DI13C.data))
34

The simplest fix is probably to extend the CISO initial conditions such that the KMT+1 level is equal to the value at level KMT instead of NaN (for all columns where KMT is not 0 (land cell) or 60 (already full depth).

 

[shiannder]

Hi Mike,

Sorry for the late reply. I hope you had a great Winter holiday.

I see there are certain values within CISO_Jint_13Ctot and CISO_Jint_14Ctot that are NaN at specific levels which throws out an error and aborts the run.

If the error is from a missing value at a specific level in the IC netcdf file, where would you suggest I make this extension to the CISO initial conditions (in the netcdf file itself)? Would I prescribe all the NaN values to a certain value?

What does KMT stand for?

 

[mlevy]

KMT is a POP variable denoting the number of active levels in a given column. Over land, KMT=0 (there are no ocean levels), otherwise it varies between 3 (the shallowest POP columns are 30m deep) and 60 (because the vertical grid only has 60 levels). So I think the problem is that there are a few points on the grid where the overflow parameterization needs to know values at KMT+1, even though that typically denotes a vertical level that is completely below the sea floor. Using the column in my previous comment as an example, if KMT=34 then POP is only going to write 34 levels of output for that column but the initial condition needs to be defined at 35 levels.

My recommendation is to replace all NaN values within the netcdf file with the bottom-most non-NaN. In python (with xarray as numpy) it might look something like

for k in range(3,60):
  ds[tracer_name].data[:,:,k] = np.where(np.isfinite(ds[tracer_name].data[:,:,k]), ds[tracer_name].data[:,:,k], ds[tracer_name].data[:,:,k-1])

This should greatly reduce the number of NaNs in the dataset; any column with active ocean will have initial conditions defined at all depths instead of just down to where POP thinks the sea floor is.