From 48d6d599a188a43be5aa6a021c407197668522d6 Mon Sep 17 00:00:00 2001 From: Mohamad Gharamti Date: Tue, 3 Sep 2024 13:21:36 -0600 Subject: [PATCH] model-mod and documentation clean up Incorporated all suggestions and comments about model-mod. Several comments have been removed, redundant loc routines are omitted, better namelist parsing, deleted unused variables, ... --- models/MARBL_column/model_mod.f90 | 232 +++++++++-------------------- models/MARBL_column/readme.rst | 19 ++- models/MARBL_column/work/input.nml | 7 +- models/README.rst | 2 +- 4 files changed, 90 insertions(+), 170 deletions(-) diff --git a/models/MARBL_column/model_mod.f90 b/models/MARBL_column/model_mod.f90 index 86b51ddd2..7c929a2a1 100644 --- a/models/MARBL_column/model_mod.f90 +++ b/models/MARBL_column/model_mod.f90 @@ -4,52 +4,47 @@ module model_mod -! 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. -use types_mod, only : r8, i8, MISSING_R8, vtablenamelength +use types_mod, only : r8, i8, MISSING_R8, vtablenamelength -use time_manager_mod, only : time_type, set_time +use time_manager_mod, only : time_type, set_time -use location_mod, only : location_type, get_close_type, & - loc_get_close_obs => get_close_obs, & - loc_get_close_state => get_close_state, & - set_location, set_location_missing, & - get_location, VERTISLEVEL +use location_mod, only : location_type, get_close_type, & + get_close_obs, get_close_state, & + set_location, set_location_missing, & + get_location, VERTISLEVEL -use utilities_mod, only : register_module, error_handler, & - E_ERR, E_MSG, & - nmlfileunit, do_output, do_nml_file, do_nml_term, & - find_namelist_in_file, check_namelist_read, & - to_upper +use utilities_mod, only : error_handler, E_ERR, E_MSG, & + nmlfileunit, do_output, do_nml_file, do_nml_term, & + find_namelist_in_file, check_namelist_read, & + to_upper -use netcdf_utilities_mod, only : nc_add_global_attribute, nc_synchronize_file, & - nc_add_global_creation_time, & - nc_begin_define_mode, nc_end_define_mode, & - NF90_MAX_NAME, nc_open_file_readonly, & - nc_get_variable, nc_get_variable_size, nc_close_file +use netcdf_utilities_mod, only : nc_add_global_attribute, nc_synchronize_file, & + nc_add_global_creation_time, & + nc_begin_define_mode, nc_end_define_mode, & + NF90_MAX_NAME, nc_open_file_readonly, & + nc_get_variable, nc_get_variable_size, nc_close_file -use state_structure_mod, only : add_domain, get_domain_size, get_model_variable_indices, & - get_varid_from_kind, get_dart_vector_index, & - get_num_domains +use state_structure_mod, only : add_domain, get_domain_size, get_model_variable_indices, & + get_varid_from_kind, get_dart_vector_index, & + get_num_domains -use obs_kind_mod, only : get_index_for_quantity, QTY_U_CURRENT_COMPONENT, & - QTY_V_CURRENT_COMPONENT, QTY_DRY_LAND, & - QTY_LAYER_THICKNESS +use obs_kind_mod, only : get_index_for_quantity, QTY_LAYER_THICKNESS -use distributed_state_mod, only: get_state +use distributed_state_mod, only : get_state -use ensemble_manager_mod, only : ensemble_type +use ensemble_manager_mod, only : ensemble_type ! These routines are passed through from default_model_mod. ! To write model specific versions of these routines ! remove the routine from this use statement and add your code to ! this the file. -use default_model_mod, only : pert_model_copies, write_model_time, & - init_time => fail_init_time, & - init_conditions => fail_init_conditions, & - convert_vertical_obs, convert_vertical_state, adv_1step +use default_model_mod, only : pert_model_copies, write_model_time, & + init_time => fail_init_time, & + init_conditions => fail_init_conditions, & + convert_vertical_obs, convert_vertical_state, adv_1step implicit none private @@ -62,9 +57,9 @@ module model_mod end_model, & static_init_model, & nc_write_model_atts, & + pert_model_copies, & get_close_obs, & get_close_state, & - pert_model_copies, & convert_vertical_obs, & convert_vertical_state, & read_model_time, & @@ -74,38 +69,40 @@ module model_mod shortest_time_between_assimilations, & write_model_time - character(len=256), parameter :: source = "model_mod.f90" -logical :: module_initialized = .false. -integer :: state_dom_id ! used to access the state structure -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 + +logical :: module_initialized = .false. +integer :: state_dom_id ! used to access the state structure +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(i8) :: model_size type(time_type) :: assimilation_time_step -real(r8) :: geolon -real(r8) :: geolat -real(r8) :: basin_depth(2,2) +real(r8) :: geolon +real(r8) :: geolat +real(r8) :: basin_depth(2,2) ! parameters to be used in specifying the DART internal state -integer, parameter :: modelvar_table_height = 13 -integer, parameter :: modelvar_table_width = 5 -integer, parameter :: modelparams_table_height = 1 -integer, parameter :: modelparams_table_width = 5 +integer, parameter :: MODELVAR_TABLE_HEIGHT = 13 +integer, parameter :: MODELVAR_TABLE_WIDTH = 5 +integer, parameter :: MODELPARAMS_TABLE_HEIGHT = 1 +integer, parameter :: MODELPARAMS_TABLE_WIDTH = 5 ! defining the variables that will be read from the namelist -character(len=256) :: template_file(2) -character(len=256) :: ocean_geometry -real(r8) :: station_location(2) -integer :: time_step_days -integer :: time_step_seconds -logical :: estimate_params -character(len=vtablenamelength) & - :: model_state_variables(modelvar_table_height * modelvar_table_width) -character(len=vtablenamelength) & - :: model_parameters(modelparams_table_height * modelparams_table_width) - -namelist /model_nml/ template_file, & +character(len=256) :: state_template_file = 'MOM.res.nc' +character(len=256) :: param_template_file = 'marbl_params.nc' +character(len=256) :: ocean_geometry = 'ocean_geometry.nc' +real(r8) :: station_location(2) = (/-64.0, 31.0/) +integer :: time_step_days = 1 +integer :: time_step_seconds = 0 +logical :: estimate_params = .false. +character(len=vtablenamelength) :: & + model_state_variables(MODELVAR_TABLE_HEIGHT * MODELVAR_TABLE_WIDTH) = '' +character(len=vtablenamelength) :: & + model_parameters(MODELPARAMS_TABLE_HEIGHT * MODELPARAMS_TABLE_WIDTH) = '' + +namelist /model_nml/ state_template_file, & + param_template_file, & ocean_geometry, & station_location, & time_step_days, & @@ -117,22 +114,19 @@ module model_mod !------------------------------------------------------------------ ! -! 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. subroutine static_init_model() integer :: iunit, io, i_dom, domain_count -character(len=vtablenamelength) :: variable_table(modelvar_table_height, modelvar_table_width), & - param_table(modelparams_table_height, modelparams_table_width) -integer :: state_qty_list(modelvar_table_height), & - param_qty_list(modelparams_table_height) -real(r8):: state_clamp_vals(modelvar_table_height, 2), & - param_clamp_vals(modelparams_table_height, 2) -logical :: update_var_list(modelvar_table_height), & - update_param_list(modelparams_table_height) +character(len=vtablenamelength) :: variable_table(MODELVAR_TABLE_HEIGHT, MODELVAR_TABLE_WIDTH), & + param_table(MODELPARAMS_TABLE_HEIGHT, MODELPARAMS_TABLE_WIDTH) +integer :: state_qty_list(MODELVAR_TABLE_HEIGHT), & + param_qty_list(MODELPARAMS_TABLE_HEIGHT) +real(r8):: state_clamp_vals(MODELVAR_TABLE_HEIGHT, 2), & + param_clamp_vals(MODELPARAMS_TABLE_HEIGHT, 2) +logical :: update_var_list(MODELVAR_TABLE_HEIGHT), & + update_param_list(MODELPARAMS_TABLE_HEIGHT) module_initialized = .true. @@ -155,6 +149,12 @@ subroutine static_init_model() assimilation_time_step = set_time(time_step_seconds, & time_step_days) +! Make sure the longitude is properly assigned +geolon = station_location(1) +geolat = station_location(2) + +if (geolon < 0.0_r8) geolon = 360.0 + geolon + ! The kind of experiments that can be performed using this code: ! 1- State estimation: ! To achieve this, you'll need to set "estimate_params" @@ -173,7 +173,7 @@ subroutine static_init_model() call verify_state_variables(model_state_variables, nfields, variable_table, & state_qty_list, state_clamp_vals, update_var_list) - state_dom_id = add_domain(template_file(1), nfields, & + state_dom_id = add_domain(state_template_file, nfields, & var_names = variable_table(1:nfields, 1), & kind_list = state_qty_list(1:nfields), & clamp_vals = state_clamp_vals, & @@ -184,7 +184,7 @@ subroutine static_init_model() call verify_state_variables(model_parameters, nfields, param_table, & param_qty_list, param_clamp_vals, update_param_list) - param_dom_id = add_domain(template_file(2), nfields, & + param_dom_id = add_domain(param_template_file, nfields, & var_names = param_table(1:nfields, 1), & kind_list = param_qty_list(1:nfields), & clamp_vals = param_clamp_vals, & @@ -215,7 +215,6 @@ function read_model_time(filename) integer :: ncid character(len=*), parameter :: routine = 'read_model_time' real(r8) :: days -type(time_type) :: mom6_time integer :: mom_base_date_in_days, mom_days mom_base_date_in_days = 139157 ! 1982 1 1 0 0 @@ -289,10 +288,6 @@ subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs ! performing the interpolation for each ensemble member individually. do ens_index = 1, ens_size - !print *, "BASIN DEPTH = ",-basin_depth(1,1) - !print *, "----------------------------------------------------" - !print *, "computing centers of layers" - !print *, "----------------------------------------------------" ! locating the layer index to be used as the upper interpolation point for this ensemble member. layer_index = nz @@ -302,26 +297,17 @@ subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs layerdepth_top = layerdepth_bottom & ! depth at top of the layer given by layer_index. + layer_thicknesses(layer_index, ens_index) - !print *, "layer = ",layer_index,", center = ",layerdepth_center,", bottom = ",layerdepth_bottom,", top = ",layerdepth_top - do while((layerdepth_center < requested_depth) .and. (layer_index > 1)) layer_index = layer_index - 1 layerdepth_bottom = layerdepth_bottom + layer_thicknesses(layer_index, ens_index) layerdepth_center = layerdepth_bottom + 0.5 * layer_thicknesses(layer_index, ens_index) layerdepth_top = layerdepth_bottom + layer_thicknesses(layer_index, ens_index) - - !print *, "layer = ",layer_index,", bottom = ",layerdepth_bottom,", center = ",layerdepth_center,", top = ",layerdepth_top - end do ! having located the index of the upper interpolation layer, we now calculate the interpolation. if((requested_depth < -basin_depth(1,1)) .or. (layerdepth_top < requested_depth)) then ! case where the requested depth is below the ocean floor, or above the ocean surface - !print *, "----------------------------------------------------" - !print *, "depth is below ocean floor or above ocean surface" - !print *, "----------------------------------------------------" - istatus(ens_index) = 1 expected_obs(ens_index) = MISSING_R8 @@ -330,25 +316,16 @@ subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs ! or the top half of the shallowest layer. In both cases, the "interpolated" value is ! simply the current value of that layer in MOM6. - !print *, "----------------------------------------------------" - !print *, "only using value from layer ", layer_index - !print *, "----------------------------------------------------" - istatus(ens_index) = 0 qty_index = get_dart_vector_index(1, 1, layer_index, state_dom_id, qty_id) state_slice = get_state(qty_index, state_handle) expected_obs(ens_index) = state_slice(ens_index) - ! print *, "final value = ",expected_obs(ens_index) - else ! case where the requested depth is above the center of some layer, and below ! the center of another. We interpolate linearly between the nearest layers above ! and below. - !print *, "----------------------------------------------------" - !print *, "interpolating between layers ", layer_index, " and ", (layer_index + 1) - istatus(ens_index) = 0 ! computing the depths at the centers of the nearest layers above and below @@ -365,13 +342,9 @@ subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs state_slice = get_state(qty_index, state_handle) val_below = state_slice(ens_index) - !print *, "corresponding to the values: ", val_above, " and ", val_below - !print *, "----------------------------------------------------" - ! linear interpolation expected_obs(ens_index) = val_above + (requested_depth - depth_above) * (val_below - val_above) / (depth_below - depth_above) - !print *, "final value = ",expected_obs(ens_index) end if end do @@ -411,12 +384,6 @@ subroutine get_state_meta_data(index_in, location, qty) call get_model_variable_indices(index_in, lon_index, lat_index, level, kind_index=local_qty) -! Make sure the longitude is properly assigned -geolon = station_location(1) -geolat = station_location(2) - -if (geolon < 0.0_r8) geolon = 360.0 + geolon - location = set_location(geolon, geolat, real(level,r8), VERTISLEVEL) if (present(qty)) then @@ -426,55 +393,6 @@ subroutine get_state_meta_data(index_in, location, qty) end subroutine get_state_meta_data -!------------------------------------------------------------------ -! Any model specific distance calcualtion can be done here -subroutine get_close_obs(gc, base_loc, base_type, locs, loc_qtys, loc_types, & - num_close, close_ind, dist, ens_handle) - -type(get_close_type), intent(in) :: gc ! handle to a get_close structure -integer, intent(in) :: base_type ! observation TYPE -type(location_type), intent(inout) :: base_loc ! location of interest -type(location_type), intent(inout) :: locs(:) ! obs locations -integer, intent(in) :: loc_qtys(:) ! QTYS for obs -integer, intent(in) :: loc_types(:) ! TYPES for obs -integer, intent(out) :: num_close ! how many are close -integer, intent(out) :: close_ind(:) ! incidies into the locs array -real(r8), optional, intent(out) :: dist(:) ! distances in radians -type(ensemble_type), optional, intent(in) :: ens_handle - -character(len=*), parameter :: routine = 'get_close_obs' - -call loc_get_close_obs(gc, base_loc, base_type, locs, loc_qtys, loc_types, & - num_close, close_ind, dist, ens_handle) - -end subroutine get_close_obs - - -!------------------------------------------------------------------ -! Any model specific distance calcualtion can be done here -subroutine get_close_state(gc, base_loc, base_type, locs, loc_qtys, loc_indx, & - num_close, close_ind, dist, ens_handle) - -type(get_close_type), intent(in) :: gc ! handle to a get_close structure -type(location_type), intent(inout) :: base_loc ! location of interest -integer, intent(in) :: base_type ! observation TYPE -type(location_type), intent(inout) :: locs(:) ! state locations -integer, intent(in) :: loc_qtys(:) ! QTYs for state -integer(i8), intent(in) :: loc_indx(:) ! indices into DART state vector -integer, intent(out) :: num_close ! how many are close -integer, intent(out) :: close_ind(:) ! indices into the locs array -real(r8), optional, intent(out) :: dist(:) ! distances in radians -type(ensemble_type), optional, intent(in) :: ens_handle - -character(len=*), parameter :: routine = 'get_close_state' - - -call loc_get_close_state(gc, base_loc, base_type, locs, loc_qtys, loc_indx, & - num_close, close_ind, dist, ens_handle) - - -end subroutine get_close_state - !------------------------------------------------------------------ ! write any additional attributes to the output and diagnostic files @@ -597,9 +515,7 @@ subroutine read_num_layers() integer :: ncid -character(len=*), parameter :: routine = 'read_num_layers' - -ncid = nc_open_file_readonly(template_file(1)) +ncid = nc_open_file_readonly(state_template_file) call nc_get_variable_size(ncid, 'Layer', nz) diff --git a/models/MARBL_column/readme.rst b/models/MARBL_column/readme.rst index 413089627..abdf04e99 100644 --- a/models/MARBL_column/readme.rst +++ b/models/MARBL_column/readme.rst @@ -45,8 +45,9 @@ The code is designed to perform 3 kinds of data assimilation (DA) experiments: 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`` + This ensures that the updated state will not be written back to the restart file. 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``. Namelist -------- @@ -55,8 +56,9 @@ The ``&model_nml`` variables and their default values are listed here: .. code-block:: fortran &model_nml - template_file = 'member_0001/RESTART/MOM.res.nc', 'member_0001/marbl_params.nc', - ocean_geometry = 'member_0001/ocean_geometry.nc', + state_template_file = 'MOM.res.nc', + param_template_file = 'marbl_params.nc', + ocean_geometry = 'ocean_geometry.nc', station_location = -64, 31 time_step_days = 1, time_step_seconds = 0, @@ -82,13 +84,14 @@ This namelist provides control over the kind of DA experiment as described abvov +-------------------------------------+--------------------+------------------------------------------------------------+ | Item | Type | Description | +=====================================+====================+============================================================+ -| ``template_file`` | character(len=256) | MARBL restart file including MARBL's prognostic variables | +| ``state_template_file`` | character(len=256) | MARBL restart file including MARBL's prognostic variables | | | | and other grid information such as the ocean layers. | | | | The state variables read from this file are listed in | | | | in the ``model_state_variables`` | -| | | | -| | | There is another template file corresponding to the | -| | | BGC parameters that we intend to estimate. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``param_template_file`` | character(len=256) | Template file corresponding to the BGC parameters that we | +| | | intend to estimate. The parameters read from this file are | +| | | listed in the ``model_parameters``. | +-------------------------------------+--------------------+------------------------------------------------------------+ | ``ocean_geometry`` | character(len=256) | The ocean geometry file is used to read ``basin_depth`` | +-------------------------------------+--------------------+------------------------------------------------------------+ diff --git a/models/MARBL_column/work/input.nml b/models/MARBL_column/work/input.nml index d76881e76..c2c1035f4 100644 --- a/models/MARBL_column/work/input.nml +++ b/models/MARBL_column/work/input.nml @@ -7,7 +7,8 @@ / &model_nml - template_file = 'ensemble/background/member_0001/RESTART/MOM.res.nc', 'ensemble/background/member_0001/marbl_params.nc', + state_template_file = 'ensemble/background/member_0001/RESTART/MOM.res.nc', + param_template_file = 'ensemble/background/member_0001/marbl_params.nc', ocean_geometry = 'ensemble/background/member_0001/ocean_geometry.nc', station_location = -64, 31 time_step_days = 1, @@ -236,8 +237,8 @@ / &model_mod_check_nml - input_state_files = 'mom_in.nc', 'param_in.nc' - output_state_files = 'mom_out.nc', 'param_out.nc' + input_state_files = 'marbl_in.nc', 'param_in.nc' + output_state_files = 'marbl_out.nc', 'param_out.nc' test1thru = 7 run_tests = 1 x_ind = 600 diff --git a/models/README.rst b/models/README.rst index 3e9f0d1a4..932d60424 100644 --- a/models/README.rst +++ b/models/README.rst @@ -27,11 +27,11 @@ DART supported models: - :doc:`lorenz_96_2scale/readme` - :doc:`lorenz_96_tracer_advection/readme` - :doc:`forced_lorenz_96/readme` +- :doc:`MARBL_column/readme` - :doc:`MITgcm_ocean/readme` - :doc:`MOM6/readme` - :doc:`mpas_atm/readme` - :doc:`mpas_ocn/readme` -- :doc:`MARBL_column/readme` - :doc:`NCOMMAS/readme` - :doc:`noah/readme` - :doc:`null_model/readme`