MARBL-DART #701

robin-armstrong · 2024-07-11T21:23:06Z

Description:

Contains model-mods for data assimilation with the Marine Biogeochemistry Library (MARBL/MOM6).

also fixes #614

Testing Datasets:

A test package for MARBL_joint_estimation can be found at /glade/u/home/rarmstrong/DART_test_package.

hkershaw-brown · 2024-07-15T14:43:01Z

…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.

…r parameter smoothing.

… 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.

… converter.

… 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.

… BATS_clim data converter.

…t, so that it works with the self-contained test package.

it is actually a netcdf file

add executables to root/.gitignore

hkershaw-brown

first batch of review comments

hkershaw-brown · 2024-07-15T17:20:29Z

assimilation_code/modules/observations/default_quantities_mod.f90

 !   QTY_MICROWAVE_BRIGHT_TEMP
+!   QTY_MICROZOOPLANKTON_CARBON


default_quantities_mod is missing many of the new ocean qtys (ocean_quantities_mod.f90)

hkershaw-brown · 2024-07-15T17:24:04Z

assimilation_code/modules/observations/ocean_quantities_mod.f90

+!   QTY_DISSOLVED_INORGANIC_SIO3
+!   QTY_MESOZOOPLANKTON_CARBON
+!   QTY_MICROZOOPLANKTON_CARBON
+!   QTY_PARAM_AUTOTROPH1_KCO2


What's the goal of having PARAM in the QTY name? Just to communicate that it is a parameter?

Do you need a qty for each parameter? Could you have QTY_PARAMETER?

In order to differentiate this param from other model parameters, we can call it QTY_BGC_PARAM so that it's only a biogeochemistry parameter. I don't think we need a quantity for every single parameter. If I am not mistaken, all MARBL parameters have the same physical properties, for instance >0.

If I remember correctly, I created multiple QTY_PARAM_... types in order to facilitate experiments where multiple parameters are being estimated at once. For experiments of this type, we need a way to distinguish one parameter from another, although they do indeed have the same physical properties. That said, we did eventually switch to only running with one parameter being estimated at a time, and for experiments of that sort, having a single QTY_BGC_PARAM would work just fine.

hkershaw-brown · 2024-07-15T17:24:59Z

models/MARBL_joint_estimation/model_mod.f90

+real(r8), parameter :: geolon = 360 - 64.0
+real(r8), parameter :: geolat = 31.0


these are hardcoded, should these be namelist options for people to put the column anywhere on Earth?

I don't see any reason to not make them namelist options, but a user will need to know that any mismatch between the location in the obs-seq file and geolat/geolon will cause DART to reject the observation. I think that was our motive for hard-coding those coordinates. We wanted to make sure that they always agreed with the coordinates written in the obs-seq file, which in turn are hard-coded into the obs converter.

hkershaw-brown · 2024-07-15T17:30:44Z

models/MARBL_joint_estimation/model_mod.f90

+! Print module information to log file and stdout.
+call register_module(source)


Suggested change

! Print module information to log file and stdout.

call register_module(source)

hkershaw-brown · 2024-07-15T17:31:52Z

models/MARBL_joint_estimation/model_mod.f90

+type(time_type) :: mom6_time
+integer :: mom_base_date_in_days, mom_days
+
+mom_base_date_in_days = 139157 ! 1982 1 1 0 0


Is this always Jan 1st 1982?

Unless changes have been made to the MARBL date/time format, yes, it's always Jan 1st 1982.

hkershaw-brown · 2024-07-15T17:39:31Z

models/MARBL_joint_estimation/model_mod.f90

+
+character(len=*), parameter :: routine = 'read_ocean_geometry'
+
+! Need nx, ny


Suggested change

! Need nx, ny

I think this comment is from the MOM6 model_mod where you need nx, ny for the size of the basin_depth variable. The column marble is (2,2) always.

hkershaw-brown · 2024-07-15T17:41:29Z

models/MARBL_joint_estimation/work/input.nml

+&preprocess_nml
+   input_obs_def_mod_file  = '/glade/derecho/scratch/rarmstrong/marbl_dart/DART/observations/forward_operators/DEFAULT_obs_def_mod.F90'
+   output_obs_def_mod_file = '/glade/derecho/scratch/rarmstrong/marbl_dart/DART/observations/forward_operators/obs_def_mod.f90'
+   input_obs_qty_mod_file  = '/glade/derecho/scratch/rarmstrong/marbl_dart/DART/assimilation_code/modules/observations/DEFAULT_obs_kind_mod.F90'
+   output_obs_qty_mod_file = '/glade/derecho/scratch/rarmstrong/marbl_dart/DART/assimilation_code/modules/observations/obs_kind_mod.f90'
+   obs_type_files          = '/glade/derecho/scratch/rarmstrong/marbl_dart/DART/observations/forward_operators/obs_def_ocean_mod.f90'
+   quantity_files          = '/glade/derecho/scratch/rarmstrong/marbl_dart/DART/assimilation_code/modules/observations/default_quantities_mod.f90', '/glade/derecho/scratch/rarmstrong/marbl_dart/DART/assimilation_code/modules/observations/ocean_quantities_mod.f90'
+   /


needs changing to relative paths (within DART) rather than pointing at Robin's directories.

hkershaw-brown · 2024-07-15T17:41:54Z

models/MARBL_joint_estimation/work/input.nml

@@ -0,0 +1,241 @@
+&perfect_model_obs_nml


will need qceff namelists in this file.

hkershaw-brown · 2024-07-15T17:42:53Z

observations/obs_converters/BATS/README.md

@@ -0,0 +1,2 @@
+This folder contains files adapted from those in DART/observations/obs\_converters/text, for the purpose of converting data from the [Bermude Atlantic Time-Series Study](https://bats.bios.asu.edu/) (BATS) into a format readable by DART.


move to .rst file for online docs

hkershaw-brown · 2024-07-15T17:43:29Z

observations/obs_converters/BATS_clim/README.md

@@ -0,0 +1,2 @@
+This folder contains files adapted from those in DART/observations/obs\_converters/text, for the purpose of converting data from the [Bermude Atlantic Time-Series Study](https://bats.bios.asu.edu/) (BATS) into a format readable by DART.


move to .rst for online docs.
Description needed for BATS_clim vs BATS

mgharamti

This is my first rundown of the code. I have not run the code but I plan to do that soon.

mgharamti · 2024-07-15T20:47:52Z

assimilation_code/modules/io/direct_netcdf_mod.f90

+   ret = nf90_inquire(ncfile_in, unlimitedDimID=unlim_dimID)
+   call nc_check(ret, 'read_variables: nf90_inquire', 'unlimitedDimID')
+
+   if (has_unlimited_dim(domain) .and. any(get_io_dim_ids(domain, i) ==  unlim_dimID )) then


Is this a marbl-only issue? Do we need this chunk of code in other models? If so, I'm wondering if we should have a separate issue for this. Perhaps this is a question to @hkershaw-brown.

yup it is an issue already #614
It affects models which have some state variables with an unlimited dimension and some without.

mgharamti · 2024-07-15T20:49:21Z

assimilation_code/modules/observations/ocean_quantities_mod.f90

+!   QTY_DISSOLVED_INORGANIC_SIO3
+!   QTY_MESOZOOPLANKTON_CARBON
+!   QTY_MICROZOOPLANKTON_CARBON
+!   QTY_PARAM_AUTOTROPH1_KCO2


In order to differentiate this param from other model parameters, we can call it QTY_BGC_PARAM so that it's only a biogeochemistry parameter. I don't think we need a quantity for every single parameter. If I am not mistaken, all MARBL parameters have the same physical properties, for instance >0.

mgharamti · 2024-07-15T20:54:15Z

models/MARBL_joint_estimation/model_mod.f90

+                              update_list = update_var_list(1:nfields))
+
+! setting up the DART parameter vector
+if(estimate_params) then


This assumes that the state is always being estimated. This is where we can generalize this and estimate:

only state

only parameters

or both

Perhaps we can visit this later, but just a note that we can come back to.

It's been a while since I was developing the parameter-only version, but as I remember, that setup still needs a separate domain for state variables. The reasoning here was that the forward operator in the parameter-only case involves running MARBL to produce climatological averages, and this operation cannot be coded into DART. Maintaining a separate domain for state variables allows us to compute those averages externally, so that DART can read them as identity observations from a "climatology" component of the state vector. That said, it might be better to rename that component within the model-mod, so that instead of state_dom_id we have something like climatology_dom_id.

mgharamti · 2024-07-15T21:03:31Z

models/MARBL_joint_estimation/model_mod.f90

+!------------------------------------------------------------------
+! Returns the number of items in the state vector as an integer. 
+
+function get_model_size()


This is not very clean, especially if the parameter domain is not active. I would count the domains and loop over them to get the model_size (which is by the way defined on line 84 but not used). I suggest a code similar to the one below:

domain_count = get_num_domains() model_size = 0 do i_dom = 1,domain_count model_size = model_size + get_domain_size(i_dom) enddo

This is related to my previous comment on active domains and the nature of the DA configuration

mgharamti · 2024-07-15T21:11:42Z

models/MARBL_joint_estimation/model_mod.f90

+    ! reading the clamp values
+
+    if (table(i, 3) /= 'NA') then
+        read(table(i,3), '(d16.8)') clamp_vals(i,1)


Is the format specification needed here? What if people added the clamps in a different format?

mgharamti · 2024-07-15T21:15:20Z

models/MARBL_param_estimation/model_mod.f90

@@ -0,0 +1,591 @@
+! DART software - Copyright UCAR. This open source software is provided


I thought we agreed that we will only do the joint state-parameter estimation one. I would stick to one marbl directory and one model mod for now. These could be consolidated into one in the future. I mean something that comes to mind now is to do no_copy_back on the state variables to turn it into a parameter estimation only.

Yes, we did agree to only do the joint state-parameter estimation model-mod. I had all three model-mods being developed on the same branch to make it easier to switch between them, which is why they're all in the same pull request; I apologize if that's bad form. In any case, MARBL_joint_estimation is the one that we agreed to pull into DART.

was pointing at Robin's directory

hkershaw-brown · 2024-07-19T13:37:59Z

note for test package, I switched to inflation not read from file

hkershaw@dec1004:/glade/work/hkershaw/DART/marbl/DART_test_package/test/DART$ diff input.nml input.nml.robin 
1,6d0
< &probit_transform_nml
< /
< 
< &algorithm_info_nml
< /
< 
73,74c67,68
<    inf_initial_from_restart    = .false.,                 .false.,
<    inf_sd_initial_from_restart = .false.,                  .false.,
---
>    inf_initial_from_restart    = .true.,                 .false.,
>    inf_sd_initial_from_restart = .true.,                  .false.,
108a103
>    filter_kind                     = 1,

mgharamti · 2024-07-23T17:09:08Z

I checked out the code and I was not able to run any of the DART executables. Here is the error that I was getting:

ERROR FROM:
  source : ensemble_manager_mod.f90
  routine: set_up_ens_distribution
  message: not enough MPI tasks for the model size, suggest at least ***** tasks

So, it turns out the length of the state vector is 7036994676854122!!
I realize the issue is in computing the variable model_size. Right now, the function get_model_size in the model_mod is computing it using local variables. When someone calls get_model_size() outside the model_mod, this fails. So, the correct approach is to compute model_size inside static_init_model as I outlined in my previous review and then get_model_size should simply look like this:

function get_model_size()
integer(i8) :: get_model_size
get_model_size = model_size 
end function get_model_size

With this change, I can now run model_mod_check and the model_size is now correctly computed as 3434.
It's still puzzling to me that this version of the code was running during Robin's internship. Could it be that we have a different version?

MARBL_joint_estimation is the model_mod we are using (see pull NCAR#701)

mgharamti

I was able to run filter and state and parameter updates. The only necessary change that is needed to run filter properly is the handling of model_size.

mgharamti · 2024-07-23T17:21:45Z

models/MARBL_joint_estimation/work/input.nml

+&filter_nml
+   single_file_in               = .false.,
+   input_state_files = '',
+   input_state_file_list = 'states_input.txt', 'params_input.txt',


If I set estimate_params = .false. in the model_nml, the code will still look for a list of parameter files. I guess that's ok given that this is a joint state-parameter estimation routine. This can be easily fixed once we work on merging the 3 different filtering strategies together.

If I set estimate_params = .false. in the model_nml, the code will still look for a list of parameter files. I guess that's ok given that this is a joint state-parameter estimation routine. This can be easily fixed once we work on merging the 3 different filtering strategies together.

merging the 3 different filtering strategies is happening as part of this pull request, right?

Helen, I'm not quite sure. I'm happy to do it but it would take some time as I'm leaving. I say, we push the joint one now after small corrections. Once I'm back, I can get to it and add flexibility to perform the 3 using one code.

There is quite a bit to do before releasing on this code. I'll go ahead an close the pull request. We can reopen it when it its ready.

hkershaw-brown · 2024-07-23T20:16:10Z

pushed robins branch to https://github.com/NCAR/DART/tree/marbl-dart

hkershaw-brown · 2024-08-30T14:00:36Z

pushed robins branch to https://github.com/NCAR/DART/tree/marbl-dart
deleting this branch, since Moha's pull request from his fork:

https://github.com/mgharamti/DART/tree/marbl-dart

hkershaw-brown self-assigned this Jul 15, 2024

robin-armstrong added 14 commits July 15, 2024 10:56

Copied in the current versions of MARBL_state_estimation (developed i…

75aa8f5

…n models/MARBL_MOM6_1D in branch marbl_mom6_1d) and MARBL_joint_estimation (developed in models/MARBL_MOM6_1D on branch param_estimation).

Added in the BATS data converter which was developed on the param_est…

77926dc

…imation branch. Also added in the observation codes for BATS variables, which were accidentally left out of the last commit.

Modified bats_to_obs.f90 so that it can generate obs_seq files fit fo…

59917c7

…r parameter smoothing.

Added QTY codes for a subset of MARBL biogeochemical parameters.

a8f0bdf

Added a draft (untested) model interface for MARBL_param_estimation.

372fb89

model_mod_check now runs without crashing for MARBL_param_estimation.…

eb03a41

… This required fixing a bug in direct_netcdf_mod, using a solution provided by Helen Kershaw.

Interpolation function for MARBL_param_estimation now runs correctly.

a5f52b1

Changed MARBL_param_estimation so that only one parameter is estimate…

6e3b47f

…d. Split the BATS data converter into two separate converters, one for sequential state or state-parameter estimation, and another for pure parameter estimation.

Cleaned up code for MARBL_param_estimation and the corresponding data…

178bf70

… converter.

Revisions to MARBL_param_estimation, which now runs data assimilation…

64bee6f

… cycles with MARBL without errors. Further testing needed to assess the quality of assimilations.

MARBL_param_estimation now estimates only a single parameter, bats_cl…

2bfc151

…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.

Made observation error variance inflation a namelist parameter in the…

6e826ea

… BATS_clim data converter.

Added QTY codes for additional MARBL parameters.

9e2d457

Changes to the joint state-parameter estimation model-mod and namelis…

b126a16

…t, so that it works with the self-contained test package.

hkershaw-brown force-pushed the marbl-dart branch from 9ed2fed to b126a16 Compare July 15, 2024 15:00

hkershaw-brown added 5 commits July 15, 2024 12:41

revert mpi utilities to main

b54a467

remove netcdf files from MARBL directories

5917450

remove binary file null

15bb370

it is actually a netcdf file

chore: remove .gitignore from various converters

c7dc847

add executables to root/.gitignore

Merge remote-tracking branch 'ncar/main' into marbl-dart

6a9ef4f

hkershaw-brown requested changes Jul 15, 2024

View reviewed changes

mgharamti reviewed Jul 15, 2024

View reviewed changes

fix: relative paths for preprocess input/output

864a61e

was pointing at Robin's directory

remove _param_estimation and _state_estimation

81c7668

MARBL_joint_estimation is the model_mod we are using (see pull NCAR#701)

mgharamti requested changes Jul 23, 2024

View reviewed changes

hkershaw-brown closed this Jul 23, 2024

mgharamti mentioned this pull request Aug 29, 2024

1D Ocean Column MARBL Interface + Obs Converter #726

Merged

15 tasks

		real(r8), parameter :: geolon = 360 - 64.0
		real(r8), parameter :: geolat = 31.0

		! Print module information to log file and stdout.
		call register_module(source)


		character(len=*), parameter :: routine = 'read_ocean_geometry'

		! Need nx, ny

		@@ -0,0 +1,2 @@
		This folder contains files adapted from those in DART/observations/obs\_converters/text, for the purpose of converting data from the [Bermude Atlantic Time-Series Study](https://bats.bios.asu.edu/) (BATS) into a format readable by DART.

		@@ -0,0 +1,591 @@
		! DART software - Copyright UCAR. This open source software is provided

MARBL-DART #701

MARBL-DART #701

Conversation

robin-armstrong commented Jul 11, 2024 • edited by hkershaw-brown Loading

Description:

Testing Datasets:

hkershaw-brown commented Jul 15, 2024 • edited Loading

hkershaw-brown left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mgharamti Jul 15, 2024 • edited Loading

Choose a reason for hiding this comment

robin-armstrong Jul 16, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mgharamti left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mgharamti Jul 15, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hkershaw-brown commented Jul 19, 2024 • edited Loading

mgharamti commented Jul 23, 2024

mgharamti left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hkershaw-brown commented Jul 23, 2024

hkershaw-brown commented Aug 30, 2024

robin-armstrong commented Jul 11, 2024 •

edited by hkershaw-brown

Loading

hkershaw-brown commented Jul 15, 2024 •

edited

Loading

mgharamti Jul 15, 2024 •

edited

Loading

robin-armstrong Jul 16, 2024 •

edited

Loading

mgharamti Jul 15, 2024 •

edited

Loading

hkershaw-brown commented Jul 19, 2024 •

edited

Loading