diff --git a/CHANGELOG.rst b/CHANGELOG.rst index dd4d2d0c6f..84bd662946 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,6 +22,10 @@ individual files. The changes are now listed with the most recent at the top. +**March 11 2024 :: SEIR model for infectious diseases. Tag v11.1.0** + +- Added SEIR model which simulates the spread of infectious diseases, for example COVID-19. + **February 13 2024 :: Fortran Standards. Tag v11.0.3** - Replace f2kcli with Fortran intrinsics for command line arguments. diff --git a/conf.py b/conf.py index 004d95ec76..f55af3ec11 100644 --- a/conf.py +++ b/conf.py @@ -21,7 +21,7 @@ author = 'Data Assimilation Research Section' # The full version, including alpha/beta/rc tags -release = '11.0.3' +release = '11.1.0' root_doc = 'index' # -- General configuration --------------------------------------------------- diff --git a/index.rst b/index.rst index 417ee8bee3..f04390488f 100644 --- a/index.rst +++ b/index.rst @@ -480,6 +480,7 @@ References models/POP/dart_pop_mod models/ROMS/readme models/rose/readme + models/seir/readme models/simple_advection/readme models/sqg/readme models/template/new_model diff --git a/models/README.rst b/models/README.rst index 09eb9bf5db..73b4dd87bd 100644 --- a/models/README.rst +++ b/models/README.rst @@ -39,6 +39,7 @@ DART supported models: - :doc:`POP/readme` - :doc:`ROMS/readme` - :doc:`rose/readme` +- :doc:`seir/readme` - :doc:`simple_advection/readme` - :doc:`sqg/readme` - :doc:`tiegcm/readme` diff --git a/models/seir/model_mod.f90 b/models/seir/model_mod.f90 new file mode 100644 index 0000000000..d71657c17c --- /dev/null +++ b/models/seir/model_mod.f90 @@ -0,0 +1,439 @@ +! 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 model_mod + +use types_mod, only : r8, i8, i4, MISSING_R8 + +use time_manager_mod, only : time_type, set_time + +use location_mod, only : location_type, set_location, get_location, & + get_close_obs, get_close_state, & + convert_vertical_obs, convert_vertical_state + +use utilities_mod, only : register_module, do_nml_file, do_nml_term, & + nmlfileunit, find_namelist_in_file, & + check_namelist_read + +use location_io_mod, only : nc_write_location_atts, nc_write_location + +use netcdf_utilities_mod, only : nc_add_global_attribute, nc_synchronize_file, & + nc_add_global_creation_time, nc_begin_define_mode, & + nc_end_define_mode + +use obs_kind_mod, only : QTY_STATE_VARIABLE + +use mpi_utilities_mod, only : my_task_id + +use random_seq_mod, only : random_seq_type, init_random_seq, random_gaussian + +use ensemble_manager_mod, only : ensemble_type, get_my_num_vars, get_my_vars + +use distributed_state_mod, only : get_state + +use state_structure_mod, only : add_domain + +use default_model_mod, only : end_model, nc_write_model_vars + +use dart_time_io_mod, only : read_model_time, write_model_time + +implicit none +private + +public :: get_model_size, & + get_state_meta_data, & + model_interpolate, & + shortest_time_between_assimilations, & + static_init_model, & + init_conditions, & + init_time, & + adv_1step, & + nc_write_model_atts + +public :: pert_model_copies, & + nc_write_model_vars, & + get_close_obs, & + get_close_state, & + end_model, & + convert_vertical_obs, & + convert_vertical_state, & + read_model_time, & + write_model_time + +character(len=*), parameter :: source = "seir/model_mod.f90" + +type(location_type) :: state_loc ! state location, compute once and store for speed +type(random_seq_type) :: random_seq +type(time_type) :: time_step + +! Input parameters +integer(i8) :: model_size = 7 +integer :: time_step_days = 0 +integer :: time_step_seconds = 3600 +real(r8) :: t_incub = 5.6_r8 ! Incubation period (days) +real(r8) :: t_infec = 3.8_r8 ! Infection time (days) +real(r8) :: t_recov = 14.0_r8 ! Recovery period (days) +real(r8) :: t_death = 7.0_r8 ! Time until death (days) +real(r8) :: alpha = 0.007_r8 ! Vaccination rate (per day) +integer(i8) :: theta = 12467 ! New birth and new residents (persons per day) +real(r8) :: mu = 0.000025_r8 ! Natural death rate (persons per day) +real(r8) :: sigma = 0.05_r8 ! Vaccination inefficacy (e.g., 95% efficient) +real(r8) :: beta = 1.36e-9_r8 ! Transmission rate (per day) +real(r8) :: kappa = 0.00308 ! Mortality rate +real(r8) :: delta_t = 0.04167_r8 ! Model time step; 1/24 := 1 hour +integer(i8) :: num_pop = 331996199 ! Population (US) +real(r8) :: pert_size = 1.0 ! Size of perturbation (lognormal pdf param) + +real(r8) :: gama, delta, lambda, rho + +! Other related model parameters +real(r8), parameter :: E0 = 1.0_r8 ! Exposed (not yet infected) +real(r8), parameter :: I0 = 1.0_r8 ! Infected (not yet quarantined) +real(r8), parameter :: Q0 = 1.0_r8 ! Quarantined (confirmed and infected) +real(r8), parameter :: R0 = 1.0_r8 ! Recovered +real(r8), parameter :: D0 = 1.0_r8 ! Dead +real(r8), parameter :: V0 = 1.0_r8 ! Vaccinated +real(r8) :: S0 ! Susceptible + +namelist /model_nml/ model_size, time_step_days, time_step_seconds, & + delta_t, num_pop, pert_size, & + t_incub, t_infec, t_recov, t_death, & + alpha, theta, beta, sigma, kappa, mu + +contains + + +!------------------------------------------------------------------ +! + +subroutine static_init_model() + +real(r8) :: x_loc +integer :: i, dom_id + +! Do any initial setup needed, including reading the namelist values +call initialize() + +! Define the locations of the model state variables +! The SEIR variables have no physical location, +! and so I'm placing all 7 variables at the same +! virtual point in space. +x_loc = 0.5_r8 +state_loc = set_location(x_loc) + +! This time is both the minimum time you can ask the model to advance +! and it sets the assimilation window. +! All observations within +/- 1/2 this interval from the current +! model time will be assimilated. +time_step = set_time(time_step_seconds, time_step_days) + +! Tell the DART I/O routines how large the model data is so they +! can read/write it. +dom_id = add_domain(model_size) + +end subroutine static_init_model + +!------------------------------------------------------------------ +! Advance the SEIR model using a four-step RK scheme + +subroutine adv_1step(x, time) + +real(r8), intent(inout) :: x(:) +type(time_type), intent(in) :: time + +real(r8), dimension(size(x)) :: xi, x1, x2, x3, x4, dx + +! 1st step +call seir_eqns(x, dx) +x1 = delta_t * dx + +! 2nd step +xi = x + 0.5_r8 * delta_t * dx +call seir_eqns(xi, dx) +x2 = delta_t * dx + +! 3rd step +xi = x + 0.5_r8 * delta_t * dx +call seir_eqns(xi, dx) +x3 = delta_t * dx + +! 4th step +xi = x + delta_t * dx +call seir_eqns(xi, dx) +x4 = delta_t * dx + +! Compute new value for x +x = x + x1/6.0_r8 + x2/3.0_r8 + x3/3.0_r8 + x4/6.0_r8 + +end subroutine adv_1step + +!------------------------------------------------------------------ +! SEIR Model Equations +! The following extended SEIR Model with Vaccination +! is adaopted from Ghostine et al. (2021): +! Ghostine, R., Gharamti, M., Hassrouny, S and Hoteit, I +! "An Extended SEIR Model with Vaccination for Forecasting +! the COVID-19 Pandemic in Saudi Arabia Using an Ensemble +! Kalman Filter" Mathematics 2021, 9, 636. +! https://dx.doi.org/10.3390/math9060636 + +subroutine seir_eqns(x, fx) + +! State: x = [S, E, I, Q, R, D, V] + +real(r8), intent(in) :: x(:) +real(r8), intent(out) :: fx(:) + +fx(1) = theta - beta * x(1) * x(3) - alpha * x(1) - mu * x(1) +fx(2) = beta * x(1) * x(3) - gama * x(2) + sigma * beta * x(7) * x(3) - mu * x(2) +fx(3) = gama * x(2) - delta * x(3) - mu * x(3) +fx(4) = delta * x(3) - (1.0_r8 - kappa) * lambda * x(4) - kappa * rho * x(4) - mu * x(4) +fx(5) = (1.0_r8 - kappa) * lambda * x(4) - mu * x(5) +fx(6) = kappa * rho * x(4) +fx(7) = alpha * x(1) - sigma * beta * x(7) * x(3) - mu * x(7) + +end subroutine seir_eqns + +!------------------------------------------------------------------ +! Perturbs a model state for generating initial ensembles. +! Returning interf_provided .true. means this code has +! added uniform small independent perturbations to a +! single ensemble member to generate the full ensemble. +subroutine pert_model_copies(state_ens_handle, ens_size, pert_amp, interf_provided) + +type(ensemble_type), intent(inout) :: state_ens_handle +integer, intent(in) :: ens_size +real(r8), intent(in) :: pert_amp +logical, intent(out) :: interf_provided + +integer :: i,j, num_my_grid_points +real(r8) :: rng + +interf_provided = .true. + +call init_random_seq(random_seq, my_task_id()+1) +! if we are running with more than 1 task, then +! we have all the ensemble members for a subset of +! the model state. which variables we have are determined +! by looking at the global index number into the state vector. + +num_my_grid_points = get_my_num_vars(state_ens_handle) + +do i=1,num_my_grid_points + + ! Lognormal Distribution + do j= 1, ens_size + rng = pert_size * random_gaussian(random_seq, 0.0_r8, 1.0_r8) + state_ens_handle%copies(j, i) = & + state_ens_handle%copies(j, i) * exp(rng) + end do +end do + + +end subroutine pert_model_copies + +!------------------------------------------------------------------ +! Returns a model state vector, x, that is some sort of appropriate +! initial condition for starting up a long integration of the model. +! At present, this is only used if the namelist parameter +! start_from_restart is set to .false. in the program perfect_model_obs. +! If this option is not to be used in perfect_model_obs, or if no +! synthetic data experiments using perfect_model_obs are planned, +! this can be a NULL INTERFACE. + +subroutine init_conditions(x) + +real(r8), dimension (model_size) :: x0 +real(r8), intent(out) :: x(:) + +x0 = (/S0, E0, I0, Q0, R0, D0, V0/) +x = x0 + +end subroutine init_conditions + +!------------------------------------------------------------------ +! Companion interface to init_conditions. Returns a time that is somehow +! appropriate for starting up a long integration of the model. +! At present, this is only used if the namelist parameter +! start_from_restart is set to .false. in the program perfect_model_obs. +! If this option is not to be used in perfect_model_obs, or if no +! synthetic data experiments using perfect_model_obs are planned, +! this can be a NULL INTERFACE. + +subroutine init_time(time) + +type(time_type), intent(out) :: time + +! for now, just set to 0 +time = set_time(0,0) + +end subroutine init_time + + +!------------------------------------------------------------------ +! Returns the number of items in the state vector as an integer. +! This interface is required for all applications. + +function get_model_size() + +integer(i8) :: get_model_size + +get_model_size = model_size + +end function get_model_size + + + +!------------------------------------------------------------------ +! Returns the smallest increment in time that the model is capable +! of advancing the state in a given implementation, or the shortest +! time you want the model to advance between assimilations. +! This interface is required for all applications. + +function shortest_time_between_assimilations() + +type(time_type) :: shortest_time_between_assimilations + +shortest_time_between_assimilations = time_step + +end function shortest_time_between_assimilations + + + +!------------------------------------------------------------------ +! Given a state handle, a location, and a model state variable quantity, +! interpolates the state variable fields to that location and returns +! the values in expected_obs. The istatus variables should be returned as +! 0 unless there is some problem in computing the interpolation in +! which case an alternate value should be returned. The itype variable +! is an integer that specifies the quantity of field (for +! instance temperature, zonal wind component, etc.). In low order +! models that have no notion of types of variables this argument can +! be ignored. For applications in which only perfect model experiments +! with identity observations (i.e. only the value of a particular +! state variable is observed), this can be a NULL INTERFACE. + +subroutine model_interpolate(state_handle, ens_size, location, iqty, expected_obs, istatus) + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: ens_size +type(location_type), intent(in) :: location +integer, intent(in) :: iqty +real(r8), intent(out) :: expected_obs(ens_size) !< array of interpolated values +integer, intent(out) :: istatus(ens_size) + +! Given the nature of the SEIR model, no interpolation +! is needed. Only identity observations are utilized. + +! This should be the result of the interpolation of a +! given quantity (iqty) of variable at the given location. +expected_obs(:) = MISSING_R8 + +! The return code for successful return should be 0. +! Any positive number is an error. +! Negative values are reserved for use by the DART framework. +! Using distinct positive values for different types of errors can be +! useful in diagnosing problems. +istatus(:) = 1 + +end subroutine model_interpolate + + + +!------------------------------------------------------------------ +! Given an integer index into the state vector structure, returns the +! associated location. A second intent(out) optional argument quantity +! (qty) can be returned if the model has more than one type of field +! (for instance temperature and zonal wind component). This interface is +! required for all filter applications as it is required for computing +! the distance between observations and state variables. + +subroutine get_state_meta_data(index_in, location, qty) + +integer(i8), intent(in) :: index_in +type(location_type), intent(out) :: location +integer, intent(out), optional :: qty + +! these should be set to the actual location and state quantity +location = state_loc !(index_in) +if (present(qty)) qty = QTY_STATE_VARIABLE + +end subroutine get_state_meta_data + + + +!------------------------------------------------------------------ +! Do any initialization/setup, including reading the +! namelist values. + +subroutine initialize() + +integer :: iunit, io + +! Print module information +call register_module(source) + +! Read the namelist +call find_namelist_in_file("input.nml", "model_nml", iunit) +read(iunit, nml = model_nml, iostat = io) +call check_namelist_read(iunit, io, "model_nml") + +! Output the namelist values if requested +if (do_nml_file()) write(nmlfileunit, nml=model_nml) +if (do_nml_term()) write( * , nml=model_nml) + +! Compute other model parameters +gama = 1.0_r8 / t_incub +delta = 1.0_r8 / t_infec +lambda = 1.0_r8 / t_recov +rho = 1.0_r8 / t_death + +! Compute initial value for S: +S0 = num_pop - E0 - I0 - Q0 - R0 - D0 + +end subroutine initialize + + +!------------------------------------------------------------------ +! Writes model-specific attributes to a netCDF file + +subroutine nc_write_model_atts(ncid, domain_id) + +integer, intent(in) :: ncid +integer, intent(in) :: domain_id +integer(i4) :: offset + +! put file into define mode. + +integer :: msize + +msize = int(model_size, i4) + +call nc_begin_define_mode(ncid) + +call nc_add_global_creation_time(ncid) + +call nc_add_global_attribute(ncid, "model_source", source) + +call nc_add_global_attribute(ncid, "model", "seir") + +call nc_write_location_atts(ncid, msize) +call nc_end_define_mode(ncid) + +! Note that location for this model isn't used +do offset = 1, msize + call nc_write_location(ncid, state_loc, offset) +enddo + +! Flush the buffer and leave netCDF file open +call nc_synchronize_file(ncid) + +end subroutine nc_write_model_atts + +!=================================================================== +! End of model_mod +!=================================================================== +end module model_mod + diff --git a/models/seir/readme.rst b/models/seir/readme.rst new file mode 100644 index 0000000000..54bb894f3b --- /dev/null +++ b/models/seir/readme.rst @@ -0,0 +1,138 @@ +SEIR +==== + +Overview +-------- + +The extended SEIR Model with Vaccination was first proposed by Ghostine et al. (2021) [1]_ +to simulate the novel coronavirus disease (COVID-19) spread. The model considers 7 +stages of infection: + + 1. Susceptible (S), + 2. Exposed (E), + 3. Infected (I), + 4. Quarantined (Q), + 5. Recovered (R), + 6. Deaths (D), + 7. Vaccinated (V). + +There are several parameters that can be changed to study different cases and regions: + + - :math:`\theta`: New births and new residents per unit of time, + - :math:`\beta`: Transmission rate divided by the population size, + - :math:`\alpha`: Vaccination rate, + - :math:`\mu`: Natural death rate, + - :math:`\gamma`: Average latent time, + - :math:`\delta`: Average quarantine time, + - :math:`\kappa`: Mortality rate, + - :math:`\lambda`: Average days until recovery, + - :math:`\rho`: Average days until death, + - :math:`\sigma`: Vaccine in-efficacy (:math:`0 \leq \sigma \leq 1`). + +Earth system models are often descritized in space. The state in these models represents +variables at different spatial locations. The variables of the SEIR model describe the +stage/phase of the disease and they do not have a physical location. To this end, +techniques such as spatial localization are not applicable in this model. DART assumes +that all 7 variables belong to the same *virtual* point in space. Any assimilated +observation will impact all 7 variables. + +The SEIR model uses identity observations. Typical observations that can be assimilated +are: + + *Recovered*, *Death* and *Vaccinated*. + +Some agencies provide data for "*Confirmed*" cases. This can be used to compute and +assimilate the number of active (which is equivelant to quarantined) cases as shown: + +*Active/Quarantined (Q) = Confirmed - Recovered (R) - Deaths (D)* + +Initial versions of the model were tested using DART_LAB. This was conducted by +**Shaniah Reece** as part of her SIParCS internship at NSF NCAR (2022). + +Namelist +-------- + +The ``&model_nml`` namelist is read from the ``input.nml`` file. Namelists +start with an ampersand ``&`` and terminate with a slash ``/``. Character +strings that contain a ``/`` must be enclosed in quotes to prevent them from +prematurely terminating the namelist. + +.. code-block:: fortran + + &model_nml + model_size = 40, + delta_t = 0.04167, + time_step_days = 0, + time_step_seconds = 3600, + num_pop = 331996199, + pert_size = 0.5, + t_incub = 5.6, + t_infec = 3.8, + t_recov = 14.0, + t_death = 7.0, + alpha = 0.000001, + theta = 12467, + mu = 0.000025, + sigma = 0.05, + beta = 0.00000000136, + kappa = 0.00308, + / + +Description of each namelist entry +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ++-------------------+----------+-------------------------------------------+ +| Item | Type | Description | ++===================+==========+===========================================+ +| model_size | integer | Number of variables in model. | ++-------------------+----------+-------------------------------------------+ +| delta_t | real(r8) | Non-dimensional timestep. This is | +| | | mapped to the dimensional timestep | +| | | specified by time_step_days and | +| | | time_step_seconds. | ++-------------------+----------+-------------------------------------------+ +| time_step_days | integer | Number of days for dimensional | +| | | timestep, mapped to delta_t. | ++-------------------+----------+-------------------------------------------+ +| time_step_seconds | integer | Number of seconds for dimensional | +| | | timestep, mapped to delta_t. | ++-------------------+----------+-------------------------------------------+ +| num_pop | integer | Population size. | ++-------------------+----------+-------------------------------------------+ +| pert_size | real(r8) | Size of perturbation used to create | +| | | an ensemble using a lognormal pdf. | ++-------------------+----------+-------------------------------------------+ +| t_incub | real(r8) | Incubation period | +| | | :math:`\equiv 1/\gamma`. | ++-------------------+----------+-------------------------------------------+ +| t_infec | real(r8) | Infection time | +| | | :math:`\equiv 1/\delta`. | ++-------------------+----------+-------------------------------------------+ +| t_recov | real(r8) | Recovery period | +| | | :math:`\equiv 1/\lambda`. | ++-------------------+----------+-------------------------------------------+ +| t_death | real(r8) | Time until death | +| | | :math:`\equiv 1/\rho`. | ++-------------------+----------+-------------------------------------------+ +| alpha | real(r8) | Vaccination rate. If study period | +| | | starts before vaccination is | +| | | available, this must be set to 0. | ++-------------------+----------+-------------------------------------------+ +| theta | integer | New birth and new residents. | ++-------------------+----------+-------------------------------------------+ +| mu | real(r8) | Natural death rate. | ++-------------------+----------+-------------------------------------------+ +| sigma | real(r8) | Vaccination inefficacy (e.g., if the | +| | | vaccine is 95% effective, then | +| | | :math:`\sigma = 1-0.95 = 0.05`). | ++-------------------+----------+-------------------------------------------+ +| beta | real(r8) | Transmission rate divided by population | +| | | size. | ++-------------------+----------+-------------------------------------------+ +| kappa | real(r8) | Mortality rate. | ++-------------------+----------+-------------------------------------------+ + +References +---------- + +.. [1] Ghostine, R.; Gharamti, M.; Hassrouny, S.; Hoteit, I. An Extended SEIR Model with Vaccination for Forecasting the COVID-19 Pandemic in Saudi Arabia Using an Ensemble Kalman Filter. Mathematics 2021, 9, 636. https://dx.doi.org/10.3390/math9060636. diff --git a/models/seir/work/input.nml b/models/seir/work/input.nml new file mode 100644 index 0000000000..f5681a265d --- /dev/null +++ b/models/seir/work/input.nml @@ -0,0 +1,245 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + +&perfect_model_obs_nml + read_input_state_from_file = .true., + single_file_in = .true. + input_state_files = "perfect_input.nc" + + write_output_state_to_file = .true., + single_file_out = .true. + output_state_files = "perfect_output.nc" + output_interval = 1, + + async = 0, + adv_ens_command = "./advance_model.csh", + + obs_seq_in_file_name = "obs_seq.in", + obs_seq_out_file_name = "obs_seq.out", + init_time_days = 0, + init_time_seconds = 0, + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + + trace_execution = .false., + output_timestamps = .false., + print_every_nth_obs = -1, + output_forward_op_errors = .false., + silence = .false., + / + +&filter_nml + single_file_in = .true., + input_state_files = '' + input_state_file_list = 'filter_input_list.txt' + + stages_to_write = 'preassim', 'analysis', 'output' + + single_file_out = .true., + output_state_files = '' + output_state_file_list = 'filter_output_list.txt' + output_interval = 1, + output_members = .true. + num_output_state_members = 20, + output_mean = .true. + output_sd = .true. + + ens_size = 3, + num_groups = 1, + perturb_from_single_instance = .true., + perturbation_amplitude = 0.2, + distributed_state = .true. + + async = 0, + adv_ens_command = "./advance_model.csh", + + obs_sequence_in_name = "obs_seq.out", + obs_sequence_out_name = "obs_seq.final", + num_output_obs_members = 20, + init_time_days = 0, + init_time_seconds = 0, + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + + inf_flavor = 5, 0, + inf_initial_from_restart = .false., .false., + inf_sd_initial_from_restart = .false., .false., + inf_deterministic = .true., .true., + inf_initial = 1.0, 1.0, + inf_lower_bound = 0.0, 1.0, + inf_upper_bound = 1000000.0, 1000000.0, + inf_damping = 0.9, 1.0, + inf_sd_initial = 0.6, 0.0, + inf_sd_lower_bound = 0.6, 0.0, + inf_sd_max_change = 1.05, 1.05, + + trace_execution = .false., + output_timestamps = .false., + output_forward_op_errors = .false., + write_obs_every_cycle = .false., + silence = .false., + / + + +&ensemble_manager_nml + / + +&assim_tools_nml + cutoff = 0.02, + sort_obs_inc = .false., + spread_restoration = .false., + sampling_error_correction = .false., + adaptive_localization_threshold = -1, + output_localization_diagnostics = .false., + localization_diagnostics_file = 'localization_diagnostics', + print_every_nth_obs = 0, + rectangular_quadrature = .true., + gaussian_likelihood_tails = .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 = .false., + read_binary_file_format = 'native' + / + +&obs_kind_nml + assimilate_these_obs_types = 'RAW_STATE_VARIABLE' + / + +&model_nml + model_size = 7, + delta_t = 0.04167, + time_step_days = 0, + time_step_seconds = 3600, + num_pop = 331996199, + pert_size = 0.5, + t_incub = 5.6, + t_infec = 3.8, + t_recov = 14.0, + t_death = 7.0, + alpha = 0.000001, + theta = 12467, + mu = 0.000025, + sigma = 0.05, + beta = 0.00000000136, + kappa = 0.00308, + / + +&utilities_nml + termlevel = 1, + module_details = .false., + logfilename = 'dart_log.out', + nmlfilename = 'dart_log.nml', + write_nml = 'file', + print_debug = .false., + / + +&mpi_utilities_nml + / + +&preprocess_nml + overwrite_output = .true. + input_obs_def_mod_file = '../../../observations/forward_operators/DEFAULT_obs_def_mod.F90' + output_obs_def_mod_file = '../../../observations/forward_operators/obs_def_mod.f90' + input_obs_qty_mod_file = '../../../assimilation_code/modules/observations/DEFAULT_obs_kind_mod.F90' + output_obs_qty_mod_file = '../../../assimilation_code/modules/observations/obs_kind_mod.f90' + obs_type_files = '../../../observations/forward_operators/obs_def_1d_state_mod.f90' + quantity_files = '../../../assimilation_code/modules/observations/oned_quantities_mod.f90' + / + +&obs_sequence_tool_nml + filename_seq = 'obs1.out', 'obs2.out', + filename_seq_list = '', + filename_out = 'obs_seq.combined', + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + print_only = .false., + gregorian_cal = .false., + / + +&obs_diag_nml + obs_sequence_name = 'obs_seq.final', + bin_width_days = -1, + bin_width_seconds = -1, + init_skip_days = 0, + init_skip_seconds = 0, + Nregions = 3, + trusted_obs = 'null', + lonlim1 = 0.00, 0.00, 0.50 + lonlim2 = 1.01, 0.50, 1.01 + reg_names = 'whole', 'lower', 'upper' + 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 = 1.0, + verbose = .true. + / + +&state_vector_io_nml + / + +&quality_control_nml + input_qc_threshold = 3.0, + outlier_threshold = -1.0, + / + +&integrate_model_nml + trace_execution = .true. + ic_file_name = 'temp_ic.nc' + ud_file_name = 'temp_uc.nc' + / + +&model_mod_check_nml + input_state_files = 'perfect_input.nc' + output_state_files = 'mmc_output.nc' + num_ens = 1 + single_file = .false. + test1thru = 7 + run_tests = 1,2,3,4,5,6 + x_ind = 4 + loc_of_interest = 0.3 + quantity_of_interest = 'QTY_STATE_VARIABLE' + interp_test_dx = 0.02 + interp_test_xrange = 0.0, 1.0 + verbose = .true. + / diff --git a/models/seir/work/obs_seq.out b/models/seir/work/obs_seq.out new file mode 100644 index 0000000000..7c6ec7e6ae --- /dev/null +++ b/models/seir/work/obs_seq.out @@ -0,0 +1,189 @@ + obs_sequence +obs_type_definitions + 0 + num_copies: 2 num_qc: 1 + num_obs: 15 max_num_obs: 15 +observations +truth +Quality Control + first: 1 last: 15 + OBS 1 + 1.5550552250651206 + 0.57764969235130070 + 0.0000000000000000 + -1 2 -1 +obdef +loc1d + 0.5000000000000000 +kind + -5 + 0 1 + 0.50000000000000000 + OBS 2 + 0.90066834931003292 + 0.30047988296891687 + 0.0000000000000000 + 1 3 -1 +obdef +loc1d + 0.5000000000000000 +kind + -6 + 0 1 + 0.29999999999999999 + OBS 3 + 332.04521140436549 + 332.12052020582382 + 0.0000000000000000 + 2 4 -1 +obdef +loc1d + 0.5000000000000000 +kind + -7 + 0 1 + 1.0000000000000000E-002 + OBS 4 + 1.2957274103345058 + 0.66709604449775439 + 0.0000000000000000 + 3 5 -1 +obdef +loc1d + 0.5000000000000000 +kind + -5 + 0 2 + 0.50000000000000000 + OBS 5 + -0.28483908888387055 + 0.30103267079449658 + 0.0000000000000000 + 4 6 -1 +obdef +loc1d + 0.5000000000000000 +kind + -6 + 0 2 + 0.29999999999999999 + OBS 6 + 664.18897402952678 + 664.13657445350088 + 0.0000000000000000 + 5 7 -1 +obdef +loc1d + 0.5000000000000000 +kind + -7 + 0 2 + 1.0000000000000000E-002 + OBS 7 + 1.2510430688864891 + 0.76733890921198655 + 0.0000000000000000 + 6 8 -1 +obdef +loc1d + 0.5000000000000000 +kind + -5 + 0 3 + 0.50000000000000000 + OBS 8 + 0.14649035752369940 + 0.30165218527115972 + 0.0000000000000000 + 7 9 -1 +obdef +loc1d + 0.5000000000000000 +kind + -6 + 0 3 + 0.29999999999999999 + OBS 9 + 996.25575565683448 + 996.14816275897965 + 0.0000000000000000 + 8 10 -1 +obdef +loc1d + 0.5000000000000000 +kind + -7 + 0 3 + 1.0000000000000000E-002 + OBS 10 + 1.5062181726711423 + 0.87797602399919394 + 0.0000000000000000 + 9 11 -1 +obdef +loc1d + 0.5000000000000000 +kind + -5 + 0 4 + 0.50000000000000000 + OBS 11 + 0.54462279403932623 + 0.30233594243716233 + 0.0000000000000000 + 10 12 -1 +obdef +loc1d + 0.5000000000000000 +kind + -6 + 0 4 + 0.29999999999999999 + OBS 12 + 1327.9926698239965 + 1328.1552851235194 + 0.0000000000000000 + 11 13 -1 +obdef +loc1d + 0.5000000000000000 +kind + -7 + 0 4 + 1.0000000000000000E-002 + OBS 13 + 1.3441050859614210 + 0.99898224420578696 + 0.0000000000000000 + 12 14 -1 +obdef +loc1d + 0.5000000000000000 +kind + -5 + 0 5 + 0.50000000000000000 + OBS 14 + 0.96724285779869490 + 0.30308378852289997 + 0.0000000000000000 + 13 15 -1 +obdef +loc1d + 0.5000000000000000 +kind + -6 + 0 5 + 0.29999999999999999 + OBS 15 + 1660.0844485052482 + 1660.1579415388080 + 0.0000000000000000 + 14 -1 -1 +obdef +loc1d + 0.5000000000000000 +kind + -7 + 0 5 + 1.0000000000000000E-002 diff --git a/models/seir/work/perfect_input.cdl b/models/seir/work/perfect_input.cdl new file mode 100644 index 0000000000..f07c245e32 --- /dev/null +++ b/models/seir/work/perfect_input.cdl @@ -0,0 +1,45 @@ +netcdf perfect_input { +dimensions: + member = 1 ; + metadatalength = 32 ; + location = 7 ; + time = UNLIMITED ; // (1 currently) +variables: + + char MemberMetadata(member, metadatalength) ; + MemberMetadata:long_name = "description of each member" ; + + double location(location) ; + location:short_name = "loc1d" ; + location:long_name = "location on a unit circle" ; + location:dimension = 1 ; + location:valid_range = 0., 1. ; + + double state(time, member, location) ; + state:long_name = "the model state" ; + + double time(time) ; + time:long_name = "valid time of the model state" ; + time:axis = "T" ; + time:cartesian_axis = "T" ; + time:calendar = "none" ; + time:units = "days" ; + +// global attributes: + :title = "true state from control" ; + :version = "$Id$" ; + :model = "SEIR" ; + :history = "NA" ; + +data: + + MemberMetadata = + "true state" ; + + location = 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5 ; + + state = 331996196, 1.0, 1.0, 1.0, 0.5, 0.3, 0.1 ; + + time = 1.0 ; + +} diff --git a/models/seir/work/quickbuild.sh b/models/seir/work/quickbuild.sh new file mode 100755 index 0000000000..a08deca17f --- /dev/null +++ b/models/seir/work/quickbuild.sh @@ -0,0 +1,59 @@ +#!/usr/bin/env bash + +# DART software - Copyright UCAR. This open source software is provided +# by UCAR, "as is", without charge, subject to all terms of use at +# http://www.image.ucar.edu/DAReS/DART/DART_download + +main() { + +export DART=$(git rev-parse --show-toplevel) +source "$DART"/build_templates/buildfunctions.sh + +MODEL=seir +LOCATION=oned + + +programs=( +closest_member_tool +filter +model_mod_check +perfect_model_obs +) + +serial_programs=( +create_fixed_network_seq +create_obs_sequence +fill_inflation_restart +integrate_model +obs_common_subset +obs_diag +obs_sequence_tool +) + +model_programs=( +) + +model_serial_programs=( +) + +# quickbuild arguments +arguments "$@" + +# clean the directory +\rm -f -- *.o *.mod Makefile .cppdefs + +# build any NetCDF files from .cdl files +cdl_to_netcdf + +# build and run preprocess before making any other DART executables +buildpreprocess + +# build +buildit + +# clean up +\rm -f -- *.o *.mod + +} + +main "$@"