diff --git a/.gitignore b/.gitignore index bac4fbf326..94d5bd9376 100644 --- a/.gitignore +++ b/.gitignore @@ -189,6 +189,7 @@ stacktest obs_rwtest test_quad_irreg_interp test_quad_reg_interp +test_table_read test_ran_unif # Directories to NOT IGNORE ... same as executable names diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 20689760e4..ff4686efe6 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,6 +22,20 @@ individual files. The changes are now listed with the most recent at the top. +**November 2 2023 :: QCEFF Input Table. Tag v11.1.0-alpha** + +- The QCEFF input table allows for the specification of QCEFF/probit + input options, per QTY, at runtime. +- This replaces the functionality of using an algorithm_info_mod specific + to the model, which meant editing algorithm_info_mod.f90 to specify + which distribution should be used for which quantity. +- The algorithm_info_mod files for the lorenz_96_tracer_advection model + examples have been replaced with set QCF tables (all_bnrhf_qcf_table.csv, + all_eakf_qcf_table.csv, state_eakf_tracer_bnrhf_qcf_table.csv, + neg_qcf_table.csv) and can be found in lorenz_96_tracer_advection/work. +- Removed the ‘global’ version of filter_kind from assim_tools_mod.f90 + and the &assim_tools_nml + **October 5 2023 :: WRF-DART tutorial diagnostic section. Tag v10.8.5** - Improvements: diff --git a/assimilation_code/modules/assimilation/algorithm_info_mod.f90 b/assimilation_code/modules/assimilation/algorithm_info_mod.f90 index 7ea474c664..aed33c0fc8 100644 --- a/assimilation_code/modules/assimilation/algorithm_info_mod.f90 +++ b/assimilation_code/modules/assimilation/algorithm_info_mod.f90 @@ -4,18 +4,21 @@ module algorithm_info_mod -use types_mod, only : r8, i8, missing_r8 +! Provides routines that give information about details of algorithms for +! observation error sampling, observation increments, and the transformations +! for regression and inflation in probit space. -use obs_def_mod, only : obs_def_type, get_obs_def_type_of_obs, get_obs_def_error_variance -use obs_kind_mod, only : get_quantity_for_type_of_obs +use types_mod, only : r8, i8, MISSING_R8, obstypelength -! Get the QTY definitions that are needed (aka kind) -use obs_kind_mod, only : QTY_STATE_VARIABLE, QTY_STATE_VAR_POWER, QTY_TRACER_CONCENTRATION, & - QTY_TRACER_SOURCE -! NOTE: Sadly, the QTY itself is not sufficient for the POWER because there is additional metadata +use obs_def_mod, only : obs_def_type, get_obs_def_type_of_obs, get_obs_def_error_variance +use obs_kind_mod, only : get_quantity_for_type_of_obs, get_name_for_quantity, get_index_for_quantity + +use utilities_mod, only : error_handler, E_ERR, E_MSG, open_file, close_file, to_upper, & + do_nml_file, do_nml_term, nmlfileunit, check_namelist_read, & + find_namelist_in_file use assim_model_mod, only : get_state_meta_data -use location_mod, only : location_type +use location_mod, only : location_type use distribution_params_mod, only : NORMAL_DISTRIBUTION, BOUNDED_NORMAL_RH_DISTRIBUTION, & GAMMA_DISTRIBUTION, BETA_DISTRIBUTION, LOG_NORMAL_DISTRIBUTION, UNIFORM_DISTRIBUTION, & @@ -24,30 +27,283 @@ module algorithm_info_mod implicit none private -! Defining parameter strings for different observation space filters -! For now, retaining backwards compatibility in assim_tools_mod requires using -! these specific integer values and there is no point in using these in assim_tools. -! That will change if backwards compatibility is removed in the future. +character(len=512) :: errstring +character(len=*), parameter :: source = 'algorithm_info_mod.f90' + +logical :: module_initialized = .false. +logical :: use_qty_defaults = .true. + +! Defining parameters for different observation space filters integer, parameter :: EAKF = 1 integer, parameter :: ENKF = 2 +integer, parameter :: KERNEL = 3 +integer, parameter :: OBS_PARTICLE = 4 integer, parameter :: UNBOUNDED_RHF = 8 integer, parameter :: GAMMA_FILTER = 11 integer, parameter :: BOUNDED_NORMAL_RHF = 101 public :: obs_error_info, probit_dist_info, obs_inc_info, & - EAKF, ENKF, BOUNDED_NORMAL_RHF, UNBOUNDED_RHF, GAMMA_FILTER - -! Provides routines that give information about details of algorithms for -! observation error sampling, observation increments, and the transformations -! for regression and inflation in probit space. -! For now, it is convenient to have these in a single module since several -! users will be developing their own problem specific versions of these -! subroutines. This will avoid constant merge conflicts as other parts of the -! assimilation code are updated. + init_algorithm_info_mod, end_algorithm_info_mod, & + EAKF, ENKF, BOUNDED_NORMAL_RHF, UNBOUNDED_RHF, & + GAMMA_FILTER, KERNEL, OBS_PARTICLE + +! type definitions for the QCF table +type obs_error_info_type + logical :: bounded_below, bounded_above + real(r8) :: lower_bound, upper_bound +end type + +type probit_inflation_type + integer :: dist_type + logical :: bounded_below, bounded_above + real(r8) :: lower_bound, upper_bound +end type + +type probit_state_type + integer :: dist_type + logical :: bounded_below, bounded_above + real(r8) :: lower_bound, upper_bound +end type + +type probit_extended_state_type + integer :: dist_type + logical :: bounded_below, bounded_above + real(r8) :: lower_bound, upper_bound +end type + +type obs_inc_info_type + integer :: filter_kind + logical :: bounded_below, bounded_above + real(r8) :: lower_bound, upper_bound +end type + +type algorithm_info_type + type(obs_error_info_type) :: obs_error_info + type(probit_inflation_type) :: probit_inflation + type(probit_state_type) :: probit_state + type(probit_extended_state_type) :: probit_extended_state + type(obs_inc_info_type) :: obs_inc_info +end type + +integer, parameter :: HEADER_LINES = 2 +character(len=129), dimension(4) :: header1 +character(len=129), dimension(25) :: header2 ! Number of table columns plus 1 + +integer, allocatable :: specified_qtys(:) +type(algorithm_info_type), allocatable :: qceff_table_data(:) + +character(len=129), allocatable :: dist_type_string_probit_inflation(:) +character(len=129), allocatable :: dist_type_string_probit_state(:) +character(len=129), allocatable :: dist_type_string_probit_extended_state(:) +character(len=129), allocatable :: filter_kind_string(:) + +! namelist +character(len=129) :: qceff_table_filename = '' + +namelist /algorithm_info_nml/ qceff_table_filename contains !------------------------------------------------------------------------- + + +subroutine init_algorithm_info_mod() + +! Gets number of lines/QTYs in the QCF table, allocates space for the table data + + +integer :: fileid +integer :: io, iunit + +integer :: numrows +integer :: nlines + +if (module_initialized) return +module_initialized = .true. + +! Read the namelist entry +call find_namelist_in_file("input.nml", "algorithm_info_nml", iunit) +read(iunit, nml = algorithm_info_nml, iostat = io) +call check_namelist_read(iunit, io, "algorithm_info_nml") + +if (do_nml_file()) write(nmlfileunit, nml=algorithm_info_nml) +if (do_nml_term()) write( * , nml=algorithm_info_nml) + + +if (qceff_table_filename == '') then + write(errstring,*) 'No QCF table file listed in namelist, using default values for all QTYs' + call error_handler(E_MSG, 'init_algorithm_info_mod:', errstring, source) + return +endif + +use_qty_defaults = .false. +fileid = open_file(trim(qceff_table_filename), 'formatted', 'read') + +! Do loop to get number of rows (or QTYs) in the table +nlines = 0 +do + read(fileid,*,iostat=io) + if(io/=0) exit + nlines = nlines + 1 +end do + +call close_file(fileid) + +numrows = nlines - HEADER_LINES + +allocate(specified_qtys(numrows)) +allocate(qceff_table_data(numrows)) +allocate(dist_type_string_probit_inflation(numrows)) +allocate(dist_type_string_probit_state(numrows)) +allocate(dist_type_string_probit_extended_state(numrows)) +allocate(filter_kind_string(numrows)) + +call read_qceff_table(qceff_table_filename) +call assert_qceff_table_version() +call verify_qceff_table_data() +call log_qceff_table_data() + +end subroutine init_algorithm_info_mod + +!------------------------------------------------------------------------ + + +subroutine read_qceff_table(qceff_table_filename) + +! Reads in the QCEFF input options from tabular data file + +character(len=129), intent(in) :: qceff_table_filename + +integer :: fileid +integer :: row +character(len=obstypelength) :: qty_string + +if (.not. module_initialized) call init_algorithm_info_mod() + +fileid = open_file(trim(qceff_table_filename), 'formatted', 'read') + +! skip the headers +read(fileid, *) header1 +read(fileid, *) header2 + +! read in table values directly to qceff_table_data type +do row = 1, size(qceff_table_data) + read(fileid, *) qty_string, qceff_table_data(row)%obs_error_info%bounded_below, qceff_table_data(row)%obs_error_info%bounded_above, & + qceff_table_data(row)%obs_error_info%lower_bound, qceff_table_data(row)%obs_error_info%upper_bound, dist_type_string_probit_inflation(row), & + qceff_table_data(row)%probit_inflation%bounded_below, qceff_table_data(row)%probit_inflation%bounded_above, & + qceff_table_data(row)%probit_inflation%lower_bound, qceff_table_data(row)%probit_inflation%upper_bound, dist_type_string_probit_state(row), & + qceff_table_data(row)%probit_state%bounded_below, qceff_table_data(row)%probit_state%bounded_above, & + qceff_table_data(row)%probit_state%lower_bound, qceff_table_data(row)%probit_state%upper_bound, dist_type_string_probit_extended_state(row), & + qceff_table_data(row)%probit_extended_state%bounded_below, qceff_table_data(row)%probit_extended_state%bounded_above, & + qceff_table_data(row)%probit_extended_state%lower_bound, qceff_table_data(row)%probit_extended_state%upper_bound, & + filter_kind_string(row), qceff_table_data(row)%obs_inc_info%bounded_below, qceff_table_data(row)%obs_inc_info%bounded_above, & + qceff_table_data(row)%obs_inc_info%lower_bound, qceff_table_data(row)%obs_inc_info%upper_bound + + call to_upper(qty_string) + specified_qtys(row) = get_index_for_quantity(qty_string) + + if(specified_qtys(row) == -1) then + write(errstring,*) trim(qty_string), ' is not a valid DART QTY' + call error_handler(E_ERR, 'read_qceff_table:', errstring, source) + endif + + ! Converting the distribution types (read in from table as a string) to its corresponding int value + call to_upper(dist_type_string_probit_inflation(row)) + + select case (trim(dist_type_string_probit_inflation(row))) + case ('NORMAL_DISTRIBUTION') + qceff_table_data(row)%probit_inflation%dist_type = NORMAL_DISTRIBUTION + case ('BOUNDED_NORMAL_RH_DISTRIBUTION') + qceff_table_data(row)%probit_inflation%dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION + case ('GAMMA_DISTRIBUTION') + qceff_table_data(row)%probit_inflation%dist_type = GAMMA_DISTRIBUTION + case ('BETA_DISTRIBUTION') + qceff_table_data(row)%probit_inflation%dist_type = BETA_DISTRIBUTION + case ('LOG_NORMAL_DISTRIBUTION') + qceff_table_data(row)%probit_inflation%dist_type = LOG_NORMAL_DISTRIBUTION + case ('UNIFORM_DISTRIBUTION') + qceff_table_data(row)%probit_inflation%dist_type = UNIFORM_DISTRIBUTION + case ('PARTICLE_FILTER_DISTRIBUTION') + qceff_table_data(row)%probit_inflation%dist_type = PARTICLE_FILTER_DISTRIBUTION + case default + write(errstring, *) 'Invalid distribution type for probit inflation: ', trim(dist_type_string_probit_inflation(row)) + call error_handler(E_ERR, 'read_qceff_table:', errstring, source) + end select + + + call to_upper(dist_type_string_probit_state(row)) + + select case (trim(dist_type_string_probit_state(row))) + case ('NORMAL_DISTRIBUTION') + qceff_table_data(row)%probit_state%dist_type = NORMAL_DISTRIBUTION + case ('BOUNDED_NORMAL_RH_DISTRIBUTION') + qceff_table_data(row)%probit_state%dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION + case ('GAMMA_DISTRIBUTION') + qceff_table_data(row)%probit_state%dist_type = GAMMA_DISTRIBUTION + case ('BETA_DISTRIBUTION') + qceff_table_data(row)%probit_state%dist_type = BETA_DISTRIBUTION + case ('LOG_NORMAL_DISTRIBUTION') + qceff_table_data(row)%probit_state%dist_type = LOG_NORMAL_DISTRIBUTION + case ('UNIFORM_DISTRIBUTION') + qceff_table_data(row)%probit_state%dist_type = UNIFORM_DISTRIBUTION + case ('PARTICLE_FILTER_DISTRIBUTION') + qceff_table_data(row)%probit_state%dist_type = PARTICLE_FILTER_DISTRIBUTION + case default + write(errstring, *) 'Invalid distribution type for probit state: ', trim(dist_type_string_probit_state(row)) + call error_handler(E_ERR, 'read_qceff_table:', errstring, source) + end select + + call to_upper(dist_type_string_probit_extended_state(row)) + + select case (trim(dist_type_string_probit_extended_state(row))) + case ('NORMAL_DISTRIBUTION') + qceff_table_data(row)%probit_extended_state%dist_type = NORMAL_DISTRIBUTION + case ('BOUNDED_NORMAL_RH_DISTRIBUTION') + qceff_table_data(row)%probit_extended_state%dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION + case ('GAMMA_DISTRIBUTION') + qceff_table_data(row)%probit_extended_state%dist_type = GAMMA_DISTRIBUTION + case ('BETA_DISTRIBUTION') + qceff_table_data(row)%probit_extended_state%dist_type = BETA_DISTRIBUTION + case ('LOG_NORMAL_DISTRIBUTION') + qceff_table_data(row)%probit_extended_state%dist_type = LOG_NORMAL_DISTRIBUTION + case ('UNIFORM_DISTRIBUTION') + qceff_table_data(row)%probit_extended_state%dist_type = UNIFORM_DISTRIBUTION + case ('PARTICLE_FILTER_DISTRIBUTION') + qceff_table_data(row)%probit_extended_state%dist_type = PARTICLE_FILTER_DISTRIBUTION + case default + write(errstring, *) 'Invalid distribution type for probit extended state: ', trim(dist_type_string_probit_extended_state(row)) + call error_handler(E_ERR, 'read_qceff_table:', errstring, source) + end select + + + ! Converting the filter kind (read in from table as a string) to its corresponding int value + call to_upper(filter_kind_string(row)) + + select case (trim(filter_kind_string(row))) + case ('EAKF') + qceff_table_data(row)%obs_inc_info%filter_kind = EAKF + case ('ENKF') + qceff_table_data(row)%obs_inc_info%filter_kind = ENKF + case ('UNBOUNDED_RHF') + qceff_table_data(row)%obs_inc_info%filter_kind = UNBOUNDED_RHF + case ('GAMMA_FILTER') + qceff_table_data(row)%obs_inc_info%filter_kind = GAMMA_FILTER + case ('BOUNDED_NORMAL_RHF') + qceff_table_data(row)%obs_inc_info%filter_kind = BOUNDED_NORMAL_RHF + case default + write(errstring, *) 'Invalid filter kind: ', trim(filter_kind_string(row)) + call error_handler(E_ERR, 'read_qceff_table:', errstring, source) + end select + +end do + +call close_file(fileid) + +end subroutine read_qceff_table + +!------------------------------------------------------------------------ + + subroutine obs_error_info(obs_def, error_variance, & bounded_below, bounded_above, lower_bound, upper_bound) @@ -58,69 +314,78 @@ subroutine obs_error_info(obs_def, error_variance, & logical, intent(out) :: bounded_below, bounded_above real(r8), intent(out) :: lower_bound, upper_bound -integer :: obs_type, obs_kind +integer :: obs_type, obs_qty integer(i8) :: state_var_index type(location_type) :: temp_loc -! Get the kind of the observation +integer :: QTY_loc(1) + +if (.not. module_initialized) call init_algorithm_info_mod() + +! Get the type of the observation obs_type = get_obs_def_type_of_obs(obs_def) ! If it is negative, it is an identity obs if(obs_type < 0) then state_var_index = -1 * obs_type - call get_state_meta_data(state_var_index, temp_loc, obs_kind) + call get_state_meta_data(state_var_index, temp_loc, obs_qty) else - obs_kind = get_quantity_for_type_of_obs(obs_type) + obs_qty = get_quantity_for_type_of_obs(obs_type) endif ! Get the default error variance error_variance = get_obs_def_error_variance(obs_def) -! Set the observation error details for each type of quantity -if(obs_kind == QTY_STATE_VARIABLE) then +!use default values if qceff_table_filename is not in namelist +if (use_qty_defaults) then bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_CONCENTRATION) then - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_SOURCE) then - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -else - write(*, *) 'Illegal obs_kind in obs_error_info' - stop + lower_bound = MISSING_R8; upper_bound = MISSING_R8 + return endif -end subroutine obs_error_info +!find location of QTY in qceff_table_data structure +QTY_loc = findloc(specified_qtys, obs_qty) + +if (QTY_loc(1) == 0) then + !use default values if QTY is not in table + bounded_below = .false.; bounded_above = .false. + lower_bound = MISSING_R8; upper_bound = MISSING_R8 + + else + bounded_below = qceff_table_data(QTY_loc(1))%obs_error_info%bounded_below + bounded_above = qceff_table_data(QTY_loc(1))%obs_error_info%bounded_above + lower_bound = qceff_table_data(QTY_loc(1))%obs_error_info%lower_bound + upper_bound = qceff_table_data(QTY_loc(1))%obs_error_info%upper_bound + +endif +end subroutine obs_error_info !------------------------------------------------------------------------- -subroutine probit_dist_info(kind, is_state, is_inflation, dist_type, & +subroutine probit_dist_info(qty, is_state, is_inflation, dist_type, & bounded_below, bounded_above, lower_bound, upper_bound) -! Computes the details of the probit transform for initial experiments -! with Molly - -integer, intent(in) :: kind +integer, intent(in) :: qty logical, intent(in) :: is_state ! True for state variable, false for obs logical, intent(in) :: is_inflation ! True for inflation transform integer, intent(out) :: dist_type logical, intent(out) :: bounded_below, bounded_above real(r8), intent(out) :: lower_bound, upper_bound -! Have input information about the kind of the state or observation being transformed +integer :: QTY_loc(1) + +! Have input information about the qty of the state or observation being transformed ! along with additional logical info that indicates whether this is an observation ! or state variable and about whether the transformation is being done for inflation -! or for regress. -! Need to select the appropriate transform. At present, options are NORMAL_PRIOR -! which does nothing or BOUNDED_NORMAL_RH_PRIOR. +! or for regress. + +! Selects the appropriate transform, which is specified in the QCF input table per QTY. +! At present, the options are NORMAL_DISTRIBUTION, BOUNDED_NORMAL_RH_DISTRIBUTION, +! GAMMA_DISTRIBUTION, BETA_DISTRIBUTION, LOG_NORMAL_DISTRIBUTION, +! UNIFORM_DISTRIBUTION, and PARTICLE_FILTER_DISTRIBUTION. ! If the BNRH is selected then information about the bounds must also be set. -! The two dimensional logical array 'bounded' is set to false for no bounds and true -! for bounded. the first element of the array is for the lower bound, the second for the upper. -! If bounded is chosen, the corresponding bound value(s) must be set in the two dimensional -! real array 'bounds'. -! For example, if my_state_kind corresponds to a sea ice fraction then an appropriate choice +! For example, if qty corresponds to a sea ice fraction then an appropriate choice ! would be: ! bounded_below = .true.; bounded_above = .true. ! lower_bound = 0.0_r8; upper_bounds = 1.0_r8 @@ -128,60 +393,51 @@ subroutine probit_dist_info(kind, is_state, is_inflation, dist_type, & ! In the long run, may not have to have separate controls for each of the input possibilities ! However, for now these are things that need to be explored for science understanding -if(is_inflation) then +if (.not. module_initialized) call init_algorithm_info_mod() + +!use default values if qceff_table_filename is not in namelist +if (use_qty_defaults) then + dist_type = NORMAL_DISTRIBUTION + bounded_below = .false.; bounded_above = .false. + lower_bound = MISSING_R8; upper_bound = MISSING_R8 + return +endif + +QTY_loc = findloc(specified_qtys, qty) + +if (QTY_loc(1) == 0) then + !use default values if QTY is not in table + dist_type = NORMAL_DISTRIBUTION + bounded_below = .false.; bounded_above = .false. + lower_bound = MISSING_R8; upper_bound = MISSING_R8 + + elseif(is_inflation) then ! Case for inflation transformation - if(kind == QTY_STATE_VARIABLE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -elseif(is_state) then + + dist_type = qceff_table_data(QTY_loc(1))%probit_inflation%dist_type + bounded_below = qceff_table_data(QTY_loc(1))%probit_inflation%bounded_below + bounded_above = qceff_table_data(QTY_loc(1))%probit_inflation%bounded_above + lower_bound = qceff_table_data(QTY_loc(1))%probit_inflation%lower_bound + upper_bound = qceff_table_data(QTY_loc(1))%probit_inflation%upper_bound + + elseif(is_state) then ! Case for state variable priors - if(kind == QTY_STATE_VARIABLE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 + + dist_type = qceff_table_data(QTY_loc(1))%probit_state%dist_type + bounded_below = qceff_table_data(QTY_loc(1))%probit_state%bounded_below + bounded_above = qceff_table_data(QTY_loc(1))%probit_state%bounded_above + lower_bound = qceff_table_data(QTY_loc(1))%probit_state%lower_bound + upper_bound = qceff_table_data(QTY_loc(1))%probit_state%upper_bound + else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -else ! This case is for observation (extended state) priors - if(kind == QTY_STATE_VARIABLE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif + + dist_type = qceff_table_data(QTY_loc(1))%probit_extended_state%dist_type + bounded_below = qceff_table_data(QTY_loc(1))%probit_extended_state%bounded_below + bounded_above = qceff_table_data(QTY_loc(1))%probit_extended_state%bounded_above + lower_bound = qceff_table_data(QTY_loc(1))%probit_extended_state%lower_bound + upper_bound = qceff_table_data(QTY_loc(1))%probit_extended_state%upper_bound + endif end subroutine probit_dist_info @@ -189,53 +445,176 @@ end subroutine probit_dist_info !------------------------------------------------------------------------ -subroutine obs_inc_info(obs_kind, filter_kind, rectangular_quadrature, gaussian_likelihood_tails, & - sort_obs_inc, spread_restoration, bounded_below, bounded_above, lower_bound, upper_bound) +subroutine obs_inc_info(obs_qty, filter_kind, bounded_below, bounded_above, & + lower_bound, upper_bound) -integer, intent(in) :: obs_kind -integer, intent(inout) :: filter_kind -logical, intent(inout) :: rectangular_quadrature, gaussian_likelihood_tails -logical, intent(inout) :: sort_obs_inc -logical, intent(inout) :: spread_restoration -logical, intent(inout) :: bounded_below, bounded_above -real(r8), intent(inout) :: lower_bound, upper_bound - -! The information arguments are all intent (inout). This means that if they are not set -! here, they retain the default values from the assim_tools_mod namelist. Bounds don't exist -! in that namelist, so default values are set in assim_tools_mod just before the call to here. +integer, intent(in) :: obs_qty +integer, intent(out) :: filter_kind +logical, intent(out) :: bounded_below, bounded_above +real(r8), intent(out) :: lower_bound, upper_bound -! Temporary approach for setting the details of how to assimilate this observation -! This example is designed to reproduce the squared forward operator results from paper +integer :: QTY_loc(1) +if (.not. module_initialized) call init_algorithm_info_mod() -! Set the observation increment details for each type of quantity -if(obs_kind == QTY_STATE_VARIABLE) then - filter_kind = BOUNDED_NORMAL_RHF +!use default values if qceff_table_filename is not in namelist +if (use_qty_defaults) then + filter_kind = EAKF bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_CONCENTRATION) then - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_SOURCE) then - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -else - write(*, *) 'Illegal obs_kind in obs_error_info' - stop + lower_bound = MISSING_R8; upper_bound = MISSING_R8 + return endif -! Default settings for now for Icepack and tracer model tests -sort_obs_inc = .false. -spread_restoration = .false. +!find location of QTY in qceff_table_data structure +QTY_loc = findloc(specified_qtys, obs_qty) -! Only need to set these two for options the original RHF implementation -!!!rectangular_quadrature = .true. -!!!gaussian_likelihood_tails = .false. +if (QTY_loc(1) == 0) then + !use default values if QTY is not in table + filter_kind = EAKF + bounded_below = .false.; bounded_above = .false. + lower_bound = MISSING_R8; upper_bound = MISSING_R8 + + else + + filter_kind = qceff_table_data(QTY_loc(1))%obs_inc_info%filter_kind + bounded_below = qceff_table_data(QTY_loc(1))%obs_inc_info%bounded_below + bounded_above = qceff_table_data(QTY_loc(1))%obs_inc_info%bounded_above + lower_bound = qceff_table_data(QTY_loc(1))%obs_inc_info%lower_bound + upper_bound = qceff_table_data(QTY_loc(1))%obs_inc_info%upper_bound + +endif end subroutine obs_inc_info !------------------------------------------------------------------------ + +subroutine assert_qceff_table_version() + +! Subroutine to ensure the correct version of the QCF table is being used + +character(1), parameter :: QCF_VERSION = '1' + +if (trim(header1(4)) /= QCF_VERSION) then + write(errstring,*) 'Using outdated/incorrect version of the QCF table' + call error_handler(E_ERR, 'assert_qceff_table_version:', errstring, source) +endif + +end subroutine assert_qceff_table_version + +!------------------------------------------------------------------------ + + +subroutine verify_qceff_table_data() + +! Subroutine to ensure that the data in the QCF table is valid + +integer :: row + +if (use_qty_defaults) return + +!Checks that all bounds are valid; currently checks that the lower bound in less than the upper +!Here we could add more specific checks if we have known limits on the bounds +do row = 1, size(qceff_table_data) + + if (qceff_table_data(row)%obs_error_info%bounded_below .and. qceff_table_data(row)%obs_error_info%bounded_above) then + if(qceff_table_data(row)%obs_error_info%lower_bound > qceff_table_data(row)%obs_error_info%upper_bound) then + write(errstring,*) 'Invalid bounds in obs_error_info' + call error_handler(E_ERR, 'verify_qceff_table_data:', errstring, source) + endif + endif + if (qceff_table_data(row)%probit_inflation%bounded_below .and. qceff_table_data(row)%probit_inflation%bounded_above) then + if(qceff_table_data(row)%probit_inflation%lower_bound > qceff_table_data(row)%probit_inflation%upper_bound) then + write(errstring,*) 'Invalid bounds in probit_inflation' + call error_handler(E_ERR, 'verify_qceff_table_data:', errstring, source) + endif + endif + if(qceff_table_data(row)%probit_state%bounded_below .and. qceff_table_data(row)%probit_state%bounded_above) then + if(qceff_table_data(row)%probit_state%lower_bound > qceff_table_data(row)%probit_state%upper_bound) then + write(errstring,*) 'Invalid bounds in probit_state' + call error_handler(E_ERR, 'verify_qceff_table_data:', errstring, source) + endif + endif + if(qceff_table_data(row)%probit_extended_state%bounded_below .and. qceff_table_data(row)%probit_extended_state%bounded_above) then + if(qceff_table_data(row)%probit_extended_state%lower_bound > qceff_table_data(row)%probit_extended_state%upper_bound) then + write(errstring,*) 'Invalid bounds in probit_extended_state' + call error_handler(E_ERR, 'verify_qceff_table_data:', errstring, source) + endif + endif + if(qceff_table_data(row)%obs_inc_info%bounded_below .and. qceff_table_data(row)%obs_inc_info%bounded_above) then + if(qceff_table_data(row)%obs_inc_info%lower_bound > qceff_table_data(row)%obs_inc_info%upper_bound) then + write(errstring,*) 'Invalid bounds in obs_inc_info' + call error_handler(E_ERR, 'verify_qceff_table_data:', errstring, source) + endif + endif +end do + + +!Ensures that there are no duplicate QTYs in the table +do row = 1, size(qceff_table_data) + if(count(specified_qtys==specified_qtys(row)) > 1) then + write(errstring,*) trim(get_name_for_quantity(specified_qtys(row))), ' has multiple entries in the table' + call error_handler(E_ERR, 'verify_qceff_table_data:', errstring, source) + endif +end do + +end subroutine verify_qceff_table_data + +!------------------------------------------------------------------------ + + +subroutine log_qceff_table_data() + +! Subroutine to write the data in QCF table to dart_log +character(len=2000) :: log_msg +integer :: row + +if (use_qty_defaults) return + +call error_handler(E_MSG, '', '', source) !Writing blank line to log +call error_handler(E_MSG, 'log_qceff_table_data:', 'Logging the data in the QCF Table', source) + +! Write the table headers to the dart_log and terminal +write(log_msg, '(A4, A6, A9, A)') header1(:) +call error_handler(E_MSG, 'log_qceff_table_data:', trim(log_msg), source) + +write(log_msg,'(3A14, 2A12, 3(A10, 2A14, 2A12), A12, 2A14, 2A12)') header2(:) +call error_handler(E_MSG, 'log_qceff_table_data:', trim(log_msg), source) + +! Write the table data to the dart_log and terminal +do row = 1, size(qceff_table_data) + write(log_msg, *) trim(get_name_for_quantity(specified_qtys(row))), qceff_table_data(row)%obs_error_info%bounded_below, qceff_table_data(row)%obs_error_info%bounded_above, & + qceff_table_data(row)%obs_error_info%lower_bound, qceff_table_data(row)%obs_error_info%upper_bound, trim(dist_type_string_probit_inflation(row)), & + qceff_table_data(row)%probit_inflation%bounded_below, qceff_table_data(row)%probit_inflation%bounded_above, & + qceff_table_data(row)%probit_inflation%lower_bound, qceff_table_data(row)%probit_inflation%upper_bound, trim(dist_type_string_probit_state(row)), & + qceff_table_data(row)%probit_state%bounded_below, qceff_table_data(row)%probit_state%bounded_above, & + qceff_table_data(row)%probit_state%lower_bound, qceff_table_data(row)%probit_state%upper_bound, trim(dist_type_string_probit_extended_state(row)), & + qceff_table_data(row)%probit_extended_state%bounded_below, qceff_table_data(row)%probit_extended_state%bounded_above, & + qceff_table_data(row)%probit_extended_state%lower_bound, qceff_table_data(row)%probit_extended_state%upper_bound, & + trim(filter_kind_string(row)), qceff_table_data(row)%obs_inc_info%bounded_below, qceff_table_data(row)%obs_inc_info%bounded_above, & + qceff_table_data(row)%obs_inc_info%lower_bound, qceff_table_data(row)%obs_inc_info%upper_bound +call error_handler(E_MSG, 'log_qceff_table_data:', trim(log_msg), source) +end do + +call error_handler(E_MSG, '', '', source) !Writing blank line to log + +end subroutine log_qceff_table_data + +!------------------------------------------------------------------------ + + +subroutine end_algorithm_info_mod() + +if (.not. module_initialized) return +module_initialized = .false. + +if (use_qty_defaults) return + +deallocate(specified_qtys) +deallocate(qceff_table_data) + +end subroutine end_algorithm_info_mod + +!---------------------------------------------------------------------- + end module algorithm_info_mod diff --git a/assimilation_code/modules/assimilation/algorithm_info_mod.nml b/assimilation_code/modules/assimilation/algorithm_info_mod.nml new file mode 100644 index 0000000000..730cc67423 --- /dev/null +++ b/assimilation_code/modules/assimilation/algorithm_info_mod.nml @@ -0,0 +1,3 @@ +&algorithm_info_mod + qceff_table_filename = '' +/ diff --git a/assimilation_code/modules/assimilation/all_eakf_algorithm_info_mod b/assimilation_code/modules/assimilation/all_eakf_algorithm_info_mod deleted file mode 100644 index 1b1f44d269..0000000000 --- a/assimilation_code/modules/assimilation/all_eakf_algorithm_info_mod +++ /dev/null @@ -1,241 +0,0 @@ -! 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 - -module algorithm_info_mod - -use types_mod, only : r8, i8, missing_r8 - -use obs_def_mod, only : obs_def_type, get_obs_def_type_of_obs, get_obs_def_error_variance -use obs_kind_mod, only : get_quantity_for_type_of_obs - -! Get the QTY definitions that are needed (aka kind) -use obs_kind_mod, only : QTY_STATE_VARIABLE, QTY_STATE_VAR_POWER, QTY_TRACER_CONCENTRATION, & - QTY_TRACER_SOURCE -! NOTE: Sadly, the QTY itself is not sufficient for the POWER because there is additional metadata - -use assim_model_mod, only : get_state_meta_data -use location_mod, only : location_type - -use distribution_params_mod, only : NORMAL_DISTRIBUTION, BOUNDED_NORMAL_RH_DISTRIBUTION, & - GAMMA_DISTRIBUTION, BETA_DISTRIBUTION, LOG_NORMAL_DISTRIBUTION, UNIFORM_DISTRIBUTION, & - PARTICLE_FILTER_DISTRIBUTION - -implicit none -private - -! Defining parameter strings for different observation space filters -! For now, retaining backwards compatibility in assim_tools_mod requires using -! these specific integer values and there is no point in using these in assim_tools. -! That will change if backwards compatibility is removed in the future. -integer, parameter :: EAKF = 1 -integer, parameter :: ENKF = 2 -integer, parameter :: UNBOUNDED_RHF = 8 -integer, parameter :: GAMMA_FILTER = 11 -integer, parameter :: BOUNDED_NORMAL_RHF = 101 - -public :: obs_error_info, probit_dist_info, obs_inc_info, & - EAKF, ENKF, BOUNDED_NORMAL_RHF, UNBOUNDED_RHF, GAMMA_FILTER - -! Provides routines that give information about details of algorithms for -! observation error sampling, observation increments, and the transformations -! for regression and inflation in probit space. -! For now, it is convenient to have these in a single module since several -! users will be developing their own problem specific versions of these -! subroutines. This will avoid constant merge conflicts as other parts of the -! assimilation code are updated. - -contains - -!------------------------------------------------------------------------- -subroutine obs_error_info(obs_def, error_variance, & - bounded_below, bounded_above, lower_bound, upper_bound) - -! Computes information needed to compute error sample for this observation -! This is called by perfect_model_obs when generating noisy obs -type(obs_def_type), intent(in) :: obs_def -real(r8), intent(out) :: error_variance -logical, intent(out) :: bounded_below, bounded_above -real(r8), intent(out) :: lower_bound, upper_bound - -integer :: obs_type, obs_kind -integer(i8) :: state_var_index -type(location_type) :: temp_loc - -! Get the kind of the observation -obs_type = get_obs_def_type_of_obs(obs_def) -! If it is negative, it is an identity obs -if(obs_type < 0) then - state_var_index = -1 * obs_type - call get_state_meta_data(state_var_index, temp_loc, obs_kind) -else - obs_kind = get_quantity_for_type_of_obs(obs_type) -endif - -! Get the default error variance -error_variance = get_obs_def_error_variance(obs_def) - -! Set the observation error details for each type of quantity -if(obs_kind == QTY_STATE_VARIABLE) then - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_CONCENTRATION) then - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_SOURCE) then - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -else - write(*, *) 'Illegal obs_kind in obs_error_info' - stop -endif - -end subroutine obs_error_info - - -!------------------------------------------------------------------------- - - -subroutine probit_dist_info(kind, is_state, is_inflation, dist_type, & - bounded_below, bounded_above, lower_bound, upper_bound) - -! Computes the details of the probit transform for initial experiments -! with Molly - -integer, intent(in) :: kind -logical, intent(in) :: is_state ! True for state variable, false for obs -logical, intent(in) :: is_inflation ! True for inflation transform -integer, intent(out) :: dist_type -logical, intent(out) :: bounded_below, bounded_above -real(r8), intent(out) :: lower_bound, upper_bound - -! Have input information about the kind of the state or observation being transformed -! along with additional logical info that indicates whether this is an observation -! or state variable and about whether the transformation is being done for inflation -! or for regress. -! Need to select the appropriate transform. At present, options are NORMAL_PRIOR -! which does nothing or BOUNDED_NORMAL_RH_PRIOR. -! If the BNRH is selected then information about the bounds must also be set. -! The two dimensional logical array 'bounded' is set to false for no bounds and true -! for bounded. the first element of the array is for the lower bound, the second for the upper. -! If bounded is chosen, the corresponding bound value(s) must be set in the two dimensional -! real array 'bounds'. -! For example, if my_state_kind corresponds to a sea ice fraction then an appropriate choice -! would be: -! bounded_below = .true.; bounded_above = .true. -! lower_bound = 0.0_r8; upper_bounds = 1.0_r8 - -! In the long run, may not have to have separate controls for each of the input possibilities -! However, for now these are things that need to be explored for science understanding - -if(is_inflation) then - ! Case for inflation transformation - if(kind == QTY_STATE_VARIABLE) then - dist_type = NORMAL_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = NORMAL_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type =NORMAL_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -elseif(is_state) then - ! Case for state variable priors - if(kind == QTY_STATE_VARIABLE) then - dist_type = NORMAL_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = NORMAL_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = NORMAL_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -else - ! This case is for observation (extended state) priors - if(kind == QTY_STATE_VARIABLE) then - dist_type = NORMAL_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = NORMAL_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = NORMAL_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -endif - -end subroutine probit_dist_info - -!------------------------------------------------------------------------ - - -subroutine obs_inc_info(obs_kind, filter_kind, rectangular_quadrature, gaussian_likelihood_tails, & - sort_obs_inc, spread_restoration, bounded_below, bounded_above, lower_bound, upper_bound) - -integer, intent(in) :: obs_kind -integer, intent(inout) :: filter_kind -logical, intent(inout) :: rectangular_quadrature, gaussian_likelihood_tails -logical, intent(inout) :: sort_obs_inc -logical, intent(inout) :: spread_restoration -logical, intent(inout) :: bounded_below, bounded_above -real(r8), intent(inout) :: lower_bound, upper_bound - -! The information arguments are all intent (inout). This means that if they are not set -! here, they retain the default values from the assim_tools_mod namelist. Bounds don't exist -! in that namelist, so default values are set in assim_tools_mod just before the call to here. - -! Temporary approach for setting the details of how to assimilate this observation -! This example is designed to reproduce the squared forward operator results from paper - - -! Set the observation increment details for each type of quantity -if(obs_kind == QTY_STATE_VARIABLE) then - filter_kind = EAKF - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_CONCENTRATION) then - filter_kind = EAKF - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_SOURCE) then - filter_kind = EAKF - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -else - write(*, *) 'Illegal obs_kind in obs_error_info' - stop -endif - -! Default settings for now for Icepack and tracer model tests -sort_obs_inc = .false. -spread_restoration = .false. - -! Only need to set these two for options the original RHF implementation -!!!rectangular_quadrature = .true. -!!!gaussian_likelihood_tails = .false. - -end subroutine obs_inc_info - -!------------------------------------------------------------------------ - -end module algorithm_info_mod diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index 72a01dc383..69da9d12dc 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -76,7 +76,9 @@ module assim_tools_mod use normal_distribution_mod, only : normal_cdf, inv_weighted_normal_cdf -use algorithm_info_mod, only : probit_dist_info, obs_inc_info +use algorithm_info_mod, only : probit_dist_info, obs_inc_info, EAKF, ENKF, & + BOUNDED_NORMAL_RHF, UNBOUNDED_RHF, GAMMA_FILTER, & + KERNEL, OBS_PARTICLE use gamma_distribution_mod, only : gamma_cdf, inv_gamma_cdf, gamma_mn_var_to_shape_scale, & gamma_gamma_prod @@ -131,22 +133,8 @@ module assim_tools_mod !---- namelist with default values -! Filter kind selects type of observation space filter -! 1 = EAKF filter -! 2 = ENKF -! 3 = Kernel filter -! 4 = particle filter -! 5 = random draw from posterior -! 6 = deterministic draw from posterior with fixed kurtosis -! 8 = Rank Histogram Filter (see Anderson 2011) -! -! special_localization_obs_types -> Special treatment for the specified observation types -! special_localization_cutoffs -> Different cutoff value for each specified obs type -! -logical :: use_algorithm_info_mod = .true. -integer :: filter_kind = 1 real(r8) :: cutoff = 0.2_r8 -logical :: sort_obs_inc = .false. +logical :: sort_obs_inc = .true. logical :: spread_restoration = .false. logical :: sampling_error_correction = .false. integer :: adaptive_localization_threshold = -1 @@ -161,7 +149,7 @@ module assim_tools_mod logical :: output_localization_diagnostics = .false. character(len = 129) :: localization_diagnostics_file = "localization_diagnostics" -! Following only relevant for filter_kind = 8 +! Following only relevant for filter_kind = UNBOUNDED_RHF logical :: rectangular_quadrature = .true. logical :: gaussian_likelihood_tails = .false. @@ -203,8 +191,7 @@ module assim_tools_mod ! compared to previous versions of this namelist item. logical :: distribute_mean = .false. -namelist / assim_tools_nml / use_algorithm_info_mod, & - filter_kind, cutoff, sort_obs_inc, & +namelist / assim_tools_nml / cutoff, sort_obs_inc, & spread_restoration, sampling_error_correction, & adaptive_localization_threshold, adaptive_cutoff_floor, & print_every_nth_obs, rectangular_quadrature, gaussian_likelihood_tails, & @@ -249,9 +236,9 @@ subroutine assim_tools_init() ! Note null_win_mod.f90 ignores distibute_mean. if (task_count() == 1) distribute_mean = .true. -! FOR NOW, can only do spread restoration with filter option 1 (need to extend this) -if(spread_restoration .and. .not. filter_kind == 1) then - write(msgstring, *) 'cannot combine spread_restoration and filter_kind ', filter_kind +if(spread_restoration) then + write(msgstring, *) 'The spread_restoration option is not supported in this version of ', & + 'DART. Contact the DAReS team if this option is needed ' call error_handler(E_ERR,'assim_tools_init:', msgstring, source) endif @@ -423,22 +410,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Need to give create_mean_window the mean copy call create_mean_window(ens_handle, ENS_MEAN_COPY, distribute_mean) -! filter kinds 1 and 8 return sorted increments, however non-deterministic -! inflation can scramble these. the sort is expensive, so help users get better -! performance by rejecting namelist combinations that do unneeded work. -if (sort_obs_inc) then - if(deterministic_inflate(inflate) .and. ((filter_kind == 1) .or. (filter_kind == 8))) then - write(msgstring, *) 'With a deterministic filter [assim_tools_nml:filter_kind = ',filter_kind,']' - write(msgstring2, *) 'and deterministic inflation [filter_nml:inf_deterministic = .TRUE.]' - write(msgstring3, *) 'assim_tools_nml:sort_obs_inc = .TRUE. is not needed and is expensive.' - call error_handler(E_MSG,'', '') ! whitespace - call error_handler(E_MSG,'WARNING filter_assim:', msgstring, source, & - text2=msgstring2,text3=msgstring3) - call error_handler(E_MSG,'', '') ! whitespace - sort_obs_inc = .FALSE. - endif -endif - ! Open the localization diagnostics file if(output_localization_diagnostics .and. my_task_id() == 0) & localization_unit = open_file(localization_diagnostics_file, action = 'append') @@ -938,11 +909,13 @@ subroutine obs_increment(ens_in, ens_size, obs, obs_var, obs_kind, obs_inc, & real(r8) :: rel_weights(ens_size) -! Declarations for bounded rank histogram filter -real(r8) :: likelihood(ens_size), like_sum +integer :: filter_kind logical :: bounded_below, bounded_above real(r8) :: lower_bound, upper_bound +! Declarations for bounded rank histogram filter +real(r8) :: likelihood(ens_size), like_sum + ! Copy the input ensemble to something that can be modified ens = ens_in @@ -973,30 +946,9 @@ subroutine obs_increment(ens_in, ens_size, obs, obs_var, obs_kind, obs_inc, & prior_var = sum((ens - prior_mean)**2) / (ens_size - 1) endif -!--------------------------begin algorithm_info control block----------------- -! More flexible abilities to control the observation space increments are -! available with this code block. It gets information about the increment method -! for the current observation is use_algorithm_info_mod is set to true in the namelist. -! This is not an extensible mechanism for doing this as the number of -! obs increments distributions and associated information goes up -! Implications for sorting increments and for spread restoration need to be examined -! further. -! Note that all but the first argument to obs_inc_info are intent(inout) so that if they -! are not set in that routine they will remain with the namelist selected values. - -! Set default values for bounds information -bounded_below = .false.; lower_bound = 0.0_r8 -bounded_above = .false.; upper_bound = 0.0_r8 - -if(use_algorithm_info_mod) & - call obs_inc_info(obs_kind, filter_kind, rectangular_quadrature, gaussian_likelihood_tails, & - sort_obs_inc, spread_restoration, bounded_below, bounded_above, lower_bound, upper_bound) - -! Could add logic to check on sort being true when not needed. -! Could also add logic to limit the use of spread_restoration to EAKF. It will fail -! in some ugly way right now. - -!----------------------------end algorithm_info control block----------------- +! Gets information about the increment method (filter_kind, bounds) for the current observation. +call obs_inc_info(obs_kind, filter_kind, bounded_below, bounded_above, & + lower_bound, upper_bound) ! The first three options in the next if block of code may be inappropriate for ! some more general filters; need to revisit @@ -1029,21 +981,21 @@ subroutine obs_increment(ens_in, ens_size, obs, obs_var, obs_kind, obs_inc, & ! note that at this point we've taken care of the cases where either the ! obs_var or the prior_var is 0, so the individual routines no longer need ! to have code to test for those cases. - if(filter_kind == 1) then + if(filter_kind == EAKF) then call obs_increment_eakf(ens, ens_size, prior_mean, prior_var, & obs, obs_var, obs_inc, net_a) - else if(filter_kind == 2) then + else if(filter_kind == ENKF) then call obs_increment_enkf(ens, ens_size, prior_var, obs, obs_var, obs_inc) - else if(filter_kind == 3) then + else if(filter_kind == KERNEL) then call obs_increment_kernel(ens, ens_size, obs, obs_var, obs_inc) - else if(filter_kind == 4) then + else if(filter_kind == OBS_PARTICLE) then call obs_increment_particle(ens, ens_size, obs, obs_var, obs_inc) - else if(filter_kind == 8) then + else if(filter_kind == UNBOUNDED_RHF) then call obs_increment_rank_histogram(ens, ens_size, prior_var, obs, obs_var, obs_inc) - else if(filter_kind == 11) then + else if(filter_kind == GAMMA_FILTER) then call obs_increment_gamma(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc) !-------------------------------------------------------------------------- - else if(filter_kind == 101) then + else if(filter_kind == BOUNDED_NORMAL_RHF) then ! Use bounded normal likelihood; Could use an arbitrary likelihood do i = 1, ens_size @@ -1077,28 +1029,13 @@ subroutine obs_increment(ens_in, ens_size, obs, obs_var, obs_kind, obs_inc, & !-------------------------------------------------------------------------- else call error_handler(E_ERR,'obs_increment', & - 'Illegal value of filter_kind in assim_tools namelist [1-8 OK]', source) + 'Illegal value of filter_kind', source) endif endif ! Add in the extra increments if doing observation space covariance inflation if(do_obs_inflate(inflate)) obs_inc = obs_inc + inflate_inc -! To minimize regression errors, may want to sort to minimize increments -! This makes sense for any of the non-deterministic algorithms -! By doing it here, can take care of both standard non-deterministic updates -! plus non-deterministic obs space covariance inflation. This is expensive, so -! don't use it if it's not needed. -if (sort_obs_inc) then - new_val = ens_in + obs_inc - ! Sorting to make increments as small as possible - call index_sort(ens_in, ens_index, ens_size) - call index_sort(new_val, new_index, ens_size) - do i = 1, ens_size - obs_inc(ens_index(i)) = new_val(new_index(i)) - ens_in(ens_index(i)) - end do -endif - ! Get the net change in spread if obs space inflation was used if(do_obs_inflate(inflate)) net_a = net_a * sqrt(my_cov_inflate) @@ -1359,7 +1296,8 @@ subroutine obs_increment_enkf(ens, ens_size, prior_var, obs, obs_var, obs_inc) real(r8) :: obs_var_inv, prior_var_inv, new_var, new_mean(ens_size) ! real(r8) :: sx, s_x2 real(r8) :: temp_mean, temp_obs(ens_size) -integer :: i +real(r8) :: new_val(ens_size) +integer :: i, ens_index(ens_size), new_index(ens_size) ! Compute mt_rinv_y (obs error normalized by variance) obs_var_inv = 1.0_r8 / obs_var @@ -1394,6 +1332,21 @@ subroutine obs_increment_enkf(ens, ens_size, prior_var, obs, obs_var, obs_inc) obs_inc(i) = new_mean(i) - ens(i) end do +! To minimize regression errors, may want to sort to minimize increments +! This makes sense for any of the non-deterministic algorithms +! By doing it here, can take care of both standard non-deterministic updates +! plus non-deterministic obs space covariance inflation. This is expensive, so +! don't use it if it's not needed. +if (sort_obs_inc) then + new_val = ens + obs_inc + ! Sorting to make increments as small as possible + call index_sort(ens, ens_index, ens_size) + call index_sort(new_val, new_index, ens_size) + do i = 1, ens_size + obs_inc(ens_index(i)) = new_val(new_index(i)) - ens(ens_index(i)) + end do +endif + ! Can also adjust mean (and) variance of final sample; works fine !sx = sum(new_mean) !s_x2 = sum(new_mean * new_mean) @@ -1570,41 +1523,6 @@ subroutine update_from_obs_inc(obs, obs_prior_mean, obs_prior_var, obs_inc, & ! Then compute the increment as product of reg_coef and observation space increment state_inc = reg_coef * obs_inc -! -! FIXME: craig schwartz has a degenerate case involving externally computed -! forward operators in which the obs prior variance is in fact exactly 0. -! adding this test allowed him to continue to use spread restoration -! without numerical problems. we don't know if this is sufficient; -! for now we'll leave the original code but it needs to be revisited. -! -! Spread restoration algorithm option. -!if(spread_restoration .and. obs_prior_var > 0.0_r8) then -! - -! Spread restoration algorithm option. -if(spread_restoration) then - ! Don't use this to reduce spread at present (should revisit this line) - net_a = min(net_a_in, 1.0_r8) - - ! Default restoration increment is 0.0 - restoration_inc = 0.0_r8 - - ! Compute the factor by which to inflate - ! These come from correl_error.f90 in system_simulation and the files ens??_pairs and - ! ens_pairs_0.5 in work under system_simulation. Assume a linear reduction from 1 - ! as a function of the net_a. Assume that the slope of this reduction is a function of - ! the reciprocal of the ensemble_size (slope = 0.80 / ens_size). These are empirical - ! for now. See also README in spread_restoration_paper documentation. - !!!factor = 1.0_r8 / (1.0_r8 + (net_a - 1.0_r8) * (0.8_r8 / ens_size)) - 1.0_r8 - factor = 1.0_r8 / (1.0_r8 + (net_a - 1.0_r8) / (-2.4711_r8 + 1.6386_r8 * ens_size)) - 1.0_r8 - !!!factor = 1.0_r8 / (1.0_r8 + (net_a**2 - 1.0_r8) * (-0.0111_r8 + .8585_r8 / ens_size)) - 1.0_r8 - - ! Variance restoration - state_mean = sum(state) / ens_size - restoration_inc = factor * (state - state_mean) - state_inc = state_inc + restoration_inc -endif - !! NOTE: if requested to be returned, correl_out is set further up in the !! code, before the sampling error correction, if enabled, is applied. !! this means it's returning a different larger value than the correl @@ -2506,33 +2424,6 @@ subroutine log_namelist_selections(num_special_cutoff, cache_override) integer :: i -select case (filter_kind) - case (1) - msgstring = 'Ensemble Adjustment Kalman Filter (EAKF)' - case (2) - msgstring = 'Ensemble Kalman Filter (ENKF)' - case (3) - msgstring = 'Kernel filter' - case (4) - msgstring = 'observation space particle filter' - case (5) - msgstring = 'random draw from posterior' - case (6) - msgstring = 'deterministic draw from posterior with fixed kurtosis' - case (7) - msgstring = 'Boxcar' - case (8) - msgstring = 'Rank Histogram Filter' - case (11) - msgstring = 'Gamma Filter' - case (101) - msgstring = 'Bounded Rank Histogram Filter' - case default - call error_handler(E_ERR, 'assim_tools_init:', 'illegal filter_kind value, valid values are 1-8', & - source) -end select -call error_handler(E_MSG, 'assim_tools_init:', 'Selected filter type is '//trim(msgstring)) - if (adjust_obs_impact) then call allocate_impact_table(obs_impact_table) call read_impact_table(obs_impact_filename, obs_impact_table, allow_any_impact_values, "allow_any_impact_values") diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.nml b/assimilation_code/modules/assimilation/assim_tools_mod.nml index d33aeead61..b0a0b76d71 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.nml +++ b/assimilation_code/modules/assimilation/assim_tools_mod.nml @@ -8,19 +8,15 @@ # enabling sampling error correction is generally beneficial # the default file is in assimilation_code/programs/gen_sampling_err_table/work -# With a deterministic filter (filter_kind == 1 or 8) -# and a deterministic inflation (filter_nml:inf_deterministic == .true.) -# sort_obs_inc is not needed and is expensive. Should be .false. +# sort_obs_inc applies to ENKF only. # specify special localization items in the same order # in both lists, the same number of items &assim_tools_nml - use_algorithm_info_mod = .true., - filter_kind = 1 cutoff = 0.2 distribute_mean = .false. - sort_obs_inc = .false. + sort_obs_inc = .true. spread_restoration = .false. sampling_error_correction = .false. adaptive_localization_threshold = -1 diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.rst b/assimilation_code/modules/assimilation/assim_tools_mod.rst index 88d0923164..9bbe135059 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.rst +++ b/assimilation_code/modules/assimilation/assim_tools_mod.rst @@ -10,25 +10,6 @@ for both mean and spread. In addition, algorithms to do a variety of flavors of particle filter, and kernel filters are included. The parallel implementation that allows each observation to update all state variables that are close to it at the same time is described in Anderson and Collins, 2007. -Filter types ------------- - -Available observation space filter types include: - -- 1 = EAKF (Ensemble Adjustment Kalman Filter, see Anderson 2001) -- 2 = ENKF (Ensemble Kalman Filter) -- 3 = Kernel filter -- 4 = Observation Space Particle filter -- 5 = Random draw from posterior (contact dart@ucar.edu before using) -- 6 = Deterministic draw from posterior with fixed kurtosis (ditto) -- 7 = Boxcar kernel filter -- 8 = Rank Histogram filter (see Anderson 2010) -- 9 = Particle filter (see Poterjoy 2016) - -We recommend using type=1, the EAKF. Note that although the algorithm is expressed in a slightly different form, the -EAKF is identical to the EnSRF (Ensemble Square Root Filter) described by Whitaker and Hamill in 2002. Highly -non-gaussian distributions may get better results from type=8, Rank Histogram filter. - Localization ------------ @@ -169,10 +150,9 @@ namelist. :: &assim_tools_nml - filter_kind = 1 cutoff = 0.2 distribute_mean = .false. - sort_obs_inc = .false. + sort_obs_inc = .true. spread_restoration = .false. sampling_error_correction = .false. adaptive_localization_threshold = -1 @@ -195,36 +175,6 @@ namelist. Description of each namelist entry ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -``filter_kind`` - *type:* integer - - Selects the variant of filter to be used. - - - 1 = EAKF (Ensemble Adjustment Kalman Filter, see Anderson 2001) - - 2 = ENKF (Ensemble Kalman Filter) - - 3 = Kernel filter - - 4 = Observation Space Particle filter - - 5 = Random draw from posterior (contact dart@ucar.edu before using) - - 6 = Deterministic draw from posterior with fixed kurtosis (ditto) - - 7 = Boxcar kernel filter - - 8 = Rank Histogram filter (see Anderson 2010) - - 9 = Particle filter (see Poterjoy 2016) - - The EAKF is the most commonly used filter. Note that although the algorithm is expressed in a slightly different - form, the EAKF is identical to the EnSRF (Ensemble Square Root Filter) described by Whitaker and Hamill in 2002. - - The Rank Histgram filter can be more successful for highly nongaussian distributions. - - Jon Poterjoy's Particle filter is included with this code release. To use, it, overwrite ``assim_tools_mod.f90`` with - ``assim_tools_mod.pf.f90`` and rebuild filter. - - :: - - - $ mv assimilation_code/modules/assimilation/assim_tools_mod.pf.f90 assimilation_code/modules/assimilation/assim_tools_mod.f90 - - There are additional namelist items in this version specific to the particle filter. Read the code for more details. - ``cutoff`` *type:* real(r8) @@ -245,10 +195,11 @@ Description of each namelist entry *type:* logical If true, the final increments from obs_increment are sorted so that the mean increment value is as small as possible. - This minimizes regression errors when non-deterministic filters or error correction algorithms are applied. HOWEVER, - when using deterministic filters (filter_kind == 1 or 8) with no inflation or a combination of a determinstic filter + Applies to ENKF only. + ``sort_ob_inc`` minimizes regression errors when non-deterministic filters or error correction algorithms are applied. HOWEVER, + when using deterministic filters with no inflation or a combination of a determinstic filter and deterministic inflation (filter_nml:inf_deterministic = .TRUE.) sorting the increments is both unnecessary and - expensive. A warning is printed to stdout and the log and the sorting is skipped. + expensive. ``spread_restoration`` *type:* logical @@ -256,6 +207,11 @@ Description of each namelist entry True turns on algorithm to restore amount of spread that would be expected to be lost if underlying obs/state variable correlation were really 0. +.. Warning:: + + ``spread_restoration`` is not supported in this version, please reach out to the DAReS team dart@ucar.edu + if you need to use spread_restoration. + ``sampling_error_correction`` *type:* logical @@ -307,12 +263,12 @@ Description of each namelist entry ``rectangular_quadrature`` *type:* logical - Only relevant for filter type 8 and recommended to leave ``.true.``. + Only relevant for filter type UNBOUNDED_RHF and recommended to leave ``.true.``. ``gaussian_likelihood_tails`` *type:* logical - Only relevant for filter type 8 and recommended to leave ``.false.``. + Only relevant for filter type UNBOUNDED_RHF and recommended to leave ``.false.``. ``close_obs_caching`` *type:* logical diff --git a/assimilation_code/modules/assimilation/filter_mod.f90 b/assimilation_code/modules/assimilation/filter_mod.f90 index 11e50c3038..d945abde24 100644 --- a/assimilation_code/modules/assimilation/filter_mod.f90 +++ b/assimilation_code/modules/assimilation/filter_mod.f90 @@ -90,9 +90,9 @@ module filter_mod use probit_transform_mod, only : transform_to_probit, transform_from_probit -use algorithm_info_mod, only : probit_dist_info +use algorithm_info_mod, only : probit_dist_info, init_algorithm_info_mod, end_algorithm_info_mod -use distribution_params_mod, only : distribution_params_type, NORMAL_DISTRIBUTION +use distribution_params_mod, only : distribution_params_type !------------------------------------------------------------------------------ @@ -166,7 +166,6 @@ module filter_mod !---------------------------------------------------------------- ! Namelist input with default values ! -logical :: use_algorithm_info_mod = .true. integer :: async = 0, ens_size = 20 integer :: tasks_per_model_advance = 1 ! if init_time_days and seconds are negative initial time is 0, 0 @@ -261,7 +260,6 @@ module filter_mod namelist /filter_nml/ async, & - use_algorithm_info_mod, & adv_ens_command, & ens_size, & tasks_per_model_advance, & @@ -361,13 +359,13 @@ subroutine filter_main() real(r8), allocatable :: prior_qc_copy(:) -call filter_initialize_modules_used() ! static_init_model called in here - ! Read the namelist entry call find_namelist_in_file("input.nml", "filter_nml", iunit) read(iunit, nml = filter_nml, iostat = io) call check_namelist_read(iunit, io, "filter_nml") +call filter_initialize_modules_used() ! static_init_model called in here + ! Record the namelist values used for the run ... if (do_nml_file()) write(nmlfileunit, nml=filter_nml) if (do_nml_term()) write( * , nml=filter_nml) @@ -1150,6 +1148,9 @@ subroutine filter_main() call end_assim_model() call trace_message('After end_model call') +! deallocate qceff_table_data structures +call end_algorithm_info_mod() + call trace_message('Before ensemble and obs memory cleanup') call end_ensemble_manager(state_ens_handle) @@ -1272,6 +1273,10 @@ subroutine filter_initialize_modules_used() call static_init_assim_model() call state_vector_io_init() call initialize_qc() + +! Initialize algorothm_info_mod and read in QCF table data +call init_algorithm_info_mod() + call trace_message('After filter_initialize_module_used call') end subroutine filter_initialize_modules_used @@ -1599,16 +1604,9 @@ subroutine filter_ensemble_inflate(ens_handle, inflate_copy, inflate, ENS_MEAN_C call get_state_meta_data(ens_handle%my_vars(j), my_state_loc, my_state_kind) ! Need to specify what kind of prior to use for each - ! Use default of untransformed if use_algorithm_info_mod is not true - if(use_algorithm_info_mod) then - call probit_dist_info(my_state_kind, .true., .true., dist_type, & - bounded_below, bounded_above, lower_bound, upper_bound) - else - ! Default is just a normal which does nothing - dist_type = NORMAL_DISTRIBUTION - bounded_below = .false. ; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = 0.0_r8 - endif + call probit_dist_info(my_state_kind, .true., .true., dist_type, & + bounded_below, bounded_above, lower_bound, upper_bound) + call transform_to_probit(grp_size, ens_handle%copies(grp_bot:grp_top, j), & dist_type, dist_params, probit_ens(1:grp_size), .false., & bounded_below, bounded_above, lower_bound, upper_bound) diff --git a/assimilation_code/modules/assimilation/filter_mod.nml b/assimilation_code/modules/assimilation/filter_mod.nml index 362138fc5f..0e3913be42 100644 --- a/assimilation_code/modules/assimilation/filter_mod.nml +++ b/assimilation_code/modules/assimilation/filter_mod.nml @@ -1,5 +1,4 @@ &filter_nml - use_algorithm_info_mod = .true., single_file_in = .false., input_state_files = '' input_state_file_list = '' diff --git a/assimilation_code/modules/assimilation/neg_algorithm_info_mod b/assimilation_code/modules/assimilation/neg_algorithm_info_mod deleted file mode 100644 index b02a2290e3..0000000000 --- a/assimilation_code/modules/assimilation/neg_algorithm_info_mod +++ /dev/null @@ -1,241 +0,0 @@ -! 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 - -module algorithm_info_mod - -use types_mod, only : r8, i8, missing_r8 - -use obs_def_mod, only : obs_def_type, get_obs_def_type_of_obs, get_obs_def_error_variance -use obs_kind_mod, only : get_quantity_for_type_of_obs - -! Get the QTY definitions that are needed (aka kind) -use obs_kind_mod, only : QTY_STATE_VARIABLE, QTY_STATE_VAR_POWER, QTY_TRACER_CONCENTRATION, & - QTY_TRACER_SOURCE -! NOTE: Sadly, the QTY itself is not sufficient for the POWER because there is additional metadata - -use assim_model_mod, only : get_state_meta_data -use location_mod, only : location_type - -use distribution_params_mod, only : NORMAL_DISTRIBUTION, BOUNDED_NORMAL_RH_DISTRIBUTION, & - GAMMA_DISTRIBUTION, BETA_DISTRIBUTION, LOG_NORMAL_DISTRIBUTION, UNIFORM_DISTRIBUTION, & - PARTICLE_FILTER_DISTRIBUTION - -implicit none -private - -! Defining parameter strings for different observation space filters -! For now, retaining backwards compatibility in assim_tools_mod requires using -! these specific integer values and there is no point in using these in assim_tools. -! That will change if backwards compatibility is removed in the future. -integer, parameter :: EAKF = 1 -integer, parameter :: ENKF = 2 -integer, parameter :: UNBOUNDED_RHF = 8 -integer, parameter :: GAMMA_FILTER = 11 -integer, parameter :: BOUNDED_NORMAL_RHF = 101 - -public :: obs_error_info, probit_dist_info, obs_inc_info, & - EAKF, ENKF, BOUNDED_NORMAL_RHF, UNBOUNDED_RHF, GAMMA_FILTER - -! Provides routines that give information about details of algorithms for -! observation error sampling, observation increments, and the transformations -! for regression and inflation in probit space. -! For now, it is convenient to have these in a single module since several -! users will be developing their own problem specific versions of these -! subroutines. This will avoid constant merge conflicts as other parts of the -! assimilation code are updated. - -contains - -!------------------------------------------------------------------------- -subroutine obs_error_info(obs_def, error_variance, & - bounded_below, bounded_above, lower_bound, upper_bound) - -! Computes information needed to compute error sample for this observation -! This is called by perfect_model_obs when generating noisy obs -type(obs_def_type), intent(in) :: obs_def -real(r8), intent(out) :: error_variance -logical, intent(out) :: bounded_below, bounded_above -real(r8), intent(out) :: lower_bound, upper_bound - -integer :: obs_type, obs_kind -integer(i8) :: state_var_index -type(location_type) :: temp_loc - -! Get the kind of the observation -obs_type = get_obs_def_type_of_obs(obs_def) -! If it is negative, it is an identity obs -if(obs_type < 0) then - state_var_index = -1 * obs_type - call get_state_meta_data(state_var_index, temp_loc, obs_kind) -else - obs_kind = get_quantity_for_type_of_obs(obs_type) -endif - -! Get the default error variance -error_variance = get_obs_def_error_variance(obs_def) - -! Set the observation error details for each type of quantity -if(obs_kind == QTY_STATE_VARIABLE) then - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_CONCENTRATION) then - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 -elseif(obs_kind == QTY_TRACER_SOURCE) then - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 -else - write(*, *) 'Illegal obs_kind in obs_error_info' - stop -endif - -end subroutine obs_error_info - - -!------------------------------------------------------------------------- - - -subroutine probit_dist_info(kind, is_state, is_inflation, dist_type, & - bounded_below, bounded_above, lower_bound, upper_bound) - -! Computes the details of the probit transform for initial experiments -! with Molly - -integer, intent(in) :: kind -logical, intent(in) :: is_state ! True for state variable, false for obs -logical, intent(in) :: is_inflation ! True for inflation transform -integer, intent(out) :: dist_type -logical, intent(out) :: bounded_below, bounded_above -real(r8), intent(out) :: lower_bound, upper_bound - -! Have input information about the kind of the state or observation being transformed -! along with additional logical info that indicates whether this is an observation -! or state variable and about whether the transformation is being done for inflation -! or for regress. -! Need to select the appropriate transform. At present, options are NORMAL_DISTRIBUTION -! which does nothing or BOUNDED_NORMAL_RH_DISTRIBUTION. -! If the BNRH is selected then information about the bounds must also be set. -! The two dimensional logical array 'bounded' is set to false for no bounds and true -! for bounded. the first element of the array is for the lower bound, the second for the upper. -! If bounded is chosen, the corresponding bound value(s) must be set in the two dimensional -! real array 'bounds'. -! For example, if my_state_kind corresponds to a sea ice fraction then an appropriate choice -! would be: -! bounded_below = .true.; bounded_above = .true. -! lower_bound = 0.0_r8; upper_bounds = 1.0_r8 - -! In the long run, may not have to have separate controls for each of the input possibilities -! However, for now these are things that need to be explored for science understanding - -if(is_inflation) then - ! Case for inflation transformation - if(kind == QTY_STATE_VARIABLE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -elseif(is_state) then - ! Case for state variable priors - if(kind == QTY_STATE_VARIABLE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -else - ! This case is for observation (extended state) priors - if(kind == QTY_STATE_VARIABLE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -endif - -end subroutine probit_dist_info - -!------------------------------------------------------------------------ - - -subroutine obs_inc_info(obs_kind, filter_kind, rectangular_quadrature, gaussian_likelihood_tails, & - sort_obs_inc, spread_restoration, bounded_below, bounded_above, lower_bound, upper_bound) - -integer, intent(in) :: obs_kind -integer, intent(inout) :: filter_kind -logical, intent(inout) :: rectangular_quadrature, gaussian_likelihood_tails -logical, intent(inout) :: sort_obs_inc -logical, intent(inout) :: spread_restoration -logical, intent(inout) :: bounded_below, bounded_above -real(r8), intent(inout) :: lower_bound, upper_bound - -! The information arguments are all intent (inout). This means that if they are not set -! here, they retain the default values from the assim_tools_mod namelist. Bounds don't exist -! in that namelist, so default values are set in assim_tools_mod just before the call to here. - -! Temporary approach for setting the details of how to assimilate this observation -! This example is designed to reproduce the squared forward operator results from paper - - -! Set the observation increment details for each type of quantity -if(obs_kind == QTY_STATE_VARIABLE) then - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_CONCENTRATION) then - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 -elseif(obs_kind == QTY_TRACER_SOURCE) then - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 -else - write(*, *) 'Illegal obs_kind in obs_error_info' - stop -endif - -! Default settings for now for Icepack and tracer model tests -sort_obs_inc = .false. -spread_restoration = .false. - -! Only need to set these two for options the original RHF implementation -!!!rectangular_quadrature = .true. -!!!gaussian_likelihood_tails = .false. - -end subroutine obs_inc_info - -!------------------------------------------------------------------------ - -end module algorithm_info_mod diff --git a/assimilation_code/modules/assimilation/one_above_algorithm_info_mod b/assimilation_code/modules/assimilation/one_above_algorithm_info_mod deleted file mode 100644 index f7e404bb3d..0000000000 --- a/assimilation_code/modules/assimilation/one_above_algorithm_info_mod +++ /dev/null @@ -1,241 +0,0 @@ -! 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 - -module algorithm_info_mod - -use types_mod, only : r8, i8, missing_r8 - -use obs_def_mod, only : obs_def_type, get_obs_def_type_of_obs, get_obs_def_error_variance -use obs_kind_mod, only : get_quantity_for_type_of_obs - -! Get the QTY definitions that are needed (aka kind) -use obs_kind_mod, only : QTY_STATE_VARIABLE, QTY_STATE_VAR_POWER, QTY_TRACER_CONCENTRATION, & - QTY_TRACER_SOURCE -! NOTE: Sadly, the QTY itself is not sufficient for the POWER because there is additional metadata - -use assim_model_mod, only : get_state_meta_data -use location_mod, only : location_type - -use distribution_params_mod, only : NORMAL_DISTRIBUTION, BOUNDED_NORMAL_RH_DISTRIBUTION, & - GAMMA_DISTRIBUTION, BETA_DISTRIBUTION, LOG_NORMAL_DISTRIBUTION, UNIFORM_DISTRIBUTION, & - PARTICLE_FILTER_DISTRIBUTION - -implicit none -private - -! Defining parameter strings for different observation space filters -! For now, retaining backwards compatibility in assim_tools_mod requires using -! these specific integer values and there is no point in using these in assim_tools. -! That will change if backwards compatibility is removed in the future. -integer, parameter :: EAKF = 1 -integer, parameter :: ENKF = 2 -integer, parameter :: UNBOUNDED_RHF = 8 -integer, parameter :: GAMMA_FILTER = 11 -integer, parameter :: BOUNDED_NORMAL_RHF = 101 - -public :: obs_error_info, probit_dist_info, obs_inc_info, & - EAKF, ENKF, BOUNDED_NORMAL_RHF, UNBOUNDED_RHF, GAMMA_FILTER - -! Provides routines that give information about details of algorithms for -! observation error sampling, observation increments, and the transformations -! for regression and inflation in probit space. -! For now, it is convenient to have these in a single module since several -! users will be developing their own problem specific versions of these -! subroutines. This will avoid constant merge conflicts as other parts of the -! assimilation code are updated. - -contains - -!------------------------------------------------------------------------- -subroutine obs_error_info(obs_def, error_variance, & - bounded_below, bounded_above, lower_bound, upper_bound) - -! Computes information needed to compute error sample for this observation -! This is called by perfect_model_obs when generating noisy obs -type(obs_def_type), intent(in) :: obs_def -real(r8), intent(out) :: error_variance -logical, intent(out) :: bounded_below, bounded_above -real(r8), intent(out) :: lower_bound, upper_bound - -integer :: obs_type, obs_kind -integer(i8) :: state_var_index -type(location_type) :: temp_loc - -! Get the kind of the observation -obs_type = get_obs_def_type_of_obs(obs_def) -! If it is negative, it is an identity obs -if(obs_type < 0) then - state_var_index = -1 * obs_type - call get_state_meta_data(state_var_index, temp_loc, obs_kind) -else - obs_kind = get_quantity_for_type_of_obs(obs_type) -endif - -! Get the default error variance -error_variance = get_obs_def_error_variance(obs_def) - -! Set the observation error details for each type of quantity -if(obs_kind == QTY_STATE_VARIABLE) then - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_CONCENTRATION) then - bounded_below = .true.; bounded_above = .true. - lower_bound = -10.0_r8; upper_bound = 1.0_r8 -elseif(obs_kind == QTY_TRACER_SOURCE) then - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 -else - write(*, *) 'Illegal obs_kind in obs_error_info' - stop -endif - -end subroutine obs_error_info - - -!------------------------------------------------------------------------- - - -subroutine probit_dist_info(kind, is_state, is_inflation, dist_type, & - bounded_below, bounded_above, lower_bound, upper_bound) - -! Computes the details of the probit transform for initial experiments -! with Molly - -integer, intent(in) :: kind -logical, intent(in) :: is_state ! True for state variable, false for obs -logical, intent(in) :: is_inflation ! True for inflation transform -integer, intent(out) :: dist_type -logical, intent(out) :: bounded_below, bounded_above -real(r8), intent(out) :: lower_bound, upper_bound - -! Have input information about the kind of the state or observation being transformed -! along with additional logical info that indicates whether this is an observation -! or state variable and about whether the transformation is being done for inflation -! or for regress. -! Need to select the appropriate transform. At present, options are NORMAL_DISTRIBUTION -! which does nothing or BOUNDED_NORMAL_RH_DISTRIBUTION. -! If the BNRH is selected then information about the bounds must also be set. -! The two dimensional logical array 'bounded' is set to false for no bounds and true -! for bounded. the first element of the array is for the lower bound, the second for the upper. -! If bounded is chosen, the corresponding bound value(s) must be set in the two dimensional -! real array 'bounds'. -! For example, if my_state_kind corresponds to a sea ice fraction then an appropriate choice -! would be: -! bounded_below = .true.; bounded_above = .true. -! lower_bound = 0.0_r8; upper_bounds = 1.0_r8 - -! In the long run, may not have to have separate controls for each of the input possibilities -! However, for now these are things that need to be explored for science understanding - -if(is_inflation) then - ! Case for inflation transformation - if(kind == QTY_STATE_VARIABLE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .true. - lower_bound = -10.0_r8; upper_bound = 1.0_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -elseif(is_state) then - ! Case for state variable priors - if(kind == QTY_STATE_VARIABLE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .true. - lower_bound = -10.0_r8; upper_bound = 1.0_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -else - ! This case is for observation (extended state) priors - if(kind == QTY_STATE_VARIABLE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .true. - lower_bound = -10.0_r8; upper_bound = 1.0_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -endif - -end subroutine probit_dist_info - -!------------------------------------------------------------------------ - - -subroutine obs_inc_info(obs_kind, filter_kind, rectangular_quadrature, gaussian_likelihood_tails, & - sort_obs_inc, spread_restoration, bounded_below, bounded_above, lower_bound, upper_bound) - -integer, intent(in) :: obs_kind -integer, intent(inout) :: filter_kind -logical, intent(inout) :: rectangular_quadrature, gaussian_likelihood_tails -logical, intent(inout) :: sort_obs_inc -logical, intent(inout) :: spread_restoration -logical, intent(inout) :: bounded_below, bounded_above -real(r8), intent(inout) :: lower_bound, upper_bound - -! The information arguments are all intent (inout). This means that if they are not set -! here, they retain the default values from the assim_tools_mod namelist. Bounds don't exist -! in that namelist, so default values are set in assim_tools_mod just before the call to here. - -! Temporary approach for setting the details of how to assimilate this observation -! This example is designed to reproduce the squared forward operator results from paper - - -! Set the observation increment details for each type of quantity -if(obs_kind == QTY_STATE_VARIABLE) then - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_CONCENTRATION) then - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .true.; bounded_above = .true. - lower_bound = -10.0_r8; upper_bound = 1.0_r8 -elseif(obs_kind == QTY_TRACER_SOURCE) then - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .false.; bounded_above = .true. - lower_bound = missing_r8; upper_bound = 0.0_r8 -else - write(*, *) 'Illegal obs_kind in obs_error_info' - stop -endif - -! Default settings for now for Icepack and tracer model tests -sort_obs_inc = .false. -spread_restoration = .false. - -! Only need to set these two for options the original RHF implementation -!!!rectangular_quadrature = .true. -!!!gaussian_likelihood_tails = .false. - -end subroutine obs_inc_info - -!------------------------------------------------------------------------ - -end module algorithm_info_mod diff --git a/assimilation_code/modules/assimilation/probit_transform_mod.f90 b/assimilation_code/modules/assimilation/probit_transform_mod.f90 index e388f2d06c..f0ea534a27 100644 --- a/assimilation_code/modules/assimilation/probit_transform_mod.f90 +++ b/assimilation_code/modules/assimilation/probit_transform_mod.f90 @@ -47,7 +47,7 @@ module probit_transform_mod ! Logical to fix bounds violations for bounded_normal_rh logical :: fix_bound_violations = .false. ! Should we use a logit transform instead of the default probit transform -logical :: use_logit_instead_of_probit = .true. +logical :: use_logit_instead_of_probit = .false. ! Set to true to do a check of the probit to/from transforms for inverse accuracy logical :: do_inverse_check = .false. diff --git a/assimilation_code/modules/assimilation/state_eakf_tracer_bnrhf_algorithm_info_mod b/assimilation_code/modules/assimilation/state_eakf_tracer_bnrhf_algorithm_info_mod deleted file mode 100644 index f01b4e5952..0000000000 --- a/assimilation_code/modules/assimilation/state_eakf_tracer_bnrhf_algorithm_info_mod +++ /dev/null @@ -1,241 +0,0 @@ -! 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 - -module algorithm_info_mod - -use types_mod, only : r8, i8, missing_r8 - -use obs_def_mod, only : obs_def_type, get_obs_def_type_of_obs, get_obs_def_error_variance -use obs_kind_mod, only : get_quantity_for_type_of_obs - -! Get the QTY definitions that are needed (aka kind) -use obs_kind_mod, only : QTY_STATE_VARIABLE, QTY_STATE_VAR_POWER, QTY_TRACER_CONCENTRATION, & - QTY_TRACER_SOURCE -! NOTE: Sadly, the QTY itself is not sufficient for the POWER because there is additional metadata - -use assim_model_mod, only : get_state_meta_data -use location_mod, only : location_type - -use distribution_params_mod, only : NORMAL_DISTRIBUTION, BOUNDED_NORMAL_RH_DISTRIBUTION, & - GAMMA_DISTRIBUTION, BETA_DISTRIBUTION, LOG_NORMAL_DISTRIBUTION, UNIFORM_DISTRIBUTION, & - PARTICLE_FILTER_DISTRIBUTION - -implicit none -private - -! Defining parameter strings for different observation space filters -! For now, retaining backwards compatibility in assim_tools_mod requires using -! these specific integer values and there is no point in using these in assim_tools. -! That will change if backwards compatibility is removed in the future. -integer, parameter :: EAKF = 1 -integer, parameter :: ENKF = 2 -integer, parameter :: UNBOUNDED_RHF = 8 -integer, parameter :: GAMMA_FILTER = 11 -integer, parameter :: BOUNDED_NORMAL_RHF = 101 - -public :: obs_error_info, probit_dist_info, obs_inc_info, & - EAKF, ENKF, BOUNDED_NORMAL_RHF, UNBOUNDED_RHF, GAMMA_FILTER - -! Provides routines that give information about details of algorithms for -! observation error sampling, observation increments, and the transformations -! for regression and inflation in probit space. -! For now, it is convenient to have these in a single module since several -! users will be developing their own problem specific versions of these -! subroutines. This will avoid constant merge conflicts as other parts of the -! assimilation code are updated. - -contains - -!------------------------------------------------------------------------- -subroutine obs_error_info(obs_def, error_variance, & - bounded_below, bounded_above, lower_bound, upper_bound) - -! Computes information needed to compute error sample for this observation -! This is called by perfect_model_obs when generating noisy obs -type(obs_def_type), intent(in) :: obs_def -real(r8), intent(out) :: error_variance -logical, intent(out) :: bounded_below, bounded_above -real(r8), intent(out) :: lower_bound, upper_bound - -integer :: obs_type, obs_kind -integer(i8) :: state_var_index -type(location_type) :: temp_loc - -! Get the kind of the observation -obs_type = get_obs_def_type_of_obs(obs_def) -! If it is negative, it is an identity obs -if(obs_type < 0) then - state_var_index = -1 * obs_type - call get_state_meta_data(state_var_index, temp_loc, obs_kind) -else - obs_kind = get_quantity_for_type_of_obs(obs_type) -endif - -! Get the default error variance -error_variance = get_obs_def_error_variance(obs_def) - -! Set the observation error details for each type of quantity -if(obs_kind == QTY_STATE_VARIABLE) then - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_CONCENTRATION) then - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_SOURCE) then - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -else - write(*, *) 'Illegal obs_kind in obs_error_info' - stop -endif - -end subroutine obs_error_info - - -!------------------------------------------------------------------------- - - -subroutine probit_dist_info(kind, is_state, is_inflation, dist_type, & - bounded_below, bounded_above, lower_bound, upper_bound) - -! Computes the details of the probit transform for initial experiments -! with Molly - -integer, intent(in) :: kind -logical, intent(in) :: is_state ! True for state variable, false for obs -logical, intent(in) :: is_inflation ! True for inflation transform -integer, intent(out) :: dist_type -logical, intent(out) :: bounded_below, bounded_above -real(r8), intent(out) :: lower_bound, upper_bound - -! Have input information about the kind of the state or observation being transformed -! along with additional logical info that indicates whether this is an observation -! or state variable and about whether the transformation is being done for inflation -! or for regress. -! Need to select the appropriate transform. At present, options are NORMAL_PRIOR -! which does nothing or BOUNDED_NORMAL_RH_PRIOR. -! If the BNRH is selected then information about the bounds must also be set. -! The two dimensional logical array 'bounded' is set to false for no bounds and true -! for bounded. the first element of the array is for the lower bound, the second for the upper. -! If bounded is chosen, the corresponding bound value(s) must be set in the two dimensional -! real array 'bounds'. -! For example, if my_state_kind corresponds to a sea ice fraction then an appropriate choice -! would be: -! bounded_below = .true.; bounded_above = .true. -! lower_bound = 0.0_r8; upper_bounds = 1.0_r8 - -! In the long run, may not have to have separate controls for each of the input possibilities -! However, for now these are things that need to be explored for science understanding - -if(is_inflation) then - ! Case for inflation transformation - if(kind == QTY_STATE_VARIABLE) then - dist_type = NORMAL_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -elseif(is_state) then - ! Case for state variable priors - if(kind == QTY_STATE_VARIABLE) then - dist_type = NORMAL_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -else - ! This case is for observation (extended state) priors - if(kind == QTY_STATE_VARIABLE) then - dist_type = NORMAL_DISTRIBUTION - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_CONCENTRATION) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - elseif(kind == QTY_TRACER_SOURCE) then - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - else - write(*, *) 'Illegal kind in obs_error_info' - stop - endif -endif - -end subroutine probit_dist_info - -!------------------------------------------------------------------------ - - -subroutine obs_inc_info(obs_kind, filter_kind, rectangular_quadrature, gaussian_likelihood_tails, & - sort_obs_inc, spread_restoration, bounded_below, bounded_above, lower_bound, upper_bound) - -integer, intent(in) :: obs_kind -integer, intent(inout) :: filter_kind -logical, intent(inout) :: rectangular_quadrature, gaussian_likelihood_tails -logical, intent(inout) :: sort_obs_inc -logical, intent(inout) :: spread_restoration -logical, intent(inout) :: bounded_below, bounded_above -real(r8), intent(inout) :: lower_bound, upper_bound - -! The information arguments are all intent (inout). This means that if they are not set -! here, they retain the default values from the assim_tools_mod namelist. Bounds don't exist -! in that namelist, so default values are set in assim_tools_mod just before the call to here. - -! Temporary approach for setting the details of how to assimilate this observation -! This example is designed to reproduce the squared forward operator results from paper - - -! Set the observation increment details for each type of quantity -if(obs_kind == QTY_STATE_VARIABLE) then - filter_kind = EAKF - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_CONCENTRATION) then - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -elseif(obs_kind == QTY_TRACER_SOURCE) then - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 -else - write(*, *) 'Illegal obs_kind in obs_error_info' - stop -endif - -! Default settings for now for Icepack and tracer model tests -sort_obs_inc = .false. -spread_restoration = .false. - -! Only need to set these two for options the original RHF implementation -!!!rectangular_quadrature = .true. -!!!gaussian_likelihood_tails = .false. - -end subroutine obs_inc_info - -!------------------------------------------------------------------------ - -end module algorithm_info_mod diff --git a/assimilation_code/programs/perfect_model_obs/perfect_model_obs.f90 b/assimilation_code/programs/perfect_model_obs/perfect_model_obs.f90 index 97593321c8..6417785e2c 100644 --- a/assimilation_code/programs/perfect_model_obs/perfect_model_obs.f90 +++ b/assimilation_code/programs/perfect_model_obs/perfect_model_obs.f90 @@ -63,7 +63,7 @@ program perfect_model_obs use mpi_utilities_mod, only : my_task_id -use algorithm_info_mod, only : obs_error_info +use algorithm_info_mod, only : init_algorithm_info_mod, obs_error_info, end_algorithm_info_mod implicit none @@ -77,7 +77,6 @@ program perfect_model_obs !----------------------------------------------------------------------------- ! Namelist with default values ! -logical :: use_algorithm_info_mod = .true. logical :: read_input_state_from_file = .false. logical :: write_output_state_to_file = .false. integer :: async = 0 @@ -112,7 +111,7 @@ program perfect_model_obs obs_seq_out_file_name = 'obs_seq.out', & adv_ens_command = './advance_model.csh' -namelist /perfect_model_obs_nml/ use_algorithm_info_mod, read_input_state_from_file,& +namelist /perfect_model_obs_nml/ read_input_state_from_file,& write_output_state_to_file, & init_time_days, init_time_seconds, async, & first_obs_days, first_obs_seconds, & @@ -552,15 +551,8 @@ subroutine perfect_main() if( qc_ens_handle%vars(i, 1) == 0 ) then ! Get the information for generating error sample for this observation - if(use_algorithm_info_mod) then - call obs_error_info(obs_def, error_variance, & - bounded_below, bounded_above, lower_bound, upper_bound) - else - ! Default is unbounded with standard error_variance - error_variance = get_obs_def_error_variance(obs_def) - bounded_below = .false. ; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = 0.0_r8 - endif + call obs_error_info(obs_def, error_variance, & + bounded_below, bounded_above, lower_bound, upper_bound) ! Capability to do a bounded normal error if(bounded_below .and. bounded_above) then @@ -666,6 +658,8 @@ subroutine perfect_main() call destroy_obs_sequence(seq) call trace_message('After ensemble and obs memory cleanup') +call end_algorithm_info_mod() + call trace_message('Perfect_model done') call timestamp_message('Perfect_model done') @@ -693,6 +687,7 @@ subroutine perfect_initialize_modules_used() call state_vector_io_init() call initialize_qc() +call init_algorithm_info_mod() end subroutine perfect_initialize_modules_used diff --git a/assimilation_code/programs/perfect_model_obs/perfect_model_obs.nml b/assimilation_code/programs/perfect_model_obs/perfect_model_obs.nml index 68f0b8ddad..37a91b74cb 100644 --- a/assimilation_code/programs/perfect_model_obs/perfect_model_obs.nml +++ b/assimilation_code/programs/perfect_model_obs/perfect_model_obs.nml @@ -1,5 +1,4 @@ &perfect_model_obs_nml - use_algorithm_info_mod = .true., read_input_state_from_file = .false., single_file_in = .false., input_state_files = "", diff --git a/conf.py b/conf.py index 87e25bad53..eea345845e 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 = '10.8.5' +release = '11.1.0-alpha' root_doc = 'index' # -- General configuration --------------------------------------------------- diff --git a/developer_tests/qceff/test_table_read.f90 b/developer_tests/qceff/test_table_read.f90 new file mode 100644 index 0000000000..9e373dbeae --- /dev/null +++ b/developer_tests/qceff/test_table_read.f90 @@ -0,0 +1,24 @@ +! 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 + +! qceff_table_filename expected as command line arguement +program test_table_read + +use algorithm_info_mod, only : init_algorithm_info_mod, end_algorithm_info_mod +use utilities_mod, only : initialize_utilities, finalize_utilities + +implicit none + +character(len=129) :: qceff_table_filename + +call initialize_utilities('test_table_read') + +call get_command_argument(1,qceff_table_filename) + +call init_algorithm_info_mod(qceff_table_filename) +call end_algorithm_info_mod() + +call finalize_utilities() + +end program test_table_read diff --git a/developer_tests/qceff/work/all_bnrhf_qceff_table.csv b/developer_tests/qceff/work/all_bnrhf_qceff_table.csv new file mode 100644 index 0000000000..3e50dfea1c --- /dev/null +++ b/developer_tests/qceff/work/all_bnrhf_qceff_table.csv @@ -0,0 +1,5 @@ +QCF table version: 1,obs_error_info,,,,probit_inflation,,,,,probit_state,,,,,probit_extended_state,,,,,obs_inc_info,,,, +QTY_TEMPLATE:,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,filter_kind,bounded_below,bounded_above,lower_bound,upper_bound +QTY_state_VARIABLE,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RHF,.false.,.false.,-888888,-888888 +QTY_TRACER_CONCENTRATION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RHF,.true.,.false.,0,-888888 +QTY_TRACER_SOURCE,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RHF,.true.,.false.,0,-888888 diff --git a/developer_tests/qceff/work/input.nml b/developer_tests/qceff/work/input.nml new file mode 100644 index 0000000000..ffa7155436 --- /dev/null +++ b/developer_tests/qceff/work/input.nml @@ -0,0 +1,28 @@ +&utilities_nml + TERMLEVEL = 1, + module_details = .false. + logfilename = 'dart_log.out' + / + +# pick a random set of inputs +&preprocess_nml + overwrite_output = .true. + 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' + 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' + obs_type_files = '../../../observations/forward_operators/obs_def_reanalysis_bufr_mod.f90', + '../../../observations/forward_operators/obs_def_radar_mod.f90', + '../../../observations/forward_operators/obs_def_metar_mod.f90', + '../../../observations/forward_operators/obs_def_dew_point_mod.f90', + '../../../observations/forward_operators/obs_def_rel_humidity_mod.f90', + '../../../observations/forward_operators/obs_def_altimeter_mod.f90', + '../../../observations/forward_operators/obs_def_gps_mod.f90', + '../../../observations/forward_operators/obs_def_vortex_mod.f90', + '../../../observations/forward_operators/obs_def_gts_mod.f90', + '../../../observations/forward_operators/obs_def_QuikSCAT_mod.f90' + quantity_files = '../../../assimilation_code/modules/observations/default_quantities_mod.f90', + / + +&obs_kind_nml +/ diff --git a/developer_tests/qceff/work/qcf_table.txt b/developer_tests/qceff/work/qcf_table.txt new file mode 100644 index 0000000000..0449b56d87 --- /dev/null +++ b/developer_tests/qceff/work/qcf_table.txt @@ -0,0 +1,3 @@ +QCF table version: 1 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind bounded_below bounded_above lower_bound upper_bound +QTY_STATE_VARIABLE .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/qcf_table_bad_qty.txt b/developer_tests/qceff/work/qcf_table_bad_qty.txt new file mode 100644 index 0000000000..3933fc478b --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_bad_qty.txt @@ -0,0 +1,3 @@ +QCF table version: 1 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind bounded_below bounded_above lower_bound upper_bound +QTY_HAIRCUT .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/qcf_table_broke.txt b/developer_tests/qceff/work/qcf_table_broke.txt new file mode 100644 index 0000000000..6585f67485 --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_broke.txt @@ -0,0 +1,3 @@ +QCF table version: 1 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind bounded_below bounded_above lower_bound upper_bound +QTY_STATE_VARIABLE .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/qcf_table_duplicates.txt b/developer_tests/qceff/work/qcf_table_duplicates.txt new file mode 100644 index 0000000000..7c4dfd9bd9 --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_duplicates.txt @@ -0,0 +1,6 @@ +QCF table version: 1 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind bounded_below bounded_above lower_bound upper_bound +QTY_VEGETATED_AREA_FRACTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_AQUIFER_WATER .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_SEAICE_SALINITY008 .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_VEGETATED_AREA_FRACTION .true. .false. 0 1000 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/qcf_table_extra_columns.txt b/developer_tests/qceff/work/qcf_table_extra_columns.txt new file mode 100644 index 0000000000..3f1236b2b6 --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_extra_columns.txt @@ -0,0 +1,3 @@ +QCF table version: 1 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind bounded_below bounded_above lower_bound upper_bound frog +QTY_STATE_VARIABLE .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 toad newt diff --git a/developer_tests/qceff/work/qcf_table_incorrect_distribution.txt b/developer_tests/qceff/work/qcf_table_incorrect_distribution.txt new file mode 100644 index 0000000000..93d10d7869 --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_incorrect_distribution.txt @@ -0,0 +1,3 @@ +QCF table version: 1 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind bounded_below bounded_above lower_bound upper_bound +QTY_STATE_VARIABLE .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 POLAR_BEAR_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/qcf_table_incorrect_filter_kind.txt b/developer_tests/qceff/work/qcf_table_incorrect_filter_kind.txt new file mode 100644 index 0000000000..38d567f833 --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_incorrect_filter_kind.txt @@ -0,0 +1,3 @@ +QCF table version: 1 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind bounded_below bounded_above lower_bound upper_bound +QTY_STATE_VARIABLE .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 PENGUIN_FILTER .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/qcf_table_lower_bound_only.txt b/developer_tests/qceff/work/qcf_table_lower_bound_only.txt new file mode 100644 index 0000000000..3916443577 --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_lower_bound_only.txt @@ -0,0 +1,6 @@ +QCF table version: 1 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind bounded_below bounded_above lower_bound upper_bound +QTY_STATE_VARIABLE .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_AQUIFER_WATER .true. .false. 0.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_SEAICE_SALINITY008 .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_VEGETATED_AREA_FRACTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/qcf_table_lower_case_dist.txt b/developer_tests/qceff/work/qcf_table_lower_case_dist.txt new file mode 100644 index 0000000000..b6750b66ec --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_lower_case_dist.txt @@ -0,0 +1,3 @@ +QCF table version: 1 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind bounded_below bounded_above lower_bound upper_bound +QTY_STATE_VARIABLE .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_normal_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/qcf_table_lower_gt_upper.txt b/developer_tests/qceff/work/qcf_table_lower_gt_upper.txt new file mode 100644 index 0000000000..4afb33e579 --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_lower_gt_upper.txt @@ -0,0 +1,6 @@ +QCF table version: 1 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind bounded_below bounded_above lower_bound upper_bound +QTY_STATE_VARIABLE .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_AQUIFER_WATER .true. .true. 10.0 0.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_SEAICE_SALINITY008 .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_VEGETATED_AREA_FRACTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/qcf_table_no_bounds_with_values.txt b/developer_tests/qceff/work/qcf_table_no_bounds_with_values.txt new file mode 100644 index 0000000000..c987d2f9e9 --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_no_bounds_with_values.txt @@ -0,0 +1,6 @@ +QCF table version: 1 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind spread_restoration bounded_below bounded_above lower_bound upper_bound +QTY_STATE_VARIABLE .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_AQUIFER_WATER .false. .false. 10.0 0.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_SEAICE_SALINITY008 .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 +QTY_VEGETATED_AREA_FRACTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/qcf_table_no_header.txt b/developer_tests/qceff/work/qcf_table_no_header.txt new file mode 100644 index 0000000000..617379ffe0 --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_no_header.txt @@ -0,0 +1,4 @@ +QTY_STATE_VARIABLE .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. .false. .false. .false. .false. -888888.0 -888888.0 +QTY_AQUIFER_WATER .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. .false. .false. .false. .false. -888888.0 -888888.0 +QTY_SEAICE_SALINITY008 .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. .false. .false. .false. .false. -888888.0 -888888.0 +QTY_VEGETATED_AREA_FRACTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. .false. .false. .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/qcf_table_v2.txt b/developer_tests/qceff/work/qcf_table_v2.txt new file mode 100644 index 0000000000..0d1b5a46df --- /dev/null +++ b/developer_tests/qceff/work/qcf_table_v2.txt @@ -0,0 +1,3 @@ +QCF table version: 2 +QTY_TEMPLATE: bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound dist_type bounded_below bounded_above lower_bound upper_bound filter_kind rectangular_quadrature gaussian_likelihood_tails sort_obs_inc spread_restoration bounded_below bounded_above lower_bound upper_bound +QTY_STATE_VARIABLE .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RH_DISTRIBUTION .false. .false. -888888.0 -888888.0 BOUNDED_NORMAL_RHF .false. .false. .false. .false. .false. .false. -888888.0 -888888.0 diff --git a/developer_tests/qceff/work/quickbuild.sh b/developer_tests/qceff/work/quickbuild.sh new file mode 100755 index 0000000000..81b1308494 --- /dev/null +++ b/developer_tests/qceff/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/buildfunctions.sh + +MODEL="none" +EXTRA="$DART"/models/template/threed_model_mod.f90 +dev_test=1 +TEST="qceff" +LOCATION="threed_sphere" + +serial_programs=( +test_table_read +) + +# quickbuild arguments +arguments "$@" + +# clean the directory +\rm -f -- *.o *.mod Makefile .cppdefs + +# build and run preprocess before making any other DART executables +buildpreprocess + +# build +buildit + +# clean up +\rm -f -- *.o *.mod + +} + +main "$@" diff --git a/developer_tests/qceff/work/runall.sh b/developer_tests/qceff/work/runall.sh new file mode 100755 index 0000000000..c354f9b4d8 --- /dev/null +++ b/developer_tests/qceff/work/runall.sh @@ -0,0 +1,56 @@ +#!/bin/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 + +# Usage: +# ./runall.sh +# ./runall.sh | grep FAIL +# ./runall.sh | grep PASS + + +should_pass () { +if [[ $? -ne 0 ]]; then + echo $1: "FAIL" +else + echo $1: "PASS" +fi +} + +should_fail () { +if [[ $? -eq 0 ]]; then + echo $1: "FAIL" +else + echo $1: "PASS" +fi +} + +./test_table_read ; should_pass "no table" + +./test_table_read qcf_table.txt ; should_pass "correct v1 table" + +./test_table_read qcf_table_v2.txt ; should_fail "detect wrong version" + +./test_table_read qcf_table_extra_columns.txt ; should_pass "extra colums" + +./test_table_read qcf_table_bad_qty.txt ; should_fail "bad qty" + +./test_table_read qcf_table_broke.txt ; should_fail "bad value" + +./test_table_read qcf_table_no_header.txt ; should_fail "no header" + +./test_table_read qcf_table_lower_gt_upper.txt ; should_fail "upper bound less than lower" + +./test_table_read qcf_table_lower_bound_only.txt ; should_pass "lower bound only" + +./test_table_read qcf_table_no_bounds_with_values.txt ; should_pass "bounds false, values for bounds" + +./test_table_read qcf_table_incorrect_filter_kind.txt ; should_fail "incorrect filter_kind" + +./test_table_read qcf_table_incorrect_distribution.txt ; should_fail "incorrect distribution" + +./test_table_read all_bnrhf_qceff_table.csv ; should_pass "lower case QTY" + +./test_table_read qcf_table_lower_case_dist.txt; should_pass "lower case dist_type" + diff --git a/guide/qceff-examples.rst b/guide/qceff-examples.rst new file mode 100644 index 0000000000..123a9a3d9c --- /dev/null +++ b/guide/qceff-examples.rst @@ -0,0 +1,184 @@ +.. _quantile tracer: + +QCEFF: Examples with the Lorenz 96 Tracer Model +=============================================== + + +The Quantile-Conserving Ensemble Filter Framework (QCEFF) tools are in the alpha release stage. +The DART development team (dart@ucar.edu) would be happy to hear about your experiences and is +anxious to build scientific collaborations using these new capabilities. + +To get started, make sure that you are on the quantile_methods branch of DART: + +.. code-block:: text + + git checkout quantile_methods + +Build the DART executables for the Lorenz 96 tracer advection model: + +.. code-block:: text + + cd DART/models/lorenz_96_tracer_advection/work + ./quickbuild.sh nompi + + +The new quantile options are set using a :ref:`qceff table ` given as a namelist +option ``qceff_table_filename`` to &algorithm_info_nml. The examples below show how to change the quantile options +using various qceff tables. You can find the .csv files for these four example in the directory +``DART/models/lorenz_96_tracer_advection/work`` + + +.. list-table:: + :header-rows: 1 + :widths: 15 60 25 + + * - example + - description + - .cvs filename + * - Example A + - boundend normal rank histogram + - all_bnrhf_qceff_table.csv + * - Example B + - Ensemble Adjustment Kalman filters + - all_eakf_qceff_table.csv + * - Example C + - EAKF for state and bounded normal rank histogram filter and priors for tracer concentration and source + - state_eakf_tracer_bnrhf_qceff_table.csv + * - Example D + - Negative tracers bounded above + - neg_qceff_table.csv + + +You can view .csv files with a text editor, or spreadsheet tool such as Google Sheets, +or Microsoft Excel. + +Example A +---------- + +Assimilating observations of state (wind) and tracer concentration using +a rank histogram observation space filter and rank histogram probit transforms for +state variable updates. This example includes adaptive inflation. + +The default model configuration has a single tracer source at gridpoint 1 along with +small uniform tracer sinks that lead to areas where the true tracer concentration is +usually 0. This is a particularly tough test for ensemble methods. + +#. Edit input.nml to set the qceff_table_filename to 'all_bnrhf_qceff_table.csv' + + .. code-block:: text + + &algorithm_info_nml + qceff_table_filename = 'all_bnrhf_qceff_table.csv' + + +#. Create a set_def.out file using create_obs_sequence, + + ``./create_obs_sequence < create_obs_sequence_input`` + +#. Create an obs_sequence.in file using create_fixed_network_seq + + ``./create_fixed_network_seq`` + + .. code:: text + + Select the default input filename , + Create a regularly repeating sequence by entering "1", + Enter "1000" for the number of observation times, + Enter "0 0" for the initial time, + Enter "0 10800" for the period, + Select the default output filename, + +#. Spin-up a model initial condition by running perfect_model_obs + + ``./perfect_model_obs`` + +#. Generate a spun-up true time series, + + ``cp perfect_output.nc perfect_input.nc`` + + + Edit input.nml to set read_input_state_from_file to .true. + + .. code:: text + + &perfect_model_obs_nml + read_input_state_from_file = .true., + + + Run ``./perfect_model_obs`` again. + +#. Run a filter assimilation, + + ``./filter`` + +#. Examine the output with your favorite tools. Looking at the analysis ensemble + for the tracer_concentration variables with indices near the source (location 1) + and far downstream from the source (location 35) is interesting. Note that the + source estimation capabilities of the model and filters are not being tested here. + + +Example B +--------- + +Using Ensemble Adjustment Kalman filters. + + +#. Edit input.nml to set the qceff_table_filename to 'all_eakf_qceff_table.csv' + + .. code-block:: text + + &algorithm_info_nml + qceff_table_filename = 'all_eakf_qceff_table.csv' + + +#. Run the filter + + ``./filter`` + +Example C +--------- + +Using Ensemble Adjustment Kalman filter for state, but bounded normal rank histogram filter and priors for tracer concentration and source. + + +#. Edit input.nml to set the qceff_table_filename to state_eakf_tracer_bnrhf_qceff_table.csv + + .. code-block:: text + + &algorithm_info_nml + qceff_table_filename = 'state_eakf_tracer_bnrhf_qceff_table.csv' + + +#. Run the filter + + ``./filter`` + +Example D +---------- + +Testing the bounded above option. Normally tracers are bounded below, but there are other quantities that may be bounded +above. There are distinct numerical challenges in implementing the quantile algorithms +for quantities that are bounded above, so flipping the sign of the tracers is a good +test. + +#. Edit input.nml to set the qceff_table_filename to neg_qceff_table.csv + + .. code-block:: text + + &algorithm_info_nml + qceff_table_filename = 'neg_qceff_table.csv' + + +#. Edit input.nml, to change the entry positive_tracer to .false. and read_input_state_from_file back to .false. + + + .. code-block:: text + + &model_nml + positive_tracer = .false., + + &perfect_model_obs_nml + read_input_state_from_file = .false., + + +#. Repeat steps 3-6 from Test A. diff --git a/guide/qceff_probit.rst b/guide/qceff_probit.rst new file mode 100644 index 0000000000..34e7478702 --- /dev/null +++ b/guide/qceff_probit.rst @@ -0,0 +1,152 @@ +.. _QCEFF: + +Quantile-Conserving Ensemble Filter Framework +============================================== + +The Quantile-Conserving Ensemble Filter Framework (QCEFF) tools are in the alpha release stage. +The DART development team (dart@ucar.edu) would be happy to hear about your experiences +and is anxious to build scientific collaborations using these new capabilities. + +The QCEFF options are set using a :ref:`qceff table ` given as a namelist option to &algorithm_info_nml. + + .. code-block:: text + + &algorithm_info_nml + qceff_table_filename = 'qceff_table.csv' + + +.. _QCEFF options: + +QCEFF options +-------------- + +QCEFF options are per quantity. For a given quantity, you specify the following +options as columns of the qceff_table: + +* Observation error information + + Used to compute sample for this observation when using perfect_model_obs + to generate noisy observations. + + * bounded_below (default .false.) + * bounded_above (default .false.) + * lower_bound + * upper_bound + + +* Probit distribution information + + Used in the computation of the probit transform. + The values needed are the bounds and the distribution type. + These can be different for all three cases (inflation, state, and extended_state priors) + + * distribution (one of :ref:`Distributions`) + * bounded_below (default .false.) + * bounded_above (default .false.) + * lower_bound + * upper_bound + + +* Observation increment information + + * filter_kind (one of :ref:`Filter kinds`) + * bounded_below (default .false.) + * bounded_above (default .false.) + * lower_bound + * upper_bound + + + +.. _qceff table: + +Creating a qceff table +----------------------- + +The table has two headers, row 1 and 2. +The first row is the version number. The second row describes the :ref:`QCEFF options` corresponding to each column of the table. +These two headers must be present in your qceff table. +The :ref:`qcf trunc table` below shows the table headers, +and an example quantity QTY_TRACER_CONCENTRATION for the first 5 columns of the table. +There is a complete table with all 25 columns in `Google Sheets `_. You can copy and edit this table as needed. + +To add a quantity, add a row to the table. +For any quantity not listed in the table, the :ref:`Default values` values will be used for all 25 options. +You only have to add rows for quantitiies that use non-default values for any of the input options. +Ensure that there are no empty rows in between the quantities listed in the spreadsheet. +Save your spreadsheet as a .csv file. + +To run filter or perfect_model_obs, put the .csv file in the directory where you are running. +Edit input.nml to set the algorithm_info_nml option qceff_table_filename, for example: + + + .. code-block:: text + + &algorithm_info_nml + qceff_table_filename = 'qceff_table.csv' + + +.. _qcf trunc table: + +.. list-table:: truncated table + :header-rows: 2 + + * - QCF table version: 1 + - + - + - + - + * - QTY + - bounded_below + - bounded_above + - lower_bound + - upper_bound + * - QTY_TRACER_CONCENTRATION + - .true. + - .false. + - 0 + - -888888 + + +.. _Filter kinds: + +Available filter kinds +----------------------- + + * EAKF (default) + * ENKF + * UNBOUNDED_RH + * GAMMA_FILTER + * BOUNDED_NORMAL_RHF + +.. _Distributions: + +Available distributions +------------------------ + + * NORMAL_DISTRIBUTION (default) + * BOUNDED_NORMAL_RH_DISTRIBUTION + * GAMMA_DISTRIBUTION + * BETA_DISTRIBUTION + * LOG_NORMAL_DISTRIBUTION + * UNIFORM_DISTRIBUTION + * PARTICLE_FILTER_DISTRIBUTION + + + +.. _Default values: + +Default values +--------------- + +If a quantity is not in the qceff table, the following default values +are used: + + * filter_kind (default EAKF) + * dist_type (default NORMAL_DISTRIBUTION) + * bounded_below (default .false.) + * bounded_above (default .false.) + * lower_bound (default -888888) + * upper_bound (default -888888) + +Note -888888 is a missing value in DART. + diff --git a/index.rst b/index.rst index 1aedcc7f87..16ff70bd10 100644 --- a/index.rst +++ b/index.rst @@ -5,10 +5,8 @@ Welcome to the Data Assimilation Research Testbed .. warning:: - Pre-release version of DART: quantile conserving and probit transform tools + Pre-release version of DART: Quantile-Conserving Ensemble Filter Framework - To get started, see the :ref:`tracer advection example` - The Data Assimilation Research Testbed (DART) is an open-source, freely available community facility for ensemble data assimilation (DA). [1]_ DART is @@ -51,6 +49,7 @@ DART includes: - Extensive documentation of its source code. - Interfaces to a variety of models and observation sets that can be used to introduce new users or graduate students to ensemble DA. +- Nonlinear and Non-Gaussian DA Capabilities DART is also designed to facilitate the combination of assimilation algorithms, models, and real or synthetic observations to allow increased @@ -62,6 +61,50 @@ These tools are intended for use by the full range of geosciencies community: beginners and experts; students and teachers; national centers and university research labs. +Nonlinear and Non-Gaussian Data Assimilation Capabilities in DART +----------------------------------------------------------------- + +The default DART algorithms assume a normal distribution to compute ensemble increments +for the observed quantity (this is the ensemble adjustment Kalman filter, or EAKF) and +then linearly regress the observation increments onto each state variable. + + +DART now implements a Quantile-Conserving Ensemble Filtering Framework :ref:`(QCEFF) `. +The QCEFF provides a very general method of computing increments for the prior ensemble of +an observed quantity by allowing the use of arbitrary distributions for the prior and the +observation error. This is especially useful for bounded quantities like tracer concentrations, +depths of things like snow or ice, and estimating model parameters that have a restricted range. +See this Monthly Weather Review article for details, +`QCEFF part1 `_. + +While the QCEFF for computing observation increments can lead to significant improvements in +analysis estimates for observed variables, those improvements can be lost when using standard +linear regression of observation increments to update other state variables. The QCEFF also +implements a capability to do regression in a probit probability integral transformed space. +Doing the regression of observation quantile increments in the transformed space guarantees +that the posterior ensembles for state variables also retain the advantages of the observation space +quantile conserving posteriors. For example, if state variables are bounded, then posterior +ensembles will respect the bounds. The posterior ensembles also respect other aspects of the +continuous prior distributions. See this Monthly Weather Review article for details, +`QCEFF part 2 `_. + +Inflation and localization, methods that improve the quality of ensemble DA, can also negate +the advantages of the QCEFF methods. For this reason, both localization and inflation can be +done in the probit-transformed quantile space as well. Combining these new methods can +significantly improve data assimilation for non-Gaussian quantities in Earth system models by +extending the capabilities of ensemble DA to general non-Gaussian and nonlinear distributions. +Transformative improvements in DA can result for many applications. The largest improvements are +found for bounded variables like tracer concentrations, snow and ice depth, soil moisture, and +similar quantities in other parts of the Earth system. Model parameters can also be estimated +with DA and large improvements can occur for bounded parameters. Variables that have distinctly +non-Gaussian prior distributions can also see large improvements. Examples can include atmospheric +quantities like moisture and cloud amount in the presence of convection, and many land surface variables. + +For instructions on how to use these tools, see :ref:`QCEFF`. + +For step-by-step examples of the QCEFF tools, you can work through +:ref:`examples with the Lorenz 96 tracer model ` + Organization of the documentation --------------------------------- @@ -230,7 +273,6 @@ References guide/downloading-dart guide/compiling-dart guide/verifying-installation - models/lorenz_96_tracer_advection/work/readme .. toctree:: :maxdepth: 2 @@ -250,6 +292,7 @@ References guide/high-level-da-workflows guide/dart-design-philosophy guide/important-capabilities-dart + guide/qceff_probit .. toctree:: :maxdepth: 2 @@ -366,6 +409,7 @@ References guide/DART_LAB/DART_LAB CLM-DART Tutorial WRF-DART Tutorial + guide/qceff-examples.rst .. toctree:: :maxdepth: 2 diff --git a/models/9var/work/input.nml b/models/9var/work/input.nml index cfe617f0ce..4ace01b342 100644 --- a/models/9var/work/input.nml +++ b/models/9var/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -89,7 +96,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 1000000.0, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/FESOM/work/input.nml b/models/FESOM/work/input.nml index 5334abff8c..f8b1d7304b 100644 --- a/models/FESOM/work/input.nml +++ b/models/FESOM/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true. input_state_files = "ENS01.2009.oce.nc" @@ -87,7 +94,6 @@ / &assim_tools_nml - filter_kind = 1 cutoff = 0.005 sort_obs_inc = .false. spread_restoration = .false. diff --git a/models/LMDZ/work/input.nml b/models/LMDZ/work/input.nml index edbf8617c5..a3d72caa37 100644 --- a/models/LMDZ/work/input.nml +++ b/models/LMDZ/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &filter_nml async = 2, adv_ens_command = "./advance_model.csh", @@ -74,7 +81,6 @@ perturbation_amplitude = 0.2 / &assim_tools_nml - filter_kind = 1, cutoff = 0.2, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/MITgcm_annulus/work/input.nml b/models/MITgcm_annulus/work/input.nml index b7301418ea..e91333406d 100644 --- a/models/MITgcm_annulus/work/input.nml +++ b/models/MITgcm_annulus/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml start_from_restart = .false., output_restart = .true., @@ -74,7 +81,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 0.15, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/MITgcm_ocean/work/input.nml b/models/MITgcm_ocean/work/input.nml index 9bc7844826..c8ddd50476 100644 --- a/models/MITgcm_ocean/work/input.nml +++ b/models/MITgcm_ocean/work/input.nml @@ -86,7 +86,6 @@ # in both lists, the same number of items &assim_tools_nml - filter_kind = 1 cutoff = 0.025 distribute_mean = .false. sort_obs_inc = .false. diff --git a/models/MOM6/work/input.nml b/models/MOM6/work/input.nml index 860e574da4..e442ab8c2a 100644 --- a/models/MOM6/work/input.nml +++ b/models/MOM6/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -86,7 +93,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 1000000.0 sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/NAAPS/work/input.nml b/models/NAAPS/work/input.nml index 13f65c076f..92c3ecf452 100644 --- a/models/NAAPS/work/input.nml +++ b/models/NAAPS/work/input.nml @@ -1,5 +1,11 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &assim_tools_nml - filter_kind = 1 cutoff = 0.03 sort_obs_inc = .false. spread_restoration = .false. diff --git a/models/NCOMMAS/work/input.nml b/models/NCOMMAS/work/input.nml index 79a7b30c8e..f8c477f70b 100644 --- a/models/NCOMMAS/work/input.nml +++ b/models/NCOMMAS/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml start_from_restart = .true., output_restart = .true., @@ -73,7 +80,6 @@ # cutoff of 0.03 (radians) is about 200km &assim_tools_nml - filter_kind = 1, cutoff = 0.20, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/POP/work/input.nml b/models/POP/work/input.nml index 8863032412..6fb54e6f2c 100644 --- a/models/POP/work/input.nml +++ b/models/POP/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true. single_file_in = .false. @@ -88,7 +95,6 @@ # if running a smaller pop case, use false to run faster. # &assim_tools_nml - filter_kind = 1 cutoff = 0.20 sort_obs_inc = .false. spread_restoration = .false. diff --git a/models/ROMS/work/input.nml b/models/ROMS/work/input.nml index 349bf35dcf..ad4ccd3aca 100644 --- a/models/ROMS/work/input.nml +++ b/models/ROMS/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true. single_file_in = .false. @@ -87,7 +94,6 @@ &assim_tools_nml - filter_kind = 1 cutoff = 0.02 sort_obs_inc = .false. spread_restoration = .false. diff --git a/models/am2/work/input.nml b/models/am2/work/input.nml index 5f453d01cd..3643d228bc 100644 --- a/models/am2/work/input.nml +++ b/models/am2/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml start_from_restart = .true., output_restart = .true., @@ -72,7 +79,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 0.2, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/bgrid_solo/work/input.nml b/models/bgrid_solo/work/input.nml index 9d65871b00..eebd179d2c 100644 --- a/models/bgrid_solo/work/input.nml +++ b/models/bgrid_solo/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -92,7 +99,6 @@ &assim_tools_nml - filter_kind = 1, cutoff = 0.20, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/cam-fv/work/algorithm_info_mod.f90 b/models/cam-fv/work/algorithm_info_mod.f90 deleted file mode 100644 index 15bfcc3bc8..0000000000 --- a/models/cam-fv/work/algorithm_info_mod.f90 +++ /dev/null @@ -1,216 +0,0 @@ -! 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 - -module algorithm_info_mod - -use types_mod, only : r8, i8, missing_r8 - -use obs_def_mod, only : obs_def_type, get_obs_def_type_of_obs, get_obs_def_error_variance -use obs_kind_mod, only : get_quantity_for_type_of_obs - -! Get the QTY definitions that are needed (aka kind) -use obs_kind_mod, only : QTY_U_WIND_COMPONENT, QTY_V_WIND_COMPONENT, QTY_SURFACE_PRESSURE, & - QTY_TEMPERATURE, QTY_SPECIFIC_HUMIDITY, QTY_CLOUD_LIQUID_WATER, & - QTY_CLOUD_ICE, QTY_GPSRO - -use assim_model_mod, only : get_state_meta_data -use location_mod, only : location_type - -use distribution_params_mod, only : NORMAL_DISTRIBUTION, BOUNDED_NORMAL_RH_DISTRIBUTION, & - GAMMA_DISTRIBUTION, BETA_DISTRIBUTION, LOG_NORMAL_DISTRIBUTION, UNIFORM_DISTRIBUTION, & - PARTICLE_FILTER_DISTRIBUTION - -implicit none -private - -! Defining parameter strings for different observation space filters -! For now, retaining backwards compatibility in assim_tools_mod requires using -! these specific integer values and there is no point in using these in assim_tools. -! That will change if backwards compatibility is removed in the future. -integer, parameter :: EAKF = 1 -integer, parameter :: ENKF = 2 -integer, parameter :: UNBOUNDED_RHF = 8 -integer, parameter :: GAMMA_FILTER = 11 -integer, parameter :: BOUNDED_NORMAL_RHF = 101 - -public :: obs_error_info, probit_dist_info, obs_inc_info, & - EAKF, ENKF, BOUNDED_NORMAL_RHF, UNBOUNDED_RHF, GAMMA_FILTER - -! Provides routines that give information about details of algorithms for -! observation error sampling, observation increments, and the transformations -! for regression and inflation in probit space. -! For now, it is convenient to have these in a single module since several -! users will be developing their own problem specific versions of these -! subroutines. This will avoid constant merge conflicts as other parts of the -! assimilation code are updated. - -contains - -!------------------------------------------------------------------------- -subroutine obs_error_info(obs_def, error_variance, & - bounded_below, bounded_above, lower_bound, upper_bound) - -! Computes information needed to compute error sample for this observation -! This is called by perfect_model_obs when generating noisy obs -type(obs_def_type), intent(in) :: obs_def -real(r8), intent(out) :: error_variance -logical, intent(out) :: bounded_below, bounded_above -real(r8), intent(out) :: lower_bound, upper_bound - -integer :: obs_type, obs_kind -integer(i8) :: state_var_index -type(location_type) :: temp_loc - -! Get the kind of the observation -obs_type = get_obs_def_type_of_obs(obs_def) -! If it is negative, it is an identity obs -if(obs_type < 0) then - state_var_index = -1 * obs_type - call get_state_meta_data(state_var_index, temp_loc, obs_kind) -else - obs_kind = get_quantity_for_type_of_obs(obs_type) -endif - -! Get the default error variance -error_variance = get_obs_def_error_variance(obs_def) - -! Set the observation error details for each type of quantity -bounded_below = .false.; bounded_above = .false. -lower_bound = missing_r8; upper_bound = missing_r8 - -end subroutine obs_error_info - - -!------------------------------------------------------------------------- - - -subroutine probit_dist_info(kind, is_state, is_inflation, dist_type, & - bounded_below, bounded_above, lower_bound, upper_bound) - -! Computes the details of the probit transform for initial experiments -! with Molly - -integer, intent(in) :: kind -logical, intent(in) :: is_state ! True for state variable, false for obs -logical, intent(in) :: is_inflation ! True for inflation transform -integer, intent(out) :: dist_type -logical, intent(out) :: bounded_below, bounded_above -real(r8), intent(out) :: lower_bound, upper_bound - -! Have input information about the kind of the state or observation being transformed -! along with additional logical info that indicates whether this is an observation -! or state variable and about whether the transformation is being done for inflation -! or for regress. -! Need to select the appropriate transform. At present, options are NORMAL_PRIOR -! which does nothing or BOUNDED_NORMAL_RH_PRIOR. -! If the BNRH is selected then information about the bounds must also be set. -! The two dimensional logical array 'bounded' is set to false for no bounds and true -! for bounded. the first element of the array is for the lower bound, the second for the upper. -! If bounded is chosen, the corresponding bound value(s) must be set in the two dimensional -! real array 'bounds'. -! For example, if my_state_kind corresponds to a sea ice fraction then an appropriate choice -! would be: -! bounded_below = .true.; bounded_above = .true. -! lower_bound = 0.0_r8; upper_bounds = 1.0_r8 - -! In the long run, may not have to have separate controls for each of the input possibilities -! However, for now these are things that need to be explored for science understanding - -select case(kind) - case(QTY_U_WIND_COMPONENT, QTY_V_WIND_COMPONENT, QTY_SURFACE_PRESSURE, QTY_TEMPERATURE) - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION -! dist_type = NORMAL_PRIOR - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - -!-------------- - case(QTY_SPECIFIC_HUMIDITY) - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION -! dist_type = NORMAL_PRIOR -! bounded_below = .false.; bounded_above = .false. - bounded_below = .true.; bounded_above = .true. - lower_bound = 0.0_r8; upper_bound = 1.0_r8 - -!-------------- - case(QTY_CLOUD_LIQUID_WATER, QTY_CLOUD_ICE) - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION -! dist_type = NORMAL_PRIOR -! bound_below = .false.; bounded_above = .false. - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - -!-------------- - case(QTY_GPSRO) - dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION -! dist_type = NORMAL_PRIOR -! bounded_below = .false.; bounded_above = .false. - bounded_below = .true.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - -!-------------- - case DEFAULT - write(*, *) 'Unexpected QTY in algorithm_info_mod ', kind - stop -end select - -end subroutine probit_dist_info - -!------------------------------------------------------------------------ - - -subroutine obs_inc_info(obs_kind, filter_kind, rectangular_quadrature, gaussian_likelihood_tails, & - sort_obs_inc, spread_restoration, bounded_below, bounded_above, lower_bound, upper_bound) - -integer, intent(in) :: obs_kind -integer, intent(inout) :: filter_kind -logical, intent(inout) :: rectangular_quadrature, gaussian_likelihood_tails -logical, intent(inout) :: sort_obs_inc -logical, intent(inout) :: spread_restoration -logical, intent(inout) :: bounded_below, bounded_above -real(r8), intent(inout) :: lower_bound, upper_bound - -! The information arguments are all intent (inout). This means that if they are not set -! here, they retain the default values from the assim_tools_mod namelist. Bounds don't exist -! in that namelist, so default values are set in assim_tools_mod just before the call to here. - -! Temporary approach for setting the details of how to assimilate this observation -! This example is designed to reproduce the squared forward operator results from paper - -select case(obs_kind) - case(QTY_U_WIND_COMPONENT, QTY_V_WIND_COMPONENT, QTY_SURFACE_PRESSURE, QTY_TEMPERATURE) - ! Set the observation increment details for each type of quantity - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .false.; bounded_above = .false. - lower_bound = missing_r8; upper_bound = missing_r8 - - case(QTY_GPSRO) - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .true.; bounded_above = .false. -! bounded_below = .false.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = missing_r8 - - case(QTY_SPECIFIC_HUMIDITY) - filter_kind = BOUNDED_NORMAL_RHF - bounded_below = .true.; bounded_above = .true. -! bounded_below = .false.; bounded_above = .false. - lower_bound = 0.0_r8; upper_bound = 1.0_r8 - - case DEFAULT - write(*, *) 'Unexpected QTY in algorithm_info_mod ', obs_kind - stop -end select - -! Default settings for now for Icepack and tracer model tests -sort_obs_inc = .false. -spread_restoration = .false. - -! Only need to set these two for options the original RHF implementation -!!!rectangular_quadrature = .true. -!!!gaussian_likelihood_tails = .false. - -end subroutine obs_inc_info - -!------------------------------------------------------------------------ - -end module algorithm_info_mod diff --git a/models/cam-fv/work/input.nml b/models/cam-fv/work/input.nml index 8c6d2abbba..8fd9fa689b 100644 --- a/models/cam-fv/work/input.nml +++ b/models/cam-fv/work/input.nml @@ -33,14 +33,14 @@ ! inf_sd_from_restart inflation restart files from the values in inf*_initial ! if needed. -&quantile_distributions_nml - fix_bound_violations = .true., - use_logit_instead_of_probit = .false. - do_inverse_check = .false. +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' / &filter_nml - use_algorithm_info_mod = .true. input_state_file_list = 'cam_init_files' input_state_files = '' single_file_in = .false. @@ -360,8 +360,6 @@ &assim_tools_nml - use_algorithm_info_mod = .true. - filter_kind = 1 cutoff = 0.15 sort_obs_inc = .false. spread_restoration = .false. diff --git a/models/cam-fv/work/new_qcf_input.nml b/models/cam-fv/work/new_qcf_input.nml deleted file mode 100644 index 92d4f8b95e..0000000000 --- a/models/cam-fv/work/new_qcf_input.nml +++ /dev/null @@ -1,436 +0,0 @@ -&probit_transform_nml - fix_bound_violations = .true., - use_logit_instead_of_probit = .false. - do_inverse_check = .false. - / - -&filter_nml - use_algorithm_info_mod = .true. - input_state_file_list = 'cam_init_files' - input_state_files = '' - single_file_in = .false. - perturb_from_single_instance = .false. - init_time_days = -1 - init_time_seconds = -1 - - stages_to_write = 'forecast','output' - - output_state_files = '' - output_state_file_list = 'cam_init_files' - output_mean = .true. - output_sd = .true. - output_members = .true. - num_output_state_members = 80 - single_file_out = .false. - write_all_stages_at_end = .false. - output_interval = 1 - - ens_size = 80 - num_groups = 1 - distributed_state = .true. - - inf_flavor = 5, 0 - inf_initial_from_restart = .true., .false. - inf_initial = 1.0, 1.0 - inf_lower_bound = 0.0, 0.0 - inf_upper_bound = 100.0, 100.0 - inf_sd_initial_from_restart = .true., .false. - inf_sd_initial = 0.6, 0.6 - inf_sd_lower_bound = 0.6, 0.6 - inf_sd_max_change = 1.05, 1.05 - inf_damping = 0.9, 0.9 - inf_deterministic = .true., .true. - - obs_sequence_in_name = 'obs_seq.out' - obs_sequence_out_name = 'obs_seq.final' - num_output_obs_members = 80 - compute_posterior = .true. - - trace_execution = .true. - output_timestamps = .true. - output_forward_op_errors = .false. - silence = .false. - / - - - first_obs_days = -1 - first_obs_seconds = -1 - last_obs_days = -1 - last_obs_seconds = -1 - obs_window_days = -1 - obs_window_seconds = -1 - adv_ens_command = 'no_CESM_advance_script' - tasks_per_model_advance = -1 Used only for models run inside filter. - write_obs_every_cycle = .false. intended for debugging when cycling inside filter. - -&perfect_model_obs_nml - read_input_state_from_file = .true. - input_state_files = "caminput.nc" - init_time_days = -1 - init_time_seconds = -1 - - write_output_state_to_file = .true. - output_state_files = "perfect_restart.nc" - - obs_seq_in_file_name = "obs_seq.in" - obs_seq_out_file_name = "obs_seq.out" - first_obs_days = -1 - first_obs_seconds = -1 - last_obs_days = -1 - last_obs_seconds = -1 - - trace_execution = .true. - output_timestamps = .true. - print_every_nth_obs = 0 - output_forward_op_errors = .false. - / - - - -&model_nml - cam_template_filename = 'caminput.nc' - cam_phis_filename = 'cam_phis.nc' - custom_routine_to_generate_ensemble = .true. - fields_to_perturb = 'QTY_TEMPERATURE', 'QTY_U_WIND_COMPONENT', 'QTY_V_WIND_COMPONENT', 'QTY_SURFACE_PRESSURE' - perturbation_amplitude = 1.0, 1.0, 1.0, 100.0 - state_variables = 'T', 'QTY_TEMPERATURE', 'NA', 'NA', 'UPDATE' - 'US', 'QTY_U_WIND_COMPONENT', 'NA', 'NA', 'UPDATE' - 'VS', 'QTY_V_WIND_COMPONENT', 'NA', 'NA', 'UPDATE' - 'Q', 'QTY_SPECIFIC_HUMIDITY', 'NA', 'NA', 'UPDATE' - 'CLDLIQ','QTY_CLOUD_LIQUID_WATER', 'NA', 'NA', 'UPDATE' - 'CLDICE','QTY_CLOUD_ICE', 'NA', 'NA', 'UPDATE' - 'PS', 'QTY_SURFACE_PRESSURE', 'NA', 'NA', 'UPDATE' - use_log_vertical_scale = .true. - use_variable_mean_mass = .false. - no_normalization_of_scale_heights = .true. - vertical_localization_coord = 'SCALEHEIGHT' - no_obs_assim_above_level = 5 - model_damping_ends_at_level = -1 - using_chemistry = .false. - assimilation_period_days = 0 - assimilation_period_seconds = 21600 - suppress_grid_info_in_output = .false. - debug_level = 0 - / - -&location_nml - horiz_dist_only = .false. - vert_normalization_pressure = 20000.0 - vert_normalization_height = 10000.0 - vert_normalization_level = 20.0 - vert_normalization_scale_height = 1.5 - approximate_distance = .true. - nlon = 283 - nlat = 144 - output_box_info = .false. - print_box_level = 0 - special_vert_normalization_obs_types = 'null' - special_vert_normalization_pressures = -888888.0 - special_vert_normalization_heights = -888888.0 - special_vert_normalization_levels = -888888.0 - special_vert_normalization_scale_heights = -888888.0 - / - - -&fill_inflation_restart_nml - write_prior_inf = .true. - prior_inf_mean = 1.01 - prior_inf_sd = 0.6 - - write_post_inf = .false. - post_inf_mean = 1.00 - post_inf_sd = 0.6 - - input_state_files = 'caminput.nc' - single_file = .false. - - verbose = .false. - / - - -&preprocess_nml - overwrite_output = .true. - 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' - 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' - obs_type_files = '../../../observations/forward_operators/obs_def_gps_mod.f90', - '../../../observations/forward_operators/obs_def_upper_atm_mod.f90', - '../../../observations/forward_operators/obs_def_reanalysis_bufr_mod.f90', - '../../../observations/forward_operators/obs_def_altimeter_mod.f90', - '../../../observations/forward_operators/obs_def_AIRS_mod.f90' - quantity_files = '../../../assimilation_code/modules/observations/atmosphere_quantities_mod.f90', - '../../../assimilation_code/modules/observations/space_quantities_mod.f90', - '../../../assimilation_code/modules/observations/chemistry_quantities_mod.f90' - '../../../assimilation_code/modules/observations/default_quantities_mod.f90' - - -&obs_kind_nml - assimilate_these_obs_types = 'RADIOSONDE_U_WIND_COMPONENT', - 'RADIOSONDE_V_WIND_COMPONENT', - 'RADIOSONDE_TEMPERATURE', - 'AIRCRAFT_U_WIND_COMPONENT', - 'AIRCRAFT_V_WIND_COMPONENT', - 'AIRCRAFT_TEMPERATURE', - 'ACARS_U_WIND_COMPONENT', - 'ACARS_V_WIND_COMPONENT', - 'ACARS_TEMPERATURE', - 'SAT_U_WIND_COMPONENT', - 'SAT_V_WIND_COMPONENT', - 'GPSRO_REFRACTIVITY', - 'AIRS_TEMPERATURE' - 'AIRS_SPECIFIC_HUMIDITY' - 'RADIOSONDE_SPECIFIC_HUMIDITY' - evaluate_these_obs_types = - 'RADIOSONDE_SURFACE_ALTIMETER', - 'MARINE_SFC_ALTIMETER', - 'LAND_SFC_ALTIMETER' - use_precomputed_FOs_these_obs_types = 'null' - / - - -&state_vector_io_nml - buffer_state_io = .false. - single_precision_output = .false. - / - - - -&ensemble_manager_nml - layout = 2 - tasks_per_node = 36 - communication_configuration = 1 - debug = .false. - / - - -&assim_tools_nml - use_algorithm_info_mod = .true. - filter_kind = 1 - cutoff = 0.15 - sort_obs_inc = .false. - spread_restoration = .false. - sampling_error_correction = .true. - adaptive_localization_threshold = -1 - adaptive_cutoff_floor = 0.0 - output_localization_diagnostics = .false. - localization_diagnostics_file = 'localization_diagnostics' - rectangular_quadrature = .true. - gaussian_likelihood_tails = .false. - close_obs_caching = .true. - adjust_obs_impact = .false. - obs_impact_filename = "" - allow_any_impact_values = .false. - convert_all_obs_verticals_first = .true. - convert_all_state_verticals_first = .true. - special_localization_obs_types = 'null' - special_localization_cutoffs = -888888.0 - print_every_nth_obs = 10000 - distribute_mean = .false. - / - - -&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 = .true. - read_binary_file_format = 'native' - / - - -&quality_control_nml - input_qc_threshold = 3.0 - outlier_threshold = 3.0 - enable_special_outlier_code = .false. - / - - -&xyz_location_nml - / - - -&utilities_nml - TERMLEVEL = 2 - module_details = .false. - logfilename = 'dart_log.out' - nmlfilename = 'dart_log.nml' - print_debug = .false. - write_nml = 'file' - / - - -&mpi_utilities_nml - reverse_task_layout = .false. - all_tasks_print = .false. - verbose = .false. - async2_verbose = .false. - async4_verbose = .false. - shell_name = '' - separate_node_sync = .false. - create_local_comm = .true. - make_copy_before_sendrecv = .false. - / - - -&obs_def_gps_nml - max_gpsro_obs = 15000000 - / - - - - - -&obs_sequence_tool_nml - num_input_files = 2 - 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 - min_lat = -90.0 - max_lat = 90.0 - min_lon = 0.0 - max_lon = 360.0 - gregorian_cal = .true. - print_only = .false. - / - - -&obs_common_subset_nml - num_to_compare_at_once = 2 - filename_seq = '' - filename_seq_list = '' - filename_out_suffix = '.common' - print_only = .false. - print_every = 10000 - calendar = 'Gregorian' - dart_qc_threshold = 3 - eval_and_assim_can_match = .false. - / - - -&obs_impact_tool_nml - input_filename = 'cross_correlations.txt' - output_filename = 'control_impact_runtime.txt' - debug = .false. - / - - -&smoother_nml - num_lags = 0 - start_from_restart = .false. - output_restart = .false. - restart_in_file_name = 'smoother_ics' - restart_out_file_name = 'smoother_restart' - / - - - - -&obs_diag_nml - obs_sequence_name = 'obs_seq.final' - obs_sequence_list = '' - first_bin_center = BOGUS_YEAR, 1, 1, 0, 0, 0 - last_bin_center = BOGUS_YEAR, 1, 2, 0, 0, 0 - bin_separation = 0, 0, 0, 6, 0, 0 - bin_width = 0, 0, 0, 6, 0, 0 - time_to_skip = 0, 0, 1, 0, 0, 0 - max_num_bins = 1000 - trusted_obs = 'null' - plevel_edges = 1036.5, 962.5, 887.5, 775, 600, 450, 350, 275, 225, 175, 125, 75, 35, 15, 2 - hlevel_edges = 200, 630, 930, 1880,3670,5680,7440,9130,10530,12290, 14650,18220,23560,29490,43000 - Nregions = 3 - reg_names = 'Northern Hemisphere', 'Tropics', 'Southern Hemisphere' - lonlim1 = 0.0, 0.0, 0.0 - lonlim2 = 360.0, 360.0, 360.0 - latlim1 = 20.0, -20.0, -90.0 - latlim2 = 90.0, 20.0, -20.0 - print_mismatched_locs = .false. - create_rank_histogram = .true. - outliers_in_histogram = .true. - use_zero_error_obs = .false. - verbose = .false. - / - - -&schedule_nml - calendar = 'Gregorian' - first_bin_start = 1601, 1, 1, 0, 0, 0 - first_bin_end = 2999, 1, 1, 0, 0, 0 - last_bin_end = 2999, 1, 1, 0, 0, 0 - bin_interval_days = 1000000 - bin_interval_seconds = 0 - max_num_bins = 1000 - print_table = .true. - / - - -&obs_seq_to_netcdf_nml - obs_sequence_name = 'obs_seq.final' - obs_sequence_list = '' - append_to_netcdf = .false. - lonlim1 = 0.0 - lonlim2 = 360.0 - latlim1 = -90.0 - latlim2 = 90.0 - verbose = .false. - / - - -&model_mod_check_nml - input_state_files = 'caminput.nc' - output_state_files = 'mmc_output.nc' - test1thru = 0 - run_tests = 1,2,3,4,5,7 - x_ind = 175001 - - quantity_of_interest = 'QTY_U_WIND_COMPONENT' - loc_of_interest = 254.727854, 39.9768545, 50000.0 - - interp_test_lonrange = 0.0, 360.0 - interp_test_dlon = 1.0 - interp_test_latrange = -90.0, 90.0 - interp_test_dlat = 1.0 - interp_test_vertrange = 10000.0, 90000.0 - interp_test_dvert = 10000.0 - interp_test_vertcoord = 'VERTISPRESSURE' - verbose = .false. - / - - - -&closest_member_tool_nml - input_restart_file_list = 'cam_in.txt' - output_file_name = 'closest_restart' - ens_size = 80 - single_restart_file_in = .false. - difference_method = 4 - use_only_qtys = '' - / - - -&perturb_single_instance_nml - ens_size = 80 - input_files = 'caminput.nc' - output_files = 'cam_pert1.nc','cam_pert2.nc','cam_pert3.nc' - output_file_list = '' - perturbation_amplitude = 0.2 - / - - -&quad_interpolate_nml - debug = 0 - / - diff --git a/models/cam-se/work/input.nml b/models/cam-se/work/input.nml index 52dc5ca5a5..621ecffcb9 100644 --- a/models/cam-se/work/input.nml +++ b/models/cam-se/work/input.nml @@ -28,6 +28,13 @@ ! inf_sd_from_restart inflation restart files from the values in inf*_initial ! if needed. +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &filter_nml input_state_files = '' input_state_file_list = 'cam_init_files' @@ -344,7 +351,6 @@ &assim_tools_nml - filter_kind = 1 cutoff = 0.15 sort_obs_inc = .false. spread_restoration = .false. diff --git a/models/cice/work/input.nml b/models/cice/work/input.nml index 9fe096838a..6b061b9954 100644 --- a/models/cice/work/input.nml +++ b/models/cice/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true. write_output_state_to_file = .true. @@ -76,7 +83,6 @@ # cutoff of 0.03 (radians) is about 200km &assim_tools_nml - filter_kind = 1 cutoff = 0.05 sort_obs_inc = .false. spread_restoration = .false. diff --git a/models/clm/work/input.nml b/models/clm/work/input.nml index daa8fb4d98..8fc480e653 100644 --- a/models/clm/work/input.nml +++ b/models/clm/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true. write_output_state_to_file = .false. @@ -114,7 +121,6 @@ # cutoff of 0.03 (radians) is about 200km &assim_tools_nml - filter_kind = 1 cutoff = 0.05 sort_obs_inc = .false. spread_restoration = .false. diff --git a/models/cm1/work/input.nml b/models/cm1/work/input.nml index 218f34f536..465a6d1d82 100644 --- a/models/cm1/work/input.nml +++ b/models/cm1/work/input.nml @@ -2,6 +2,13 @@ ! For high-resolution models with large DART states, ! use 'distributed_state = .true.' +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &filter_nml async = 2 adv_ens_command = 'advance_model.csh' @@ -250,7 +257,6 @@ &assim_tools_nml adaptive_localization_threshold = -1 cutoff = 15000.0 - filter_kind = 1 print_every_nth_obs = 100 rectangular_quadrature = .true. sampling_error_correction = .false. diff --git a/models/coamps_nest/templates/EXPERIMENT_EXAMPLE/input.nml b/models/coamps_nest/templates/EXPERIMENT_EXAMPLE/input.nml index 5636b43894..2bc0b72bcf 100644 --- a/models/coamps_nest/templates/EXPERIMENT_EXAMPLE/input.nml +++ b/models/coamps_nest/templates/EXPERIMENT_EXAMPLE/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml start_from_restart = .true., output_restart = .true., @@ -60,7 +67,6 @@ &assim_tools_nml - filter_kind = 1, cutoff = 0.125, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/dynamo/work/input.nml b/models/dynamo/work/input.nml index c3c0921059..290276fdaf 100644 --- a/models/dynamo/work/input.nml +++ b/models/dynamo/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml start_from_restart = .false., output_restart = .true., @@ -74,7 +81,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 0.00001, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/forced_barot/work/input.nml b/models/forced_barot/work/input.nml index 2ad26c71b6..9f17e2bb1f 100644 --- a/models/forced_barot/work/input.nml +++ b/models/forced_barot/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml start_from_restart = .true., output_restart = .true., @@ -76,7 +83,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 0.02, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/forced_lorenz_96/work/input.nml b/models/forced_lorenz_96/work/input.nml index e11a177c5e..57fa73ea4a 100644 --- a/models/forced_lorenz_96/work/input.nml +++ b/models/forced_lorenz_96/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -87,7 +94,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 1000000.0, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/gitm/work/input.nml b/models/gitm/work/input.nml index f3afac2137..5e5d47f4ef 100644 --- a/models/gitm/work/input.nml +++ b/models/gitm/work/input.nml @@ -1,3 +1,9 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / &filter_nml input_state_files = '' @@ -84,7 +90,6 @@ # cutoff of 0.03 (radians) is about 200km &assim_tools_nml - filter_kind = 1, cutoff = 0.60, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/ikeda/work/input.nml b/models/ikeda/work/input.nml index 57f57246c3..2f65e7217f 100644 --- a/models/ikeda/work/input.nml +++ b/models/ikeda/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .false., single_file_in = .true. @@ -86,7 +93,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 1000000.0, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/lorenz_04/work/input.nml b/models/lorenz_04/work/input.nml index 19b372df2a..2385b45b64 100644 --- a/models/lorenz_04/work/input.nml +++ b/models/lorenz_04/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -84,7 +91,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 1000000.0, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/lorenz_63/work/input.nml b/models/lorenz_63/work/input.nml index 90504677f9..2a020e875b 100644 --- a/models/lorenz_63/work/input.nml +++ b/models/lorenz_63/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -86,7 +93,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 1000000.0 sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/lorenz_84/work/input.nml b/models/lorenz_84/work/input.nml index 6634c45d1f..2c8baa0680 100644 --- a/models/lorenz_84/work/input.nml +++ b/models/lorenz_84/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -89,7 +96,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 1000000.0, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/lorenz_96/work/input.nml b/models/lorenz_96/work/input.nml index 39bbef5f02..c1ba2dca4e 100644 --- a/models/lorenz_96/work/input.nml +++ b/models/lorenz_96/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -86,7 +93,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 0.02, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/lorenz_96_2scale/work/input.nml b/models/lorenz_96_2scale/work/input.nml index e7edc3fa92..024038aeca 100644 --- a/models/lorenz_96_2scale/work/input.nml +++ b/models/lorenz_96_2scale/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -85,7 +92,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 1000000.0, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/lorenz_96_tracer_advection/work/all_bnrhf_qceff_table.csv b/models/lorenz_96_tracer_advection/work/all_bnrhf_qceff_table.csv new file mode 100644 index 0000000000..890ecb1983 --- /dev/null +++ b/models/lorenz_96_tracer_advection/work/all_bnrhf_qceff_table.csv @@ -0,0 +1,5 @@ +QCF table version: 1,obs_error_info,,,,probit_inflation,,,,,probit_state,,,,,probit_extended_state,,,,,obs_inc_info,,,, +QTY_TEMPLATE:,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,filter_kind,bounded_below,bounded_above,lower_bound,upper_bound +QTY_STATE_VARIABLE,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RHF,.false.,.false.,-888888,-888888 +QTY_TRACER_CONCENTRATION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RHF,.true.,.false.,0,-888888 +QTY_TRACER_SOURCE,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RHF,.true.,.false.,0,-888888 \ No newline at end of file diff --git a/models/lorenz_96_tracer_advection/work/all_eakf_qceff_table.csv b/models/lorenz_96_tracer_advection/work/all_eakf_qceff_table.csv new file mode 100644 index 0000000000..5b6286816c --- /dev/null +++ b/models/lorenz_96_tracer_advection/work/all_eakf_qceff_table.csv @@ -0,0 +1,5 @@ +QCF table version: 1,obs_error_info,,,,probit_inflation,,,,,probit_state,,,,,probit_extended_state,,,,,obs_inc_info,,,, +QTY_TEMPLATE:,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,filter_kind,bounded_below,bounded_above,lower_bound,upper_bound +QTY_STATE_VARIABLE,.false.,.false.,-888888,-888888,NORMAL_DISTRIBUTION,.false.,.false.,-888888,-888888,NORMAL_DISTRIBUTION,.false.,.false.,-888888,-888888,NORMAL_DISTRIBUTION,.false.,.false.,-888888,-888888,EAKF,.false.,.false.,-888888,-888888 +QTY_TRACER_CONCENTRATION,.true.,.false.,0,-888888,NORMAL_DISTRIBUTION,.true.,.false.,0,-888888,NORMAL_DISTRIBUTION,.true.,.false.,0,-888888,NORMAL_DISTRIBUTION,.true.,.false.,0,-888888,EAKF,.true.,.false.,0,-888888 +QTY_TRACER_SOURCE,.true.,.false.,0,-888888,NORMAL_DISTRIBUTION,.true.,.false.,0,-888888,NORMAL_DISTRIBUTION,.true.,.false.,0,-888888,NORMAL_DISTRIBUTION,.true.,.false.,0,-888888,EAKF,.true.,.false.,0,-888888 \ No newline at end of file diff --git a/models/lorenz_96_tracer_advection/work/input.nml b/models/lorenz_96_tracer_advection/work/input.nml index bbbe929578..955e328a72 100644 --- a/models/lorenz_96_tracer_advection/work/input.nml +++ b/models/lorenz_96_tracer_advection/work/input.nml @@ -1,5 +1,11 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml - use_algorithm_info_mod = .true., read_input_state_from_file = .false., single_file_in = .true. input_state_files = "perfect_input.nc" @@ -29,7 +35,6 @@ / &filter_nml - use_algorithm_info_mod = .true., single_file_in = .true., input_state_files = 'perfect_input.nc' input_state_file_list = '' @@ -84,18 +89,10 @@ silence = .false., / -&probit_transform_nml - fix_bound_violations = .false., - use_logit_instead_of_probit = .false. - do_inverse_check = .true. - / - &ensemble_manager_nml / &assim_tools_nml - use_algorithm_info_mod = .true., - filter_kind = 1, cutoff = 1000000.0 sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/lorenz_96_tracer_advection/work/neg_qceff_table.csv b/models/lorenz_96_tracer_advection/work/neg_qceff_table.csv new file mode 100644 index 0000000000..01d64ad744 --- /dev/null +++ b/models/lorenz_96_tracer_advection/work/neg_qceff_table.csv @@ -0,0 +1,5 @@ +QCF table version: 1,obs_error_info,,,,probit_inflation,,,,,probit_state,,,,,probit_extended_state,,,,,obs_inc_info,,,, +QTY_TEMPLATE:,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,filter_kind,bounded_below,bounded_above,lower_bound,upper_bound +QTY_STATE_VARIABLE,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.false.,-888888,-888888,BOUNDED_NORMAL_RHF,.false.,.false.,-888888,-888888 +QTY_TRACER_CONCENTRATION,.false.,.true.,-888888,0,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.true.,-888888,0,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.true.,-888888,0,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.true.,-888888,0,BOUNDED_NORMAL_RHF,.false.,.true.,-888888,0 +QTY_TRACER_SOURCE,.false.,.true.,-888888,0,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.true.,-888888,0,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.true.,-888888,0,BOUNDED_NORMAL_RH_DISTRIBUTION,.false.,.true.,-888888,0,BOUNDED_NORMAL_RHF,.false.,.true.,-888888,0 \ No newline at end of file diff --git a/models/lorenz_96_tracer_advection/work/readme.rst b/models/lorenz_96_tracer_advection/work/readme.rst deleted file mode 100644 index 8ae7d72681..0000000000 --- a/models/lorenz_96_tracer_advection/work/readme.rst +++ /dev/null @@ -1,118 +0,0 @@ -.. _quantile tracer: - - -Quantile conserving and probit transform tools -============================================== - -This file contains instructions for using the lorenz_96_tracer model with DART -quantile conserving and probit transform filtering tools. These tools are still -being refined, but are working for the examples described. The DART development -team (dart@ucar.edu) would be happy to hear about your experiences and is -anxious to build scientific collaborations using these new capabilities. - - -Steps for reproducing basic tests: - -Test A: Assimilating observations of state (wind) and tracer concentration using -a rank histogram observation space filter and rank histogram probit transforms for -state variable updates. Example includes adaptive inflation. -The default model configuration has a single tracer source at gridpoint 1 along with -small uniform tracer sinks that lead to areas where the true tracer concentration is -usually 0. This is a particularly tough test for ensemble methods. - -#. Build all executables, - - ``./quickbuild.sh nompi`` -#. Create a set_def.out file using create_obs_sequence: - - ``./create_obs_sequence < create_obs_sequence_input`` - -#. Create an obs_sequence.in file using ``./create_fixed_network_seq`` - - .. code:: text - - ./create_fixed_network_seq - Select the default input filename , - Create a regularly repeating sequence by entering "1", - Enter "1000" for the number of observation times, - Enter "0 0" for the initial time, - Enter "0 10800" for the period, - Select the default output filename, - -#. Spin-up a model initial condition by running perfect_model_obs - - ``./perfect_model_obs`` - -#. Generate a spun-up true time series, - - .. code:: text - - cp perfect_output.nc perfect_input.nc - Use a text editor to change read_input_state_from_file to .true. in the file input.nml - Run "./perfect_model_obs" again - -#. Run a filter assimilation, - - ``./filter`` - -#. Examine the output with your favorite tools. Looking at the analysis ensemble - for the tracer_concentration variables with indices near the source (location 1) - and far downstream from the source (location 35) is interesting. Note that the - source estimation capabilities of the model and filters are not being tested here. - - -Test B: Using default ensemble adjustment Kalman filters. - -The new quantile options are controlled by Fortran code in the module -algorithm_info_mod.f90 in the assimilation_code/modules/assimilation directory. -More information about the control can be found in that module. The tests below -replace the default version of that module with others that change certain options. -Doing a diff between these modules shows how the control is being changed for the -following tests in that module. The tests below -replace the default version of that module with others that change certain options. - -#. In directory assimilation_code/modules/assimilation, - - ``cp all_eakf_algorithm_info_mod algorithm_info_mod.f90`` - -#. Recompile all programs in this directory, - - ``./quickbiuld.sh nompi`` - -#. Run the filter - ``./filter`` - -Test C: Using default ensemble adjustment Kalman filter for state, but bounded normal rank histogram filter and priors for tracer concentration and source. - -#. In directory assimilation_code/modules/assimilation, - - ``cp state_eakf_tracer_bnrhf_algorithm_info_mod algorithm_info_mod.f90`` - -#. Recompile all programs in this directory, - - ``./quickbiuld.sh nompi`` - -#. Run the filter - ``./filter`` - -Test D: Testing bounded above option - -Normally tracers are bounded below, but there are other quantities that may be bounded -above. There are distinct numerical challenges in implementing the quantile algorithms -for quantities that are bounded above, so flipping the sign of the tracers is a good -test. - -#. In directory assimilation_code/modules/assimilation, - - ``cp neg_algorithm_info_mod algorithm_info_mod.f90`` - -#. Recompile all programs in this directory, - - ``./quickbiuld.sh nompi`` - -#. In the file input.nml, change the entry positive_tracer to .false. Also, change the - entry read_input_state_from_file back to .false. - -#. Repeat steps 3-6 from Test A. - - diff --git a/models/lorenz_96_tracer_advection/work/state_eakf_tracer_bnrhf_qceff_table.csv b/models/lorenz_96_tracer_advection/work/state_eakf_tracer_bnrhf_qceff_table.csv new file mode 100644 index 0000000000..3ae9561edd --- /dev/null +++ b/models/lorenz_96_tracer_advection/work/state_eakf_tracer_bnrhf_qceff_table.csv @@ -0,0 +1,5 @@ +QCF table version: 1,obs_error_info,,,,probit_inflation,,,,,probit_state,,,,,probit_extended_state,,,,,obs_inc_info,,,, +QTY_NAME:,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,dist_type,bounded_below,bounded_above,lower_bound,upper_bound,filter_kind,bounded_below,bounded_above,lower_bound,upper_bound +QTY_STATE_VARIABLE,.false.,.false.,-888888,-888888,NORMAL_DISTRIBUTION,.false.,.false.,-888888,-888888,NORMAL_DISTRIBUTION,.false.,.false.,-888888,-888888,NORMAL_DISTRIBUTION,.false.,.false.,-888888,-888888,EAKF,.false.,.false.,-888888,-888888 +QTY_TRACER_CONCENTRATION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RHF,.true.,.false.,0,-888888 +QTY_TRACER_SOURCE,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RH_DISTRIBUTION,.true.,.false.,0,-888888,BOUNDED_NORMAL_RHF,.true.,.false.,0,-888888 \ No newline at end of file diff --git a/models/mpas_atm/work/input.nml b/models/mpas_atm/work/input.nml index 67869b6174..10a2171162 100644 --- a/models/mpas_atm/work/input.nml +++ b/models/mpas_atm/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true. single_file_in = .false. @@ -92,7 +99,6 @@ / &assim_tools_nml - filter_kind = 1 cutoff = 0.10 distribute_mean = .false. convert_all_obs_verticals_first = .true. diff --git a/models/mpas_ocn/work/input.nml b/models/mpas_ocn/work/input.nml index a3b35c0cfb..1aa03fe07f 100644 --- a/models/mpas_ocn/work/input.nml +++ b/models/mpas_ocn/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml start_from_restart = .true., output_restart = .true., @@ -74,7 +81,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 1000000.0, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/noah/work/input.nml b/models/noah/work/input.nml index 86f921a95a..0f52583682 100644 --- a/models/noah/work/input.nml +++ b/models/noah/work/input.nml @@ -1,5 +1,12 @@ # This namelist is for both NOAH and NOAH-MP +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &model_nml lsm_model_choice = 'noahMP_36' domain_shapefiles = 'RESTART.2003051600_DOMAIN1_01' @@ -115,7 +122,6 @@ # cutoff = 0.06 (radians) is about 400 km &assim_tools_nml - filter_kind = 1 cutoff = 0.05 sort_obs_inc = .false. spread_restoration = .false. diff --git a/models/null_model/work/input.nml b/models/null_model/work/input.nml index c3e8905428..cbb0e06326 100644 --- a/models/null_model/work/input.nml +++ b/models/null_model/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -97,7 +104,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 1000000.0, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/pe2lyr/work/input.nml b/models/pe2lyr/work/input.nml index ddfc37e57b..e25edd4337 100644 --- a/models/pe2lyr/work/input.nml +++ b/models/pe2lyr/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml start_from_restart = .true., output_restart = .true., @@ -70,7 +77,6 @@ perturbation_amplitude = 0.2 / &assim_tools_nml - filter_kind = 1, cutoff = 0.02, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/rose/work/input.nml b/models/rose/work/input.nml index 786486552b..8f52f2fb5e 100644 --- a/models/rose/work/input.nml +++ b/models/rose/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml start_from_restart = .false., output_restart = .true., @@ -73,7 +80,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 0.2, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/simple_advection/work/input.nml b/models/simple_advection/work/input.nml index 13f84b1220..cd436c8f4a 100644 --- a/models/simple_advection/work/input.nml +++ b/models/simple_advection/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true. single_file_in = .true. @@ -91,7 +98,6 @@ # large for the model to converge. to test that the model is # doing a successful assimilation, change cutoff to 0.02 and rerun. &assim_tools_nml - filter_kind = 1 cutoff = 100000000.0 sort_obs_inc = .false. spread_restoration = .false. diff --git a/models/sqg/work/input.nml b/models/sqg/work/input.nml index 62b0acaee3..5981b6ca67 100644 --- a/models/sqg/work/input.nml +++ b/models/sqg/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml start_from_restart = .true., output_restart = .true., @@ -79,7 +86,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 100000.0, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/template/work/oned_input.nml b/models/template/work/oned_input.nml index 66b5a07bb8..65ee6dcb9d 100644 --- a/models/template/work/oned_input.nml +++ b/models/template/work/oned_input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -86,7 +93,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 0.02, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/template/work/threed_input.nml b/models/template/work/threed_input.nml index 90608b10f6..8b50853c78 100644 --- a/models/template/work/threed_input.nml +++ b/models/template/work/threed_input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true., single_file_in = .true. @@ -86,7 +93,6 @@ / &assim_tools_nml - filter_kind = 1, cutoff = 1000000.0 sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/tiegcm/work/input.nml b/models/tiegcm/work/input.nml index c92d64e746..7420562fa9 100644 --- a/models/tiegcm/work/input.nml +++ b/models/tiegcm/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &quality_control_nml / @@ -104,7 +111,6 @@ / &assim_tools_nml - filter_kind = 1 cutoff = 0.2 sort_obs_inc = .false. spread_restoration = .false. diff --git a/models/wrf/work/input.nml b/models/wrf/work/input.nml index 30b6f0cad0..aec2653ea6 100644 --- a/models/wrf/work/input.nml +++ b/models/wrf/work/input.nml @@ -1,3 +1,10 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &perfect_model_obs_nml read_input_state_from_file = .true. single_file_in = .false. @@ -96,7 +103,6 @@ # localization radius, where the influence of an obs decreases # to ~half at 300 km, and ~0 at the edges of the area. &assim_tools_nml - filter_kind = 1, cutoff = 0.05, sort_obs_inc = .false., spread_restoration = .false., diff --git a/models/wrf_hydro/work/input.nml b/models/wrf_hydro/work/input.nml index 8b22e40e68..7ab439a0ca 100644 --- a/models/wrf_hydro/work/input.nml +++ b/models/wrf_hydro/work/input.nml @@ -7,6 +7,13 @@ # domain_order = 'hydro', 'parameters' # domain_shapefiles = 'restart.hydro.nc', 'parameters.nc' +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + &model_nml assimilation_period_days = 0 assimilation_period_seconds = 3600 @@ -122,7 +129,6 @@ / -# filter_kind=9 => JPoterjoy particle filter # cutoff is in radians; for the earth, 0.05 is about 300 km. # cutoff is defined to be the half-width of the localization radius # so 0.05 radians for cutoff is about an 600 km effective @@ -134,7 +140,6 @@ # max_link_distance = 2000m = cutoff = 0.000157 radians &assim_tools_nml - filter_kind = 1 cutoff = 0.000157 sort_obs_inc = .false. spread_restoration = .false.