Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

1D Ocean Column MARBL Interface + Obs Converter #726

Merged
merged 38 commits into from
Sep 10, 2024

Conversation

mgharamti
Copy link
Contributor

@mgharamti mgharamti commented Aug 29, 2024

Description:

A new model interface to the biogeochemistry model MARBL. The work done here is contributed by Robin Armstrong. Thanks Robin!

Fixes issue

This address all comments, suggestions, and questions raised on the original PR MARBL-DART #701. It also fixes issue #614

Types of changes

  • Bug fix (non-breaking change which fixes an issue)
  • New feature (non-breaking change which adds functionality)
  • Breaking change (fix or feature that would cause existing functionality to not work as expected)
  • Documentation update

Documentation changes needed?

  • My change requires a change to the documentation.
    • I have updated the documentation accordingly.

Tests

Ran several DART program to assess the behavior: model_mod_check, filter, and bats_to_obs
My tests are here:

  • Model Interface: /glade/work/gharamti/DART/bgc_work/models/MARBL_column/work
  • Obs Converter: /glade/work/gharamti/DART/bgc_work/observations/obs_converters/BATS

Checklist for merging

  • Updated changelog entry
  • Documentation updated
  • Update conf.py

Checklist for release

  • Merge into main
  • Create release from the main branch with appropriate tag
  • Delete feature-branch

Testing Datasets

  • Dataset needed for testing available upon request
  • Dataset download instructions included
  • No dataset needed

robin-armstrong and others added 30 commits July 15, 2024 10:56
…n models/MARBL_MOM6_1D in branch marbl_mom6_1d) and MARBL_joint_estimation (developed in models/MARBL_MOM6_1D on branch param_estimation).
…imation branch. Also added in the observation codes for BATS variables, which were accidentally left out of the last commit.
… This required fixing a bug in direct_netcdf_mod, using a solution provided by Helen Kershaw.
…d. Split the BATS data converter into two separate converters, one for sequential state or state-parameter estimation, and another for pure parameter estimation.
… cycles with MARBL without errors. Further testing needed to assess the quality of assimilations.
…imatology.py now records observation error in units of standard deviation, and the script also generates a surrogate netCDF file containing the BATS empirical climatology.
…t, so that it works with the self-contained test package.
it is actually a netcdf file
add executables to root/.gitignore
was pointing at Robin's directory
MARBL_joint_estimation is the model_mod we are using (see pull NCAR#701)
Replaced model directory 'MARBL_joint_estimation' with
'MARBL_column' and compiled to make sure the code
works
Removed the hard-coded location of BATS station
and added it as `station_location` in the input.nml
This allows the users to run the code for any location
on earth (i.e., it doesn't have to be Bermuda station).

I also cleaned up the input.nml and stripped out the
old path names to Robin's directories.
The model size is now computed based on the number
of active domains. Also, the accessor function
get_model_size can now be used from outside the
model_mod without the need for local variables.

The code as it stands should be able to run:
- state estimation only
- parameter estimation only
- state and parameters estimation
Removed unused routines and unnecessary formatting
I removed the QTYs associated with specific parameters
and added a single definition: 'QTY_BGC_PARAM'
The changes are now reflected in the
default_quantities_mod

Did some cleaning in the model-mod and changed
some configuration in the input.nml file
to be able to run the interpolation test
I went through all model-mod checks and made sure all
tests are successful. The commit also includes some
more cleaning for the model-mod and the input.nml
Updated the BATS station location in the obs_converter.
Previously it assumed a lon of 64 while it should be
64W or 296.
Also, the BATS text data file has been modified and
that required changing the input parameters for
reading the quantities properly.

filter has been run several times and obs_seq.final
was assessed with no apparent issues.
Added information in the readme.rst file
Removed the old readme.md file and added
a new rst readme file describing the converter
and its namelist options.
Added the scripts and code to generate climatology data
to the BATS converter. This includes a python script and
fortran code. The namelists were also consolidated.
The BATS_clim converter has been removed.

Documentation of both the model and the converter
was also modified for clarity.
Copy link
Member

@hkershaw-brown hkershaw-brown left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Moha,
Thanks for this. I haven't looked at the code yet.
There are some warnings in the docs build.

(py-docs) [hkershaw:DOCS]() > make html
Running Sphinx v6.2.1
/Users/hkershaw/DART/DOCS/py-docs/lib/python3.9/site-packages/urllib3/__init__.py:35: NotOpenSSLWarning: urllib3 v2 only supports OpenSSL 1.1.1+, currently the 'ssl' module is compiled with 'LibreSSL 2.8.3'. See: https://github.com/urllib3/urllib3/issues/3020
  warnings.warn(
building [mo]: targets for 0 po files that are out of date
writing output... 
building [html]: targets for 237 source files that are out of date
updating environment: [new config] 237 added, 0 changed, 0 removed
reading sources... [100%] theory/readme                                                                                                                       
/Users/hkershaw/DART/pull_requests/pull_726/models/MARBL_column/readme.rst:28: ERROR: Unexpected indentation.
/Users/hkershaw/DART/pull_requests/pull_726/models/MARBL_column/readme.rst:34: ERROR: Unexpected indentation.
/Users/hkershaw/DART/pull_requests/pull_726/models/MARBL_column/readme.rst:45: ERROR: Unexpected indentation.
/Users/hkershaw/DART/pull_requests/pull_726/models/MARBL_column/readme.rst:46: WARNING: Block quote ends without a blank line; unexpected unindent.
/Users/hkershaw/DART/pull_requests/pull_726/observations/obs_converters/BATS/readme.rst:20: ERROR: Unexpected indentation.
/Users/hkershaw/DART/pull_requests/pull_726/observations/obs_converters/BATS/readme.rst:56: ERROR: Unexpected indentation.
/Users/hkershaw/DART/pull_requests/pull_726/observations/obs_converters/BATS/readme.rst:79: WARNING: Inline literal start-string without end-string.
/Users/hkershaw/DART/pull_requests/pull_726/observations/obs_converters/BATS/readme.rst:97: WARNING: Title underline too short.

Climatology:
~~~~~~~~~~~
looking for now-outdated files... none found
pickling environment... done
checking consistency... /Users/hkershaw/DART/pull_requests/pull_726/models/MARBL_column/readme.rst: WARNING: document isn't included in any toctree
/Users/hkershaw/DART/pull_requests/pull_726/observations/obs_converters/BATS/readme.rst: WARNING: document isn't included in any toctree
done
preparing documents... done
writing output... [100%] theory/readme                                                                                                                        
generating indices... genindex done
writing additional pages... search done
copying images... [100%] guide/images/WRF_tutorial_evolution2.png                                                                                             
copying downloadable files... [100%] ../DOCS/guide/DART_LAB/presentation/DART_LAB_Section05.pdf                                                               
copying static files... done
copying extra files... done
dumping search index in English (code: en)... done
dumping object inventory... done
build succeeded, 10 warnings.

The HTML pages are in dart-docs/html.

https://readthedocs.org/api/v2/build/25474391.txt

In addition to adding to the toc:

The converter should be added to 'converter programs' https://docs.dart.ucar.edu/en/latest/observations/obs_converters/README.html#converter-programs

The model should be added to 'supported models'
https://docs.dart.ucar.edu/en/latest/models/README.html

Improved text in the documentation.
New model interface and obs converter are now
added to the toctree and can be accessed on docs.dart.edu
@hkershaw-brown
Copy link
Member

a question/note on the parameter estimation only.

This is the docs:

"Parameters Estimation only: where only the parameters are constrained using the data. DART will still need to read in the state to construct ensemble covariances and compute innovations. The only difference is that the ensemble increments are only regressed onto the unknown parameters. To achive this goal, you’ll need to set estimate_params = .true. and turn the update status for all state variable to NO_COPY_BACK

This ensures that the state will not be updated. The NO_COPY_BACK option is added as the 5th entry in the state table (after the variable name, its associated quantity and its physical bounds) within &model_nml"

NO_COPY_BACK means the state is not written out to a restart file. The state is still updated by the assimilation, and inflated.

Copy link
Member

@hkershaw-brown hkershaw-brown left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Next batch of review comments for model_mod.

Couple of definite changes:

  • namelist items need default values
  • param file and state file template files, since you have separate namelist options for turning these on and off, I think it would be better to have a param_template_file and a state_template_file
  • model_size is i8 in dart
  • for get_close_obs|state this is just passing through to the location mod, so you can use the location mod code directly.

Alphabetical order / style notes:

A question on the interpolation, there is no horizontal interpolation so just confirming that that is correct?

The model_time is set for days from 1982. Is this true for all MARBL runs?

There are some other bits and bobs of suggestions for unused variables and old print statements.

Comment on lines 7 to 9
! This is a template showing the interfaces required for a model to be compliant
! with the DART data assimilation infrastructure. Do not change the arguments
! for the public routines.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not the template model

Either remove this comment or edit it, e.g.

Suggested change
! This is a template showing the interfaces required for a model to be compliant
! with the DART data assimilation infrastructure. Do not change the arguments
! for the public routines.
! This is the MARBL_columun model_mod for the DART data assimilation infrastructure.
! Do not change the arguments for the public routines.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

Comment on lines 91 to 94
integer, parameter :: modelvar_table_height = 13
integer, parameter :: modelvar_table_width = 5
integer, parameter :: modelparams_table_height = 1
integer, parameter :: modelparams_table_width = 5
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

parameters (constants) CAPITAL LETTERS for the names

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

Comment on lines 120 to 123
! Called to do one time initialization of the model. As examples,
! might define information about the model size or model timestep.
! In models that require pre-computed static data, for instance
! spherical harmonic weights, these would also be computed here.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
! Called to do one time initialization of the model. As examples,
! might define information about the model size or model timestep.
! In models that require pre-computed static data, for instance
! spherical harmonic weights, these would also be computed here.
! Called to do one time initialization of the model.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

type(time_type) :: mom6_time
integer :: mom_base_date_in_days, mom_days

mom_base_date_in_days = 139157 ! 1982 1 1 0 0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this date always 1982 1 1 0 0 for MARBL?
For MOM6 the data formate is in days from year 1

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I believe this the time format for MARBL. I think it has to do with dates related to obs availability. I recall Robin spending some time on this before figuring out the correct format.

Comment on lines 292 to 295
!print *, "BASIN DEPTH = ",-basin_depth(1,1)
!print *, "----------------------------------------------------"
!print *, "computing centers of layers"
!print *, "----------------------------------------------------"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
!print *, "BASIN DEPTH = ",-basin_depth(1,1)
!print *, "----------------------------------------------------"
!print *, "computing centers of layers"
!print *, "----------------------------------------------------"

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

integer :: param_dom_id ! used to access MARBL internal parameters
integer :: nfields ! number of fields in the state or parameter vector
integer :: nz ! the number of vertical layers
integer :: model_size
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

model_size is an i8 integer in dart.

Suggested change
integer :: model_size
integer(i8) :: model_size

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

end do


! performing the interpolation for each ensemble member individually.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a question/confirmation for my own understanding. There is no horizontal interpolation, in MARBL_column.
All four gridpoints have identical values for a given level so you only do vertical interpolation, and you just pick point (1,1). You could pick (2,2) and the interpolation would give the same answer.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, you're right. All four corners are identical columns so technically you can select any corner. This also means that no horizontal interpolation is needed.

integer, parameter :: modelparams_table_width = 5

! defining the variables that will be read from the namelist
character(len=256) :: template_file(2)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because you have this namelist as 2 values, you're requiring that people give two filenames. So even if they

If you only give one file, but estimate_params = .true., you have you get errors like this:

***************** RUNNING    TEST 0    ***********************
 -- Reading the model_mod namelist and implicitly running static_init_model
**************************************************************
 
 Assimilate_these_obs_types:
    BATS_OXYGEN
    BATS_ALKALINITY
    BATS_INORGANIC_CARBON
    BATS_SILICATE
    BATS_PHOSPHATE
    BATS_NITRATE
 Evaluate_these_obs_types:
    BATS_ORGANIC_CARBON
    BATS_NITROGEN
 Use the precomputed Prior Forward Operators for these obs types:
    none
 
 ERROR FROM:
  source : netcdf_utilities_mod.f90
  routine: load_state_variable_info nf90_open
  message: 
 
 
 : NetCDF: Malformed URL
  message: ...

It would be better to have two namelist options, one for each of state, parameter:
state_template_file = 'state.nc'
parameter_template_file = 'param.nc

In the same way that you hae separate namelist entries for 'model_state_variables' and 'model_parameters'

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

Comment on lines 98 to 102
character(len=256) :: ocean_geometry
real(r8) :: station_location(2)
integer :: time_step_days
integer :: time_step_seconds
logical :: estimate_params
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

all the namelist variables need default values.

otherwise you have the value set to whatever junk is in memory if the namelist value is not set in input.nml

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

+-------------------------------------+--------------------+------------------------------------------------------------+
| Item | Type | Description |
+=====================================+====================+============================================================+
| ``template_file`` | character(len=256) | MARBL restart file including MARBL's prognostic variables |
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

template file is a an array of length 2
see other comments, about making this two namelist options.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

Copy link
Member

@hkershaw-brown hkershaw-brown left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see comment before.

Incorporated all suggestions and comments about
model-mod. Several comments have been removed,
redundant loc routines are omitted, better
namelist parsing, deleted unused variables, ...
@mgharamti
Copy link
Contributor Author

Helen, thanks for the useful comments. I think I addressed all your comments/suggestions and incorporated them in the code. Thanks for clarifying the no-copy-back option. I've edited the text in the documentation to make it more clear. Please let me know if you have any questions.

@hkershaw-brown
Copy link
Member

cool thanks Moha, I'll move on to the obs converter

Copy link
Member

@hkershaw-brown hkershaw-brown left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

next batch of comments on obs_converters.

! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is a relic from svn, no need to have it anymore

Suggested change
! $Id$

Comment on lines 43 to 44
real(r8), parameter :: bats_lon = 360.0_r8 - 64.0_r8
real(r8), parameter :: bats_lat = 31.0_r8
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

CAPITAL LETTERS for parameters

Suggested change
real(r8), parameter :: bats_lon = 360.0_r8 - 64.0_r8
real(r8), parameter :: bats_lat = 31.0_r8
real(r8), parameter :: BATS_LON = 360.0_r8 - 64.0_r8
real(r8), parameter :: BATS_LAT = 31.0_r8

Comment on lines 47 to 50
character(len=256) :: text_input_file, obs_out_dir
integer :: max_lines
real(r8) :: obs_err_var_inflation
logical :: debug
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

need default values for namelist variables

real(r8) :: qc, obs_err
real(r8) :: lat, lon, vert, ovalue

! the uncertainties corresponding to the observations above
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not sure what this comment is referencing

! adapted for more generic use 11 Mar 2010, nancy collins, ncar/image
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine create_3d_obs(lat, lon, vval, vkind, obsv, otype, oerr, day, sec, qc, obs)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there reason to have a create_3d_obs in this converter rather than using the create_3d_obs in
obs_utilities_mod?

https://github.com/NCAR/DART/blob/main/observations/obs_converters/utilities/obs_utilities_mod.f90

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think there is any special reason for that, I've removed the ones in the converter and used the tools in obs_utilities_mod and got the same answer. So, I'm cleaning this up as recommended.


! unit corrections from BATS to MARBL
if((OTYPE_ORDERING(otype_index) == BATS_ALKALINITY) .or. (OTYPE_ORDERING(otype_index) == BATS_INORGANIC_CARBON)) then
ovalue = ovalue*1.026
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

just a comment, 1.026 is a wild unit conversion.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

agreed :)

Comment on lines 93 to 94
os.system("rm -f "+output_txt_path)
output_txt = open(output_txt_path, "w")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this rm -f seems dangerous.

Doesn't open(output_txt_path, "w") overwrite the file if it exists?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed

Comment on lines 98 to 99
os.system("rm -f "+output_nc_path)
output_nc = nc.Dataset(output_nc_path, "w")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same rm -f comment

bats_data = open(input_file_path, "r")
line_number = 0

for line in bats_data.readlines():
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure how much python reviewing you want on this,

with open is recommended rather than open and close.
You can slice with readlines to skip the first nlines.

with open(input_file_path, 'r') as bats_data:
    for line in bats_data.readlines()[first_data_line:]

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like with open more. The code will be modified accordingly

@@ -0,0 +1,156 @@
import os
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

#!/usr/bin/env python3

! QTY_DISSOLVED_ORGANIC_P
! QTY_DISSOLVED_INORGANIC_IRON
! QTY_SURFACE_CHLOROPHYLL
! QTY_LAYER_THICKNESS
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

qty_layer_thickness is in twice in this file

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

edit: there are a few duplicates in here:

Duplicate lines containing 'QTY' found:
!   QTY_ALKALINITY (appears 2 times)
!   QTY_DISSOLVED_INORGANIC_CARBON (appears 2 times)
!   QTY_DISSOLVED_INORGANIC_IRON (appears 2 times)
!   QTY_DISSOLVED_INORGANIC_SIO3 (appears 2 times)
!   QTY_DISSOLVED_ORGANIC_CARBON (appears 2 times)
!   QTY_DISSOLVED_ORGANIC_NITROGEN (appears 2 times)
!   QTY_DISSOLVED_ORGANIC_P (appears 2 times)
!   QTY_DISSOLVED_OXYGEN (appears 2 times)
!   QTY_LAYER_THICKNESS (appears 2 times)
!   QTY_MESOZOOPLANKTON_CARBON (appears 2 times)
!   QTY_MICROZOOPLANKTON_CARBON (appears 2 times)
!   QTY_NITRATE_CONCENTRATION (appears 2 times)
!   QTY_PHOSPHATE_CONCENTRATION (appears 2 times)
!   QTY_PHYTOPLANKTON_BIOMASS (appears 2 times)
!   QTY_SURFACE_CHLOROPHYLL (appears 2 times)

@hkershaw-brown
Copy link
Member

@hkershaw-brown todo: test netcdf IO change

Did a cleanup of the obs converter:
- Removed comments
- Removed unused variables
- Improved style
- Better use of DART's modules
- Improved initialization and functions
- Ran and tested the converter after the changes were made
@mgharamti
Copy link
Contributor Author

@hkershaw-brown I incorporated all of your comments and suggestions from the obs converter.

@hkershaw-brown
Copy link
Member

@mgharamti do you have a wrfhydro case (or other model) you can run with this branch?
It would be good to get a sanity check on the netcdf changes.

@mgharamti
Copy link
Contributor Author

@hkershaw-brown this is to confirm that the netcdf changes are good. I've ran filter with different wrf_hydro cases and all look good.

@hkershaw-brown
Copy link
Member

nice thanks @mgharamti

@hkershaw-brown
Copy link
Member

@mgharamti this is good to go, but let's release it next week (paranoia for netcdf changes on a Friday)

@mgharamti
Copy link
Contributor Author

sounds good to me. thanks for your help reviewing this.

@hkershaw-brown hkershaw-brown merged commit 71d70e1 into NCAR:main Sep 10, 2024
4 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants