diff --git a/.gitignore b/.gitignore index b8f113f92f..7874132504 100644 --- a/.gitignore +++ b/.gitignore @@ -162,6 +162,8 @@ ssec_satwnd gts_to_dart littler_tf_dart rad_3dvar_to_dart +bats_to_clim_obs +bats_to_obs L1_AMSUA_to_netcdf convert_airs_L2 convert_amsu_L1 diff --git a/CHANGELOG.rst b/CHANGELOG.rst index ea49e40a59..d31552f944 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,6 +22,22 @@ individual files. The changes are now listed with the most recent at the top. +**September 10 2024:: MARBL_column. Tag 11.8.0** + +- Interface for MARBL_column for DART: + + - state estimation + - state and parameter estimation + - parameter estimation only + +- BATS observation converter, and BATS climatology scripting + +*contributed by Robin Armstrong* + +Bugfix: + +- fix for IO for NetCDF files when only some variables have the unlimited dimension + **August 29 2024 :: Bug fixes for shortest_time_between_assimilations and get_close_init. Tag 11.7.1** Bug fixes: diff --git a/assimilation_code/modules/io/direct_netcdf_mod.f90 b/assimilation_code/modules/io/direct_netcdf_mod.f90 index 9419e19c73..aea2d8dae3 100644 --- a/assimilation_code/modules/io/direct_netcdf_mod.f90 +++ b/assimilation_code/modules/io/direct_netcdf_mod.f90 @@ -93,7 +93,7 @@ module direct_netcdf_mod get_index_start, get_index_end , get_num_dims, & create_diagnostic_structure, & end_diagnostic_structure, & - has_unlimited_dim + has_unlimited_dim, get_io_dim_ids use io_filenames_mod, only : get_restart_filename, inherit_copy_units, & stage_metadata_type, get_file_description, & @@ -858,16 +858,17 @@ subroutine read_variables(ncfile_in, var_block, start_var, end_var, domain) counts(:) = 1 slice_start(:) = 1 ! default to read all dimensions start at 1 - - if (has_unlimited_dim(domain)) then + + 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 counts(num_dims) = 1 ! one slice of unlimited dimesion - counts(1:num_dims-1) = get_dim_lengths(domain, i) ! the state + counts(1:get_num_dims(domain, i)) = get_dim_lengths(domain, i) ! the state ! read latest time slice - hack to get started with tiegcm ! not sure if it will always be the last time slice - ret = nf90_inquire(ncfile_in, unlimitedDimID=unlim_dimID) - call nc_check(ret, 'read_variables: nf90_inquire', 'unlimitedDimID') if (unlim_dimID /= -1) then ! unlimited dimension exists ret = nf90_inquire_dimension(ncfile_in, unlim_dimID, len=slice_start(num_dims)) @@ -1588,15 +1589,17 @@ subroutine write_variables(ncid, var_block, start_var, end_var, domain, & slice_start(:) = 1 ! default to read all dimensions starting at 1 counts(:) = 1 - if (has_unlimited_dim(domain)) then + ret = nf90_inquire(ncid, 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 counts(num_dims) = 1 ! one slice of unlimited dimesion counts(1:get_num_dims(domain, i)) = get_dim_lengths(domain, i) ! write the latest time slice - HK hack to get started with tiegcm ! not sure if it will always be the last time slice - ret = nf90_inquire(ncid, unlimitedDimID=unlim_dimID) - call nc_check(ret, 'write_variables: nf90_inquire', 'unlimitedDimID') + if (unlim_dimID /= -1) then ! unlimited dimension exists ret = nf90_inquire_dimension(ncid, unlim_dimID, len=slice_start(num_dims)) call nc_check(ret, 'write_variables: nf90_inquire dimension', 'unlimitedDim length') diff --git a/assimilation_code/modules/observations/default_quantities_mod.f90 b/assimilation_code/modules/observations/default_quantities_mod.f90 index 84fc0cea8a..ebc2d068e9 100644 --- a/assimilation_code/modules/observations/default_quantities_mod.f90 +++ b/assimilation_code/modules/observations/default_quantities_mod.f90 @@ -72,6 +72,8 @@ ! QTY_CWP_PATH_ZERO ! QTY_DISSOLVED_INORGANIC_CARBON ! QTY_DISSOLVED_INORGANIC_IRON +! QTY_DISSOLVED_INORGANIC_SIO3 +! QTY_DISSOLVED_ORGANIC_CARBON ! QTY_DISSOLVED_ORGANIC_NITROGEN ! QTY_DISSOLVED_ORGANIC_P ! QTY_DISSOLVED_OXYGEN @@ -158,6 +160,7 @@ ! QTY_LANDMASK ! QTY_LARGE_SCALE_STATE desc="state varies with large time/space scale" ! QTY_LATENT_HEAT_FLUX +! QTY_LAYER_THICKNESS ! QTY_LEAF_AREA_INDEX ! QTY_LEAF_CARBON ! QTY_LEAF_NITROGEN @@ -167,7 +170,9 @@ ! QTY_LIVE_STEM_CARBON ! QTY_LIVE_STEM_NITROGEN ! QTY_MEAN_SOURCE +! QTY_MESOZOOPLANKTON_CARBON ! QTY_MICROWAVE_BRIGHT_TEMP +! QTY_MICROZOOPLANKTON_CARBON ! QTY_MINERAL_ACCUM ! QTY_MINERAL_COARSE ! QTY_MINERAL_NUCLEUS @@ -382,6 +387,8 @@ ! QTY_WIND_TURBINE_POWER ! QTY_W_CURRENT_COMPONENT ! QTY_X_LAMBDA +! QTY_COLUMN_DEPTH +! QTY_BGC_PARAM ! ! END DART PREPROCESS QUANTITY DEFINITIONS diff --git a/assimilation_code/modules/observations/ocean_quantities_mod.f90 b/assimilation_code/modules/observations/ocean_quantities_mod.f90 index fb4d0a3bc5..26e04124fb 100644 --- a/assimilation_code/modules/observations/ocean_quantities_mod.f90 +++ b/assimilation_code/modules/observations/ocean_quantities_mod.f90 @@ -32,11 +32,17 @@ ! QTY_PHYTOPLANKTON_BIOMASS ! QTY_ALKALINITY ! QTY_DISSOLVED_INORGANIC_CARBON +! QTY_DISSOLVED_ORGANIC_CARBON ! QTY_DISSOLVED_ORGANIC_NITROGEN ! QTY_DISSOLVED_ORGANIC_P ! QTY_DISSOLVED_INORGANIC_IRON ! QTY_SURFACE_CHLOROPHYLL ! QTY_LAYER_THICKNESS +! QTY_COLUMN_DEPTH +! QTY_DISSOLVED_INORGANIC_SIO3 +! QTY_MESOZOOPLANKTON_CARBON +! QTY_MICROZOOPLANKTON_CARBON +! QTY_BGC_PARAM ! ! ! fixme - these units are hardcoded in obs_diag and shouldn't be ! QTY_U_WIND_COMPONENT diff --git a/conf.py b/conf.py index 7d4d9a8f4f..833c43425d 100644 --- a/conf.py +++ b/conf.py @@ -21,7 +21,7 @@ author = 'Data Assimilation Research Section' # The full version, including alpha/beta/rc tags -release = '11.7.1' +release = '11.8.0' root_doc = 'index' # -- General configuration --------------------------------------------------- diff --git a/guide/available-observation-converters.rst b/guide/available-observation-converters.rst index 6b84da9cae..aa3b1cdc76 100644 --- a/guide/available-observation-converters.rst +++ b/guide/available-observation-converters.rst @@ -11,6 +11,7 @@ Each directory has at least one converter: - ``AURA``: See ``DART/observations/obs_converters/AURA`` - ``Aviso+/CMEMS``: :doc:`../observations/obs_converters/AVISO/AVISO` - ``Ameriflux``: :doc:`../observations/obs_converters/Ameriflux/level4_to_obs` +- ``BATS``: :doc:`../observations/obs_converters/BATS/readme` - ``CHAMP``: :doc:`../observations/obs_converters/CHAMP/work/README` - ``cice``: :doc:`../observations/obs_converters/cice/cice_to_obs` - ``CNOFS``: See ``DART/observations/obs_converters/CNOFS`` diff --git a/index.rst b/index.rst index 8b067e83fe..06b4686699 100644 --- a/index.rst +++ b/index.rst @@ -361,6 +361,7 @@ References observations/obs_converters/Ameriflux/fluxnetfull_to_obs observations/obs_converters/Ameriflux/level4_to_obs observations/obs_converters/CHAMP/work/README + observations/obs_converters/BATS/readme observations/obs_converters/cice/cice_to_obs observations/obs_converters/CONAGUA/README observations/obs_converters/COSMOS/COSMOS_to_obs @@ -466,6 +467,7 @@ References models/lorenz_96_tracer_advection/readme models/forced_lorenz_96/readme models/MITgcm_ocean/readme + models/MARBL_column/readme models/MOM6/readme models/mpas_atm/readme models/mpas_atm/mpas_dart_obs_preprocess diff --git a/models/MARBL_column/model_mod.f90 b/models/MARBL_column/model_mod.f90 new file mode 100644 index 0000000000..7c929a2a17 --- /dev/null +++ b/models/MARBL_column/model_mod.f90 @@ -0,0 +1,554 @@ +! by UCAR, "as is", without charge, subject to all terms of use at +! http://www.image.ucar.edu/DAReS/DART/DART_download +! + +module model_mod + +! 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 time_manager_mod, only : time_type, set_time + +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 : 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 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_LAYER_THICKNESS + +use distributed_state_mod, only : get_state + +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 + +implicit none +private + +! routines required by DART code - will be called from filter and other +! DART executables. +public :: get_model_size, & + get_state_meta_data, & + model_interpolate, & + end_model, & + static_init_model, & + nc_write_model_atts, & + pert_model_copies, & + get_close_obs, & + get_close_state, & + convert_vertical_obs, & + convert_vertical_state, & + read_model_time, & + adv_1step, & + init_time, & + init_conditions, & + 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(i8) :: model_size +type(time_type) :: assimilation_time_step +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 + +! defining the variables that will be read from the namelist +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, & + time_step_seconds, & + model_state_variables, & + estimate_params, & + model_parameters +contains + +!------------------------------------------------------------------ +! +! 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) + +module_initialized = .true. + +! Read values from the namelist + +call find_namelist_in_file("input.nml", "model_nml", iunit) +read(iunit, nml = model_nml, iostat = io) + +call check_namelist_read(iunit, io, "model_nml") + +! Record the namelist values used for the run +if (do_nml_file()) write(nmlfileunit, nml=model_nml) +if (do_nml_term()) write( * , nml=model_nml) + +! This time is both the minimum time you can ask the model to advance +! (for models that can be advanced by filter) and it sets the assimilation +! window. All observations within +/- 1/2 this interval from the current +! model time will be assimilated. If this is not settable at runtime +! feel free to hardcode it and remove from the namelist. +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" +! in the model_nml to ".false." +! 2- State and Parameters estimation: +! To achieve this, set "estimate_params" to ".true." +! 3- Parameters estimation only: +! To achieve this, we still need to read in the state +! in order to compute covariances and innovations. +! First, set "estimate_params" to ".true." and +! Next, set the 5th column in the state table as +! "NO_COPY_BACK". The parameters update status should +! be kept as "Update" + +! setting up the DART state vector +call verify_state_variables(model_state_variables, nfields, variable_table, & + state_qty_list, state_clamp_vals, update_var_list) + + 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, & + update_list = update_var_list(1:nfields)) + +! setting up the DART parameter vector +if(estimate_params) then + call verify_state_variables(model_parameters, nfields, param_table, & + param_qty_list, param_clamp_vals, update_param_list) + + 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, & + update_list = update_param_list(1:nfields)) +end if + +call read_num_layers ! setting the value of nz +call read_ocean_geometry ! determining the basin depth + +! Determine the domain and model size +domain_count = get_num_domains() + +model_size = 0 +do i_dom = 1,domain_count + model_size = model_size + get_domain_size(i_dom) +enddo + +end subroutine static_init_model + +!------------------------------------------------------------------ +! Reads the simulation length from a netCDF file. + +function read_model_time(filename) + +character(len=*), intent(in) :: filename +type(time_type) :: read_model_time + +integer :: ncid +character(len=*), parameter :: routine = 'read_model_time' +real(r8) :: days +integer :: mom_base_date_in_days, mom_days + +mom_base_date_in_days = 139157 ! 1982 1 1 0 0 +ncid = nc_open_file_readonly(filename, routine) + +call nc_get_variable(ncid, 'Time', days, routine) + +call nc_close_file(ncid, routine) + +mom_days = int(days) + mom_base_date_in_days + +read_model_time = set_time(0,mom_days) + +end function read_model_time + +!------------------------------------------------------------------ +! Returns the number of items in the state vector as an integer. + +function get_model_size() + +integer(i8) :: get_model_size + +if ( .not. module_initialized ) call static_init_model + +get_model_size = model_size + +end function get_model_size + +!------------------------------------------------------------------ +! Given a state handle, a location, and a state quantity, +! interpolates the state variable fields to that location and returns +! the values in expected_obs. The istatus variables should be returned as +! 0 unless there is some problem in computing the interpolation in +! which case a positive istatus should be returned. + +subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs, istatus) + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: ens_size +type(location_type), intent(in) :: location +integer, intent(in) :: qty +real(r8), intent(out) :: expected_obs(ens_size) !< array of interpolated values +integer, intent(out) :: istatus(ens_size) + +integer :: qty_id, thickness_id, layer_index +integer(i8) :: qty_index, thickness_index, ens_index +real(8) :: requested_depth, layerdepth_bottom, layerdepth_center, & + layerdepth_top, depth_below, depth_above, val_below, val_above +real(8) :: layer_thicknesses(nz, ens_size) +real(8) :: state_slice(ens_size) +real(8) :: loc_temp(3) + +if ( .not. module_initialized ) call static_init_model + +qty_id = get_varid_from_kind(state_dom_id, qty) + +! extracting the requested depth value from `location`. +loc_temp = get_location(location) +requested_depth = loc_temp(3) + +! extracting the current layer thicknesses for each ensemble member. +thickness_id = get_varid_from_kind(state_dom_id, QTY_LAYER_THICKNESS) + +do layer_index = 1, nz + ! layer thicknesses are always extracted from the grid cell at index (1, 1), since the four grid + ! cells are supposed to represent identical columns. + thickness_index = get_dart_vector_index(1, 1, layer_index, state_dom_id, thickness_id) + layer_thicknesses(layer_index, :) = get_state(thickness_index, state_handle) +end do + + +! performing the interpolation for each ensemble member individually. +do ens_index = 1, ens_size + + ! locating the layer index to be used as the upper interpolation point for this ensemble member. + layer_index = nz + layerdepth_bottom = -basin_depth(1,1) ! depth at bottom of the layer given by layer_index. + layerdepth_center = layerdepth_bottom & ! depth at center of the layer given by layer_index. + + .5*layer_thicknesses(layer_index, ens_index) + layerdepth_top = layerdepth_bottom & ! depth at top of the layer given by layer_index. + + layer_thicknesses(layer_index, ens_index) + + 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) + 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 + + istatus(ens_index) = 1 + expected_obs(ens_index) = MISSING_R8 + + else if((layer_index == nz) .or. (layerdepth_center < requested_depth)) then + ! case where the requested depth is either in the bottom half of the deepest layer, + ! or the top half of the shallowest layer. In both cases, the "interpolated" value is + ! simply the current value of that layer in MOM6. + + 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) + + 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. + + istatus(ens_index) = 0 + + ! computing the depths at the centers of the nearest layers above and below + depth_above = layerdepth_center + depth_below = layerdepth_top - layer_thicknesses(layer_index, ens_index) & + - .5*layer_thicknesses(layer_index + 1, ens_index) + + ! extracting the quantity values at the layers above and below + qty_index = get_dart_vector_index(1, 1, layer_index, state_dom_id, qty_id) + state_slice = get_state(qty_index, state_handle) + val_above = state_slice(ens_index) + + qty_index = get_dart_vector_index(1, 1, layer_index + 1, state_dom_id, qty_id) + state_slice = get_state(qty_index, state_handle) + val_below = state_slice(ens_index) + + ! linear interpolation + expected_obs(ens_index) = val_above + (requested_depth - depth_above) * (val_below - val_above) / (depth_below - depth_above) + + end if +end do + +end subroutine model_interpolate + + + +!------------------------------------------------------------------ +! Returns the smallest increment in time that the model is capable +! of advancing the state in a given implementation, or the shortest +! time you want the model to advance between assimilations. + +function shortest_time_between_assimilations() + +type(time_type) :: shortest_time_between_assimilations + +if ( .not. module_initialized ) call static_init_model + +shortest_time_between_assimilations = assimilation_time_step + +end function shortest_time_between_assimilations + + +!------------------------------------------------------------------ +! Given an integer index into the state vector, returns the +! associated location and optionally the physical quantity. + +subroutine get_state_meta_data(index_in, location, qty) + +integer(i8), intent(in) :: index_in +type(location_type), intent(out) :: location +integer, intent(out), optional :: qty + +integer :: lon_index, lat_index, level, local_qty + +if ( .not. module_initialized ) call static_init_model + +call get_model_variable_indices(index_in, lon_index, lat_index, level, kind_index=local_qty) + +location = set_location(geolon, geolat, real(level,r8), VERTISLEVEL) + +if (present(qty)) then + qty = local_qty +endif + + +end subroutine get_state_meta_data + +!------------------------------------------------------------------ +! write any additional attributes to the output and diagnostic files + +subroutine nc_write_model_atts(ncid, domain_id) + +integer, intent(in) :: ncid ! netCDF file identifier +integer, intent(in) :: domain_id + +if ( .not. module_initialized ) call static_init_model + +! put file into define mode. + +call nc_begin_define_mode(ncid) + +call nc_add_global_creation_time(ncid) + +call nc_add_global_attribute(ncid, "model_source", source ) +call nc_add_global_attribute(ncid, "model", "MARBL_column") + +call nc_end_define_mode(ncid) + +! Flush the buffer and leave netCDF file open +call nc_synchronize_file(ncid) + +end subroutine nc_write_model_atts + +!------------------------------------------------------------------ +! Verify that the namelist was filled in correctly, and check +! that there are valid entries for the dart_kind. +! Returns a table with columns: +! +! netcdf_variable_name ; dart_qty_string ; lowerbound ; upperbound ; update_string + +subroutine verify_state_variables(state_variables, ngood, table, qty_list, clamp_vals, update_var) + +character(len=*), intent(inout) :: state_variables(:) +integer, intent(out) :: ngood +character(len=*), intent(out) :: table(:,:) +integer, intent(out) :: qty_list(:) ! kind number +real(r8), intent(out) :: clamp_vals(:,:) +logical, intent(out) :: update_var(:) ! logical update + +integer :: nrows, i +character(len=NF90_MAX_NAME) :: varname, dartstr, lowerbound, upperbound, update +character(len=256) :: string1, string2 + +if ( .not. module_initialized ) call static_init_model + +nrows = size(table,1) + +ngood = 0 + +MyLoop : do i = 1, nrows + + varname = trim(state_variables(5*i -4)) + dartstr = trim(state_variables(5*i -3)) + lowerbound = trim(state_variables(5*i -2)) + upperbound = trim(state_variables(5*i -1)) + update = trim(state_variables(5*i )) + + call to_upper(update) + + table(i,1) = trim(varname) + table(i,2) = trim(dartstr) + table(i,3) = trim(lowerbound) + table(i,4) = trim(upperbound) + table(i,5) = trim(update) + + if ( table(i,1) == ' ' .and. table(i,2) == ' ' .and. table(i,3) == ' ' .and. table(i,4) == ' ' .and. table(i,5) == ' ') exit MyLoop ! Found end of list. + + if ( table(i,1) == ' ' .or. table(i,2) == ' ' .or. table(i,3) == ' ' .or. table(i,4) == ' ' .or. table(i,5) == ' ') then + string1 = 'model_nml:model_state_variables not fully specified' + call error_handler(E_ERR,'verify_state_variables',string1) + endif + + ! Make sure DART qty is valid + + qty_list(i) = get_index_for_quantity(dartstr) + if( qty_list(i) < 0 ) then + write(string1,'(''there is no obs_kind <'',a,''> in obs_kind_mod.f90'')') trim(dartstr) + call error_handler(E_ERR,'verify_state_variables',string1) + endif + + ! Make sure the update variable has a valid name + + select case (update) + case ('UPDATE') + update_var(i) = .true. + case ('NO_COPY_BACK') + update_var(i) = .false. + case default + write(string1,'(A)') 'only UPDATE or NO_COPY_BACK supported in model_state_variable namelist' + write(string2,'(6A)') 'you provided : ', trim(varname), ', ', trim(dartstr), ', ', trim(update) + call error_handler(E_ERR,'verify_state_variables',string1, text2=string2) + end select + + ! reading the clamp values + + if (table(i, 3) /= 'NA') then + read(table(i,3), *) clamp_vals(i,1) + else + clamp_vals(i,1) = MISSING_R8 + endif + + if (table(i,4) /= 'NA') then + read(table(i,4), *) clamp_vals(i,2) + else + clamp_vals(i,2) = MISSING_R8 + endif + + ngood = ngood + 1 +enddo MyLoop + + +end subroutine verify_state_variables + +!------------------------------------------------------------ +! Read number of vertical layers from mom6 template file +subroutine read_num_layers() + +integer :: ncid + +ncid = nc_open_file_readonly(state_template_file) + +call nc_get_variable_size(ncid, 'Layer', nz) + +call nc_close_file(ncid) + +end subroutine read_num_layers + +!------------------------------------------------------------ +subroutine read_ocean_geometry() + +integer :: ncid + +character(len=*), parameter :: routine = 'read_ocean_geometry' + +if ( .not. module_initialized ) call static_init_model + +ncid = nc_open_file_readonly(ocean_geometry) + +call nc_get_variable(ncid, 'D', basin_depth, routine) + +call nc_close_file(ncid) + +end subroutine read_ocean_geometry + +!------------------------------------------------------------------ + +subroutine end_model() + + +end subroutine end_model + +!=================================================================== +! End of model_mod +!=================================================================== +end module model_mod + diff --git a/models/MARBL_column/readme.rst b/models/MARBL_column/readme.rst new file mode 100644 index 0000000000..abdf04e99f --- /dev/null +++ b/models/MARBL_column/readme.rst @@ -0,0 +1,136 @@ +MARBL_column +============ + +**MARBL** stands for the Marine Biogeochemistry Library; it's a modular biogeochemical modeling suite for next-generatioon models. +It simulates marine ecosystem dynamics and the coupled cycles of carbon, nitrogen, phosphorus, iron, silicon, and oxygen. +It is a component of the Community Earth System Model (`CESM `_) and is often coupled to +other physical ocean models such as the Modular Ocean Model (`MOM6 `_). + +For a detailed description of the model, the reader is refered to Long et al., 2021 [1]_. + +The MARBL source code can be downloaded from https://github.com/marbl-ecosys/MARBL. You may also want to read +MARBL's documentation `here `_. + +This MARBL-DART interface was developed by `Robin Armstrong `_. Thanks Robin! + +Overview +-------- +**MARBL_column** provides an interface between DART and a 1D (column) configuration of MARBL. +It's designed to be used at locations where in-situ data is abundant such as +`Bermuda Atlantic Time-series Study (BATS) station `_, +`Weather Station Mike `_, +`Hawaii Ocean Time-series (HOT) `_ ... to name a few. + +The code is designed to perform 3 kinds of data assimilation (DA) experiments: + +#. **State Estimation:** where the prognostic state variables of MARBL such as nitrate concerntration are updated. + To achive this, you'll need to set + + ``estimate_params = .false.`` within ``&model_nml`` in the namelist file ``input.nml`` + +#. **State and Parameters Estimation:** where both the state and a set of model parameters are updated. + MARBL has a long list of uncertain model parameters that can be constrained alongside the state. + This usually improves the prediction skill of the model and alleviates some of its biases. + To achieve this DA exercise, you'll need to set + + ``estimate_params = .true.`` + + The combined DART state will be of the form :math:`Z_k = \left[ \mathbf{x}_k, \boldsymbol{\theta} \right]^T` + where :math:`Z_k` is the joint state, :math:`\mathbf{x}` and parameters, :math:`\boldsymbol{\theta}` + vector at time :math:`t_k` + +#. **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 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 +-------- +The ``&model_nml`` variables and their default values are listed here: + +.. code-block:: fortran + + &model_nml + 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, + model_state_variables = 'NO3 ', 'QTY_NITRATE_CONCENTRATION ', '0.0', 'NA', 'UPDATE ', + 'SiO3 ', 'QTY_DISSOLVED_INORGANIC_SIO3 ', '0.0', 'NA', 'UPDATE ', + 'PO4 ', 'QTY_PHOSPHATE_CONCENTRATION ', '0.0', 'NA', 'UPDATE ', + 'Fe ', 'QTY_DISSOLVED_INORGANIC_IRON ', '0.0', 'NA', 'UPDATE ', + 'DIC ', 'QTY_DISSOLVED_INORGANIC_CARBON', '0.0', 'NA', 'UPDATE ', + 'O2 ', 'QTY_DISSOLVED_OXYGEN ', '0.0', 'NA', 'UPDATE ', + 'DOC ', 'QTY_DISSOLVED_ORGANIC_CARBON ', '0.0', 'NA', 'UPDATE ', + 'DON ', 'QTY_DISSOLVED_ORGANIC_NITROGEN', '0.0', 'NA', 'UPDATE ', + 'DOP ', 'QTY_DISSOLVED_ORGANIC_P ', '0.0', 'NA', 'UPDATE ', + 'ALK ', 'QTY_ALKALINITY ', '0.0', 'NA', 'UPDATE ', + 'microzooC', 'QTY_MICROZOOPLANKTON_CARBON ', '0.0', 'NA', 'UPDATE ', + 'mesozooC ', 'QTY_MESOZOOPLANKTON_CARBON ', '0.0', 'NA', 'UPDATE ', + 'h ', 'QTY_LAYER_THICKNESS ', '0.0', 'NA', 'NO_COPY_BACK' + estimate_params = .true. + model_parameters = 'autotroph_settings(1)%kDOP ', 'QTY_BGC_PARAM', '0.0', 'NA', 'UPDATE' + / + +This namelist provides control over the kind of DA experiment as described abvove. + ++-------------------------------------+--------------------+------------------------------------------------------------+ +| Item | Type | Description | ++=====================================+====================+============================================================+ +| ``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`` | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``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`` | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``station_location`` | real(2) | Longitude and latitude of the ocean column location. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``time_step_days`` | integer | The number of days to advance the model for each | +| | | assimilation. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``time_step_seconds`` | integer | In addition to ``time_step_days``, the number | +| | | of seconds to advance the model for each assimilation. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``model_state_variables`` | character(:,5) | Strings that associate MARBL variables with a DART | +| | | quantity. They also describe their physical bounds and | +| | | whether or not to write the updated values to the restart | +| | | files. These variables will be read from the MARBL restart | +| | | file and modified by the assimilation. Some (perhaps all) | +| | | will be used by the forward observation operators. If the | +| | | 5th column is ``UPDATE``, the output files will have the | +| | | modified (assimilated,posterior) values. If the 5th | +| | | column is ``NO_COPY_BACK``, that variable will not be | +| | | written to the restart files. **The DART diagnostic files | +| | | will always have the (modified) posterior values.** | +| | | Diagnostic variables that are useful for the calculation | +| | | of the forward observation operator but have no impact on | +| | | the forecast trajectory of the model could have a value of | +| | | ``NO_COPY_BACK``. The 3rd and 4th column list the minimum | +| | | and maximum allowed values for the updated variables. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``estimate_params`` | logical | A switch to turn on/off parameter estimation. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``model_parameters`` | character(:,5) | Similar to ``model_state_variables``, this is a list of | +| | | parameters that will take part in the DART state and | +| | | would possibly get updated. | ++-------------------------------------+--------------------+------------------------------------------------------------+ + + +References +---------- +.. [1] Long, Matthew C., J. Keith Moore, Keith Lindsay, Michael Levy, Scott C. Doney, + Jessica Y. Luo, Kristen M. Krumhardt, Robert T. Letscher, Maxwell Grover, and Zephyr T. Sylvester. + "Simulations with the marine biogeochemistry library (MARBL)." + Journal of Advances in Modeling Earth Systems 13, no. 12 (2021): e2021MS002647. diff --git a/models/MARBL_column/work/input.nml b/models/MARBL_column/work/input.nml new file mode 100644 index 0000000000..c2c1035f47 --- /dev/null +++ b/models/MARBL_column/work/input.nml @@ -0,0 +1,255 @@ + +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + +&model_nml + 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, + time_step_seconds = 0, + model_state_variables = 'NO3 ', 'QTY_NITRATE_CONCENTRATION ', '0.0', 'NA', 'UPDATE ', + 'SiO3 ', 'QTY_DISSOLVED_INORGANIC_SIO3 ', '0.0', 'NA', 'UPDATE ', + 'PO4 ', 'QTY_PHOSPHATE_CONCENTRATION ', '0.0', 'NA', 'UPDATE ', + 'Fe ', 'QTY_DISSOLVED_INORGANIC_IRON ', '0.0', 'NA', 'UPDATE ', + 'DIC ', 'QTY_DISSOLVED_INORGANIC_CARBON', '0.0', 'NA', 'UPDATE ', + 'O2 ', 'QTY_DISSOLVED_OXYGEN ', '0.0', 'NA', 'UPDATE ', + 'DOC ', 'QTY_DISSOLVED_ORGANIC_CARBON ', '0.0', 'NA', 'UPDATE ', + 'DON ', 'QTY_DISSOLVED_ORGANIC_NITROGEN', '0.0', 'NA', 'UPDATE ', + 'DOP ', 'QTY_DISSOLVED_ORGANIC_P ', '0.0', 'NA', 'UPDATE ', + 'ALK ', 'QTY_ALKALINITY ', '0.0', 'NA', 'UPDATE ', + 'microzooC', 'QTY_MICROZOOPLANKTON_CARBON ', '0.0', 'NA', 'UPDATE ', + 'mesozooC ', 'QTY_MESOZOOPLANKTON_CARBON ', '0.0', 'NA', 'UPDATE ', + 'h ', 'QTY_LAYER_THICKNESS ', '0.0', 'NA', 'NO_COPY_BACK' + estimate_params = .true. + model_parameters = 'autotroph_settings(1)%kDOP ', 'QTY_BGC_PARAM', '0.0', 'NA', 'UPDATE' + / + +&filter_nml + single_file_in = .false., + input_state_files = '', + input_state_file_list = 'states_input.txt', 'params_input.txt' + + stages_to_write = 'preassim', 'analysis', 'output', + + single_file_out = .false., + output_state_files = '', + output_state_file_list = 'states_output.txt', 'params_output.txt' + output_interval = 1, + output_members = .true., + num_output_state_members = 3, + output_mean = .true., + output_sd = .true., + write_all_stages_at_end = .true., + + ens_size = 3, + num_groups = 1, + perturb_from_single_instance = .false., + perturbation_amplitude = 1e-2, + distributed_state = .true., + + async = 0, + adv_ens_command = "./advance_model.csh", + + obs_sequence_in_name = 'BATS_147612.out', + obs_sequence_out_name = 'obs_seq.final', + num_output_obs_members = 3, + init_time_days = -1, + init_time_seconds = -1, + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + + inf_flavor = 5, 0, + inf_initial_from_restart = .true., .false., + inf_sd_initial_from_restart = .true., .false., + inf_deterministic = .true., .true., + inf_initial = 1.0, 1.0, + inf_lower_bound = 0.0, 1.0, + inf_upper_bound = 100.0, 1000000.0, + inf_damping = 0.9, 1.0, + inf_sd_initial = 0.6, 0.0, + inf_sd_lower_bound = 0.6, 0.0, + inf_sd_max_change = 1.05, 1.05, + + trace_execution = .true., + output_timestamps = .false., + output_forward_op_errors = .false., + silence = .false., + / + +&fill_inflation_restart_nml + write_prior_inf = .true., + prior_inf_mean = 1.0, + prior_inf_sd = 0.6, + + write_post_inf = .false., + post_inf_mean = 1.00, + post_inf_sd = 0.6, + + input_state_files = 'ensemble/background/member_0001/RESTART/MOM.res.nc', 'ensemble/background/member_0001/marbl_params.nc', + single_file = .false., + verbose = .false. + / + +&quality_control_nml + input_qc_threshold = 3.0, + outlier_threshold = 3.0, +/ + +&location_nml + / + +&perfect_model_obs_nml + read_input_state_from_file = .true., + single_file_in = .false., + input_state_files = "MOM.res.nc", + + write_output_state_to_file = .true., + single_file_out = .false., + output_state_files = "MOM.res.output.nc", + output_interval = 1, + + async = 0, + adv_ens_command = "./advance_model.csh", + + obs_seq_in_file_name = "obs_seq.in", + obs_seq_out_file_name = "obs_seq.out", + init_time_days = 0, + init_time_seconds = 0, + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + + trace_execution = .false., + output_timestamps = .false., + print_every_nth_obs = -1, + output_forward_op_errors = .false., + silence = .false., + / + +&ensemble_manager_nml + / + +&assim_tools_nml + cutoff = 1000000.0, + sort_obs_inc = .false., + spread_restoration = .false., + sampling_error_correction = .false., + adaptive_localization_threshold = -1, + distribute_mean = .false. + output_localization_diagnostics = .false., + localization_diagnostics_file = 'localization_diagnostics', + print_every_nth_obs = 0 + / + +&cov_cutoff_nml + select_localization = 1 + / + +®_factor_nml + select_regression = 1, + input_reg_file = "time_mean_reg", + save_reg_diagnostics = .false., + reg_diagnostics_file = "reg_diagnostics" + / + +&obs_sequence_nml + write_binary_obs_sequence = .false. + / + +&obs_kind_nml + 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' + / + +&utilities_nml + TERMLEVEL = 1, + module_details = .false., + logfilename = 'dart_log.out', + nmlfilename = 'dart_log.nml', + write_nml = 'none' + / + +&preprocess_nml + input_obs_def_mod_file = '../../../observations/forward_operators/DEFAULT_obs_def_mod.F90' + output_obs_def_mod_file = '../../../observations/forward_operators/obs_def_mod.f90' + input_obs_qty_mod_file = '../../../assimilation_code/modules/observations/DEFAULT_obs_kind_mod.F90' + output_obs_qty_mod_file = '../../../assimilation_code/modules/observations/obs_kind_mod.f90' + obs_type_files = '../../../observations/forward_operators/obs_def_ocean_mod.f90' + quantity_files = '../../../assimilation_code/modules/observations/default_quantities_mod.f90', + '../../../assimilation_code/modules/observations/ocean_quantities_mod.f90' + / + +&obs_sequence_tool_nml + filename_seq = 'obs_seq.one', 'obs_seq.two', + filename_out = 'obs_seq.processed', + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + print_only = .false., + gregorian_cal = .false. + / + +&obs_diag_nml + obs_sequence_name = '' + obs_sequence_list = 'obs_seq_list.txt' + first_bin_center = 2005, 2, 24, 0, 0, 0 + last_bin_center = 2010, 2, 28, 0, 0, 0 + bin_separation = 0, 0, 1, 0, 0, 0 + bin_width = 0, 0, 1, 0, 0, 0 + time_to_skip = 0, 0, 0, 0, 0, 0 + max_num_bins = 1000 + plevel = -888888.0 + hlevel = -1, -5, -10, -50, -100, -250, -500, -1000, -2000, -3000, -4000, -5000, -6000 + mlevel = -888888 + plevel_edges = -888888.0 + hlevel_edges = -888888.0 + mlevel_edges = -888888 + Nregions = 1 + lonlim1 = 0 + lonlim2 = 360 + latlim1 = -90 + latlim2 = 90 + reg_names = 'null' + trusted_obs = 'null' + create_rank_histogram = .true. + outliers_in_histogram = .false. + use_zero_error_obs = .false. + verbose = .false. + / + +&state_vector_io_nml + / + +&model_mod_check_nml + 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 + loc_of_interest = 296., 31., -500.0 + quantity_of_interest = 'QTY_DISSOLVED_OXYGEN' + interp_test_vertcoord = 'VERTISHEIGHT' + interp_test_lonrange = 295., 297. + interp_test_dlon = 0.10 + interp_test_latrange = 30., 32. + interp_test_dlat = 0.10 + interp_test_vertrange = 0.0, -100.00 + interp_test_dvert = -10.0 + verbose = .true. + / diff --git a/models/MARBL_column/work/quickbuild.sh b/models/MARBL_column/work/quickbuild.sh new file mode 100755 index 0000000000..0e847b5a03 --- /dev/null +++ b/models/MARBL_column/work/quickbuild.sh @@ -0,0 +1,59 @@ +#!/usr/bin/env bash + +# DART software - Copyright UCAR. This open source software is provided +# by UCAR, "as is", without charge, subject to all terms of use at +# http://www.image.ucar.edu/DAReS/DART/DART_download + +main() { + +export DART=$(git rev-parse --show-toplevel) +source "$DART"/build_templates/buildfunctions.sh + +MODEL=MARBL_column +LOCATION=threed_sphere + + +programs=( +closest_member_tool +filter +model_mod_check +perfect_model_obs +) + +serial_programs=( +create_fixed_network_seq +create_obs_sequence +fill_inflation_restart +integrate_model +obs_common_subset +obs_diag +obs_sequence_tool +) + +model_programs=( +) + +model_serial_programs=( +) + +# quickbuild arguments +arguments "$@" + +# clean the directory +\rm -f -- *.o *.mod Makefile .cppdefs + +# build any NetCDF files from .cdl files +cdl_to_netcdf + +# build and run preprocess before making any other DART executables +buildpreprocess + +# build +buildit + +# clean up +\rm -f -- *.o *.mod + +} + +main "$@" diff --git a/models/README.rst b/models/README.rst index 73b4dd87bd..932d60424e 100644 --- a/models/README.rst +++ b/models/README.rst @@ -27,6 +27,7 @@ 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` diff --git a/observations/forward_operators/obs_def_ocean_mod.f90 b/observations/forward_operators/obs_def_ocean_mod.f90 index daf95fc0b3..cc60436583 100644 --- a/observations/forward_operators/obs_def_ocean_mod.f90 +++ b/observations/forward_operators/obs_def_ocean_mod.f90 @@ -70,6 +70,14 @@ ! FERRYBOX_SALINITY, QTY_SALINITY, COMMON_CODE ! FERRYBOX_TEMPERATURE, QTY_TEMPERATURE, COMMON_CODE ! OCEAN_COLOR, QTY_SURFACE_CHLOROPHYLL, COMMON_CODE +! BATS_OXYGEN, QTY_DISSOLVED_OXYGEN, COMMON_CODE +! BATS_ALKALINITY, QTY_ALKALINITY, COMMON_CODE +! BATS_ORGANIC_CARBON, QTY_DISSOLVED_ORGANIC_CARBON, COMMON_CODE +! BATS_INORGANIC_CARBON, QTY_DISSOLVED_INORGANIC_CARBON, COMMON_CODE +! BATS_NITROGEN, QTY_DISSOLVED_ORGANIC_NITROGEN, COMMON_CODE +! BATS_NITRATE, QTY_NITRATE_CONCENTRATION, COMMON_CODE +! BATS_SILICATE, QTY_DISSOLVED_INORGANIC_SIO3, COMMON_CODE +! BATS_PHOSPHATE, QTY_PHOSPHATE_CONCENTRATION, COMMON_CODE ! END DART PREPROCESS TYPE DEFINITIONS ! From Ibrahim - 19 May 2009 diff --git a/observations/obs_converters/BATS/bats_climatology.py b/observations/obs_converters/BATS/bats_climatology.py new file mode 100644 index 0000000000..a184869133 --- /dev/null +++ b/observations/obs_converters/BATS/bats_climatology.py @@ -0,0 +1,142 @@ +#!/usr/bin/env python3 +import numpy as np +import netCDF4 as nc +import matplotlib.pyplot as plt + +################################################################# +####################### SCRIPT PARAMETERS ####################### +################################################################# + +# input and output datafiles +input_file_path = 'bats_bottle.txt' +output_txt_path = 'bats_climatology.txt' +output_nc_path = 'bats_climatology.nc' + +first_data_line = 61 # file is read starting at this line (so that header information is ignored) +month_columns = [18, 19] # first and last columns in the data file containing month values +depth_columns = [64, 69] # first and last columns in the data file containing depth values +num_obs_types = 6 # number of variables being read from the datafile +max_num_data = 16660 # upper bound on the number of data points to be processed + +# first and last columns in the data file containing values for each variable +obs_val_columns = [[113, 119], # O2 + [137, 143], # CO2 + [145, 151], # Alk + [153, 159], # NO31 + [170, 176], # PO41 + [178, 184]] # Si1 + +# layer interface depths for data in prog_z.nc output file from MARBL, run "ncdump -v "z_i" prog_z.nc" to generate. +marbl_depths = [0, 5, 15, 25, 40, 62.5, 87.5, 112.5, 137.5, 175, 225, 275, 350, 450, + 550, 650, 750, 850, 950, 1050, 1150, 1250, 1350, 1450, 1625, 1875, 2250, + 2750, 3250, 3750, 4250, 4750, 5250, 5750, 6250] + +# set to 'True' in order to produce a histogram of BATS observation vs depth for each month +plot_sample_hist = False + +# human-readable names of the variables being read from the file +obs_names = ["O2", "CO2", "Alk", "NO31", "PO41", "Si1"] + +################################################################# +####################### MAIN PROGRAM ############################ +################################################################# + +clim_means = np.zeros((num_obs_types, 12, len(marbl_depths))) +sq_means = np.zeros((num_obs_types, 12, len(marbl_depths))) +sample_counts = np.zeros((num_obs_types, 12, len(marbl_depths)), dtype=int) + +with open(input_file_path, "r") as bats_data: + for line in bats_data.readlines()[first_data_line-1:]: + + # extracting the depth and month where this observation was taken + depth_val = float(line[depth_columns[0]-1:depth_columns[1]]) + month_index = int(line[month_columns[0]-1:month_columns[1]]) - 1 + + # assigning this observation to one of the MARBL layers + depth_index = len(marbl_depths) - 1 + + while((depth_val < marbl_depths[depth_index]) and (depth_index > 0)): + depth_index -= 1 + + # adding this observation into the climatology + for obs_type_index in range(num_obs_types): + obs_val = float(line[obs_val_columns[obs_type_index][0]-1:obs_val_columns[obs_type_index][1]]) + + if(abs(obs_val + 999.0) > 1e-8): # this filters out 'missing' values + prev_mean = clim_means[obs_type_index, month_index, depth_index] + prev_sqmean = sq_means[obs_type_index, month_index, depth_index] + prev_samples = sample_counts[obs_type_index, month_index, depth_index] + + new_mean = (prev_samples*prev_mean + obs_val)/(prev_samples + 1) + new_sqmean = (prev_samples*prev_sqmean + obs_val**2)/(prev_samples + 1) + + clim_means[obs_type_index, month_index, depth_index] = new_mean + sq_means[obs_type_index, month_index, depth_index] = new_sqmean + + sample_counts[obs_type_index, month_index, depth_index] += 1 + + print("month =",month_index,", type =",obs_type_index,", depth =",marbl_depths[depth_index],", value =",obs_val) + + +# setting up a CSV file to write the climatology into; this will be read by the DART data converter. +output_txt = open(output_txt_path, "w") + +# setting up a netCDF file to write the climatology into; this will be useful for creating MARBL-DART diagnostics. +output_nc = nc.Dataset(output_nc_path, "w") + +output_nc.createDimension("Month") +output_nc.createDimension("Layer") +output_nc.createVariable("Depth", "double", ("Layer",)) +output_nc["Depth"][:] = marbl_depths + +for obs_type_index in range(num_obs_types): + output_nc.createVariable(obs_names[obs_type_index]+"_value", "double", ("Month", "Layer",)) + output_nc.createVariable(obs_names[obs_type_index]+"_error_sd", "double", ("Month", "Layer",)) + output_nc.createVariable(obs_names[obs_type_index]+"_samples", "int", ("Month", "Layer",)) + +for month_index in range(12): + for obs_type_index in range(num_obs_types): + for depth_index in range(len(marbl_depths)): + samples = sample_counts[obs_type_index, month_index, depth_index] + output_nc[obs_names[obs_type_index]+"_samples"][month_index, depth_index] = samples + + if(samples > 0): + depth = marbl_depths[depth_index] + mean = clim_means[obs_type_index, month_index, depth_index] + sq_mean = sq_means[obs_type_index, month_index, depth_index] + + variability_err = (sq_mean - mean**2)/samples + obs_error = (.1*mean)**2 + uncertainty = np.sqrt(variability_err + obs_error) + + output_txt.write(str(month_index+1)+","+str(obs_type_index+1)+","+str(depth)+","+str(mean)+","+str(uncertainty)+"\n") + + output_nc[obs_names[obs_type_index]+"_value"][month_index, depth_index] = mean + output_nc[obs_names[obs_type_index]+"_error_sd"][month_index, depth_index] = uncertainty + +output_txt.close() +output_nc.close() + +if(plot_sample_hist): + month_names = ["January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December"] + month_plot = 12*[None] + + fig = plt.figure(figsize = (6, 70)) + + for month_index in range(12): + month_plot[month_index] = fig.add_subplot(12, 1, month_index + 1) + month_plot[month_index].set_xlabel("Depth (meters)") + month_plot[month_index].set_xscale("log") + month_plot[month_index].set_ylabel("Number of Data Points ("+month_names[month_index]+")") + month_plot[month_index].set_yscale("log") + + (month_index == 0) and month_plot[month_index].set_title("Depth Distribution of BATS Data") + + for obs_type_index in range(num_obs_types): + month_plot[month_index].plot(marbl_depths, sample_counts[obs_type_index, month_index, :], marker = "o", markerfacecolor = "none", label = obs_names[obs_type_index]) + + month_plot[month_index].legend() + + plt.legend() + plt.savefig("bats_depth_hist.pdf") + plt.close(fig) diff --git a/observations/obs_converters/BATS/bats_to_clim_obs.f90 b/observations/obs_converters/BATS/bats_to_clim_obs.f90 new file mode 100644 index 0000000000..60d1ad2b42 --- /dev/null +++ b/observations/obs_converters/BATS/bats_to_clim_obs.f90 @@ -0,0 +1,247 @@ +! DART software - Copyright UCAR. This open source software is provided +! by UCAR, "as is", without charge, subject to all terms of use at +! http://www.image.ucar.edu/DAReS/DART/DART_download +! + +! This file is meant to read a text file containing bottle data from the +! Bermuda Atlantic Time-Series Study (https://bats.bios.asu.edu/), which +! is converted to climatological estimates with monthly resolution. + +program bats_to_clim_obs + +use types_mod, only : r8 +use utilities_mod, only : initialize_utilities, finalize_utilities, & + open_file, close_file, & + find_namelist_in_file, check_namelist_read, & + error_handler, E_ERR, E_MSG, nmlfileunit, & + do_nml_file, do_nml_term +use time_manager_mod, only : time_type, set_calendar_type, set_date, & + operator(>=), increment_time, get_time, & + operator(-), GREGORIAN, operator(+), & + print_date +use location_mod, only : VERTISHEIGHT, VERTISPRESSURE +use obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq, & + static_init_obs_sequence, init_obs, write_obs_seq, & + init_obs_sequence, get_num_obs, set_copy_meta_data, & + set_qc_meta_data, destroy_obs_sequence +use obs_utilities_mod, only : create_3d_obs, add_obs_to_seq +use obs_kind_mod, only : BATS_OXYGEN, BATS_INORGANIC_CARBON, BATS_ALKALINITY, & + BATS_NITRATE, BATS_PHOSPHATE, BATS_SILICATE + +implicit none + +integer, parameter :: NUM_SCALAR_OBS = 6 ! maximum number of scalar observation variables that will + ! be assimilated at each observation. + +! this array defines the order in which observations are read from the file +integer, parameter :: OTYPE_ORDERING(NUM_SCALAR_OBS) & + = (/BATS_OXYGEN, BATS_INORGANIC_CARBON, BATS_ALKALINITY, BATS_NITRATE, & + BATS_PHOSPHATE, BATS_SILICATE/) + +real(r8), parameter :: MIN_OBS_ERROR = 0.1_r8 +real(r8), parameter :: BATS_LON = 360.0_r8 - 64.0_r8 +real(r8), parameter :: BATS_LAT = 31.0_r8 + +! namelist variables, changeable at runtime +character(len=256) :: text_input_file = 'bats_climatology.txt' +character(len=256) :: obs_out_dir = 'obs_seq_files' +integer :: max_lines = 3000 +real(r8) :: obs_err_var_inflation = 10.0 +logical :: debug = .true. + +namelist /bats_to_clim_obs_nml/ text_input_file, max_lines, obs_err_var_inflation, obs_out_dir, debug + +! local variables +character (len=294) :: input_line, obs_out_file +character (len=6) :: month_str + +integer :: oday, osec, rcio, iunit, char_index, comma_index, line_number, otype_index, & + month, month_old, num_copies, num_qc, max_obs +integer :: comma_locs(4) +logical :: first_obs, new_obs_seq + +real(r8) :: qc, obs_err +real(r8) :: vert, ovalue + +type(obs_sequence_type) :: obs_seq +type(obs_type) :: obs, prev_obs +type(time_type) :: time_obs, prev_time + + +! start of executable code + +call initialize_utilities('bats_to_clim_obs') + +! Read the namelist entries +call find_namelist_in_file("input.nml", "bats_to_clim_obs_nml", iunit) +read(iunit, nml = bats_to_clim_obs_nml, iostat = rcio) +call check_namelist_read(iunit, rcio, "bats_to_clim_obs_nml") + +! Record the namelist values used for the run +if (do_nml_file()) write(nmlfileunit, nml=bats_to_clim_obs_nml) +if (do_nml_term()) write( * , nml=bats_to_clim_obs_nml) + +! time setup +call set_calendar_type(GREGORIAN) + +! open input text file + +iunit = open_file(text_input_file, 'formatted', 'read') +if (debug) print *, 'opened input file ' // trim(text_input_file) + +max_obs = NUM_SCALAR_OBS*max_lines +num_copies = 1 +num_qc = 1 + +! call the initialization code, and initialize two empty observation types +call static_init_obs_sequence() +call init_obs(obs, num_copies, num_qc) +call init_obs(prev_obs, num_copies, num_qc) +first_obs = .true. + +! Set the DART data quality control. +qc = 0.0_r8 + +! used later to manage splitting of data into obs-sequence files +month_old = 0 + +line_number = 0 ! counts the number of lines that have been read so far + +obsloop: do ! no end limit - have the loop break when input ends + ! read in entire text line into a buffer + read(iunit, "(A)", iostat=rcio) input_line + line_number = line_number + 1 + + if (rcio /= 0) then + if (debug) print *, 'got bad read code from input file at line ',line_number,', rcio = ', rcio + exit obsloop + endif + + ! finding the locations of commas within the string + + comma_index = 0 + char_index = 0 + + do while(comma_index < 4) + char_index = char_index + 1 + + if(input_line(char_index:char_index) == ',') then + comma_index = comma_index + 1 + comma_locs(comma_index) = char_index + end if + end do + + ! reading the month value of this line + + read(input_line(1:(comma_locs(1) - 1)), *, iostat=rcio) month + if(rcio /= 0) then + if(debug) print *, "got bad read code getting month at line ",line_number,", rcio = ",rcio + exit obsloop + end if + + ! reading the observation type + + read(input_line((comma_locs(1) + 1):(comma_locs(2) - 1)), *, iostat=rcio) otype_index + if(rcio /= 0) then + if(debug) print *, "got bad read code getting obs type at line ",line_number,", rcio = ",rcio + exit obsloop + end if + + ! reading the observation depth + + read(input_line((comma_locs(2) + 1):(comma_locs(3) - 1)), *, iostat=rcio) vert + if(rcio /= 0) then + if(debug) print *, "got bad read code getting vert at line ",line_number,", rcio = ",rcio + exit obsloop + end if + + ! reading the observtaion value + + read(input_line((comma_locs(3) + 1):(comma_locs(4) - 1)), *, iostat=rcio) ovalue + if(rcio /= 0) then + if(debug) print *, "got bad read code getting observation value at line ",line_number,", rcio = ",rcio + exit obsloop + end if + + ! reading the observtaion uncertainty (standard deviation) + + read(input_line((comma_locs(4) + 1):len(input_line)), *, iostat=rcio) obs_err + if(rcio /= 0) then + if(debug) print *, "got bad read code getting observation uncertainty at line ",line_number,", rcio = ",rcio + exit obsloop + end if + + ! inflating the observation error according to the number of cycles in MARBL-DART + + obs_err = obs_err*sqrt(obs_err_var_inflation) + + ! 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 + obs_err = obs_err*1.026 + end if + + ! For ensemble smoothing, each observation is considered to be a measurement + ! of the same climatological "year", which is arbitrarily labeled as year 1601 + ! (the first year in DART's internal calendar). Observations are given identical + ! day/hour/minute timestamps to avoid time-ordering errors in the obs-seq file. + + time_obs = set_date(1601, 1, 1, hours=0, minutes=0) + call get_time(time_obs, osec, days=oday) + + if(debug) then + print *, "adding climatological mean for month ",month + print *, " \__ type: ",OTYPE_ORDERING(otype_index) + print *, " vert: ",vert + print *, " value: ",ovalue + print *, " uncertainty: ",obs_err + print *, " oday: ",oday + print *, " osec: ",osec + end if + + new_obs_seq = (first_obs .or. (month /= month_old)) + month_old = month + + ! if necessary, saving the old observation sequence and beginning a new one. + if(new_obs_seq) then + ! dumping the observations so far into their own file + if ( (.not. first_obs) .and. (get_num_obs(obs_seq) > 0) ) then + if (debug) print *, 'writing obs_seq, obs_count = ', get_num_obs(obs_seq) + call write_obs_seq(obs_seq, obs_out_file) + call destroy_obs_sequence(obs_seq) + first_obs = .true. + endif + + write(month_str, "(I0.2)") month + obs_out_file = trim(obs_out_dir)//"/BATS_clim_"//trim(month_str)//".out" + + ! create a new, empty obs_seq file. + call init_obs_sequence(obs_seq, num_copies, num_qc, max_obs) + + ! the first one needs to contain the string 'observation' and the + ! second needs the string 'QC'. + call set_copy_meta_data(obs_seq, 1, 'observation') + call set_qc_meta_data(obs_seq, 1, 'Data QC') + end if + + ! adding the observation + call create_3d_obs(BATS_LAT, BATS_LON, vert, VERTISHEIGHT, & + ovalue, OTYPE_ORDERING(otype_index), obs_err, & + oday, osec, qc, obs) + + call add_obs_to_seq(obs_seq, obs, time_obs, prev_obs, prev_time, first_obs) +end do obsloop + +! putting any remaining observations into an obs sequence file + +if ( (.not. first_obs) .and. (get_num_obs(obs_seq) > 0) ) then + if (debug) print *, 'writing obs_seq, obs_count = ', get_num_obs(obs_seq) + call write_obs_seq(obs_seq, obs_out_file) + call destroy_obs_sequence(obs_seq) +endif + +! end of main program +call finalize_utilities() + +end program bats_to_clim_obs diff --git a/observations/obs_converters/BATS/bats_to_obs.f90 b/observations/obs_converters/BATS/bats_to_obs.f90 new file mode 100644 index 0000000000..e6e0059206 --- /dev/null +++ b/observations/obs_converters/BATS/bats_to_obs.f90 @@ -0,0 +1,323 @@ +! DART software - Copyright UCAR. This open source software is provided +! by UCAR, "as is", without charge, subject to all terms of use at +! http://www.image.ucar.edu/DAReS/DART/DART_download +! + +! This file is meant to read a text file containing bottle data from the +! Bermuda Atlantic Time-Series Study (https://bats.bios.asu.edu/). + +program bats_to_obs + +use types_mod, only : r8 +use utilities_mod, only : initialize_utilities, finalize_utilities, & + open_file, close_file, & + find_namelist_in_file, check_namelist_read, & + error_handler, E_ERR, E_MSG, nmlfileunit, & + do_nml_file, do_nml_term +use time_manager_mod, only : time_type, set_calendar_type, set_date, & + operator(>=), increment_time, get_time, & + operator(-), GREGORIAN, operator(+), & + print_date +use location_mod, only : VERTISHEIGHT, VERTISPRESSURE +use obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq, & + static_init_obs_sequence, init_obs, write_obs_seq, & + init_obs_sequence, get_num_obs, set_copy_meta_data, & + set_qc_meta_data, destroy_obs_sequence +use obs_utilities_mod, only : create_3d_obs, add_obs_to_seq +use obs_kind_mod, only : BATS_OXYGEN, BATS_INORGANIC_CARBON, BATS_ALKALINITY, & + BATS_NITRATE, BATS_PHOSPHATE, BATS_SILICATE + +implicit none + +integer, parameter :: NUM_SCALAR_OBS = 6 ! maximum number of scalar observation variables that will + ! be assimilated at each observation. + +! this array defines the order in which observations are read from the file +integer, parameter :: OTYPE_ORDERING(NUM_SCALAR_OBS) & + = (/BATS_OXYGEN, BATS_INORGANIC_CARBON, BATS_ALKALINITY, BATS_NITRATE, & + BATS_PHOSPHATE, BATS_SILICATE/) + +real(r8), parameter :: MIN_OBS_ERROR = 0.1_r8 +real(r8), parameter :: BATS_LON = 360.0_r8 - 64.0_r8 +real(r8), parameter :: BATS_LAT = 31.0_r8 + +! namelist variables, changeable at runtime +character(len=256) :: text_input_file = 'bats_bottle.txt' +character(len=256) :: obs_out_dir = 'obs_seq_files' +integer :: max_lines = 68000 +integer :: read_starting_at_line = 61 +integer :: date_firstcol = 14 +integer :: hourminute_firstcol = 35 +integer :: lat_cols(2) = (/42, 47/) +integer :: lon_cols(2) = (/51, 56/) +integer :: vert_cols(2) = (/64, 69/) +real(r8) :: obs_uncertainties(NUM_SCALAR_OBS) = 0.2_r8 +logical :: debug = .true. +integer :: scalar_obs_cols(2, NUM_SCALAR_OBS) = reshape( (/ & + 113, 137, 145, 153, 170, 178, & + 119, 143, 151, 159, 176, 184 /), & + shape(scalar_obs_cols), order=(/2,1/) ) + +namelist /bats_to_obs_nml/ text_input_file, max_lines, read_starting_at_line, date_firstcol, & + hourminute_firstcol, lat_cols, lon_cols, vert_cols, scalar_obs_cols, & + obs_uncertainties, obs_out_dir, debug + +! local variables +character (len=294) :: input_line, obs_out_file +character (len=6) :: daystr + +integer :: oday, day_bin, day_bin_old, osec, rcio, iunit, line_number, otype_index +integer :: year, month, day, hour, minute, hourminute_raw, date_raw +integer :: num_copies, num_qc, max_obs +integer :: num_processed(NUM_SCALAR_OBS) + +logical :: first_obs, new_obs_seq + +real(r8) :: qc, obs_err +real(r8) :: vert, ovalue +real(r8) :: running_sum(NUM_SCALAR_OBS), running_sqsum(NUM_SCALAR_OBS), & + maxvals(NUM_SCALAR_OBS), minvals(NUM_SCALAR_OBS) + +type(obs_sequence_type) :: obs_seq +type(obs_type) :: obs, prev_obs +type(time_type) :: time_obs, prev_time + + +! start of executable code + +call initialize_utilities('bats_to_obs') + +! Read the namelist entries +call find_namelist_in_file("input.nml", "bats_to_obs_nml", iunit) +read(iunit, nml = bats_to_obs_nml, iostat = rcio) +call check_namelist_read(iunit, rcio, "bats_to_obs_nml") + +! Record the namelist values used for the run +if (do_nml_file()) write(nmlfileunit, nml=bats_to_obs_nml) +if (do_nml_term()) write( * , nml=bats_to_obs_nml) + +! time setup +call set_calendar_type(GREGORIAN) + +! open input text file + +iunit = open_file(text_input_file, 'formatted', 'read') +if (debug) print *, 'opened input file ' // trim(text_input_file) + +max_obs = NUM_SCALAR_OBS*max_lines +num_copies = 1 +num_qc = 1 + +! call the initialization code, and initialize two empty observation types +call static_init_obs_sequence() +call init_obs(obs, num_copies, num_qc) +call init_obs(prev_obs, num_copies, num_qc) +first_obs = .true. + +! Set the DART data quality control. +qc = 0.0_r8 + +! used later to manage splitting of data into obs-sequence files +day_bin_old = 0 + +line_number = 0 ! counts the number of lines that have been read so far + +do otype_index = 1, NUM_SCALAR_OBS + running_sum(otype_index) = 0.0_r8 ! these arrays will be used to compute means, + running_sqsum(otype_index) = 0.0_r8 ! variances, mins, and max's of observation values. + num_processed(otype_index) = 0 + maxvals(otype_index) = 0.0_r8 + minvals(otype_index) = 100000.0_r8 +end do + +obsloop: do ! no end limit - have the loop break when input ends + ! read in entire text line into a buffer + read(iunit, "(A)", iostat=rcio) input_line + line_number = line_number + 1 + + if(line_number < read_starting_at_line) then + cycle obsloop + end if + + if (rcio /= 0) then + if (debug) print *, 'got bad read code from input file at line ',line_number,', rcio = ', rcio + exit obsloop + endif + + ! extracting the date when the observation was taken + + read(input_line(date_firstcol:(date_firstcol + 7)), *, iostat=rcio) date_raw + if(rcio /= 0) then + if(debug) print *, "got bad read code getting date at line ",line_number,", rcio = ",rcio + exit obsloop + + else if(date_raw == -999) then + cycle obsloop ! missing date + + end if + + read(input_line(date_firstcol:(date_firstcol + 3)), *, iostat=rcio) year + if(rcio /= 0) then + if(debug) print *, "got bad read code getting year at line ",line_number,", rcio = ",rcio + exit obsloop + end if + + read(input_line((date_firstcol + 4):(date_firstcol + 5)), *, iostat=rcio) month + if(rcio /= 0) then + if(debug) print *, "got bad read code getting month at line ",line_number,", rcio = ",rcio + exit obsloop + end if + + read(input_line((date_firstcol + 6):(date_firstcol + 7)), *, iostat=rcio) day + if(rcio /= 0) then + if(debug) print *, "got bad read code getting day at line ",line_number,", rcio = ",rcio + exit obsloop + end if + + if((month == 02) .and. (day == 29)) then ! we do not assimilate data taken on leap days + cycle obsloop + end if + + ! extracting the time of day when the observation was taken + + read(input_line(hourminute_firstcol:(hourminute_firstcol + 3)), *, iostat=rcio) hourminute_raw + if(rcio /= 0) then + if(debug) print *, "got bad read code parsing raw hour-minute at line ",line_number,", rcio = ",rcio + exit obsloop + + else if(hourminute_raw == -999) then + cycle obsloop ! missing timestamp + + else if(hourminute_raw < 60) then + hour = 0 + + else + read(input_line(hourminute_firstcol:(hourminute_firstcol + 1)), *, iostat=rcio) hour + + if(rcio /= 0) then + if(debug) print *, "got bad read code getting hour at line ",line_number,", rcio = ",rcio + exit obsloop + + end if + end if + + read(input_line((hourminute_firstcol + 2):(hourminute_firstcol + 3)), *, iostat=rcio) minute + if(rcio /= 0) then + if(debug) print *, "got bad read code getting minute at line ",line_number,", rcio = ",rcio + exit obsloop + end if + + ! extracting the observation location + + read(input_line(vert_cols(1):vert_cols(2)), *, iostat=rcio) vert + if(rcio /= 0) then + if(debug) print *, "got bad read code getting vert at line ",line_number,", rcio = ",rcio + exit obsloop + + else if(abs(vert + 999.0) < 1.0d-8) then + cycle obsloop ! missing vertical coordinate + + else + vert = -vert ! MOM6 depth coordinate becomes negative as you increase depth + end if + + ! extracting the observation values + + otype_loop : do otype_index = 1, NUM_SCALAR_OBS + read(input_line(scalar_obs_cols(1, otype_index):scalar_obs_cols(2, otype_index)), *, iostat=rcio) ovalue + + if(rcio /= 0) then + if(debug) print *, "got bad read code getting observation type ",otype_index," at line ",line_number,", rcio = ",rcio + exit obsloop + + else if(abs(ovalue + 999.0) < 1.0d-8) then + cycle otype_loop ! missing observation + + end if + + ! 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 + end if + + time_obs = set_date(year, month, day, hours=hour, minutes=minute) + call get_time(time_obs, osec, days=oday) + + if(debug) then + call print_date(time_obs, "adding observation taken on") + print *, " \__ observation type: ",OTYPE_ORDERING(otype_index) + print *, " vert: ",vert + print *, " observation value: ",ovalue + end if + + num_processed(otype_index) = num_processed(otype_index) + 1 + running_sum(otype_index) = running_sum(otype_index) + ovalue + running_sqsum(otype_index) = running_sqsum(otype_index) + ovalue**2 + maxvals(otype_index) = max(maxvals(otype_index), ovalue) + minvals(otype_index) = min(minvals(otype_index), ovalue) + + ! determining which obs-sequence file to put this observation into + day_bin = oday + if(osec > 43200) day_bin = day_bin + 1 + new_obs_seq = (first_obs .or. (day_bin /= day_bin_old)) + day_bin_old = day_bin + + ! if necessary, saving the old observation sequence and beginning a new one. + if(new_obs_seq) then + ! dumping the observations so far into their own file + if ( (.not. first_obs) .and. (get_num_obs(obs_seq) > 0) ) then + if (debug) print *, 'writing obs_seq, obs_count = ', get_num_obs(obs_seq) + call write_obs_seq(obs_seq, obs_out_file) + call destroy_obs_sequence(obs_seq) + first_obs = .true. + endif + + write(daystr, "(I6)") day_bin + obs_out_file = trim(obs_out_dir)//"/BATS_"//daystr//".out" + + ! create a new, empty obs_seq file. + call init_obs_sequence(obs_seq, num_copies, num_qc, max_obs) + + ! the first one needs to contain the string 'observation' and the + ! second needs the string 'QC'. + call set_copy_meta_data(obs_seq, 1, 'observation') + call set_qc_meta_data(obs_seq, 1, 'Data QC') + end if + + obs_err = max(obs_uncertainties(otype_index)*ovalue, MIN_OBS_ERROR) + + call create_3d_obs(BATS_LAT, BATS_LON, vert, VERTISHEIGHT, & + ovalue, OTYPE_ORDERING(otype_index), obs_err, & + oday, osec, qc, obs) + + call add_obs_to_seq(obs_seq, obs, time_obs, prev_obs, prev_time, first_obs) + end do otype_loop +end do obsloop + +! putting any remaining observations into an obs sequence file + +if ( (.not. first_obs) .and. (get_num_obs(obs_seq) > 0) ) then + if (debug) print *, 'writing obs_seq, obs_count = ', get_num_obs(obs_seq) + call write_obs_seq(obs_seq, obs_out_file) + call destroy_obs_sequence(obs_seq) +endif + +print *, "" +print *, "SUMMARY:" +print *, "" + +do otype_index = 1, NUM_SCALAR_OBS + print *, "observation type: ",OTYPE_ORDERING(otype_index) + print *, "\__ total observations: ",num_processed(otype_index) + print *, " mean: ",running_sum(otype_index)/num_processed(otype_index) + print *, " variance: ",running_sqsum(otype_index)/num_processed(otype_index) & + - (running_sum(otype_index)/num_processed(otype_index))**2 + print *, " maximum value: ",maxvals(otype_index) + print *, " minimum value: ",minvals(otype_index) + print *, "" +end do + +! end of main program +call finalize_utilities() + +end program bats_to_obs diff --git a/observations/obs_converters/BATS/readme.rst b/observations/obs_converters/BATS/readme.rst new file mode 100644 index 0000000000..d72d34f9a3 --- /dev/null +++ b/observations/obs_converters/BATS/readme.rst @@ -0,0 +1,108 @@ +BATS +==== + +Overview: +--------- +BATS stands for the `Bermuda Atlantic Time-series Study `_. +BATS has collected data on the physical, biological, and chemical properties of +the ocean every month since 1988. BATS was established to uncover mysteries of the +deep ocean by analyzing important hydrographic and biological parameters +throughout the water column. + +Data Source: +------------ +Water column data from BATS (roughly 31N, 64W) can be downloaded from https://bats.bios.asu.edu/bats-data/ +The data is stored in in different forms. This converter operates on the ASCII formatted file, +often named ``bats_bottle.txt`` + +The data is huge extending from Oct 1988 to present day. It consists of a long list of information +such as +`Depth, Oxygen, CO2, Nitrate, Phosphate, Silicate, Alkalinity, Organic Carbon, Bacteria, ...` + +Observation Converter: +---------------------- +The obs converter is a program called ``bats_to_obs`` and has a namelist by the name ``&bats_to_obs_nml`` +Namelists start with an ampersand '&' and terminate with a slash '/'. + + .. code-block:: fortran + + &bats_to_obs_nml + text_input_file = "../bats_bottle.txt" + max_lines = 68000 + read_starting_at_line = 61 + date_firstcol = 14 + hourminute_firstcol = 35 + lat_cols = 42, 47 + lon_cols = 51, 56 + vert_cols = 64, 69 + scalar_obs_cols = 113, 119, + 137, 143, + 145, 151, + 153, 159, + 170, 176, + 178, 184 + obs_uncertainties = 0.2, + 0.2, + 0.2, + 0.2, + 0.2, + 0.2 + obs_out_dir = '../obs_seq_files', + debug = .true. + / + +This namelist provides control over the kind of observations to extract from the file in addition to their uncertainties. +In its current form, the observations that are extracted from the data file are: + +``BATS_OXYGEN``, ``BATS_INORGANIC_CARBON``, ``BATS_ALKALINITY``, ``BATS_NITRATE``, ``BATS_PHOSPHATE``, ``BATS_SILICATE`` + ++-------------------------------------+--------------------+------------------------------------------------------------+ +| Contents | Type | Description | ++=====================================+====================+============================================================+ +| ``text_input_file`` | character(len=256) | Pathname to the data file: ``bats_bottle.txt`` | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``max_lines`` | integer | Upper bound on the number of lines in the file that record | +| | | observations. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``read_starting_at_line`` | integer | Skip the information in the header of the file. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``date_firstcol`` | integer | First column of the YYYYMMDD date code at each line. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``hourminute_firstcol`` | integer | First column of the HHMM time stamp at each line. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``lat_cols`` | integer(2) | First and last columns where latitude is recorded. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``lon_cols`` | integer(2) | First and last columns where longitude is recorded. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``vert_cols`` | integer(2) | First and last columns where depth is recorded. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``scalar_obs_cols`` | integer(:, 2) | ith row of this table should list the first and last | +| | | columns where the value of the ith observation variable | +| | | is recorded. Ordering of observation variables is defined | +| | | by the OTYPE_ORDERING parameter in bats_to_obs.f90. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``obs_uncertainties`` | real(:) | ith entry of this list gives the uncertainty associated | +| | | with the ith observation variable. | +| | | | +| | | The observation error variance is defind as the square of | +| | | the product of the ``obs_uncertainties`` and the | +| | | observation value. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``obs_out_dir`` | character(len=256) | Pathname to obs_seq files resulting from the converter. | ++-------------------------------------+--------------------+------------------------------------------------------------+ +| ``debug`` | logical | A switch that makes the converetr prints useful | +| | | information as it runs. | ++-------------------------------------+--------------------+------------------------------------------------------------+ + +Climatology: +~~~~~~~~~~~~ +On top of assimilating real-time data, we often observe the quasi-cyclostationary behavior of the biogeochemical system +over the period of one year, and we update MARBL parameters by comparing this observed climatology to a climatology +predicted by MARBL. This usually involves running different forms of the ensmeble smoother where the model is re-run +using the updated parameters over long periods of times. + +To access the observed climatology at BATS, the script ``bats_climatology.py`` can be used to generate the climatology +by averaging the data over time. The program ``bats_to_clim_obs`` can then be executed to generate DART-style +observation sequence files using the climatological data. This code also supports +`Multiple Data Assimilation (MDA) `_ in +which the observations are assimilated multiple times with inflated observation error variance. diff --git a/observations/obs_converters/BATS/work/input.nml b/observations/obs_converters/BATS/work/input.nml new file mode 100644 index 0000000000..165b306376 --- /dev/null +++ b/observations/obs_converters/BATS/work/input.nml @@ -0,0 +1,76 @@ + +&bats_to_obs_nml + text_input_file = "../bats_bottle.txt" + max_lines = 68000 ! upper bound on the number of lines in the file that record observations + read_starting_at_line = 61 + date_firstcol = 14 ! first column of the YYYYMMDD date code at each line + hourminute_firstcol = 35 ! first column of the HHMM time stamp at each line + lat_cols = 42, 47 ! first and last columns where latitude, longitude, and vertical are recorded + lon_cols = 51, 56 + vert_cols = 64, 69 + scalar_obs_cols = 113, 119, ! i^th row of this table should list the first and last + 137, 143, ! columns where the value of the i^th observation variable + 145, 151, ! is recorded. Ordering of observation variables is defined + 153, 159, ! by the OTYPE_ORDERING parameter in bats_to_obs.f90. + 170, 176, + 178, 184 + obs_uncertainties = 0.2, ! i^th entry of this list gives the uncertainty associated + 0.2, ! with the i^th observation variable. Ordering of observation + 0.2, ! variables is defined by the OTYPE_ORDERING parameter in + 0.2, ! bats_to_obs.f90. + 0.2, ! + 0.2 + obs_out_dir = '../obs_seq_files', + debug = .true. + / + +&bats_to_clim_obs_nml + text_input_file = "../bats_climatology.txt" + max_lines = 3000 + obs_err_var_inflation = 10 + obs_out_dir = '../obs_seq_files', + debug = .true. + / + +&preprocess_nml + input_obs_def_mod_file = '../../../../observations/forward_operators/DEFAULT_obs_def_mod.F90' + output_obs_def_mod_file = '../../../../observations/forward_operators/obs_def_mod.f90' + input_obs_qty_mod_file = '../../../../assimilation_code/modules/observations/DEFAULT_obs_kind_mod.F90' + output_obs_qty_mod_file = '../../../../assimilation_code/modules/observations/obs_kind_mod.f90' + obs_type_files = '../../../../observations/forward_operators/obs_def_ocean_mod.f90' + quantity_files = '../../../../assimilation_code/modules/observations/default_quantities_mod.f90', + '../../../../assimilation_code/modules/observations/ocean_quantities_mod.f90' + / + +&obs_kind_nml + assimilate_these_obs_types = 'BATS_OXYGEN', + 'BATS_INORGANIC_CARBON', + 'BATS_ALKALINITY', + 'BATS_NITRATE', + 'BATS_PHOSPHATE', + 'BATS_SILICATE' + / + +&location_nml + / + +&utilities_nml + module_details = .false. + / + +&obs_sequence_nml + write_binary_obs_sequence = .false. + / + +&obs_sequence_tool_nml + filename_seq = 'obs_seq.out' + filename_seq_list = '' + filename_out = 'obs_seq.copy' + print_only = .false. + gregorian_cal = .true. + first_obs_days = -1 + first_obs_seconds = -1 + last_obs_days = -1 + last_obs_seconds = -1 + / + diff --git a/observations/obs_converters/BATS/work/quickbuild.sh b/observations/obs_converters/BATS/work/quickbuild.sh new file mode 100755 index 0000000000..ba99a14b5b --- /dev/null +++ b/observations/obs_converters/BATS/work/quickbuild.sh @@ -0,0 +1,40 @@ +#!/usr/bin/env bash + +# DART software - Copyright UCAR. This open source software is provided +# by UCAR, "as is", without charge, subject to all terms of use at +# http://www.image.ucar.edu/DAReS/DART/DART_download + +main() { + +export DART=$(git rev-parse --show-toplevel) +source "$DART"/build_templates/buildconvfunctions.sh + +CONVERTER=BATS +LOCATION=threed_sphere + +programs=( +bats_to_obs +bats_to_clim_obs +obs_sequence_tool +advance_time +) + +# build arguments +arguments "$@" + +# clean the directory +\rm -f -- *.o *.mod Makefile .cppdefs + +# build and run preprocess before making any other DART executables +buildpreprocess + +# build +buildconv + + +# clean up +\rm -f -- *.o *.mod + +} + +main "$@" diff --git a/observations/obs_converters/README.rst b/observations/obs_converters/README.rst index 0425020a93..151c326641 100644 --- a/observations/obs_converters/README.rst +++ b/observations/obs_converters/README.rst @@ -361,6 +361,7 @@ converters) include: - ``AURA``: See ``./AURA`` - ``Aviso+/CMEMS``: :doc:`./AVISO/AVISO` - ``Ameriflux``: :doc:`./Ameriflux/level4_to_obs` +- ``BATS``: :doc:`./BATS/readme` - ``CHAMP``: :doc:`./CHAMP/work/README` - ``cice``: :doc:`./cice/cice_to_obs` - ``CNOFS``: See ``./CNOFS``