NOAA-CEFI-Regional-Ocean-Modeling / ocean_BGC

3 stars 4 forks source link

Add support for a COBALT_param_doc file #59

Closed theresa-morrison closed 1 month ago

theresa-morrison commented 1 month ago

Using the file parser from MOM6, add the option to read parameters from COBALT_input and COBALT_override files and write parameters to COBALT_param_doc files.

With these changes, COBALT_input and COBALT_override files must be listed in the input namelist. As a result all existing experiments will need their namelists updated and COBALT_input files generated. The files can be blank and the default values will be used.

&cobalt_input_nml
             parameter_filename = 'INPUT/COBALT_input','INPUT/COBALT_override'
/

Descriptions of the parameters have not been added with this PR and some parameters lack units. For the description, the parameter name has been used but should be updated to include a note about what the parameter is used for and a recommended range (or that the parameter not be adjusted).

As @amoebaliz suggested, here are the needed changes to make to an xml to use the parameter doc file.

  1. In the <input> section for each experiment or for the primary experiment, add the new namelist option:
    <namelist name = cobalt_input_nml >
    parameter_filename = 'INPUT/COBALT_input','INPUT/COBALT_override'
    </namelist>
  2. You can create a COBALT_input file in the same location as your MOM_input file and link to is in a similar way, for example:
    <dataSource site="ncrc">$(ESMG_H)/ESMG-configs/Arctic6/mid_12/MOM_input</dataSource>
    <dataSource site="ncrc">$(ESMG_H)/ESMG-configs/Arctic6/mid_12/COBALT_input</dataSource>
  3. Alternatively you can write a COBALT_input file directly from the xml. Toward the end of the input section, where the MOM_override file is defined in you can add:
    
    <csh type='always'>
    <![CDATA[
    truncate -s 0 $work/INPUT/MOM_override
    cat > $work/INPUT/MOM_override << MOM_OVERRIDE_EOF
    MOM_OVERRIDE_EOF

truncate -s 0 $work/INPUT/COBALT_input cat > $work/INPUT/COBALT_input << COBALT_INPUT_EOF ! Put the contents of your COBALT_input file here COBALT_INPUT_EOF ]]>

This is how and where you can add a COBALT_override file

truncate -s 0 $work/INPUT/COBALT_override cat > $work/INPUT/COBALT_override << COBALT_OVERRIDE_EOF ! Put the contents of your COBALT_override file here. ! For parameters which are not specified in the COBALT_input file, just set the new value, ex: ! do_case2_mod = True ! For parameters which are specified in the COBALT_input file, specify an override, ex: ! #override do_case2_mod = True COBALT_OVERRIDE_EOF


5. Finally, you can update the post processing section to include the COBALT_parameter_doc files.
<postProcess>
  <csh><![CDATA[
     cd $work
     #Make a directory to trick FRE to pick up and archive in ascii
     mkdir extra.results
     mv *velocity_truncations MOM_parameter_doc* SIS_parameter_doc* COBALT_parameter_doc* seaice.stats* ocean.stats* extra.results/
     cp diag_table extra.results/.
     cp $scriptName extra.results/runscript.csh
     #When the ocean uses a mask_table the ocean_static.nc file produced by the model run has holes in coordinates (geolon,geolat)
     #This causes problems for ferret and python tools for analysis.
     #Copy a non-masked version of ocean_static.nc to be saved as a history file to be used by the analysis scipts .
     #cp $work/INPUT/ocean_static_no_mask_table.nc $work/
  ]]></csh>
  <xi:include xpointer="xpointer(//freInclude[@name='MOM6_postprocess']/postProcess/node())"/>
</postProcess>
yichengt900 commented 1 month ago

Thanks @theresa-morrison . I have updated the github actions testing workflow to accommodate this PR. Please update your feature branch, and it should be able to pass the RT.

I can confirm that this new functionality works properly. For example, if I set remin_ramp_scale = 100.0 in COBALT_input and use COBALT_override to revert it to the default value, the model results are identical to the reference.

@charliestock and anyone else, it would be great if you could double-check the changes in generic_tracers/generic_COBALT.F90 and ensure you are comfortable with these updates.

theresa-morrison commented 1 month ago

As @amoebaliz suggested, here are the needed changes to make to an xml to use the parameter doc file.

  1. In the <input> section for each experiment or for the primary experiment, add the new namelist option:
    <namelist name = cobalt_input_nml >
    parameter_filename = 'INPUT/COBALT_input','INPUT/COBALT_override'
    </namelist>
  2. You can create a COBALT_input file in the same location as your MOM_input file and link to is in a similar way, for example:
    <dataSource site="ncrc">$(ESMG_H)/ESMG-configs/Arctic6/mid_12/MOM_input</dataSource>
    <dataSource site="ncrc">$(ESMG_H)/ESMG-configs/Arctic6/mid_12/COBALT_input</dataSource>
  3. Alternatively you can write a COBALT_input file directly from the xml. Toward the end of the input section, where the MOM_override file is defined in you can add:
    
    <csh type='always'>
    <![CDATA[
    truncate -s 0 $work/INPUT/MOM_override
    cat > $work/INPUT/MOM_override << MOM_OVERRIDE_EOF
    MOM_OVERRIDE_EOF

truncate -s 0 $work/INPUT/COBALT_input cat > $work/INPUT/COBALT_input << COBALT_INPUT_EOF ! Put the contents of your COBALT_input file here COBALT_INPUT_EOF ]]>

This is how and where you can add a COBALT_override file

truncate -s 0 $work/INPUT/COBALT_override cat > $work/INPUT/COBALT_override << COBALT_OVERRIDE_EOF ! Put the contents of your COBALT_override file here. ! For parameters which are not specified in the COBALT_input file, just set the new value, ex: ! do_case2_mod = True ! For parameters which are specified in the COBALT_input file, specify an override, ex: ! #override do_case2_mod = True COBALT_OVERRIDE_EOF


5. Finally, you can update the post processing section to include the COBALT_parameter_doc files.
<postProcess>
  <csh><![CDATA[
     cd $work
     #Make a directory to trick FRE to pick up and archive in ascii
     mkdir extra.results
     mv *velocity_truncations MOM_parameter_doc* SIS_parameter_doc* COBALT_parameter_doc* seaice.stats* ocean.stats* extra.results/
     cp diag_table extra.results/.
     cp $scriptName extra.results/runscript.csh
     #When the ocean uses a mask_table the ocean_static.nc file produced by the model run has holes in coordinates (geolon,geolat)
     #This causes problems for ferret and python tools for analysis.
     #Copy a non-masked version of ocean_static.nc to be saved as a history file to be used by the analysis scipts .
     #cp $work/INPUT/ocean_static_no_mask_table.nc $work/
  ]]></csh>
  <xi:include xpointer="xpointer(//freInclude[@name='MOM6_postprocess']/postProcess/node())"/>
</postProcess>
theresa-morrison commented 1 month ago

With the latest commit, parameters with default values that included a conversion from per-day or per-year ito per-second in the default have been updated so that there default is the per-day or per-year value.

These value are better understood in the community which will make it easier for users to adjust them.

In addition, avoiding division in calculation of the default is important for producing accurate values in the parameter doc files. Otherwise, when the rounded values from the parameter doc file are specified rather than the default the difference changes the model results.

This commit did not impact the parameters with calculations in the default that appeared to be related to stoichiometry. These parameters should be clearly labeled in their descriptions possibly with a warning that users should not adjust these values.

theresa-morrison commented 1 month ago

There are two relevant options for get_param that we can use here: do_not_read : If present and true, do not read a value for this parameter, although it might be logged. scale : A scaling factor that the parameter is multiplied by before it is returned.

do_not read would be good for the stoichiometry parameters that users may want to make sure are consistent, but should not be able to easily (or accidentally) change.

scale is similar to what Xiao suggested in the documentation meeting yesterday, and would be a good option for the time unit conversions that some of the parameters have.

charliestock commented 1 month ago

Thanks again for the great work on this pull request Theresa. It is an exciting development!

The flexibility offered by "do_not_read" and "scale" sounds great. For "scale" it seems that we have a limited number of primary conversions that may make this approach tractable:

1) Time: My inclination is to require all inputs in "days" or "days-1" and then convert to "seconds" or "seconds-1" using sperd = 86400. There are some parameters that we have converted from "year" or "year-1" in the past, but I don't think there are many of these and it may be easiest to limit the number. Thoughts?

2) Concentrations: We could have the input unit be micromoles kg-1 and use a uniform conversion micromol2mol = 1e6 to convert?

3) Nutrient ratios: There are a number of ratios that are commonly reported in the literature relative to carbon (i.e., iron to carbon) that are converted into ratios relative to nitrogen for the model. A conversion c2n = 106.0/16.0 would do the trick for these.

4) Light: Most of the photosynthesis literature reports parameters, such as the initial slope of photosynthesis-irradiance curve or the compensation irradiance in "micromoles quanta" of photosynthetically available radiation. We use a conversion from Smith and Morel (1974) to convert this to watts: under typical conditions 2.77e18 quanta sec-1 = a watt. So there are 2.77e18/6.022e17 = 4.6 micromole quanta per sec-1 per Watt (note that 6.022e17 is a micromole). How about something like micromolQpersec2W = 2.77e18/6.022e17? Any suggestions on how to shorten?

I may have missed a few, but I think this would cover most of them.

Lastly, I'm wondering when the best time to accept this initial pull request is? It seems that we'd want to accept the foundation of the pull request ASAP so that we can begin work on filling in parameter details, etc. Do folks think we should pull things over before "scaling" settings and then work to refine and rescale together? Or is it worth implementing/testing the rescaling first?

Have we confirmed that this doesn't change answers for a more complete run like NWA12?

Again - excellent work Theresa and all who have contributed!

Charlie

yichengt900 commented 1 month ago

Thanks, @theresa-morrison, for this great work, and thank you, @Charlie, for the detailed and thoughtful comments. I agree with your points and appreciate the clarity on the primary conversions.

I have a minor comment regarding the zooplankton swimming parameters that could help us maintain answers (it's not activate at this moment but the get_param call will lock the parameter so it's better to use different parameter names for each groups) in the NWA12 case.

Regarding the do_case2_mod=false setting, do we prefer this as the default?

I also support accepting this PR sooner rather than later. We can revisit the rescaling in a separate PR. This would allow other developers to benefit from the enhancements and begin working with the new functionality immediately.

Great work, @theresa-morrison, and thanks to everyone involved!

theresa-morrison commented 1 month ago
1. Time: My inclination is to require all inputs in "days" or "days-1" and then convert to "seconds" or "seconds-1" using sperd = 86400.  There are some parameters that we have converted from "year" or "year-1" in the past, but I don't think there are many of these and it may be easiest to limit the number.  Thoughts?

I agree with this, there are about 8 spery conversions. I think it would be fine to keep those around for now.

2. Concentrations: We could have the input unit be micromoles kg-1 and use a uniform conversion micromol2mol = 1e6 to convert?

I like this because it will make it very clear what the 1e6 is for.

3. Nutrient ratios: There are a number of ratios that are commonly reported in the literature relative to carbon (i.e., iron to carbon) that are converted into ratios relative to nitrogen for the model.  A conversion c2n = 106.0/16.0 would do the trick for these.

I can do this, I will just search for "106.0/16.0" which should get most instances.

4. Light: Most of the photosynthesis literature reports parameters, such as the initial slope of photosynthesis-irradiance curve or the compensation irradiance in "micromoles quanta" of photosynthetically available radiation.  We use a conversion from Smith and Morel (1974) to convert this to watts: under typical conditions 2.77e18 quanta sec-1 = a watt.  So there are 2.77e18/6.022e17 = 4.6 micromole quanta per sec-1 per Watt (note that 6.022e17 is a micromole).  How about something like micromolQpersec2W = 2.77e18/6.022e17?  Any suggestions on how to shorten?

I think micromolQpersec2W is clear and I can include a reference to Smith and Morel (1974) and some of the description here where the parameter is defined. These scaling factors will only be in the parameter read statements so it is okay if they are long.

Lastly, I'm wondering when the best time to accept this initial pull request is? It seems that we'd want to accept the foundation of the pull request ASAP so that we can begin work on filling in parameter details, etc. Do folks think we should pull things over before "scaling" settings and then work to refine and rescale together? Or is it worth implementing/testing the rescaling first?

Have we confirmed that this doesn't change answers for a more complete run like NWA12?

I can make the changes and do some basic testing on my workstation today. I have done short tests with the NWA12 but I would recommend we run ~1 year in case there are parameters that aren't active in the first month or few days of the year. I have only been looking at the ocean stats file, but we should look at the COBALT values in the standard out file as well.

I also agree with @yichengt900 that we can put the rescaling in a separate PR to expedite this. And I will change do_case2_mod=true back to the default.

charliestock commented 1 month ago

OK, I like the plan of including the rescaling as a second line of code in the initial pull request, and seeing if we can take these out in a subsequent one. We can discuss whether we want to the case2 modifications to be the default today. I actually lean a bit toward making the default "false" since this is the case that is most consistent with the physics.

yichengt900 commented 1 month ago

@theresa-morrison, unfortunately, your last commit broke the automated testing again. I will have to take a look to determine why the tests are not happy with the latest changes.

theresa-morrison commented 1 month ago

@yichengt900 it looks like the variables that use "do_not_read" are being set to 0.0, I think that is probably causing the issue with the testing.

yichengt900 commented 1 month ago

@theresa-morrison, thanks. I think this makes sense because I keep getting FATAL: NaN in input field of reproducing_EFP_sum(_2d) in the automated testing with your last commit. Also, I noticed that once you introduced do_not_read, it changed the baseline answers. I suggest we revert the changes related to do_not_read and test it again.