Adding a new SEIR model that simulates Infectious Diseases

[mgharamti]

Description:

Added a new model to DART. The model simulates the spread of infectious diseases for example COVID-19. The model has 7 state variables and 10 other parameters. This version of the code can be used for state estimation only. Future versions may also include parameter estimation capabilities. The model makes use of identity observations such as: number of people recovered, number of those who died or the number of people who got vaccinated.

Fixes issue

NA

Types of changes

• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)
• Documentation update

Documentation changes needed?

• My change requires a change to the documentation.
  • I have updated the documentation accordingly.

Tests

Ran model_mod_check, perfect_model_obs, filter with different obs_seq files.

Checklist for merging

• Updated changelog entry
• Documentation updated
• Update conf.py

Checklist for release

• Merge into main
• Create release from the main branch with appropriate tag
• Delete feature-branch

Testing Datasets

• Dataset needed for testing available upon request
• Dataset example is included
• No dataset needed

[hkershaw-brown]

Hi Moha, this looks really good. Quick question before I review this:
At the moment this pull request is into NCAR:seir from mgharamti:seir. Do you want to do the pull request into NCAR:main and release this model_mod as part of dart?

[mgharamti]

Opps, didn't realize that the PR wasn't made into main. Yes, the intention is to push this model to the main public dart code. Do you want me to cancel and do the PR again?

[hkershaw-brown]

no problem I switched to be be into main. Cheers, Helen

[hkershaw-brown]

d86fdc5 swapped :raw-math: for :math: since readthedocs not rendering the :raw-math:

[hkershaw-brown]

Hi Moha, this looks great.

It looks like a lot of suggestions in the review, but it is not really.

I'd like to have the comments, specific to SEIR rather than keeping in the advice from the template model_mod. This (I think) makes reading the code easier for someone using the SEIR model.

The other two changes are

• the perturb - I don't think you are perturbing based on the qty or location, so you can simplify the perturb call (see suggestion)
• For the location, since everything is 0.5_r8 you don't need to have an array for locations.

Lastly, I'm in two minds about this: there is no check that the module has been initialized in the public routines. e.g.
if ( .not. module_initialized ) call static_init_model

This is a hot take, but I'm ok not having this check in there.

Just for a quick check I compiled and ran the model with ifort/gfortran (with & without debugging) on my mac. I haven't checked the equations - I'm assuming you're happy with the science.

Cheers,
Helen

[mgharamti]

Helen, thanks for all the great suggestions especially for the location and perturbation.
For the check to initialized modules, I'm also ok either way. I just felt it's a simple model and wanted to keep the code volume minimal but happy to change if needed.
As far as the science go, I think it should be ok. I compared the equations to the original work I did with my collaborators and the code that Shaniah left us with. It should be ok.

[nancycollins]

if you make the state_loc array only a single value i don't think this routine is going to do what it should. it compiles and runs but i think it will write only a single location to the output netcdf file in the last index location because this form assumes the last argument is an index not a count, leaving the rest of the location array uninitialized. you could either put state_loc back as a single-valued array, or make a loop from i=1 to msize around this call and pass in i as the last arg here.

(full disclosure - i didn't download and run this to be sure. could someone delete any existing netcdf output files, run it, and dump the 'locations' array to be sure?)

[mgharamti]

Maybe this confirms your suspesions Nancy. Here is the output from pmo

gharamti@cisl-reyno: seir/work (mgharamti-seir)$ ncdump -v location perfect_output.nc
netcdf perfect_output {
dimensions:
	metadatalength = 64 ;
	member = 1 ;
	time = UNLIMITED ; // (1 currently)
	NMLlinelen = 129 ;
	NMLnlines = 245 ;
	location = 7 ;
variables:
	char MemberMetadata(member, metadatalength) ;
		MemberMetadata:long_name = "Metadata for each copy/member" ;
	char inputnml(NMLnlines, NMLlinelen) ;
		inputnml:long_name = "input.nml contents" ;
	double time(time) ;
		time:long_name = "time" ;
		time:calendar = "none" ;
		time:units = "days" ;
	double location(location) ;
		location:description = "location coordinates" ;
		location:location_type = "loc1d" ;
		location:long_name = "location on unit circle" ;
		location:storage_order = "X" ;
		location:units = "none" ;
	double state(time, member, location) ;
		state:units = "none" ;

// global attributes:
		:title = "perfect_output.nc" ;
		:assim_model_source = "direct_netcdf_mod.f90" ;
		:creation_date = "YYYY MM DD HH MM SS = 2024 02 28 15 53 57" ;
		:model_source = "seir/model_mod.f90" ;
		:model = "seir" ;
		:DART_note = " pmo restart" ;
data:

 location = _, _, _, _, _, _, 0.5 ;
}

[mgharamti]

In any case, location is really meaningless here .. so whatever that array has I don't think it impacts the solution of the state unless I'm missing something.

[nancycollins]

since location really doesn't mean anything here you could remove that call to nc_write_location() completely. see if that upsets any of the diagnostic tools or utilities. i suspect not. i think that's cleaner than only writing part of a variable.

p.s. if it does, then i'd vote for the loop solution:

do i = 1, msize
call nc_write_location(ncid, state_loc, i)
enddo

[mgharamti]

I removed the nc_write_location and nc_write_location_atts and the code compiled and ran with no issues. Here is a filter output:

ncdump -h filter_output.nc                    
netcdf filter_output {
dimensions:
	metadatalength = 64 ;
	member = 3 ;
	time = UNLIMITED ; // (1 currently)
	NMLlinelen = 129 ;
	NMLnlines = 245 ;
	location = 7 ;
variables:
	char MemberMetadata(member, metadatalength) ;
		MemberMetadata:long_name = "Metadata for each copy/member" ;
	char inputnml(NMLnlines, NMLlinelen) ;
		inputnml:long_name = "input.nml contents" ;
	double time(time) ;
		time:long_name = "time" ;
		time:calendar = "none" ;
		time:units = "days" ;
	double state(time, member, location) ;
		state:units = "none" ;
	double state_mean(time, location) ;
		state_mean:units = "none" ;
	double state_sd(time, location) ;
		state_sd:units = "none" ;
	double state_priorinf_mean(time, location) ;
		state_priorinf_mean:units = "unitless" ;
	double state_priorinf_sd(time, location) ;
		state_priorinf_sd:units = "unitless" ;

// global attributes:
		:title = "filter_output.nc" ;
		:assim_model_source = "direct_netcdf_mod.f90" ;
		:creation_date = "YYYY MM DD HH MM SS = 2024 02 28 17 18 37" ;
		:model_source = "seir/model_mod.f90" ;
		:model = "seir" ;
		:DART_note = "output prior inflation sd" ;

[mgharamti]

Hi Helen, I just did that. Hopefully, I didn't break anything.

[mgharamti]

To follow up, I did change my cdl file into this:

ncdump perfect_input.nc 
netcdf perfect_input {
dimensions:
	member = 1 ;
	phase = 7 ;
	metadatalength = 32 ;
	time = UNLIMITED ; // (1 currently)
variables:
	char PhaseMetadata(phase, metadatalength) ;
		PhaseMetadata:long_name = "description of each phase" ;
	double state(time, phase) ;
		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:

 PhaseMetadata =
  "Disease phase: Susceptible (S)",
  "Disease phase: Exposed     (E)",
  "Disease phase: Infected    (I)",
  "Disease phase: Quarantined (Q)",
  "Disease phase: Recovered   (R)",
  "Disease phase: Dead        (D)",
  "Disease phase: Vaccinated  (V)" ;

 state =
  331996196, 1, 1, 1, 0.5, 0.3, 0.1 ;

 time = 1 ;
}

Added phase dimension and described its metadata. I had to keep the member dimension to be able to run pmo. After running pmo and filter, the location shows up as 7 in the state dimension. Ideally, this model would need a dummy location module (i.e., not a spatial one). I am concerned this might be confusing for users. I am more inclined to returning the location dimension and adding some more documentation about it in the readme file. What are your thoughts?

[hkershaw-brown]

Hi Moha,
this location dimension is hardcoded in state_structure_mod, it is not even the location_mod involved with this.

DART/assimilation_code/modules/io/state_structure_mod.f90

Lines 446 to 448 in 77bb8c2

	state%domain(dom_id)%unique_dim_names(1) = 'location'
	state%domain(dom_id)%unique_dim_names(2) = 'member'
	state%domain(dom_id)%unique_dim_names(3) = 'time'

its a bad assumption in state_structure (I think it was put in assuming every model using a simple vector would look like lorenz_96- location on a unit circle.).
I'll add this problem to the list of state structure & IO issues.
Sorry about messing you around on this pull request. I agree with you for this pull request, keep location and note in the documentation.

[nancycollins]

if it's got to be there, then to avoid uninitialized data in the netcdf array
i'd put a simple loop around the write and do it msize times:

integer(i8) :: offset

do offset=1, msize
  call nc_write_location(ncid, state_loc, offset)
enddo

! and maybe add a comment near this code saying that for 
! this model the location isn't currently used.

[mgharamti]

Done

[nancycollins]

completely optional comment - now that source is a fixed string and not expanded when it's checked out (like subversion used to do), this can be (len=*) and avoid allocating more space than is needed.

[mgharamti]

Done

[hkershaw-brown]

Nice work Moha! approved