fluidity_options.rnc

include "spud_base.rnc"

include "adaptivity_options.rnc"
include "diagnostic_algorithms.rnc"
include "input_output.rnc"
include "solvers.rnc"
include "stabilisation.rnc"
include "mesh_options.rnc"
include "physical_parameters.rnc"
include "prognostic_field_options.rnc"
include "prescribed_field_options.rnc"
include "spatial_discretisation.rnc"
include "temporal_discretisation.rnc"
include "flredecomp.rnc"
include "multiphase_interaction.rnc"
include "equation_of_state.rnc"

start =
   (
      ## The root node of the options dictionary.
      element fluidity_options {
         comment,
         ## Model output files are named according to the simulation
         ## name, e.g. [simulation_name]_0.vtu. Non-standard
         ## characters in the simulation name should be avoided.
         element simulation_name {
            anystring
         },
         ## Option problem_type does not change the tree.  It is just used for options checking.
         element problem_type {
            element string_value {
               # Lines is a hint to the gui about the size of the text box.
               # It is not an enforced limit on string length.
               attribute lines { "1" },
               ( "fluids" | "oceans" | "multimaterial" | "stokes" | "large_scale_ocean_options" | "foams" | "multiphase" )
            },
            comment
         },
         geometry,
         ## Input/output options
         element io {
            ## Format for dump files. Only vtk for now.
            element dump_format {
               element string_value{
                  "vtk"
               }
            },
            (
               ## Period between dumps in time units.
               ##
               ## Specifies the period between each dump of the solution to disk.
               ## A value of 0.0 indicates that there would be a dump at every timestep.
               element dump_period {
                (
                     element constant {
                     real
                   }|
                   ## Python function prescribing real input. Functions should be of the form:
                   ##
                   ##  def val(t):
                   ##     # Function code
                   ##     return # Return value
                   ##
                   ## 
                   element python {
                     python_code
                   }
                )  
               }|                 
               ## Dump period, in timesteps.
               ## 
               ## Specifies the number of timesteps between each dump of the solution to disk.
               ## A value of 0 indicates a dump at every timestep.
               element dump_period_in_timesteps {
                 (
                     element constant {
                     integer
                   }|
                   ## Python function prescribing integer input. Functions should be of the form:
                   ##
                   ##  def val(t):
                   ##     # Function code
                   ##     return # Return value
                   ##
                   ## 
                   element python {
                     python_code
                   }
                )   
               }
            ),
            ## Disable dump at simulation start
            element disable_dump_at_start {
               comment
            }?,
            ## Disable dump at simulation end
            element disable_dump_at_end {
               comment
            }?,
            # every CPUDUM seconds write results to disc.
            ## This is usually disabled.
            element cpu_dump_period {
               real
            }?,
            ## The period between dumps in walltime seconds. This is usually disabled.
            element wall_time_dump_period {
               real
            }?,
            ## Number of dumps before we overwrite previous dumps.
            element max_dump_file_count {
               integer
            }?,
            (
               ## The mesh on to which all the fields will be
               ## interpolated for VTK output.
               element output_mesh {
                  attribute name { "VelocityMesh" }
               }|
               ## The mesh on to which all the fields will be
               ## interpolated for VTK output.
               element output_mesh {
                  attribute name { "PressureMesh" }
               }|
               ## The mesh on to which all the fields will be
               ## interpolated for VTK output.
               element output_mesh {
                  attribute name { "CoordinateMesh" }
               }|
               ## The mesh on to which all the fields will be
               ## interpolated for VTK output.
               element output_mesh {
                  attribute name { xsd:string }
               }
            ),
            ## Options for convergence analysis.
            element convergence {
               ## Whether to enable the creation of a convergence
               ## file, giving details of the convergence of each
               ## field over the global nonlinear iteration loop.
               ## The .convergence file is in the same format as the .stat file.
               element convergence_file {
                  comment
               }?,
               ## Write state to a vtu on every iteration.
               ## 
               ## This is a useful debugging tool if things are not converging.
               ## To prevent an excessive number of files being accumulated previous
               ## timestep files will be overwritten hence it is best to use
               ## in conjunction with /timestepping/nonlinear_iterations/terminate_if_not_converged
               element convergence_vtus {
                  comment
               }?
            }?,
            ## Whether to enable dumping of checkpointing output.
            ##
            ## See http://amcg.ese.ic.ac.uk/index.php?title=Local:Checkpointing_from_new_options
            element checkpointing {
               ## Checkpointing period, in dumps. Non-negative value
               ## required. A value of zero indicates that checkpoints
               ## should be created at every dump. If
               ## /io/max_dumpfile_count is exceeded then earlier
               ## checkpoints may be overwritten.
               element checkpoint_period_in_dumps {
                  integer
               },
               ## Enable to checkpoint at simulation start.
               element checkpoint_at_start {
                  comment
               }?,
               ## Enable to force a checkpoint at simulation end.
               element checkpoint_at_end {
                  comment
               }?,
               comment
            }?,
            ## Diagnostic output (.stat file) options
            element stat {
               ## Enable to write diagnostic output at simulation start
               element output_at_start {
                  comment
               }?,
               ## Enable to write diagnostic output immediately before mesh adapts
               element output_before_adapts {
                  comment
               }?,
               ## Enable to write diagnostic output immediately after mesh adapts
               element output_after_adapts {
                  comment
               }?,
               comment
            },
            ## Specification of detectors. Note that when running in parallel the detector output is in binary format even if binary_output is not enabled. When running in serial, although the output is in principle still generated in ascii format if binary_output is not enabled, it is not certain that it is working well. Hence, it is recommended to enable binary_output and work with binary files. 
            element detectors {
               (
                  ## A single static detector
                  element static_detector {
                     attribute name { xsd:string },
                     (
                        element location {
                           real_dim_vector
                     }|
                        ## File containing the detectors positions in binary form
                        element from_checkpoint_file {
                           attribute file_name { xsd:string },
                           ## The format of the input file containing field data.
                           element format {
                             element string_value {
                                 "binary"
                             }
                           }
                        }
                     )
                  }|
                  ## A single lagrangian detector
                  element lagrangian_detector {
                     attribute name { xsd:string },
                     (
                     ## This is the initial location of a detector that moves with the fluid velocity.
                        element location {
                           real_dim_vector
                     }|
                        ## File containing the detectors positions in binary form
                        element from_checkpoint_file {
                           attribute file_name { xsd:string },
                           ## The format of the input file containing field data.
                           element format {
                             element string_value {
                                 "binary"
                             }
                           }
                        }
                     )
                  }|
                  ## Detectors with their locations specified via a python function or from a file. Allows detector arrays to be added.
                  element detector_array {
                     attribute name { xsd:string },
                     ## The number of detectors prescribed by the python function.
                     element number_of_detectors {
                        integer
                     },
                     (
                        ## Create fixed detectors.
                        element static {
                           empty
                        }|
                        ## Create detectors which move with the fluid velocity.
                        element lagrangian {
                           empty
                        }
                     ),
                     (
                        ## Python function prescribing dimensional vector input. Functions should be of the form:
                        ##
                        ##  def val(t):
                        ##     # Function code
                        ##     return # Return value
                        ##
                        ## The return value must have length number_of_detectors.
                        ##
                        ## *** IMPORTANT NOTE ***
                        ##
                        ## The t argument is for future use only - currently detector locations are only set at simulation start.
                        element python {
                          python_code
                        }|
                        ## File containing the detectors positions in binary form
                        element from_file {
                           attribute file_name { xsd:string },
                           ## The format of the input file containing field data.
                           element format {
                             element string_value {
                                 "binary"
                             }
                           }
                        }|
                        ## File containing the detectors positions in binary form
                        element from_checkpoint_file {
                           attribute file_name { xsd:string },
                           ## The format of the input file containing field data.
                           element format {
                              element string_value {
                                 "binary"
                              }
                           }
                        }
                     )
                  }
               )*,
               (
                  ## By default Fluidity will fail if a detector has left the domain.
                  element fail_outside_domain {
                     empty
                  }|
                  ## Enable to write NaN values to detector output if a detector has left the domain.
                  element write_nan_outside_domain {
                     empty
                  }
               ),
               ## Enable to write detector output in ascii format. This output format will not work in parallel, if not enabled format defaults to binary.
               element ascii_output {
                  comment
               }?,
               ## Enable to let detectors move with the domain if mesh_movement is enabled.
               element move_with_mesh {
                  empty
               }?,
               lagrangian_timestepping?
            }?, 
            ## Options to create even more output in the logs:
            ##
            ## Note that the main option to control the log output is given on the command line:
            ##
            ## -v0  only output error and warnings
            ##
            ## -v1  also give "navigational information", to indicate where in the code we currently are
            ##
            ## -v2  also give any additional information (mins and maxes of fields, etc.)
            ##
            element log_output {
               ## Log all allocates and deallocates done for meshes, fields, sparsities and matrices.
               ##
               ## NOTE: Requires -v2
               element memory_diagnostics {
                  empty
               }?
            }?,
            ## include time as VTK FieldData with the name "TimeValue".
            ## Paraview versions >= 5.5 will automatically use this to
            ## set the time level.
            element time_as_metadata {
               empty
            }?
         },
         ## Particle options. Particles, like Lagrangian detectors, move with a specified velocity field.
         ## In addition, particles can carry attributes, being properties of the particles. These attributes
         ## can be constant, or based on a python expression that depends on other attribute values, at the end
         ## of the previous timestep, or field values interpolated at the particle location. Planned future
         ## functionality includes the ability to spawn and delete particles to ensure sufficient particles are
         ## present in every region, and the ability to interpolate from particles back to the (Eulerian) mesh.
         element particles {
            (
               ## Particles groups containing multiple particle arrays or subgroups. The particle group specifies
               ## which velocity field is used, as well as attributes carried and (eventually) spawning
               ## and deleting options.
               element particle_group {
                  attribute name { xsd:string },
                  element velocity_field {
                     attribute name { xsd:string }
                  },
                  ## Particle subgroups arrays specified within a particle group. All subgroups utilize the same velocity field
                  ## and carry the same attributes. Subgroups can have different initial positions specified via a python
                  ## function of from a file.
                  element particle_subgroup {
                     attribute name { xsd:string },
                     ## The number of particles prescribed by the python function.
                     element number_of_particles {
                        integer
                     },
                     (
                        ##Initial position of particles
                        element initial_position {
                           ## Python function prescribing dimensional vector input. Functions should be of the form:
                           ##
                           ##  def val(t):
                           ##     # Function code
                           ##     return # Return value
                           ##
                           ## The return value must have length number_of_particles.
                           ##
                           ## *** IMPORTANT NOTE ***
                           ##
                           ## The t argument is for future use only - currently particle locations are only set at simulation start.
                           element python {
                              python_code
                           }|
                           ## File containing the particle positions in binary form
                           element from_file {
                              attribute file_name { xsd:string },
                              ## The format of the input file containing field data.
                              element format {
                                 element string_value {
                                    "binary"
                                 }
                              }
                           }
                        } 
                     ),
                     (
                        ##Define particle attributes
                        element attributes {
                           element attribute {
                              attribute name { xsd:string },
                              generic_particle_attribute_choice
                           }+
                        }?
                     )
                  }+
               }+
            )*,
            (
               ## By default Fluidity will fail if a particle has left the domain.
               element fail_outside_domain {
                  empty
               }|
               ## Enable to write NaN values to particle output if a particle has left the domain.
               element write_nan_outside_domain {
                  empty
               }
            ),
            ## Enable to write particle output in ascii format
            element ascii_output {
               comment
            }?,
            ## Enable to let particles move with the domain if mesh_movement is enabled.
            element move_with_mesh {
               empty
            }?,
            lagrangian_timestepping
         }?,
         ## Options dealing with time discretisation
         element timestepping {
            ## Current simulation time. At the start of the simulation this
            ## is the start time.
            element current_time {
               real,
               ## The following excerpt from the Udunits
               ## documentation explains the time unit encoding by
               ## example:
               ##
               ## The specification:
               ##
               ## seconds since 1992-10-8 15:15:42.5 -6:00
               ##
               ## indicates seconds since October 8th, 1992 at 3
               ## hours, 15 minutes and 42.5 seconds in the afternoon
               ## in the time zone which is six hours to the west of
               ## Coordinated Universal Time (i.e.  Mountain Daylight
               ## Time). The time zone specification can also be
               ## written without a colon using one or two-digits
               ## (indicating hours) or three or four digits
               ## (indicating hours and minutes).
               ##
               ## Time units are particularly required in situations
               ## where the problem (time-varying) boundary conditions
               ## and/ initial conditions are a function of time as
               ## defined by a calendar.  Examples include atmospheric
               ## forcing and climatology. The current time, specified
               ## above, is zero at the reference data/time.
               element time_units{attribute date { xsd:string }}?
            },
            ## The time step size. If adaptive time stepping is used
            ## then this is the initial time step size.
            element timestep {
               real
            },
            ## Simulation time at which the simulation should end.
            element finish_time {
               real
            },
            ## Timestep after which the simulation should end.
            element final_timestep {
               integer
            }?,
            ## Maximum CPU time (in seconds) before the simulation terminates
            element cpu_time_limit {
               real
            }?,
            ## Maximum wall time (secs) taken up before
            ## simulation terminates writing results to disc.
            ## 
            ## This is usually disabled.
            element wall_time_limit {
               real
            }?,
            ## maximum number of non-linear iterations.
            ## 
            ## Manual suggests 2
            element nonlinear_iterations {
               integer,
               ## tolerance for non-linear iteration.
               ## Manual suggests 1.0E-12
               element tolerance {
                  real,
                  (
                      ## Select the norm with which you want the tolerance to be tested.
                      ##
                      ## The infinity norm.
                      element infinity_norm {
                        empty
                      }|
                      ## Select the norm with which you want the tolerance to be tested.
                      ##
                      ## The l2 norm.
                      element l2_norm {
                        empty
                      }|
                      ## Select the norm with which you want the tolerance to be tested.
                      ##
                      ## The l2 norm evaluated on a control volume mesh.
                      element cv_l2_norm {
                        empty
                      }
                  )
               }?,
               ## Terminate the simulation if the number of
               ## nonlinear_iterations is reached
               ## and the tolerance criterion is not met.
               ## This is mostly useful as a debugging option if you
               ## suspect the solution is not converging.
               element terminate_if_not_converged {
                  empty
               }?,
               ## Number of non-linear iterations for the first time step after adapting the mesh.
               ## This option will work only if the mesh_adaptivity is switched on.
               element nonlinear_iterations_at_adapt {
                  integer
               }?
            }?,
            ## Vary the timestep according to the courant number.
            element adaptive_timestep {
               ## The timestep will be adjusted (within the tolerance
               ## and bounds specified) to target this courant
               ## number. Timestep adapts occur at the end of each
               ## timestep and after a mesh adapt.
               element requested_cfl {
                  real
               },
               timestep_cfl_number_options,
               ## Minimum time step size.
               ## Manual suggests 0.0
               element minimum_timestep {
                  ## If enabled, signals model termination if a timestep less
                  ## than or equal to the minimum_timestep is requested. The
                  ## model will stop at the end of the timestep in order to
                  ## allow for the latest output to be written. 
                  element terminate_if_reached {
                     comment
                  }?,
                  real
               }?,
               ## Maximum time step size.
               ## Manual suggests 1.E+10
               element maximum_timestep {
                  real
               }?,
               ## The maximum ratio by which the timestep is allowed
               ## to increase in a timestep adapt. e.g., a value of
               ## 1.1 indicates that the timestep may be increased by
               ## at most 10%.
               element increase_tolerance {
                  real
               }?,
               ## Specify whether you want to calculate a new timestep
               ## at the first timestep or not.
               element at_first_timestep {
                  empty
               }?
            }?,
            ## Activate if you want to terminate the simulation once a
            ## steady state is reached.
            ## 
            ## Enable/disable fields in this check under each field in
            ## steady_state options.
            # Preprocessor legacy: STEDER = 0. is equivalent to inactive.
            element steady_state {
               ## Enter the tolerance to which you want a steady state to be judged.
               element tolerance {
                  real,
                  (
                      ## Select the norm with which you want the tolerance to be tested.
                      ##
                      ## The infinity norm.
                      element infinity_norm {
                        empty
                      }|
                      ## Select the norm with which you want the tolerance to be tested.
                      ##
                      ## The l2 norm.
                      element l2_norm {
                        empty
                      }|
                      ## Select the norm with which you want the tolerance to be tested.
                      ##
                      ## The l2 norm evaluated on a control volume mesh.
                      element cv_l2_norm {
                        empty
                      }
                  )
               },
               ## If activated compare the above tolerance to the rate
               ## of change of the fields. Otherwise compare it
               ## directly to the change in the field.
               element acceleration_form {
                  empty
               }?,
               ## Write out the changes in the tested fields to a .steady_state
               ## file.
               element steady_state_file {
                  (
                     ## Write steady state output in binary format
                     element binary_output {
                        comment
                     }|
                     ## Write steady state output in plain text format
                     element plain_text_output {
                        comment
                     }
                  ),
                  comment
               }?
            }?
         },
         physical_parameter_options?,
         ## The material or phase options
         element material_phase {
            attribute name { xsd:string },
            equation_of_state_options?,
            ## Subgridscale parameterisations
            element subgridscale_parameterisations {
               ## Lagrangian-averaged Navier-Stokes equations 
               element LANS {
                  (
                     ## smoothing length specified as isotropic homogeneous
                     element alpha_isotropic_homogeneous {
                        real
                     }|
                     element alpha_anisotropic_homogeneous_cartesian {
                        real_dim_symmetric_tensor
                     }
                  )
                  #               (
                  #                  element leray {
                  #                     empty
                  #                  }|
                  #                  element LANS_momentum_form {
                  #                     empty
                  #                  }|
                  #                  element LANS_stress_form {
                  #                     empty
                  #                  }
                  #               )
               }?,
               ## This is the generic length scale (vertical turbulence mixed layer) model, 
               ## based here on Warner et al 2005, Ocean Modelling 8:81-113,
               ## which is itself based on the works of Umlauf and Burchard 2003.
               ##
               ## The GLS model encodes four individual turbulence closure models which can be chosen via 'option' below.
               ##
               ## You will need to make sure that DistanceToTop and DistanceToBottom fields (under geometry/ocean_boundaries) 
               ## are switched on, as well as PerturbationDensity.
               element GLS {
                  ## The base GLS option:
                  ## 1. k-kl which is a variant of Mellor-Yamada 2.5
                  ## 2. k-epsilon
                  ## 3. k-omega
                  ## 4. 'gen' from Warner et al 2005
                  ##
                  ## k-epsilon is recommended.
                  element option {
                     element string_value {
                     "k-kl"|"k-epsilon"|"k-omega"|"gen"
                     }
                  },
                  ## The stability function choice:
                  ## 1. KanthaClayson-94 which corresponds to Kantha and Clayson 1994
                  ## 2. Galperin-88 which corresponds to Galperin et al 1988
                  ## 3. Canuto-01-A which corresponds to choice A from Canuto et al 2001
                  ## 4. Canuto-01-B which corresponds to choice B from Canuto et al 2001
                  ##
                  ## Canuto A or B are recommended.
                  element stability_function {
                     element string_value {
                     "KanthaClayson-94"|"GibsonLaunder-78"|"Canuto-01-A"|"Canuto-01-B"
                     }
                  },                  
                  ## The wall function choice:
                  ## 1. None - pick this unless you're using k-kl
                  ## 2. Mellor and Yamada (1980) - parabolic shape
                  ## 3. Burchard (1998) - symmetric linear shape
                  ## 4. Burchard (2001) - Used for infinitely deep basins
                  ## 5. Blumberg et al (1992) - open channel flow
                  ##
                  ## See manual for equations.
                  element wall_function {
                     element string_value {
                     "none"|"MellorYamda"|"Burchard98"|"Burchard01"|"Blumberg"
                     }
                  }?,
                  ## Smooth buoyancy frequency before using it?
                  element smooth_buoyancy{
                    empty
                  }?,
                  ## Smooth velocity shear before using it?
                  element smooth_shear{
                    empty
                  }?,
                  ## Do you want the boundary conditions to be set automatically?
                  ## Make sure the ocean geometry settings are enabled for this option
                  element calculate_boundaries {
                     element string_value {
                     "neumann"|"dirichlet"
                     },
                     element top_surface_ids {
                        integer_vector
                     },
                     element bottom_surface_ids {
                        integer_vector
                     }
                  }?,
                  ## Perform relaxation of the diffusivity and viscosity in the GLS model.
                  ## Value should be >=0 and < 1. 0 is no relaxation (i.e. always use the
                  ## most up-to-date value) and 0.9 would represent making the current value
                  ## be 0.9*old_value + 0.1*new_value. Default is 0.0
                  ## If being used with adaptivity, ensure you switch on interpolation of the 
                  ## GLSVerticalDiffusivity and GLSVerticalViscosity fields. You will seg fault
                  ## otherwise.
                  element relax_diffusivity {
                     real
                  }?,
                  ## Add extra parameterisation for internal wave breaking at base of MLD. This
                  ## is based on the same parameterisation in NEMO and smooths the TKE down 
                  ## the water column based on an exponential function
                  ## Set the two parameters: % of TKE to smooth down and length scale to do this over
                  element ocean_parameterisation {
                     element lengthscale {
                        real
                     },
                     element percentage {
                        real
                     }
                  }?,
                  ## Turbulent kinetic energy. Make sure that the Diffusivity tensor field in here is set to diagnostic/internal.
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSTurbulentKineticEnergy" },
                     (
                        element prognostic {
                           velocity_mesh_choice,
                           prognostic_scalar_field,
                           ## Minimum value of TKE in m2s2. A typical value will be around 1e-6
                           element minimum_value {
                              real
                           }
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  },
                  ## Generic second equation used in GLS. 
                  ## Make sure that the Diffusivity tensor field in here is set to diagnostic/internal.
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSGenericSecondQuantity" },
                     (
                        element prognostic {
                           velocity_mesh_choice,
                           prognostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                      )
                  },
                  ## Background viscosity
                  element tensor_field {
                     attribute rank { "2" },
                     attribute name { "GLSBackgroundViscosity" },
                     (
                        element prescribed {
                           mesh_choice,
                           prescribed_tensor_field
                        }|
                        element diagnostic {
                           mesh_choice,
                           (
                              sediment_concentration_dependent_viscosity_algorithm
                           ),
                           diagnostic_tensor_field
                        }
                     )
                  },                  
                  ## Background diffusivity
                  element tensor_field {
                     attribute rank { "2" },
                     attribute name { "GLSBackgroundDiffusivity" },
                     (
                        element prescribed {
                           mesh_choice,
                           prescribed_tensor_field
                        }
                     )
                  },
                  ## Eddy viscosity K_M
                  element tensor_field {
                     attribute rank { "2" },
                     attribute name { "GLSEddyViscosityKM" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  },
                  ## Eddy diffusivity K_H
                  element tensor_field {
                     attribute rank { "2" },
                     attribute name { "GLSEddyDiffusivityKH" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  },
                  ## Length scale (a diagnostic with GLS)
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSLengthScale" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,                  
                  ## Unedited TKE. The TKE filed has the upper and lower surfaces
                  ## altered with Dirichlet conditions for ouput. This is the
                  ## unedited surface.
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSTurbulentKineticEnergyOriginal" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Buoyancy frequency
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSBuoyancyFrequency" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Velocity shear
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSVelocityShear" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Shear production P
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSShearProduction" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Buoyancy production B
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSBuoyancyProduction" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Dissipation epsilon
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSDissipationEpsilon" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Stability function S_M
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSStabilityFunctionSM" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Stability function S_H
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSStabilityFunctionSH" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Source1 - TKE source term
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSSource1" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Source2 - Second Quantity source term
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSSource2" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Absorption1 - TKE absorption term
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSAbsorption1" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Absorption2 - Second Quantity absorption term
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSAbsorption2" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## GLS Wall function
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSWallFunction" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Vertical eddy viscosity
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSVerticalViscosity" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,                  
                  ## Vertical eddy diffusivity
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "GLSVerticalDiffusivity" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?
               }?,
               element Mellor_Yamada {
                  ## Kinetic Energy
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "KineticEnergy" },
                     (
                        element prognostic {
                           velocity_mesh_choice,
                           prognostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  },
                  ## Turbulent Length Scale x Kinetic Energy
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "TurbulentLengthScalexKineticEnergy" },
                     (
                        element prognostic {
                           velocity_mesh_choice,
                           prognostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                      )
                  },
                  ## Vertical Viscosity
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "VerticalViscosity" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?,
                  ## Vertical Diffusivity of Temperature
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "VerticalDiffusivity" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element aliased {
                           generic_aliased_field
                        }
                     )
                  }?
               }?,
               ## Trivial case in which the user supplies the diffusivity.
               element prescribed_diffusivity {
                  element tensor_field {
                     attribute rank { "2" },
                     attribute name { "PrescribedDiffusivity" },
                     element prescribed {
                        mesh_choice,
                        prescribed_tensor_field_no_adapt
                     }
                  }
               }?,
               ## Standard k-epsilon turbulence model (see e.g. Ferziger and Peric(2002) p.295).
               ## Solves 2 additional equations for TurbulentDissipation (epsilon) and TurbulentKineticEnergy (k),
               ## in order to close momentum equations.
               ## Generates an isotropic 'eddy viscosity', added to molecular viscosity field, that
               ## carries the influence of turbulence onto the velocity field. See the manual for details.
               element k-epsilon {
                  ## Turbulent kinetic energy (k).
                  ## 1. 'k_esilon'-type boundary conditions are recommended for this field.
                  ## 2. Turn on diffusivity, source and absorption diagnostic fields.
                  ## 3. Use the same mesh as Velocity.
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "TurbulentKineticEnergy" },
                     (
                        element prognostic {
                           velocity_mesh_choice,
                           prognostic_scalar_field
                        }|
                        element prescribed {
                           velocity_mesh_choice,
                           prescribed_scalar_field
                        }
                     )
                  },
                  ## TurbulentDissipation (epsilon).
                  ## 1. 'k_epsilon'-type boundary conditions are recommended for this field.
                  ## 2. Turn on diffusivity, source and absorption diagnostic fields.
                  ## 3. Use the same mesh as Velocity.
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "TurbulentDissipation" },
                     (
                        element prognostic {
                           velocity_mesh_choice,
                           prognostic_scalar_field
                        }|
                        element prescribed {
                           velocity_mesh_choice,
                           prescribed_scalar_field
                        }
                     )
                  },
                  ## Set the value of the background (laminar) viscosity field here.
                  ## Make it small to see the influence of the eddy viscosity.
                  ## IMPORTANT: make velocity/viscosity diagnostic.
                  element tensor_field {
                     attribute rank { "2" },
                     attribute name { "BackgroundViscosity" },
                     (
                        element prescribed {
                           velocity_mesh_choice,
                           prescribed_tensor_field
                        }|
                        element diagnostic {
                           velocity_mesh_choice,
                           (
                              sediment_concentration_dependent_viscosity_algorithm|
                              tensor_python_diagnostic_algorithm
                           ),
                           diagnostic_tensor_field
                        }
                     )
                  },
                  ## Set the value of the background (laminar) diffusivity field here.
                  ## This is used to calculate scalar field diffusivities.
                  element tensor_field {
                     attribute rank { "2" },
                     attribute name { "BackgroundDiffusivity" },
                     (
                        element prescribed {
                           velocity_mesh_choice,
                           prescribed_tensor_field
                        }
                     )
                  }?,
                  ## Eddy viscosity (turbulent diffusion of velocity).
                  ## This is a fictitious isotropic viscosity, added to normal viscosity field, that
                  ## carries the influence of turbulence onto the velocity field.
                  element tensor_field {
                     attribute rank { "2" },
                     attribute name { "EddyViscosity" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_tensor_field
                        }
                     )
                  },
                  ## Scalar component of the eddy viscosity tensor.
                  ## This will appear in the stat file.
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "ScalarEddyViscosity" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }|
                        element prescribed {
                           velocity_mesh_choice,
                           prescribed_scalar_field
                        }
                     )
                  },
                  ## Integral length scale of the turbulence (diagnostic).
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "LengthScale" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }
                     )
                  },
                  ## f_1 damping coefficient for low_Re k-epsilon model
                  ## Required for low_Re boundaries. If no low_Re boundaries are present
                  ## this will be set to 1.0 throughout the domain and will not
                  ## affect the result
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "f_1" },
                     element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field                        
                     }
                  },
                  ## f_2 damping coefficient for low_Re k-epsilon model
                  ## Required for low_Re boundaries. If no low_Re boundaries are present
                  ## this will be set to 1.0 throughout the domain and will not
                  ## affect the result
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "f_2" },
                     element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field                        
                     }
                  },
                  ## f_mu damping coefficient for low_Re k-epsilon model
                  ## Required for low_Re boundaries. If no low_Re boundaries are present
                  ## this will be set to 1.0 throughout the domain and will not
                  ## affect the result
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "f_mu" },
                     element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field                        
                     }
                  },
                  ## Limit the maximum value of damping function values. This helps stability when using low_Re number
                  ## boundary conditions. If not using low_Re boundaries this will have no effect. 
                  ## Recommended value (default): 10.0
                  element max_damping_value {real},
                  ## Describes distance to nearest solid wall.
                  ## Required for low_Re boundaries.
                  ## For simple geometries the simplest method of providing this information is to use a python function.
                  ## For complex geometries where this is not possible precursive Eikonal equation or Poisson equation
                  ## simulations must be run using Fluidity to determine the values for this field.
                  ## Details of how this is done can be found in:
                  ## Tucker, P 2011: "Hybrid Hamilton/Jacobi/Poisson wall distance function model"
                  ## Elias et al 2007: "Simple finite element-based computation of distance functions in unstructured grids"
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "DistanceToWall" },
                     element prescribed {
                        velocity_mesh_choice,
                        prescribed_scalar_field
                     }
                  }?,
                  ## VLES Filter function (diagnostic).
                  ## Run a VLES by switching this option on.
                  ## If so then filtering is applied to the eddy viscosity.
                  element scalar_field {
                     attribute rank { "0" },
                     attribute name { "VLESFilter" },
                     (
                        element diagnostic {
                           internal_algorithm,
                           velocity_mesh_choice,
                           diagnostic_scalar_field
                        }
                     )
                  }?,
                  ## Eddy-viscosity coefficient: nu_T = density * C_mu * k**2 / epsilon.
                  ## Recommended value (default): 0.09.
                  element C_mu {real},
                  ## TurbulentDissipation production coefficient.
                  ## Recommended value (default): 1.44.
                  element C_eps_1 {real},
                  ## TurbulentDissipation destruction coefficient.
                  ## Recommended value (default): 1.92.
                  element C_eps_2 {real},
                  ## Turbulent Schmidt number (eddy viscosity coefficient from k equation).
                  ## This is also the ratio of eddy viscosity to eddy diffusivity for use in other scalar fields.
                  ## Recommended value (default): 1.0.
                  element sigma_k {real},
                  ## Turbulent Schmidt number (eddy-viscosity coefficient from epsilon equation).
                  ## Recommended value (default): 1.3.
                  element sigma_eps {real},
                  ## The Schmidt number (ratio of viscous diffusion rate to momentum diffusion rate) for 
                  ## massive scalar fields, or Prandtl number (ratio of viscous diffusion rate to thermal diffusion rate) 
                  ## for thermal fields. This is used to calculate the turbulent buoyancy term.
                  ## Recommended value (default): 1.0.
                  element sigma_p {real},
                  ## Time discretisation of the source terms in the k and epsilon equations and also
                  ## of the eddy viscosity
                  element time_discretisation {  
                     ## Implicit/explicit control (THETA) of source terms and eddy viscosity in
                     ## the k-epsilon model
                     ##  = 0.  -- explicit
                     ##  = 0.5 -- Crank-Nicolson
                     ##  = 1.  -- implicit
                     element theta {real},
                     ## Each term in the k and epsilon equations can be implemented as either a 
                     ## source or as an absorbtion term. See the manual under 
                     ## 'Parameterisations/Turbulent flow modelling and simulation/Reynolds Averaged Navier Stokes (RANS) Modelling/Standard k − ε Turbulence Model/Time Discretisation and Coupling' 
                     ## for more information.
                     element source_term_implementation { 
                        element production_term { 'source'|'absorbtion' },
                        element destruction_term { 'source'|'absorbtion' },
                        element buoyancy_term { 'source'|'absorbtion' }
                     }
                  },
                  ## The source terms in the k-epsilon model are calculated before the field is 
                  ## solved. The process requires inversion of a mass matrix. This element contains
                  ## options for defining how this is done. For P1 meshes mass lumping can be used,
                  ## for other discretisations the mass matrix can not be easily inverted and so a 
                  ## solve must be carried out, with solver options specified.
                  element mass_terms {  
                     (
                        element lump_mass { empty }
                        |
                        element use_consistent_mass_matrix { 
                           element solver {
                              linear_solver_options_asym_scalar
                           }
                        }                     
                     )
                  },
                  ## produce vtu output of individual kk and eps source terms and set prescribed 
                  ## source terms for k and epsilon fields (for mms tests)
                  element debugging_options{
                     ## output fields for each of the source terms in the k-epsilon model
                     element source_term_output_fields{
                        element scalar_field {
                           attribute rank { "0" },
                           attribute name { "TurbulentKineticEnergy_production_term" },
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_scalar_field                        
                           }
                        }?,
                        element scalar_field {
                           attribute rank { "0" },
                           attribute name { "TurbulentKineticEnergy_destruction_term" },
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_scalar_field                        
                           }
                        }?,
                        element scalar_field {
                           attribute rank { "0" },
                           attribute name { "TurbulentKineticEnergy_buoyancy_term" },
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_scalar_field                        
                           }
                        }?,
                        element scalar_field {
                           attribute rank { "0" },
                           attribute name { "TurbulentDissipation_production_term" },
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_scalar_field                        
                           }
                        }?,
                        element scalar_field {
                           attribute rank { "0" },
                           attribute name { "TurbulentDissipation_destruction_term" },
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_scalar_field                        
                           }
                        }?,
                        element scalar_field {
                           attribute rank { "0" },
                           attribute name { "TurbulentDissipation_buoyancy_term" },
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_scalar_field                        
                           }
                        }?
                     },
                     ## Enable to apply prescribed source terms into the k and epsilon equations.
                     ## This is useful for MMS tests
                     element prescribed_source_terms{
                        element scalar_field {
                           attribute rank { "0" },
                           attribute name { "TurbulentKineticEnergyPrescribedSource" },
                           element prescribed {
                              velocity_mesh_choice,
                              prescribed_scalar_field                        
                           }
                        }?,
                        element scalar_field {
                           attribute rank { "0" },
                           attribute name { "TurbulentDissipationPrescribedSource" },
                           element prescribed {
                              velocity_mesh_choice,
                              prescribed_scalar_field                        
                           }
                        }?
                     },
                     ## enabling this option disables the production term in the k and epsilon equations
                     element disable_production_term{ empty }?,
                     ## enabling this option disables the destruction term in the k and epsilon equation
                     element disable_destruction_term{ empty }?,
                     ## enabling this option disables the bouyancy term in the k and epsilon equation
                     element disable_buoyancy_term{ empty }?,
                     ## enabling this option enables the low-Re number damping functions regardless of whether there is
                     ## a low-Re boundary condition
                     element enable_lowRe_damping{ empty }?,
                     ## Enabling this option disables feedback from the k-epsilon model back into 
                     ## the rest of the model. Reynolds stress tensor is set to zero by zeroing 
                     ## EddyViscosity and the added source term, based on k, in the momentum equation.
                     ## Hence, Viscosity will always stay as the BackgroundViscosity,
                     ## diffusivities will remain at the relevant BackgroundDiffusivity, and the 
                     ## momentum equation will be as if there were no turbulence model present. 
                     ## 
                     ## ScalarEddyViscosity is still calculated as normal.
                     element zero_reynolds_stress_tensor{ empty }?
                  }?
               }?
            }?,
            ## Pressure
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Pressure" },
               ## Field type
               (
                  element prognostic {
                     # mesh choice with PressureMesh as first option
                     pressure_mesh_choice,
                     prognostic_pressure_field
                  }|
                  element prescribed {
                     # mesh choice with PressureMesh as first option
                     pressure_mesh_choice,
                     prescribed_scalar_field
                  }|
                  ## Compute pressure from Density and InternalEnergy
                  ## via a compressible equation of state.
                  element diagnostic {
                     # mesh choice with PressureMesh as first option
                     pressure_mesh_choice,
                     internal_algorithm,
                     diagnostic_scalar_field
                  }|
                  element aliased {
                     attribute material_phase_name { xsd:string },
                     attribute field_name {"Pressure" }
                  }
               )
            }?,
            ## Density
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Density" },
               ## Field type
               (
                  element diagnostic {
                     internal_algorithm,
                     velocity_mesh_choice,
                     diagnostic_scalar_field
                  }|
                  element prognostic {
                     pressure_mesh_choice,
                     prognostic_density_field
                  }|
                  element aliased {
                     attribute material_phase_name { xsd:string },
                     attribute field_name {"Density" }
                  }
               )
            }?,
            ## Velocity vector and momentum options
            element vector_field {
               attribute rank { "1" },
               attribute name { "Velocity" },
               ## Field type
               (
                  element prognostic {
                     velocity_mesh_choice,
                     prognostic_velocity_field
                  }|
                  element prescribed {
                     velocity_mesh_choice,
                     prescribed_vector_field
                  }|
                  element diagnostic {
                     velocity_mesh_choice,
                     vector_python_diagnostic_algorithm,
                     diagnostic_vector_field
                  }|                  
                  element aliased {
                     attribute material_phase_name { xsd:string },
                     attribute field_name {"Velocity" }
                  }
               )
            },
            scalar_field_choice*,
            vector_field_choice*,
            tensor_field_choice*,
            population_balance*,
            ## Phase interaction options for Fluidity's multiphase flow model
            element multiphase_properties {
               comment,
               (element particle_diameter { real }|
               element particle_diameter_use_scalar_field { anystring })?,
               element apply_diameter_cap {
                  # default value 1.0E-12
                  element lower_cap { real }
               }?,
               ## If this is the fluid phase, 
               ## the effective conductivity is required 
               ## for the inter-phase heat transfer term.
               element effective_conductivity { 
                     real 
               }?,
               ## The specific heat (at constant volume) is required 
               ## for the inter-phase heat transfer term.
               element specific_heat { 
                     real
               }?
            }?,
            sediment?
         }+,
         mesh_adaptivity_options?,
         ## turbine model
         element turbine_model {
           ## specifies a turbine
           element turbine {
             attribute name {xsd:string},
             ## Model type
               (
                  ## Enforces the turbine model by using the dirichlet condition.
                  element dirichlet {
                     ## Name of the upstream (weak or strong) dirichlet boundary condition where the turbine model shall be applied to.
                     element boundary_condition_name_1 {
                       attribute name {xsd:string}
                     },
                     ## Name of the downstream (weak or strong) dirichlet boundary condition where the turbine model shall be applied to.
                     element boundary_condition_name_2 {
                       attribute name {xsd:string}
                     },
                     ## These coordinates define the upstream point where the free surface will be evaluated in order to calculate the head.
                     element free_surface_point_1 {
                       real_dim_vector
                     },
                     ## These coordinate define the downstream point where the free surface will be evaluated in order to calculate the head.
                     element free_surface_point_2 {
                       real_dim_vector
                     },
                     ## Python function val(h) takes the head h (=free_surface_point_1-free_surface_point_2) and returns the total flux through turbine (in m^3/s). A return value 
                     ## of 0.0 will deactivate the turbine model which results in the dirichlet conditions prescribed at these boundaries.
                     ## A flux value will be outflow for boundary_condition_name_1 and inflow for boundary_condition_name_2 and therefore sign(h)==sign(val(h)) should be true to achieve stability.
                     element head_flux {
                       python_code
                     }
                  }|
                  ## Implements the turbine model by using the fluxes from the discontinues galerkin model. Works only if dg is used. 
                  element flux {
                       ## Model type
                       (
                         element penalty {
                              ## Python function val(h) takes the pressure jump h and returns the flux factor between 0 and 1. A return value 
                              ## of 0.0 will deactivate the turbine model which results in the dirichlet conditions prescribed at these boundaries 
                              ## while a factor of +INF represents a connected domain.
                              element factor {
                                  python_code
                              }
                           }|
                         element dg {
                              ## Python function val(h) takes the pressure jump h and returns the flux factor between 0 and 1. A return value 
                              ## of 0.0 will deactivate the turbine model which results in the dirichlet conditions prescribed at these boundaries 
                              ## while a factor of 1 represents a connected domain.
                              element factor {
                                  python_code
                              }

                           }
                       ),
                       ## Name of the upstream (weak or strong) dirichlet boundary condition where the turbine model shall be applied to.
                       element boundary_condition_name_1 {
                         attribute name {xsd:string}
                       },
                       ## Name of the downstream (weak or strong) dirichlet boundary condition where the turbine model shall be applied to.
                       element boundary_condition_name_2 {
                         attribute name {xsd:string}
                       },
                       (
                         ## The turbine is always running. This is mainly used for debugging purposes.
                         element always_on {
                              empty
                           }|
                         ## The turbine is never running. This is mainly used for debugging purposes.
                         element always_off {
                              empty
                           }
                       )
                  }
               )
           }*
         }?, 
         biology?,
         forcing?,
         flredecomp?,
         multiphase_interaction?,
         element population_balance_continuous_phase_name{ string }?
      }
   )      

population_balance = 
   (
      ##
      element population_balance {
         (
            ## Determine initial weighted_abscissa and weights using the PD algorithm and the 
            ## initial conditions set in the diagnostic moment fields.
            ## 
            ## In this case it does not matter what is set in the weighted_abscissa and weights
            ## fields initial conditions option
            element calculate_initial_conditions_from_moments {
               empty
            }|
            ## Use initial conditions set in the weighted_abscissa and weights prognostic
            ## fields
            ## 
            ## In this case it does not matter what is set in the moment field initial conditions 
            ## options
            element use_prognostic_field_initial_conditions {
               empty
            }
         ),
         attribute name { xsd:string },
         element abscissa {
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Abscissa_0" },
               ## Field type
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            },
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Abscissa_1" },
               ## Field type
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            },
            element scalar_field {
               attribute rank { "0" },
               attribute name { xsd:string },
               ## Field type
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            }*
         },
         element weights {
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Weight_0" },
               ## Field type
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }
            },
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Weight_1" },
               ## Field type
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }
            },
            element scalar_field {
               attribute rank { "0" },
               attribute name { xsd:string },
               ## Field type
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }
            }*
         },
         element weighted_abscissa {
            element scalar_field {
               attribute rank { "0" },
               attribute name { "WeightedAbscissa_0" },
               ## Field type
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }
            },
            element scalar_field {
               attribute rank { "0" },
               attribute name { "WeightedAbscissa_1" },
               ## Field type
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }
            },
            element scalar_field {
               attribute rank { "0" },
               attribute name { xsd:string },
               ## Field type
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }
            }*
         },
         ## This option only works if you are using a continuous mesh for the population balance scalar fields.
         ## A DG mesh will use inverse mass calculated at element level.
         ## In the population balance DQMOM algorithm all sources for the 
         ## advection-diffusion equations are calculated at the quadrature points first 
         ## and then interpolated back to the nodal points. 
         ## For a continuous mesh, this interpolation is made by inverting the mass matrix. 
         ## Mass lumping is fairly quick and can be used with 1st order continuous meshes. 
         ## But the error is high for higher order continuous meshes and full mass matrix should be used.
         element adv_diff_source_term_interpolation {
            ## The mass matrix is lumped to obtain source fields at nodal points.
            element use_mass_lumping { empty } |
            ## Full mass marix is used and the linear system is solved iteratively to obtain source fields at nodal points.
            element use_full_mass_matrix {
               ## This block specifies the solver options for solving the linear system to interpolate advection-diffusion source terms at nodal points.
               element solver {
                  linear_solver_options_asym_scalar
               }
            }
         },
         element moments {
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Moment_0" },
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  (
                     ## Initial condition for WholeMesh
                     ##
                     ## Only specify one condition if not using mesh regions.
                     ## Otherwise select other initial_condition option, specify region_ids
                     ## and distinct names.  Then add extra intial conditions for other regions.
                     element initial_condition {
                        attribute name { "WholeMesh" },
                        input_choice_initial_condition_real
                     }|
                     ## Multiple initial_conditions are allowed if specifying
                     ## different values in different
                     ## regions of the mesh (defined by region_ids).  In this case
                     ## each initial_condition
                     ## requires a distinct name for the options dictionary.
                     element initial_condition {
                        attribute name { string },
                        region_ids?,
                        input_choice_initial_condition_real
                     }
                  )+,
                  diagnostic_scalar_field
               }
            },
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Moment_1" },
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  (
                     ## Initial condition for WholeMesh
                     ##
                     ## Only specify one condition if not using mesh regions.
                     ## Otherwise select other initial_condition option, specify region_ids
                     ## and distinct names.  Then add extra intial conditions for other regions.
                     element initial_condition {
                        attribute name { "WholeMesh" },
                        input_choice_initial_condition_real
                     }|
                     ## Multiple initial_conditions are allowed if specifying
                     ## different values in different
                     ## regions of the mesh (defined by region_ids).  In this case
                     ## each initial_condition
                     ## requires a distinct name for the options dictionary.
                     element initial_condition {
                        attribute name { string },
                        region_ids?,
                        input_choice_initial_condition_real
                     }
                  )+,
                  diagnostic_scalar_field
               }
            },
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Moment_2" },
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  (
                     ## Initial condition for WholeMesh
                     ##
                     ## Only specify one condition if not using mesh regions.
                     ## Otherwise select other initial_condition option, specify region_ids
                     ## and distinct names.  Then add extra intial conditions for other regions.
                     element initial_condition {
                        attribute name { "WholeMesh" },
                        input_choice_initial_condition_real
                     }|
                     ## Multiple initial_conditions are allowed if specifying
                     ## different values in different
                     ## regions of the mesh (defined by region_ids).  In this case
                     ## each initial_condition
                     ## requires a distinct name for the options dictionary.
                     element initial_condition {
                        attribute name { string },
                        region_ids?,
                        input_choice_initial_condition_real
                     }
                  )+,
                  diagnostic_scalar_field
               }
            },
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Moment_3" },
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  (
                     ## Initial condition for WholeMesh
                     ##
                     ## Only specify one condition if not using mesh regions.
                     ## Otherwise select other initial_condition option, specify region_ids
                     ## and distinct names.  Then add extra intial conditions for other regions.
                     element initial_condition {
                        attribute name { "WholeMesh" },
                        input_choice_initial_condition_real
                     }|
                     ## Multiple initial_conditions are allowed if specifying
                     ## different values in different
                     ## regions of the mesh (defined by region_ids).  In this case
                     ## each initial_condition
                     ## requires a distinct name for the options dictionary.
                     element initial_condition {
                        attribute name { string },
                        region_ids?,
                        input_choice_initial_condition_real
                     }
                  )+,
                  diagnostic_scalar_field
               }
            },
            element scalar_field {
               attribute rank { "0" },
               attribute name { xsd:string },
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  (
                     ## Initial condition for WholeMesh
                     ##
                     ## Only specify one condition if not using mesh regions.
                     ## Otherwise select other initial_condition option, specify region_ids
                     ## and distinct names.  Then add extra intial conditions for other regions.
                     element initial_condition {
                        attribute name { "WholeMesh" },
                        input_choice_initial_condition_real
                     }|
                     ## Multiple initial_conditions are allowed if specifying
                     ## different values in different
                     ## regions of the mesh (defined by region_ids).  In this case
                     ## each initial_condition
                     ## requires a distinct name for the options dictionary.
                     element initial_condition {
                        attribute name { string },
                        region_ids?,
                        input_choice_initial_condition_real
                     }
                  )+,
                  diagnostic_scalar_field
               }
            }*
         },
         element statistics {
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Mean" },
               ## Field type
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            }?,
            element scalar_field {
               attribute rank { "0" },
               attribute name { "StandardDeviation" },
               ## Field type
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            }?,
            element scalar_field {
               attribute rank { "0" },
               attribute name { "Skew" },
               ## Field type
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            }?,
            element scalar_field {
               attribute rank { "0" },
               attribute name { "SauterMeanDia" },
               ## Field type
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            }?,
            element scalar_field {
               attribute rank { "0" },
               attribute name { "MeanDia10" },
               ## Field type
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            }?
         },   

         ## Source terms for the population balance equation like growth, 
         ## breakage, aggegation and internal dispersion
         element population_balance_source_terms {
            ## Growth term in the internal coordinate space
            element growth{ 
               ## Enter the growth degree (r) for the power-law growth G(\xi)=\xi^{r}
               element power_law_growth{ real }
            }?,
            ## Dispersion coefficient D. Currently it only works for constant dispersion and can't handle dispersion varying with internal coordinate
            element internal_dispersion{ real }?,
            ##
            element aggregation{
               ##
               element aggregation_frequency{
                  ##
                  element constant_aggregation{ real }|
                  ## aggregation frequency is given by \beta_{alpha gamma} = \xi_{alpha}^3 + \xi_{gamma}^3
                  element hydrodynamic_aggregation{ empty }|
                  ## aggregation frequency is given by \beta_{alpha gamma} = \beta_o (\xi_{alpha} + \xi_{gamma}). Enter the value of constant \beta_o below
                  element sum_aggregation{ real }|
                  ## 
                  element laakkonen_2007_aggregation{ 
                     ## enter the value of the coefficient C5 in Laakkonen's aggregation rate expression. default: 0.88  
                     element C5{ real }? 
                  }
               }
            }?,
            ##
            element breakage{
               ##
               element breakage_frequency{
                  ##
                  element constant_breakage{ real }|
                  ## breakage frequency is given by a_{alpha} = p * \xi_{alpha}^r
                  element power_law_breakage{
                     ## coefficient p in breakage frequency a_{alpha} = p * \xi_{alpha}^r
                     element coefficient{ real },
                     ## degree r in breakage frequency a_{alpha} = p * \xi_{alpha}^r
                     element degree{ real }
                  }|
                  element laakkonen_breakage{ empty }
               },
               element distribution_function{
                  ##
                  element symmetric_fragmentation{ empty }|
                  ## McCoy and Madras (2003) used P(v|v') = 2/v' for v<v' and 0 otherwise, where v stands for volume. 
                  ## When written in terms of particle diameter, this converts to P(x|x')=(6x^2)/x'^3 for x<x' and 0 otherwise.
                  element mcCoy_madras_2003{ empty } |
                  element laakkonen_2007{ empty }
               }
            }?
         },

         ## At each vertex a linear system is solved to determine source terms
         ## for the abscissa and weight transport equations. This requires inversion
         ## of a matrix that has terms that are a function of the weighted abscissa.
         ## If any two abscissa happen to lie very close to one another this can cause
         ## this matrix to be ill-conditioned which can cause large errors or solver
         ## failures when attempting to invert the matrix. There several options for 
         ## evading this problem:
         ##
         ## (a) Set the abscissa and weight source terms to zero at the offending vertex.
         ## (b) Perturbate the weighted absciss used to construct the matrix until the 
         ## matrix is no longer ill conditioned and can be inverted.
         ## (c) Use the source terms of neighbouring cells (not implemented)
         element ill_conditioned_matrices {
            ## The required condition number of the matrix. Recommended value of 1e-12
            element required_condition_number{ real },
            (
               ## This is a better option when much of the domain has a minimum 0'th moment
               element set_source_to_zero { empty }|
               element perturbate{
                  ## The amount to perturbate the weighted abscissa. The values are 
                  ## iteratively perturbated by rand(-1,1)*perturbation until the matrix
                  ## exceeds the required condition number. Too low a value can cause
                  ## very long run times but a higher value will lead to larger errors.
                  ## After every ten perturbations at any one grid location a message is 
                  ## printed to the log file noting the difficulties in solving the linear 
                  ## system.
                  ## Recommended value of 1e-7
                  element perturbation{ real }
               }|
               ## Do not check for ill-conditioning of matrix A
               element do_nothing{ empty }
            )
         },
         ## weights of zero can cause instabilities. It is recommended that a minimum
         ## weight of approx 1e-10 is used.
         element minimum_weight{ real },
         ## In case the user wants to apply a minimum cap for weighted abscissas. Not Recommended. 
         element minimum_weighted_abscissa{ real }?,
         ## Scaling factor to be used for the Sauter Mean Diameter / Mean Diameter calculated in the 
         ## statistics. It is useful in a non-dimensional implementation of the PBE.
         element scaling_factor_Dia{ real }?,
         ## This is the submerged specific gravity, R, of this population balance
         ## as a function of internal coordinate, Z. 
         ## It will be used in the equation of state.
         ## 
         ## R(Z) = (rho_c(Z) - rho_a)/(rho_a)
         ## 
         ## Where: rho_c(Z) is the density of the dispersed phase and rho_a is 
         ## the ambient fluid density
         ##
         ## Must be given as a python function in the form:
         ## 
         ##  def R(Z): 
         ##     # Function code
         ##     return # Return value of R
         element submerged_specific_gravity {
            python_code
         }?,
         element output_abscissa_and_weights{ empty }?,
         element output_sources{ empty }?,
         ## This will apply all source terms for the weights and weighted_abscissa scalar 
         ## fields as absorption terms.
         element apply_source_as_absorption{ 
            ## For including a sponge region for weights and weighted_abscissas (which are 
            ## applied as absorption terms) a sponge field is defined which will be analogous to 
            ## an absorption term. This field is zero everywhere except where a large absorption
            ## is present. 
            element include_sponge_region{
               ## Name of the scalar field which is the sponge scalar field. Note: This sponge field must
               ## be defined the the same material phase.
               element sponge_scalar_field_name{ string }
            }?         
         }?,
         ## This will apply selected source terms for the weights and weighted_abscissa scalar 
         ## fields as absorption terms wherever the source is negative.
         element apply_source_as_absorption_for_negative_source_only{
            ## For including a sponge region for weights and weighted_abscissas (which are 
            ## applied as absorption terms) a sponge field is defined which will be analogous to 
            ## an absorption term. This field is zero everywhere except where a large absorption
            ## is present. 
            element include_sponge_region{
               ## Name of the scalar field which is the sponge scalar field. Note: This sponge field must
               ## be defined the the same material phase.
               element sponge_scalar_field_name{ string }
            }?
         }?
      }
   )

sediment = 
   (
      ## A sediment model. See the manual for details on how to use this model.
      element sediment {
         ## A single sediment field with discrete characteristics
         ## 
         ## Notes: 
         ##
         ## - a sinking velocity is required for sediment fields, this can be set for each
         ## sediment field under prognostic/SinkingVelocity
         ##
         ## - continuity must be the same for all fields. i.e. use the same mesh for all fields
         ##
         ## - rentrainment bc's can be used to allow sediment to be picked up from the bed
         ## due to turbulence. This must be set on the same boundaries as the SedimentDepositon
         ## field
         ## 
         ## - see the manual for more details on how to use this model 
         element scalar_field {
            attribute rank { "0" },
            attribute name { string },
            ## Field type
            element prognostic {
               velocity_mesh_choice,
               prognostic_scalar_field,
               ## Sediment bedload diagnostic which records the sediment 
               ## deposited through the prescribed boundary.
               ##
               ## The units are in unit distances as a depth of sediment.
               ##
               ## - The equation type must be set to SedimentBedload
               ##
               ## - Spatial discretisation should be as the parent sediment class
               ##
               ## - Initial conditions can be set as normal
               ##
               ## - Most other settings are irrelevant
               ##     - No solver settings are required as the sediment model only works
               ##       for p1 discretisations where the mass matrix can be easily inverted.
               ##     - The bedload is always calculated completely explicitly
               ##     - You cannot apply Source or Absorbtion terms to the bedload
               ##     - There is no diffusion of the bedload
               element scalar_field {
                  attribute rank { "0" },
                  attribute name { "Bedload" },
                  ## Field type
                  (
                     element prognostic {
                        ## Surface ids over which to calculate the bedload:
                        element surface_ids {
                           integer_vector
                        },
                        prognostic_scalar_field,
                        ## Disables calculation of change to bedload due to 
                        ## erosion/deposition (generally only used for mms testing):
                        element disable_calculation { empty }?
                     }
                  )
               },
               ## Sediment bedload diagnostic which records the sediment 
               ## deposited through the prescribed boundary.
               ##
               ## The units are in unit distances/second.
               element scalar_field {
                  attribute rank { "0" },
                  attribute name { "BedloadVolumeFraction" },
                  ## Field type
                  element diagnostic {
                     diagnostic_scalar_field
                  }
               },
               ## Sediment bedload deposit rate diagnostic.
               ##
               ## The units are in unit distances/second.
               element scalar_field {
                  attribute rank { "0" },
                  attribute name { "BedloadDepositRate" },
                  ## Field type
                  element diagnostic {
                     diagnostic_scalar_field
                  }
               }?,
               ## Sediment bedload erosion rate diagnostic.
               ##
               ## The units are in unit distances as a depth of sediment.
               element scalar_field {
                  attribute rank { "0" },
                  attribute name { "BedloadErosionRate" },
                  ## Field type
                  element diagnostic {
                     diagnostic_scalar_field
                  }
               }?,
               ## This is the unhindered sinking velocity of the sediment particles. 
               ## 
               ## With the sinking velocity field set to diagnostic, this field is used, along with the sediment 
               ## concentration, to calculate the hindered sinking velocity of sediment particles. 
               ##
               ## This uses the equation developed by Richardson and Zaki [1954]. 
               ##
               ## If the 'SinkingVelocity' is not set to diagnostic, this field will not be used.
               ##
               ## As with the sinking velocity, this velocity is in the direction of gravity so if the substance
               ## floats upwards, this field should be negative.
               ##
               ## This should be on the same mesh as the associated SinkingVelocity field.
               element scalar_field {
                  attribute name { "UnhinderedSinkingVelocity" },
                  attribute rank { "0" },
                  element prescribed {
                     velocity_mesh_choice,
                     prescribed_scalar_field
                  }
               }?,
               ## This is the submerged specific gravity, R, of this sediment. 
               ## It will be used with the sediment concentration in the equation
               ## of state, as well as in erosion algorithms.
               ## 
               ## R = (rho_s - rho_a)/(rho_a)
               ##
               ## Where: rho_s is the sediment density and rho_a is the ambient fluid density
               element submerged_specific_gravity {
                  real
               },
               ## This is the diameter of the grain.
               ## Diameter must be specified here or under each sediment class.
               element diameter {
                  real
               }?,
               ## Set porosity to use. Default is 0.3. 1.0 is a solid material.
               element bed_porosity {
                  real
               }?,
               ## Erodability of the sediment grain. A value of one means that only the
               ## critical shear stress is used to determine if a grain can be put into 
               ## suspension. A value of zero means these grains can never be resuspended.
               ## Default is 1.
               element erodability {
                  real
               }?,
               ## Critical shear stress of a grain. If not switched on a value
               ## will be computed according to the particle Reynolds number and Shield's
               ## criterion.
               element critical_shear_stress {
                  real
               }?
            }
         }*,
         ## The median sediment diameter in the active layer of the bed. 
         ## Required for some reentrainment algorithms.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "SedimentBedActiveLayerD50" },
            (
               element diagnostic {
                  ## Surface ids over which to calculate the bedload:
                  element surface_ids {
                     integer_vector
                  },
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            )
         },
         ## The standard deviation of sediment in the active layer of the
         ## bed. Required for some reentrainment algorithms
         element scalar_field {
            attribute rank { "0" },
            attribute name { "SedimentBedActiveLayerSigma" },
            (
               element diagnostic {
                  ## Surface ids over which to calculate the bedload:
                  element surface_ids {
                     integer_vector
                  },
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            )
         },
         ## ZeroSedimentConcentrationViscosity field:
         ##
         ## Field for the viscosity of this a fluid with zero 
         ## concentration of sediment.
         ## Required if using a diagnostic viscosity
         ## in a sediment problem where sediment concentration
         ## dependent viscosity is required.
         element tensor_field {
            attribute rank { "2" },
            attribute name { "ZeroSedimentConcentrationViscosity" },
            (
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_tensor_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }?
      }
   )

prognostic_density_field =
   (
      scalar_equation_choice?,
      ## Spatial discretisation options
      element spatial_discretisation {
         (
            ## Continuous Galerkin formulation.
            element continuous_galerkin {
               advection_stabilisation_options,
               ## Discretisation options for the advection terms.
               element advection_terms {
                  ## Integrate the advection terms of the momentum equation by
                  ## parts.
                  element integrate_advection_by_parts {
                     comment
                  }?
               },
               ## Discretisation options for the mass terms.
               element mass_terms {
                  ## Lump the mass matrix
                  element lump_mass_matrix {
                     empty
                  }?
               },
               comment
            }|
            ## Discontinuous galerkin formulation
            discontinuous_galerkin_options|
            ## Use a control volume discretisation.
            element control_volumes {
              spatial_control_volume_options
            }
         ),
         ## Conservative discretisation of field advection equation
         ##  TBETA=1. -- conservative (divergence form)
         ##  TBETA=0. -- non-conservative
         ##  0. < TBETA < 1.
         element conservative_advection {
           real
         }
      },
      element temporal_discretisation {
         ## Implicit/explicit control (TTHETA)
         ##  =0.  -- explicit
         ##  =0.5 -- Crank-Nicholson
         ##  =1.  -- implicit
         element theta {
            real
         },
         temporal_control_volume_options?
      },
      (
         ## Solver
         element solver {
            linear_solver_options_asym_scalar
         }|
         ## Assume this field is being solved explicitly and skip the solver.
         ##
         ## ONLY AVAILABLE FOR PURE CONTROL VOLUME SPATIAL DISCRETISATIONS.
         ##
         ## Assumes lhs matrix only has diagonal lumped mass (times
         ## density if appropriate for equation)
         ## and divides the rhs by this.
         element explicit {
            empty
         }
      )?,
      (
         ## Initial condition for WholeMesh
         ##
         ## Only specify one condition if not using mesh regions.
         ## Otherwise select other initial_condition option, specify region_ids
         ## and distinct names.  Then add extra intial conditions for other regions.
         element initial_condition {
            attribute name { "WholeMesh" },
            input_choice_initial_condition_real
         }|
         ## Multiple initial_conditions are allowed if specifying
         ## different values in different
         ## regions of the mesh (defined by region_ids).  In this case
         ## each initial_condition
         ## requires a distinct name for the options dictionary.
         element initial_condition {
            attribute name { string },
            region_ids?,
            input_choice_initial_condition_real
         }
      )+,
      ## Boundary conditions
      element boundary_conditions {
         attribute name { string },
         ## Surface id:
         element surface_ids {
            integer_vector
         },
         ## Type
         (
            element type {
               attribute name { "dirichlet" },
               ## Apply the dirichlet bc weakly.  Available
               ## automatically with discontinuous_galerkin,
               ## and control_volume spatial_discretisations.
               ## If not selected boundary conditions are applied strongly.
               element apply_weakly {
                  ## If the initial condition and boundary conditions
                  ## differ, setting this option will cause the initial
                  ## condition on the boundary to be overwritten with
                  ## the boundary condition. Since you are applying the
                  ## boundary condition weakly, you probably do *not*
                  ## want this.
                  element boundary_overwrites_initial_condition {
                     empty
                  }?
               }?,
               input_choice_real
            }
         )
      }*,
      ## source term
      element scalar_field {
         attribute name { "Source" },
         attribute rank { "0" },
         (
            element prescribed {
               prescribed_scalar_field_no_adapt,
               recalculation_options?
            }|
            element diagnostic {
              internal_algorithm,
              diagnostic_scalar_field_no_adapt
            } 
         )
      }?,
      ## Absorption term
      element scalar_field {
         attribute name { "Absorption" },
         attribute rank { "0" },
         element prescribed {
            prescribed_scalar_field_no_adapt
         }
      }?,
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options,
      scalar_convergence_options,
      prognostic_detector_options,
      prognostic_particle_options,
      scalar_steady_state_options,
      adaptivity_options_prognostic_scalar_field,
      interpolation_algorithm_scalar,
      discrete_properties_algorithm_scalar?,
      ## Set the priority of this field
      ## This determines the order in which scalar_fields are solved for:
      ##  - higher numbers have the highest priority
      ##  - lower numbers (including negative) have the lowest priority
      ##  - default if not set is 0
      element priority {
         integer
      }?
   )

geometry = 
   (
      ## Options dealing with the specification of geometry
      element geometry {
         ## Dimension of the problem.
         ## <b>This can only be set once</b>
         element dimension {
            element integer_value {
               attribute rank {"0"},
               ("3"|"2"|"1")
            }
         },
         ## The position mesh
         element mesh {
            attribute name { "CoordinateMesh" },
            mesh_info
         },
         ## The velocity mesh
         element mesh {
            attribute name { "VelocityMesh" },
            mesh_info
         }?,
         ## The pressure mesh
         element mesh {
            attribute name { "PressureMesh" },
            mesh_info
         }?,
         element mesh {
            attribute name { xsd:string },
            mesh_info,
            element exclude_from_mesh_adaptivity{empty}?
         }*,
         ## Quadrature
         element quadrature {
            ## Quadrature degree
            ## 
            ## note: this specifies the degree of quadrature,
            ## not the number of gauss points
            element degree {
               integer
            },
            ## Surface quadrature degree
            ## 
            ## note: this specifies the degree of surface
            ## quadrature not the number of surface gauss points
            element surface_degree {
               integer
            }?,
            ## Sets the degree of quadrature on each quadrilateral
            ## face of the control volume. 
            ##
            ## Defaults to 1 if
            ## unselected which is the same as pre-new options
            ## behaviour.
            element controlvolume_surface_degree {
               integer
            }?,
            ## Select which family of quadrature rules to use.
            ## The default is family_cools.
            ## family_wandzura allows for degree up to 30
            ## on triangular meshes.
            ## family_grundmann_moeller allows for degree up to
            ## 29 on simplicial meshes in arbitrary dimension.
            element quadrature_family {
               element string_value {
                  ( "family_cools" | "family_grundmann_moeller" | "family_wandzura" ) 
               }
            }?
         },
         ## This causes the change of variables associated with each element in 
         ## the mesh to be stored rather than calculated every time it is used. 
         ## This should speed up computations at a cost of some memory. 
         ## 
         ## The cache is automatically regenerated after mesh movement or 
         ## adaptivity and is automatically disabled for non-linear positions fields.
         element disable_geometric_data_cache {
            empty
         }?,
         ## Options specifying that the problem is on the surface of the n-sphere.
         element spherical_earth {
            ## Enabling this option approximates the curvature of the Earth as an
            ## nth degree polynomial, where n is the polynomial degree of the coordinate mesh.
            element superparametric_mapping {
               empty
            }?,
            ## The simplicial elements (as defined by the CoordinateMesh) are bent on the n-sphere
            ## using an analytical expression (as opposed to superparametric where it's approximated
            ## by a higher order polynomial), the derivatives are calculated by taking the analytical
            ## derivative on the gauss points.
            element analytical_mapping {
               empty
            }?
         }?,
         ## Options specifying the top surface and bottom of the domain
         ## used in various ocean calculations.
         element ocean_boundaries {
            ## Specify the surface ids that make up the top of the domain,
            ## i.e. the free surface or rigid lid.
            element top_surface_ids {
               integer_vector
            },
            ## Specify the surface ids that make up the bottom.
            element bottom_surface_ids {
               integer_vector
            },
            ## Diagnostic field giving the distance to the top surface.
            element scalar_field {
               attribute rank { "0" },
               attribute name { "DistanceToTop" },
               element diagnostic {
                  internal_algorithm,
                  element mesh {
                     attribute name {  "CoordinateMesh" }
                  },
                  diagnostic_scalar_field
               }
            },
            ## Diagnostic field giving the distance to ocean bottom.
            element scalar_field {
               attribute rank { "0" },
               attribute name { "DistanceToBottom" },
               element diagnostic {
                  internal_algorithm,
                  element mesh {
                     attribute name {  "CoordinateMesh" }
                  },
                  diagnostic_scalar_field
               }
            }
         }?
      }
   )

lagrangian_timestepping = 
   (
      element lagrangian_timestepping {
         ## Number of subdivisions of the timestep
         ## increase this if you are not happy with your 
         ## detector trajectory accuracy, or if particles
         ## are jumping out of the domain a lot
         element subcycles {
            integer
         },
         ## Tolerance for deciding if detector is in a given
         ## element. Recommended value 1.0e-10.
         element search_tolerance {
            real
         },
         (
            ## Use explicit runge kutta method with
            ## guided search particle tracking
            element explicit_runge_kutta_guided_search {
               ## Number of RK stages
               ## For the RK4 method, it should be 4.
               element n_stages {
                  integer
               },
               ## ERK stage array. This is an array
               ## containing the lower-triangular
               ## part of the Butcher weight matrix
               ## A that explains how to compute the
               ## RK stages.  See
               ## http://en.wikipedia.org/wiki/Runge–Kutta_methods#Explicit_Runge.E2.80.93Kutta_methods
               ## for notation.  The array is stored
               ## in the following order:
               ## [a_{21},a_{31},a_{32},...,a_{s1},a_{s2},a_{s(s-1)}]
               ## and so the array has size s(s-1)/2
               ## where s is the number of stages.
               ## For the RK4 method, it should be
               ## [0.5,0,0.5,0,0,1]
               element stage_weights {
                  real_vector
               },
               ## ERK timestep weights. This is the
               ## b vector that explains how to
               ## compute the timestep from the RK
               ## stages.  See
               ## http://en.wikipedia.org/wiki/Runge–Kutta_methods#Explicit_Runge.E2.80.93Kutta_methods
               ## for notation.  It should have size
               ## s where s is the number of stages.
               ## For the RK4 method, it should be
               ## [1/6,1/3,1/3,1/6]
               element timestep_weights {
                  real_vector
               }
            }|
            ## Use explicit Forward Euler method with
            ## guided search particle tracking
            element forward_euler_guided_search {
               empty
            }|
            ## Use classical Runge-Kutta method with
            ## guided search particle tracking
            element rk4_guided_search {
               empty
            }
         )
      }
   )

# Default child of diagnostic scalar field
diagnostic_scalar_field =
   (
      diagnostic_output_options,
      diagnostic_scalar_stat_options,
      scalar_convergence_options,
      diagnostic_detector_options,
      diagnostic_particle_options,
      scalar_steady_state_options,
      adaptivity_options_scalar_field,
      recalculation_options?,
      interpolation_algorithm_scalar?
   )
   
# Default child of diagnostic scalar field without adaptivity options
diagnostic_scalar_field_no_adapt =
   (
      diagnostic_output_options,
      diagnostic_scalar_stat_options,
      diagnostic_detector_options,
      diagnostic_particle_options
   )

# Default child of diagnostic scalar field
diagnostic_scalar_field_tidal_range =
   (
      diagnostic_output_options,
      diagnostic_scalar_stat_options,
      diagnostic_detector_options,
      adaptivity_options_scalar_field,
      (
          element spin_up_time {
             real
          }
      )       
   )

# Default child of diagnostic vector field
# Currently, this is empty, but in future this might include
# options that are general to all diagnostic vector fields
diagnostic_vector_field =
   (
      diagnostic_output_options,
      diagnostic_vector_stat_options,
      vector_convergence_options,      
      diagnostic_detector_options,
      vector_steady_state_options,      
      adaptivity_options_vector_field,
      recalculation_options?,
      interpolation_algorithm_vector?
   )

diagnostic_vector_field_bed_shear_stress =
   (
      element density {
         real
      }, 
      (
         ## Three options for calculation of bed shear stress are available:
         element calculation_method {
            ## density*drag_coeff*|u|*u 
            ## This is calculated at all boundaries
            element drag_coefficient {
               real
            }|
            ## nu * (grad(u) + grad(u).T)
            ## This uses the velocity gradient at the boundaries to calculate the 
            ## bed shear stress. Valid for non parameterised velocity boundaries.
            element velocity_gradient {    
               empty
            }
         }
      ),
      diagnostic_output_options,
      diagnostic_vector_stat_options,
      diagnostic_detector_options,
      adaptivity_options_vector_field
   )


# Default child of diagnostic tensor field
# Currently, this is empty, but in future this might include
# options that are general to all diagnostic tensor fields
diagnostic_tensor_field =
   (
      diagnostic_output_options,
      diagnostic_tensor_stat_options,
      adaptivity_options_tensor_field,
      interpolation_algorithm_vector?
   )
   
# Richardson number field. This is a normal diagnostic scalar field, but with
# Richardson number metric options added
adaptivity_options_richardson_number_field.adaptivity_options =
   (
      ## Do not use an interpolation error driven metric for this field
      element no_interpolation_measure {
        comment
      }
   )
adaptivity_options_richardson_number_field.adaptivity_options |= adaptivity_options_scalar_field.adaptivity_options
adaptivity_options_richardson_number_field =
   (
      element adaptivity_options {
         adaptivity_options_richardson_number_field.adaptivity_options,
         ## An isotropic metric formulation based on the Richardson number. Uses
         ## the logic that wherever the Richardson number is small, we expect
         ## to need resolution. Generates edge lengths using:
         ##
         ##   Edge length = min_edge_length if Ri <= min_ri
         ##                 max_edge_length if Ri >= max_ri
         ##                 a linear fit between min_edge_length and max_edge_length otherwise
         element richardson_number_metric {
            ## Richardson number at which we have minimum edge length (default 0.0)
            element min_ri {
              real
            }?,
            ## Richardson number at which we have maximum edge length
            element max_ri {
              real
            },
            ## Minimum edge length that can be requested by the Richardson number
            ## metric
            element min_edge_length {
              real
            },
            ## Maximum edge length that can be requested by the Richardson number
            ## metric
            element max_edge_length {
              real
            },
            ## Enable to preserve anisotropy when merging with other metric
            ## formulations
            element anisotropy_preserving_merge {
               comment
            }?,
            comment
         }?,
         adaptivity_preprocessing
      }?
   )
diagnostic_richardson_number_field = diagnostic_scalar_field_no_adapt
diagnostic_richardson_number_field &= adaptivity_options_richardson_number_field

diagnostic_cv_gradient_vector_field =
   (
      ## Choose whether the mass matrix is lumped or not
      element lump_mass_matrix {
            empty
      }?,
      ## Solver options are necessary if you're not lumping your mass or if you're field isn't dg
      element solver {
         linear_solver_options_sym
      }?,
      ## Normalise the gradient by its magnitude
      element normalise {
        empty
      }?,
      diagnostic_output_options,
      diagnostic_vector_stat_options,
      vector_convergence_options,
      diagnostic_detector_options,
      vector_steady_state_options,
      adaptivity_options_vector_field,
      recalculation_options?
   )

diagnostic_gradient_vector_field =
   (
      ## Solver
      element solver {
         linear_solver_options_sym
      },
      diagnostic_output_options,
      diagnostic_vector_stat_options,
      vector_convergence_options,
      diagnostic_detector_options,
      vector_steady_state_options,
      adaptivity_options_vector_field,
      recalculation_options?
   )

diagnostic_cv_divergence_scalar_field =
   (
      # No solver options because it can be solved directly!
      diagnostic_output_options,
      diagnostic_scalar_stat_options,
      scalar_convergence_options,
      diagnostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_scalar_field,
      recalculation_options?
   )

diagnostic_fe_divergence_scalar_field =
   (
      ## Solver
      element solver {
         linear_solver_options_sym
      },
      diagnostic_output_options,
      diagnostic_scalar_stat_options,
      scalar_convergence_options,
      diagnostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_scalar_field,
      recalculation_options?
   )

# three optional input vectors for user-specified rotation matrix
rotation_matrix_components =
   (
      ## Select if you want to specify the unit normal direction
      ## of the rotation matrix.
      ## If off then fluidity computes the normal.  
      element normal_direction {
         input_choice_real_dim_vector
      }?,
      ## Select if you want to specify the first unit tangent direction
      ## of the rotation matrix (in dim > 1 simulations).
      ## If off then fluidity computes the tangent.  
      element tangent_direction_1 {
         input_choice_real_dim_vector
      }?,
      ## Select if you want to specify the second unit tangent direction
      ## of the rotation matrix (in dim > 2 simulations).
      ## If off then fluidity computes the tangent.  
      element tangent_direction_2 {
         input_choice_real_dim_vector
      }?
   )

velocity_components_choice =
   (
      (
         element align_bc_with_cartesian {
            element x_component {
               input_choice_real_bc_component,
               bc_smoothing
            }?,
            element y_component {
               input_choice_real_bc_component,
               bc_smoothing
            }?,
            element z_component {
               input_choice_real_bc_component,
               bc_smoothing
            }?
         }|
         element align_bc_with_surface {
            element normal_component {
               input_choice_real_plus_field,
               bc_smoothing
            }?,
            element tangent_component_1 {
               input_choice_real_plus_field,
               bc_smoothing
            }?,
            element tangent_component_2 {
               input_choice_real_plus_field,
               bc_smoothing
            }?,
            rotation_matrix_components,
            ## this will calculate the determinant of the
            ## rotation matrix for every boundary node
            ## and dump a vtu with the node 
            ## normals and tangenials 1/2
            element debugging_mode{empty}?
         }
      )
   )

# Mainly tidal harmonic diagnostics
free_surface_diagnostic_options = 
   (
      (
         ## Activate this for tidal harmonic analysis 
          free_surface_history_algorithm?,
         ## Activate this for tidal harmonic analysis 
          tidal_harmonic_algorithm?
      )
   )

velocity_boundary_conditions =
   (
      (
         element type {
            attribute name { "dirichlet" },
            ## Apply the dirichlet bc weakly.  Available automatically
            ## with a discontinuous_galerkin Velocity
            ## spatial_discretisation.  Available if you
            ## integrate_continuity_by_parts with a
            ## continuous_galerkin Pressure or use a control_volume
            ## Pressure spatial_discretisation and/or
            ## integrate_advection_by_parts under Velocity
            ## spatial_discretisation with continuous_galerkin.
            ##
            ## If not selected boundary conditions are applied strongly.
            element apply_weakly {
               ## If the initial condition and boundary conditions
               ## differ, setting this option will cause the initial
               ## condition on the boundary to be overwritten with
               ## the boundary condition. Since you are applying the
               ## boundary condition weakly, you probably do *not*
               ## want this.
               element boundary_overwrites_initial_condition {
                  empty
               }?
            }?,
            velocity_components_choice
         }|
         element type {
            attribute name { "neumann" },
            velocity_components_choice
         }|
         ## Add a bulk formulae boundary condition. Only makes sense
         ## on the Velocity field.
         element type {
            attribute name { "bulk_formulae" },
               empty
         }|
         element type {
           attribute name { "free_surface" },
           (
              ## The default free_surface implements the physical boundary
              ## condition p=0 (this is full pressure without subtracting the hydrostatic component)
              ## For viscous fluids the correct boundary condition is a no normal
              ## stress condition, which includes a viscosity term. When using
              ## this option a prognostic FreeSurface field is required.
              element no_normal_stress {
                 ## Treat the free surface evolution equation separately from the momentum equation.  
                 ## i.e. stagger the solves.
                 ##
                 ## Requires a prognostic FreeSurface field.
                 element explicit {
                   empty
                 }?,
                 ## Force the normals used in the free surface boundary condition to be purely in the
                 ## radial direction (up or down), instead of using the facet normals. This option is
                 ## recommended in idealised cylindrical and spherical domains (only) when a surface-aligned
                 ## strong Dirichlet boundary condition is used to prescribe the tangential velocities. In
                 ## this case we need the normals to be orthogonal to the prescribed tangential direction.
                 element radial_normals {
                   empty
                 }?
              }?,
              ## When selected the density at the surface is calculated using the Density field.  Therefore
              ## this requires a Density field!
              element variable_density {
                 empty
              }?,
              ## When selected a hydrostatic pressure distribution on the outside is assumed
              ## using the specified density, i.e. instead of applying p=0 (resp.
              ## normal_stress=0) we apply p=p_external (resp. normal_stress=p_external) where p_external is determined by dp_external/dz=external_density*g
              element external_density {
                 input_choice_real_plus_field
              }?,
              ## This options adds a surface stabilisation term to the free surface. Works only for cg velocity so far. 
              ## Note: Once activated, the stabilisation term will occur in all free surface areas in the domain. 
              ## IN DEVELOPMENT
              element surface_stabilisation {
                  ## Scale factor for the surface stabilisation.
                  element scale_factor{
                    real
                  }
              }?
           )
         }|
         ## Apply quadratic drag. Specify drag coefficient. If you
         ## want to exactly replicate results from using the OCEDRA
         ## option, set this to 0.003 and remember to apply to both
         ## bottom and sides.
         element type {
            attribute name { "drag" },
            input_choice_real,
            (
              ## Use a quadratic drag.
              ##
              ## This means that the drag coefficient is nondimensional.
              element quadratic_drag {
                ## Use the Manning-Strickler formulation:
                ## n^2*g*|u|*u/H^(1/3)
                ## where n is the Manning coefficient, g is gravity, u is the velocity vector and H is the water heigth at that point.
                ##
                ## The coefficient given above defines the Manning coefficient [s/m^(1/3)] (a typical value for sand is 0.02)
                element manning-strickler {
                   empty
                }?
              }|
              ## Use a linear drag (basically just a surface absorption term).
              ##
              ## This means that the drag coefficient has units of momentum.
              element linear_drag {
                empty
              }
            )
         }|


         ## Apply wind forcing specified by stress or wind velocity.
         ## Replaces windy.dat and windy.py
         element type {
            attribute name { "wind_forcing" },
            (
               ## Wind forcing with user specified wind stress
               ##
               ## <b> Note that the stress needs to be specified
               ## using the same density units as the reference_density 
               ## under equation of state.</b>So if you use the recommended
               ## non-dimensional value of 1.0 for reference_density and
               ## your calculated stress is in kg m^-1s^-2 and the dimensional
               ## reference_density is 1000.0 kg m^-3, you need to divide
               ## the calculated stress in SI units by 1000.0.
               element wind_stress {
                  input_choice_real_dim_minus_one_vector|
                  element from_netcdf {
                     ## The format of this file should conform to NetCDF CF 1.x
                     ## (http://cf-pcmdi.llnl.gov/).
                     attribute file_name { xsd:string },
                     attribute east_west { xsd:string },
                     attribute north_south { xsd:string },
                     comment
                  }
               }|
               ## Wind forcing with user specified 10m wind velocity
               element wind_velocity {
                  ## Specify wind drag coefficient (dimensionless)
                  ## Suggested value: 4.0e-4
                  element wind_drag_coefficient {
                     input_choice_real
                  },
                  ## Density of air. 
                  ##
                  ## <b>Note that you have to specify
                  ## this density in the same units as the 
                  ## reference_density under equation of state.</b>
                  ## So with a typicial value of rho_air=1.3 kgm^-3
                  ## and rho_water=1000 kgm^-3, if you fill in the 
                  ## recommended (non-dimensional) value of 1.0 for 
                  ## reference_density, this field needs to be 1.3e-3.
                  element density_air {
                     real
                  },
                  ## Specify wind velocity
                  element wind_velocity {
                     input_choice_real_dim_minus_one_vector|
                     element from_netcdf {
                        ## The format of this file should conform to NetCDF CF 1.x
                        ## (http://cf-pcmdi.llnl.gov/)
                        attribute file_name { xsd:string },
                        attribute east_west { xsd:string },
                        attribute north_south { xsd:string },
                        comment
                     }
                  }
               }
            )
         }|

         ## When using control_volumes under Pressure
         ## spatial_discretisation or when using
         ## integrate_continuity_by_parts with continuous_galerkin
         ## Pressure and continuous_galerkin Velocity this
         ## boundary condition type imposes a weak no normal flow
         ## boundary condition on the surface_ids specified.
         element type {
            attribute name { "no_normal_flow" },
            empty
         }|
         ## Implements a prescribed normal flow. Works for DG only. This
         ## implements a weakly imposed bounadry condition normal to the surface.
         ## Note: Positive points out of the domain, negative into the domain.
         element type {
            attribute name { "prescribed_normal_flow" },
               input_choice_real_bc_component
         }|
         ## Flux boundary condition. 
         ##
         ## Weakly enforces a flux across a boundary
         element type {
            attribute name { "flux" },
            element align_bc_with_cartesian {
               element x_component {
                  input_choice_real_bc_component
               }?,
               element y_component {
                  input_choice_real_bc_component
               }?,
               element z_component {
                  input_choice_real_bc_component
               }?
            }
         }
      )
   )
   
# Output options for prognostic fields
prognostic_scalar_output_options =
   (
      ## Specify what is written to vtu dump files.
      element output {
         ## Exclude this field from dump files.
         element exclude_from_vtu {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## on the previous timestep.
         ## (included under the name: Old<field_name> )
         element include_previous_time_step {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## used in the nonlinear iteration.
         ## (included under the name: Nonlinear<field_name> )
         element include_nonlinear_field {
            empty
         }?,
         ## Output a file details the convergence (or otherwise) of
         ## this field with every advective nonlinear
         ## iteration.
         ## ONLY WORKS FOR PURE CONTROL VOLUME DISCRETISATIONS.
         element convergence_file {
            comment
         }?
      }
   )

# Output options for pressure (can't have a convergence file)
pressure_output_options =
   (
      ## Specify what is written to vtu dump files.
      element output {
         ## Exclude this field from dump files.
         element exclude_from_vtu {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## on the previous timestep.
         ## (included under the name: Old<field_name> )
         element include_previous_time_step {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## used in the nonlinear iteration.
         ## (included under the name: Nonlinear<field_name> )
         element include_nonlinear_field {
            empty
         }?,
         ## Write out some extra debugging vtu files that can be used
         ## to analyse what goes on in the pressure projection steps.
         ## WARNING: this may create a huge amount of vtu files, as 
         ## multiple files are written per nonlinear iteration.
         element debugging_vtus {
            empty
         }?
      }
   )

# Output options for prognostic fields
prognostic_vector_output_options =
   (
      ## Specify what is written to dump files.
      element output {
         ## Exclude this field from dump files.
         element exclude_from_vtu {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## on the previous timestep.
         ## (included under the name: Old<field_name> )
         element include_previous_time_step {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## used in the nonlinear iteration.
         ## (included under the name: Nonlinear<field_name> )
         element include_nonlinear_field {
            empty
         }?
      }
   )
   
# Field output options for all other fields
field_output_options =
   (
      ## Specify what is written to dump files.
      element output {
         ## Exclude this field from dump files.
         element exclude_from_vtu {
            comment
         }?,
         ## Field will be checkpointed (this is not performed by default).
         element checkpoint {
            empty
         }?
      }
   )
field_output_options_disabled =
   (
      ## Specify what is written to vtu dump files.
      element output {
         (
            ## Exclude this field from dump files.
            element exclude_from_vtu {
               comment
            }|
            ## Include this field in dump files.
            element include_in_vtu {
               comment
            }
         )
      }
   )
   
diagnostic_output_options = field_output_options
prescribed_output_options = field_output_options

# Options for inclusion/exclusion of standard field statistics from the .stat
# file
include_stat =
   (
      ## Include this field in the .stat file (magnitude and components)
      element include_in_stat {
         comment
      }
   )
exclude_components_from_stat =
   (
      ## Include just the magnitude of this field in the .stat file
      ## (excluding the components)
      element exclude_components_from_stat {
         comment
      }
   )
exclude_stat =
   (
      ## Exclude this field from the .stat file.
      element exclude_from_stat {
         comment
      }
   )

# Diagnostic statistics options for prognostic scalar fields
prognostic_scalar_stat_options = 
   (
      ## Specify what is added to .stat files
      element stat {
        prognostic_scalar_stat_options.stat
      }
   )
  
# Diagnostic statistics for all other scalar fields
prognostic_scalar_stat_no_old_or_nonlinear_options =
   (
      ## Specify what is added to .stat files
      element stat {
         prognostic_scalar_stat_no_old_or_nonlinear_options.stat
         
      }
   )

diagnostic_scalar_stat_options = prognostic_scalar_stat_no_old_or_nonlinear_options
prescribed_scalar_stat_options = prognostic_scalar_stat_no_old_or_nonlinear_options

# Combining of stat elements for scalar fields
prognostic_scalar_stat_options.stat = prognostic_scalar_stat_no_old_or_nonlinear_options.stat
prognostic_scalar_stat_options.stat &=
   (
      ## Enable to include the previous timestep value of this field in the .stat file.
      element include_previous_time_step {
         comment
      }?,
      ## Enable to include the values of this field in the nonlinear
      ## iteration in the .stat file.
      element include_nonlinear_field {
         comment
      }?
   )
prognostic_scalar_stat_no_old_or_nonlinear_options.stat = 
   (
      exclude_stat?,
      cv_stats?,
      surface_integral_stats_scalar*,
      mixing_stats*
   )   
   
# Diagnostic statistics options for vector fields, with enabled by default
vector_field_stat_options_enabled_default = include_stat
vector_field_stat_options_enabled_default |= exclude_components_from_stat
vector_field_stat_options_enabled_default |= exclude_stat

# Diagnostic statistics options for vector fields, with enabled by default
vector_field_stat_options_disabled_default = exclude_stat
vector_field_stat_options_disabled_default |= exclude_components_from_stat
vector_field_stat_options_disabled_default |= include_stat

# Diagnostic statistics options for tensor fields, with enabled by default
tensor_field_stat_options_enabled_default = include_stat
tensor_field_stat_options_enabled_default |= exclude_components_from_stat
tensor_field_stat_options_enabled_default |= exclude_stat

# Diagnostic statistics options for tensor fields, with enabled by default
tensor_field_stat_options_disabled_default = exclude_stat
tensor_field_stat_options_disabled_default |= exclude_components_from_stat
tensor_field_stat_options_disabled_default |= include_stat

# Diagnostic statistics for prognostic vector fields
prognostic_velocity_stat_options =
   (
      ## Specify what is added to .stat files
      element stat {
         (
            prognostic_velocity_stat_options.stat
         )
      }      
   )

# Diagnostic statistics for all other vector fields
prognostic_vector_stat_no_old_or_nonlinear_options =
   (
      ## Specify what is added to .stat files
      element stat {
         prognostic_vector_stat_no_old_or_nonlinear_options.stat
      }
   )
diagnostic_vector_stat_options = prognostic_vector_stat_no_old_or_nonlinear_options
prescribed_vector_stat_options = prognostic_vector_stat_no_old_or_nonlinear_options

# Diagnostic statistics for tensor fields
prognostic_tensor_stat_no_old_or_nonlinear_options =
   (
      ## Specify what is added to .stat files
      element stat {
         prognostic_tensor_stat_no_old_or_nonlinear_options.stat
      }
   )
diagnostic_tensor_stat_options = prognostic_tensor_stat_no_old_or_nonlinear_options
prescribed_tensor_stat_options = prognostic_tensor_stat_no_old_or_nonlinear_options

# Combining of stat elements for vector fields
prognostic_velocity_stat_options.stat = prognostic_vector_stat_no_old_or_nonlinear_options.stat
prognostic_velocity_stat_options.stat &=
   (
      ## Specify how the previous timestep value of this field is added to the .stat file.
      element previous_time_step {
         vector_field_stat_options_disabled_default
      },
      ## Specify how the values of this field used in the nonlinear iteration are added to the .stat file.
      element nonlinear_field {
         vector_field_stat_options_disabled_default
      },
      element compute_body_forces_on_surfaces {
        attribute name { xsd:string },
        (
          ## What surface IDs do you want to do the calculation over?
          element surface_ids {
            integer_vector
          },
          ## Enable to output the pressure and viscous terms separately (as well
          ## as the total force)
          element output_terms {
             comment
          }?
        )
      }*,
      ## Compute the divergence of this field at the Gauss points
      ## and return its stats.  This is a direct measure of the
      ## divergence at the gauss points rather than a discrete measure
      ## at the nodes (provided by several other diagnostic fields).
      element divergence_stats {
        empty
      }?,
      ## Calculate the error in the conservation of momentum
      ## IN PROGRESS - Does not include all terms!
      element calculate_momentum_conservation_error {
         empty
      }?
   )
prognostic_vector_stat_no_old_or_nonlinear_options.stat =
   (
      vector_field_stat_options_enabled_default,
      surface_integral_stats_vector*
   )

# Combining of stat elements for tensor fields
prognostic_tensor_stat_no_old_or_nonlinear_options.stat =
   (
      tensor_field_stat_options_enabled_default
   )

# Convergence options for prognostic scalar fields
scalar_convergence_options =
   (
      ## Decide whether this field is tested for convergence
      ## during nonlinear iterations
      ## (if /timestepping/nonlinear_iterations and
      ## /timestepping/nonlinear_iterations/tolerance are
      ## enabled).
      ## Also specifies whether the field is added to the 
      ## convergence file (if /io/convergence_file is enabled).
      element convergence {
         (
            ## Include this field in convergence testing
            ## (if /timestepping/nonlinear_iterations and
            ## /timestepping/nonlinear_iterations/tolerance are
            ## enabled) and file (if /io/convergence_file is enabled)
            element include_in_convergence {
               comment
            }|
            ## Exclude this field from convergence testing and file 
            element exclude_from_convergence {
               comment
            }
         )
      }
   )

# Convergence statistics options for prognostic vector fields (velocity)
vector_convergence_options =
   (
      ## Decide whether this field is tested for convergence
      ## during nonlinear iterations
      ## (if /timestepping/nonlinear_iterations and
      ## /timestepping/nonlinear_iterations/tolerance are
      ## enabled).
      ## Also specifies whether the field is added to the 
      ## convergence file (if /io/convergence_file is enabled).
      element convergence {
         (
            ## Include this field (magnitude and components)
            ## in convergence testing
            ## (if /timestepping/nonlinear_iterations and
            ## /timestepping/nonlinear_iterations/tolerance are
            ## enabled) and file (if /io/convergence_file is enabled)
            element include_in_convergence {
               comment
            }|
            ## Include just the magnitude of this field 
            ## in convergence testing
            ## (if /timestepping/nonlinear_iterations and
            ## /timestepping/nonlinear_iterations/tolerance are
            ## enabled) and file (if /io/convergence_file is enabled)
            ## i.e. excluding the components
            element exclude_components_from_convergence {
               comment
            }|
            ## Exclude this field entirely from convergence testing and file 
            element exclude_from_convergence {
               comment
            }
         )
      }
   )

# Steady state options for prognostic scalar fields
scalar_steady_state_options =
   (
      ## Decide whether this field is tested for a steady state
      ## between timesteps
      ## (if /timestepping/steady_state is
      ## enabled).
      element steady_state {
         (
            ## Include this field in steady state testing
            ## (if /timestepping/steady_state is
            ## enabled)
            element include_in_steady_state {
               comment
            }|
            ## Exclude this field from steady state testing
            element exclude_from_steady_state {
               comment
            }
         )
      }
   )

# Steady state statistics options for prognostic vector fields (velocity)
vector_steady_state_options =
   (
      ## Decide whether this field is tested for a steady state
      ## between timesteps
      ## (if /timestepping/steady_state is
      ## enabled).
      element steady_state {
         (
            ## Include this field (magnitude and components)
            ## in steady state testing
            ## (if /timestepping/steady_state is enabled)
            element include_in_steady_state {
               comment
            }|
            ## Include just the magnitude of this field 
            ## in steady state testing
            ## (if /timestepping/steady_state is
            ## enabled)
            ## i.e. excluding the components
            element exclude_components_from_steady_state {
               comment
            }|
            ## Exclude this field entirely from convergence testing and file 
            element exclude_from_steady_state {
               comment
            }
         )
      }
   )

# Options for whether a field is to be included in detector output.
detector_options_enabled_default = 
   (
      ## Specify what is added to detector files
      element detectors {
         (
            ## This field is output at each detector location.
            element include_in_detectors {
               comment
            }|
            ## This field is not output at detector locations.
            element exclude_from_detectors {
               comment
            }
         )
      }
   )

# Options for whether a field is to be included in detector output.
detector_options_disabled_default = 
   (
      ## Specify what is added to detector files
      element detectors {
         (
            ## This field is not output at detector locations.
            element exclude_from_detectors {
               comment
            }|
            ## This field is output at each detector location.
            element include_in_detectors {
               comment
            }
         )
      }
   )

# Detector output defaults on for prognostic and diagnostic fields, 
# off for prescribed.
prognostic_detector_options = detector_options_enabled_default
diagnostic_detector_options = detector_options_enabled_default
prescribed_detector_options = detector_options_disabled_default

# Specify whether or not values for this field are interpolated to particle locations..
particle_options_disabled_default = 
   (
      ## Specify what is added to particles
      element particles {
         (
            ## Field value is not interpolated to particle locations.
            element exclude_from_particles {
               comment
            }|
            ## Field value is interpolated to each particle location.
            element include_in_particles {
               ## Store values of field at previous timestep/non-linear iteration, for use when updating particle attributes.
               element store_old_field {
                  comment
               }?
            }
         )
      }
   )

# Particle defaults off for prognostic, diagnostic and prescribed fields
prognostic_particle_options = particle_options_disabled_default
diagnostic_particle_options = particle_options_disabled_default
prescribed_particle_options = particle_options_disabled_default

generic_aliased_field =
   (
      attribute material_phase_name { xsd:string },
      attribute field_name { xsd:string }
   )

generic_particle_attribute_choice =
   (
      ## Constant value
      element constant {
         real
      }|
      ## Python function prescribing real input. Functions should be of the form:
      ##
      ##  def val(X, t):
      ##     # Function code
      ##     return # Return value
      ##
      ## where X is a tuple of length geometry dimension.
      element python {
         python_code
      }|
      ## Python function prescribing real input. Functions should be of the form:
      ##
      ##  def val (X, t, fields):
      ##     # Function code
      ##     return # Return value
      ##
      ##  where X is a tuple of length geometry dimension.
      ##
      ##  E.g. if Temperature field and x position are required
      ##  def val (X, t, fields):
      ##    x = X[0]
      ##    y = x + fields["Temperature"]
      ##    return y
      element python_fields {
         python_code,
         ## Store attribute values from previous timestep on particles.
         element store_old_attribute {
            comment
         }?
      }|
      ## File containing the particle positions in binary form
      element from_checkpoint_file {
         attribute file_name { xsd:string },
         ## The format of the input file containing field data.
         element format {
            element string_value {
               "binary"
            }
         }
      }
   )

# This is the choice of additional scalar field to be solved for
scalar_field_choice =
   (
# The first is a generic field, which may be used for any user-defined field
# that FLUIDITY knows nothing about, or a generic diagnostic
      (
         element scalar_field {
            attribute rank { "0" },
            attribute name { xsd:string },
            ## Field type
            (
               element prognostic {
                  velocity_mesh_choice,   
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element diagnostic {
                  scalar_diagnostic_algorithms,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Prognostic scalar fields below this
         element ___Prognostic_Fields_Below___ {
            empty
         }|

# This is the long list of fields that FLUIDITY knows about
# -- First is a list of fields that are primarily prognostic,
#    but can be set to prescribed, or aliased...
# -- The list is in order of most frequently used.

         ## Salinity
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Salinity" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Temperature
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Temperature" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Background Temperature
         element scalar_field {
            attribute rank { "0" },
            attribute name { "BackgroundTemperature" },
            (
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|
         ## Passive Tracer
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Tracer" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Free Surface
         ## Use in combination with a free_surface boundary condition on velocity
         ## Without the no_normal_stress option, this field should be diagnostic
         ## and is optional. With the no_normal_stress option this field is required
         ## and should be prognostic.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "FreeSurface" },
            (
               ## Free Surface elevation, to be used in combination with:
               ##
               ## * a free_surface boundary condition applied to the Velocity
               ## field, without the no_normal_stress_option.
               ## It gives you a full dimensional field (constant over the vertical)
               ## of the free surface elevation.
               ##
               ## * equation type ShallowWater for Velocity. It simply divides 
               ## Pressure by g, the gravitational constant.
               element diagnostic {
                  internal_algorithm,
                  ## Must be on the same mesh as Pressure
                  pressure_mesh_choice,
                  diagnostic_scalar_field
               }|
               ## Free Surface
               ## NOTE: the prognostic FreeSurface field only works in combination
               ## with the free_surface boundary condition applied to the Velocity
               ## field in combination with the no_normal_stress_option.
               ## It gives you a full dimensional field (constant over the vertical)
               ## of the free surface elevation. Initial and boundary conditions
               ## for free surface should be set on this field.
               element prognostic {
                  ## Must be on the same mesh as Pressure
                  pressure_mesh_choice,
                  (
                     ## Initial condition for WholeMesh
                     ##
                     ## Only specify one condition if not using mesh regions.
                     ## Otherwise select other initial_condition option, specify region_ids
                     ## and distinct names.  Then add extra intial conditions for other regions.
                     element initial_condition {
                        attribute name { "WholeMesh" },
                        input_choice_initial_condition_real
                     }|
                     ## Multiple initial_conditions are allowed if specifying
                     ## different values in different
                     ## regions of the mesh (defined by region_ids).  In this case
                     ## each initial_condition
                     ## requires a distinct name for the options dictionary.
                     element initial_condition {
                        attribute name { string },
                        region_ids?,
                        input_choice_initial_condition_real
                     }
                  )+,
                  ## Solver - only necessary in combination with an explicit no_normal_stress free_surface bc under Velocity.
                  element solver {
                     linear_solver_options_sym
                  },
                  prognostic_scalar_output_options,
                  prognostic_scalar_stat_options,
                  scalar_convergence_options,
                  prognostic_detector_options,
                  prognostic_particle_options,
                  scalar_steady_state_options,
                  adaptivity_options_prognostic_scalar_field,
                  interpolation_algorithm_scalar
               }
               
            )
         }|
         ## Wetting and drying alpha coefficient. Alpha is 1 in dry and 0 in wet regions.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "WettingDryingAlpha" },
            (
               ## Free Surface
               ## NOTE: the diagnostic WettingDryingAlpha only works in combination
               ## with the free_surface boundary condition applied to the Velocity
               ## field. It gives you a 3D field (constant over the vertical)
               ## of the wetting and drying alpha coefficient.
               element diagnostic {
                  internal_algorithm,
                  # this is hard-coded on the PressureMesh as long as the Pressure is
                  # if this is no longer true, it should be option-checked to be on the
                  # same mesh as Pressure
                  ## Must be on the same mesh as Pressure
                  pressure_mesh_choice,
                  diagnostic_scalar_field
               }
               
            )
         }|
         ## Second Fluid
         element scalar_field {
            attribute rank { "0" },
            attribute name { "SecondFluid" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Diffuse Interface
         element scalar_field {
            attribute rank { "0" },
            attribute name { "DiffuseInterface" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## If enabled, decomposes Pressure by solving for the balanced part of 
         ## Pressure using a "geopressure" solver:
         ## 
         ##   f = - grad p_gp + g
         ##
         ## By choosing an appropriate mesh (typically velocity mesh order + 1)
         ## for the balanced part of pressure, physical balance can be
         ## represented to a higher degree of accuracy.
         ##
         ## If buoyancy is included and HydrostaticPressure or
         ## HydrostaticPressureGradient are enabled, the "geopressure" solver
         ## is modified:
         ##
         ##   f = - grad p_gp + g - grad p_hp
         element scalar_field {
            attribute rank { "0" },
            attribute name { "GeostrophicPressure" },
            (
               element prognostic {
                  ## The GeostrophicPressure mesh. Must be continuous.
                  ## 
                  ## <b>WARNING: It is usual for this to be a higher degree
                  ## mesh than the velocity mesh</b>
                  element mesh {
                     attribute name { xsd:string },
                     comment
                  },
                  prognostic_geostrophic_pressure_field
               }
            )
         }|
         ## If enabled, decomposes Pressure so that HydrostaticPressure balances against the density.
         ## Solves:
         ##
         ## grav_direction dot grad hp = grav*density
         ##
         ## for hp.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HydrostaticPressure" },
            (
               element prognostic {
                  mesh_choice,
                  prognostic_hydrostatic_pressure_field
               }
            )
         }|
         ## **UNDER TESTING**
         ##
         ## If enabled, decomposes Pressure so that VerticalBalancePressure balances against the density.
         ## Solves:
         ##
         ## (grav_direction dot grad) (grav_direction dot grad) vbp = (grav_direction dot grad) grav*density
         ##
         ## for vbp.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "VerticalBalancePressure" },
            (
               element prognostic {
                  ## This needs to be a continuous field!
                  mesh_choice,
                  prognostic_vertical_balance_pressure_field
               }
            )
         }|
         ## FoamVelocityPotential field:
         ##
         ## Required in foam flow simulations where the foam velocity is not a
         ## prescribed field. It solves Laplace's equation for a scalar potential, 
         ## (the foam velocity is then obtained as the gradient of the potential).
         element scalar_field {
            attribute rank { "0" },
            attribute name { "FoamVelocityPotential" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_foam_velocity_potential_field
               }
            )
         }|
         ## MaterialVolumeFraction field:
         ##
         ## Volume fraction of this material.
         ## Required in multimaterial simulations.
         ##  - if prognostic solves for the volume fraction
         ##  - if prescribed uses a specified volume fraction
         ##  - if diagnostic solves for the final material volume fraction
         ## Only 1 diagnostic MaterialVolumeFraction field allowed per
         ## simulation or solves for all the volume fractions based on
         ## the SumMaterialVolumeFractions field.
         ## 
         ## A diagnostic MaterialVolumeFraction field is currently required for
         ## compressible multimaterial simulations (even if only 1 material).
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialVolumeFraction" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field,
                  cap_option?,
                  surface_tension_option?
               }|
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field,
                  cap_option?
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field,
                  cap_option?
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## MaterialDensity field:
         ##
         ## Field for the density of this material.
         ## Required in compressible multimaterial simulations.
         ## Can be:
         ##  - diagnostic if using a linear equation of state
         ##  - prognostic if a compressible simulation
         ## (note that if you set a multimaterial
         ## equation of state and this field is
         ## prognostic then its initial condition
         ## will be overwritten by the density that
         ## satisfies the initial pressure and
         ## the equation of state)
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialDensity" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## MaterialInternalEnergy field:
         ##
         ## Field for the internal energy of this material.
         ## Required in multimaterial compressible simulations
         ## with full stiffened_gas (perfect gas) eos.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialInternalEnergy" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## InternalEnergy field:
         ##
         ## Field for the internal energy of a single material.
         ## Required in compressible simulations
         ## with full stiffened_gas (perfect gas) eos.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "InternalEnergy" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## SumMaterialVolumeFractions field:
         ##
         ## Sums the prognostic MaterialVolumeFraction fields.
         ## - diagnostic: sums all the volume fractions in the other
         ##   material phases
         element scalar_field {
            attribute rank { "0" },
            attribute name { "SumMaterialVolumeFractions" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## CopiedField - This field copies the previous timesteps
         ## values from another (specified) field at every iteration
         ## and then solves the field using different (again, specified)
         ## scheme and solution options.
         ## For instance, this field can be used to create a diffused
         ## field to adapt to.
         ## Unless someone requests otherwise this is only currently possible
         ## for fields within the same material_phase.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "CopiedField" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  attribute copy_from_field { string },
                  prognostic_scalar_field
               }
            )
         }|
         ## Calculate the stream function of 2D incompressible flow. Note 
         ## that this *only* makes sense for proper 2D (not pseudo-2D) simulations.
         ## Requires a continuous mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "StreamFunction" },
            (
               element prognostic {
                  mesh_choice,
                  prognostic_stream_function_field
               }
            )
         }|
         ## ***UNDER TESTING***
         ## Calculate the stream function of 2D incompressible flow for multiply connected domains.
         ## Note that this *only* makes sense for proper 2D (not pseudo-2D) simulations.
         ## Requires a continuous mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MultiplyConnectedStreamFunction" },
            attribute depends { "Velocity" },
            (
               element prognostic {
                  mesh_choice,
                  prognostic_multipath_stream_function_field
               }
            )
         }|
         ## Phytoplankton
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Phytoplankton" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|
         ## Zooplankton
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Zooplankton" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|
         ## Nutrient
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Nutrient" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|
         ## Detritus
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Detritus" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|
         ## Chlorophyll
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Chlorophyll" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|
         ## Ammonium
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Ammonium" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|

         ## PhaseVolumeFraction
         ## Required multiphase problem types
         element scalar_field {
            attribute rank { "0" },
            attribute name { "PhaseVolumeFraction" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field,
                  cap_option?
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element diagnostic {
                  velocity_mesh_choice,
                  internal_algorithm,
                  diagnostic_scalar_field_no_adapt,
                  cap_option?
               }
            )
         }|

         # Insert new prognostic scalar fields here using the template:
         #        element scalar_field {
         #            attribute rank { "0" },
         #            attribute name { "NewFieldName" },
         #            (
         #               element prognostic {
         #                  velocity_mesh_choice,
         #                  prognostic_scalar_field
         #               }|
         #               element prescribed {
         #                  velocity_mesh_choice,
         #                  prescribed_scalar_field
         #               }|
         #               element aliased {
         #                  generic_aliased_field
         #               }
         #            )
         #        }
         
# -- Second is a list of fields that are primarily prescribed,
#    but can be aliased. An example is wind velocity.
# -- The list is in order of most frequently used.

         ## Prescribed scalar fields below this
         element ___Prescribed_fields_below___ {
            empty
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "DistanceToSideBoundaries" },
            (
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## This is field should only be used with FEMDEM and adapt_at_first_timestep.
         ## It is used for adapting the mesh at the first time step
         ## and then it is removed from the state.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "FirstAdaptDummy" },
            (
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## This field should only be used with foam simulations.
         ## It is the local Plateau border length per unite volume.
         ## lambda = 1.71/(bubble_radius)^2
         element scalar_field {
            attribute rank { "0" },
            attribute name { "DrainageLambda" },
            (
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
#
# Insert new prescribed scalar fields here using the template:
#        element scalar_field {
#            attribute rank { "0" },
#            attribute name { "NewFieldName" },
#            (
#               element prescribed {
#                  velocity_mesh_choice,
#                  prescribed_scalar_field
#               }|
#               element aliased {
#                  generic_aliased_field
#               }
#            )
#        }
#
# -- Last is a list of fields that are primarily diagnostic,
#    but can be aliased. An example is Tidal Range.
# -- The list is in order of most frequently used.
#
         ## Diagnostic scalar fields below this
         element ___Diagnostic_Fields_Below___ {
            empty
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "PerturbationDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## ControlVolumeDivergence:
         ##
         ## div field
         ##
         ## Divergence of the velocity field where
         ## the divergence operator is defined using
         ## the control volume C^T matrix.
         ## This assumes that the test space is discontinuous
         ## control volumes.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "ControlVolumeDivergence" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name { string },
                  velocity_mesh_choice,
                  diagnostic_cv_divergence_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## FiniteElementDivergence:
         ##
         ## div field
         ##
         ## Divergence of the velocity field where
         ## the divergence operator is defined using
         ## the finite element C^T matrix.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "FiniteElementDivergence" },
            (
               element diagnostic {
                  legacy_internal_algorithm,
                  attribute field_name { string },
                  velocity_mesh_choice,
                  element integrate_divergence_by_parts {
                     empty
                  }?,
                  diagnostic_fe_divergence_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Diffusive dissipation
         ##
         ## Only coded for 2D
         ##
         ##   -g*drho/dy
         ##
         ## where rho is the Density scalar field
         ##
         ## BE AWARE OF WHAT THE DENSITY SCALAR FIELD IS
         ##
         ##   Note the actual diffusive dissipation = -2*g*kappa*drho/dy so need
         ##   to multiply by kappa in post-processing (2 subject to definition) where
         ##   drho/dt + u.grad(rho) = kappa grad^2(rho)
         ##   This assumes kappa is isotropic and constant
         ##   For a linear equation of state kappa = diffusivity (thermal or haline)
         ##
         ## c.f. Winters 1995, Journal of Fluid Mechanics
         element scalar_field {
            attribute rank { "0" },
            attribute name { "DiffusiveDissipation" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Richardson Number:
         ##
         ##  Ri = \frac{N^2}{(\frac{\partial u}{\partial z})^2 + (\frac{\partial u}{\partial z})^2}
         ## with 
         ##  N^2 = -\frac{g}{\rho_0}\frac{\partial \rho}{\partial z}
         ##
         ## Limitations:
         ##  - Gravity must be constant.
         ##  - Assumes gravity is in -ve final coordinate direction.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "RichardsonNumber" },
            attribute depends { "Velocity,PerturbationDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_richardson_number_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## CFLNumber
         ##
         ## See http://amcg.ese.ic.ac.uk/index.php?title=Local:Diagnostics#CFL_Number
         ##
         ## Adapting to this field is not recommended
         element scalar_field {
            attribute rank { "0" },
            attribute name { "CFLNumber" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## ControlVolumeCFLNumber
         ##
         ## Courant Number as defined on a control volume mesh
         ##
         ## Adapting to this field is not recommended
         element scalar_field {
            attribute rank { "0" },
            attribute name { "ControlVolumeCFLNumber" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## DG_CourantNumber
         ##
         ## Courant Number as defined on a DG mesh
         ##
         ## Adapting to this field is not recommended
         element scalar_field {
            attribute rank { "0" },
            attribute name { "DG_CourantNumber" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Buoyancy Adjustment diffusivity
         ##
         ## The element-wise diffusion applied
         ## in the buoyancy adjustment by vertical mixing
         ## scheme available for tracer fields.
         ##
         ## This is an element-wise P0 field, so it is best
         ## calculated in a discontinuous function space.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "BuoyancyAdjustmentDiffusivity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## The equilibrium pressure perturbation 
         ## 
         ## Further details to be added
         ## Currently a test field
         element scalar_field {
            attribute rank { "0" },
            attribute name { "EquilibriumPressure" },
            (
               element diagnostic {
                  internal_algorithm,
                  pressure_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## The equilibrium pressure perturbation 
         ## 
         ## Further details to be added
         ## Currently a test field
         element scalar_field {
            attribute rank { "0" },
            attribute name { "PressureRelativeToExternal" },
            (
               element diagnostic {
                  internal_algorithm,
                  pressure_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Discontinuity detector
         ##
         ## takes value 1 where detector is triggered
         ##
         ## Adapting to this field is not recommended
         element scalar_field {
            attribute rank { "0" },
            attribute name { "DiscontinuityDetector" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## CVMaterialDensityCFLNumber
         ##
         ## Courant Number as defined on a control volume mesh and
         ## incorporating the MaterialDensity.
         ## Requires a MaterialDensity field!
         ##
         ## Adapting to this field is not recommended
         element scalar_field {
            attribute rank { "0" },
            attribute name { "CVMaterialDensityCFLNumber" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## This scalar field is meant to replace DENTRAF.
         ## Basically, if you use new options, DENTRAF is no longer needed
         ## No repointing is done from this field to DENTRAF.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "CopyofDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            )
         }|
         ## add a MaterialVolume scalar_field to calculate the spatially varying 
         ## volume of a material (requires a MaterialVolumeFraction)
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialVolume" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## add a MaterialMass scalar_field to calculate the spatially varying 
         ## mass of a material (requires a MaterialVolumeFraction and either a
         ## MaterialDensity field or a reference_density in a linear eos)
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialMass" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Calculates the material density based on the bulk Pressure
         ## (and MaterialInternalEnergy if appropriate) for the equation
         ## of state of this material.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialEOSDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Calculates the MaterialPressure based on the MaterialDensity
         ## (and MaterialInternalEnergy if appropriate) for the equation
         ## of state of this material.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialPressure" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Calculates the BulkMaterialPressure based on the MaterialDensity
         ## and MaterialVolumeFraction (and MaterialInternalEnergy if appropriate) 
         ## for the equation of state of all materials.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "BulkMaterialPressure" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Buoyancy adjustment (mixing by diffusion) coefficient
         element scalar_field {
            attribute rank { "0" },
            attribute name { "BuoyancyAdjustment" },
            (element diagnostic {
               internal_algorithm,
               mesh_choice,
               diagnostic_scalar_field
            }
            | element aliased { generic_aliased_field })
         }|
         ## Grid Reynolds number
         element scalar_field {
            attribute rank { "0" },
            attribute name { "GridReynoldsNumber" },
            (element diagnostic {
               internal_algorithm,
               velocity_mesh_choice,
               ## Include the density field in the Grid Reynolds number calculation.
               ## A Density field must exist in the material_phase for this.
               ## This assumes that the viscosity used is the dynamic viscosity.
               element include_density_field {
                  comment
               }?,
               diagnostic_scalar_field
            }
            | element aliased { generic_aliased_field })
         }|
         ## GridPecletNumber
         ##
         ## Peclet Number Pe = U*dx/2*diffusivity
         ##
         ## Also see the test case 'grid_peclet_number'
         ## if you wish to see the effect of changing the 
         ## diffusivity on a 1D, cg-discretised tracer-field
         ##
         ## Adapting to this field is not recommended
         element scalar_field {
            attribute rank { "0" },
            attribute name { "GridPecletNumber" },
            (
               element diagnostic {
                  internal_algorithm,
                  ## Mesh on which to calculate dx
                  mesh_choice,
                  ## This is the name of the scalar field
                  ## to calculate the Peclet number for
                  ## Note this field needs to have a diffusivity
                  element field_name { string },
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Horizontal velocity divergence:
         ##
         ## div_H velocity
         ##
         ## Uses the gravity field direction to determine the horizontal plane.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HorizontalVelocityDivergence" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## Velocity divergence:
         ##
         ## div velocity
         ##
         element scalar_field {
            attribute rank { "0" },
            attribute name { "VelocityDivergence" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         
         ## Kinetic energy density:
         ##
         ##  1/2 rho*|u|^2
         ##
         ## where rho is the Density scalar field
         ##
         ## BE AWARE OF WHAT IS IN THE DENSITY SCALAR FIELD
         ##
         ## In the Boussinesq approximation 
         ## rho -> rho_0 (reference density)
         ## to get this use the square of the 
         ## L2norm of the velocity field
         ## from the stat file and multiply by
         ## 1/2*rho_0 in post-processing
         element scalar_field {
            attribute rank { "0" },
            attribute name { "KineticEnergyDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         
         ## Gravitational potential energy density:
         ##
         ## -rho*(g dot (r - r_0))
         ##
         ## where rho is the Density scalar field
         ## r_0 is the potential energy zero point
         ## and g is the gravity vector
         ##
         ## BE AWARE OF WHAT IS IN THE DENSITY SCALAR FIELD
         ##
         ## Limitations:
         ##  - Requires a constant gravity direction.
         ##  - The Density and GravitationalPotentialEnergyDensity fields must be on the same mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "GravitationalPotentialEnergyDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field,
                  ## Coordinate of a point with a potential energy of zero.
                  element zero_point {
                     real_dim_vector
                  }
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Isopycnal coordinate
         ## Not parallelised
         ##
         ##  z_star(x,t) = 1/A int_V' H(rho(x',t)-rho(x,t)) dV'
         ##
         ## where rho is the Density scalar field, A is the width/area of the domain
         ##
         ## BE AWARE OF WHAT IS IN THE DENSITY SCALAR FIELD
         ##
         ## Limitations:
         ##  - You need to specify a (fine) mesh to redistribute the Density onto
         ##  - Requires a constant gravity direction.
         ##  - The Density and GravitationalPotentialEnergyDensity fields must be on the same mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "IsopycnalCoordinate" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  ## This is the mesh onto which we redistribute the PerturbationDensity
                  element fine_mesh {
                     attribute name { string }
                  },
                  diagnostic_scalar_field                 
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Background potential energy density:
         ##
         ##   PE_b = rho*g*z_star
         ##
         ## where rho is the Density scalar field and
         ## z_star is the isopycnal coordinate. 
         ##
         ##
         ## BE AWARE OF WHAT IS IN THE DENSITY SCALAR FIELD
         ##
         ##
         ## Limitations:
         ##  - Requires a constant gravity direction.
         ##  - Requires a Density scalar field.
         ##  - Requires the IsopycnalCoordinate diagnostic field.
         ##  - The Density and GravitationalPotentialEnergyDensity
         ##  fields must be on the same mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "BackgroundPotentialEnergyDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field                 
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         
         ## Ertel potential vorticity:
         ##
         ##  (f + curl u) dot grad rho'
         ##
         ## Limitations:
         ##  - Requires a geometry dimension of 3.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "PotentialVorticity" },
            attribute depends { "Velocity,PerturbationDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Relative potential vorticity:
         ##
         ##   curl u dot grad rho'
         element scalar_field {
            attribute rank { "0" },
            attribute name { "RelativePotentialVorticity" },
            attribute depends { "Velocity,PerturbationDensity" },
            (
              element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Calculate the horizontal stream function psi where:
         ##   \partial_x \psi = -v
         ##   \partial_y \psi = u
         ## where u and v are perpendicular to the gravity direction. Applies a
         ## strong Dirichlet boundary condition of 0 on all boundaries.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HorizontalStreamFunction" },
            attribute depends { "Velocity" },
            (
              element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  ## Solver
                  element solver {
                     linear_solver_options_sym
                  },
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Speed:
         ##
         ##  |u|
         ##
         ## Limitations:
         ##  - The Speed and Velocity fields must be on the same mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Speed" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|


         ## Both fields must be in this material_phase.
         ## Assumes both fields are on the same mesh as the AbsoluteDifference field.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "AbsoluteDifference" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name_a { string },
                  attribute field_name_b { string },
                  mesh_choice,
                  diagnostic_scalar_field,
                  ## Evaluate the absolute difference once the average difference has been removed?
                  element relative_to_average {
                    empty
                  }?,
                  ## Ignore boundary nodes (i.e. zero them when calculating the difference)
                  element ignore_boundaries {
                    empty
                  }?
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## Absolute Difference between two scalar fields.
         ##
         ## Both fields must be in this material_phase.
         ## Assumes both fields are on the same mesh as the AbsoluteDifference field.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "ScalarAbsoluteDifference" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name_a { string },
                  attribute field_name_b { string },
                  mesh_choice,
                  diagnostic_scalar_field,
                  ## Evaluate the absolute difference once the average difference has been removed?
                  element relative_to_average {
                    empty
                  }?,
                  ## Ignore boundary nodes (i.e. zero them when calculating the difference)
                  element ignore_boundaries {
                    empty
                  }?
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Galerkin projection of one field onto another mesh.
         ##
         ## The field must be in this material_phase.
         ## 
         ## NOTE: you need the solver options if the mesh
         ## of this field is continuous.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "GalerkinProjection" },
            (
               element diagnostic {
                  legacy_internal_algorithm,
                  element source_field_name { string },
                  mesh_choice,
                  ## Lump the mass matrix of the galerkin projection
                  ## less accurate but faster and might give smoother result.                  
                  element lump_mass {
                     empty
                  }?,
                  element solver {
                    linear_solver_options_sym
                  }?,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Primary production of Phytoplankton. This is calculated by
         ## the ocean biology module and will not be calculated unless
         ## ocean biology is being simulated.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "PrimaryProduction" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Grazing of Phytoplankton by Zooplankton. This is calculated by
         ## the ocean biology module and will not be calculated unless
         ## ocean biology is being simulated.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "PhytoplanktonGrazing" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )

         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "TidalRange" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field_tidal_range
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "FreeSurfaceHistory" },
            (
               element diagnostic {
                  free_surface_history_algorithm,
                  velocity_mesh_choice
               }
            )
         }|
         ## Output the universal numbering of the mesh on which this field is based.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "UniversalNumber" },
            (
               element diagnostic {
                  legacy_internal_algorithm,
                  mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Output the processors which own the nodes of the mesh on which this field is based.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "NodeOwner" },
            (
               element diagnostic {
                  legacy_internal_algorithm,
                  mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## An estimate of the edge wieghts whilst adapting using Zoltan
         ## Note: you *must* turn on inteprolation otherwise this field
         ## will be emptied by an adapt.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaxEdgeWeightOnNodes" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## SumVelocityDivergence
         ##
         ## Diagnostic field used in multiphase simulations.
         ## Sums up the divergence of each phase's apparent velocity, i.e. \sum{ div(vfrac*u) }
         element scalar_field {
            attribute rank { "0" },
            attribute name { "SumVelocityDivergence" },
            (
               element diagnostic {
                  velocity_mesh_choice,
                  internal_algorithm,
                  element integrate_divergence_by_parts {
                     empty
                  }?,
                  ## Test the divergence operator with the CV dual function space. 
                  ## This is useful if the pressure is solved with a CV discretisation 
                  ## or with a CG but with the continuity tested with the CV dual. 
                  element test_with_cv_dual {
                     comment
                  }?,
                  diagnostic_scalar_field_no_adapt,
                  element solver {
                     linear_solver_options_sym
                  }
               }
            )
         }|
         
         ## CompressibleContinuityResidual
         ##
         ## Computes the residual of the continuity equation used in compressible multiphase flow simulations
         ## i.e. vfrac_c*d(rho_c)/dt + div(rho_c*vfrac_c*u_c) + \sum_i{ rho_c*div(vfrac_i*u_i) }
         ## where _c and _i denote the compressible and incompressible phases respectively.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "CompressibleContinuityResidual" },
            (
               element diagnostic {
                  velocity_mesh_choice,
                  internal_algorithm,
                  diagnostic_scalar_field_no_adapt,
                  element solver {
                     linear_solver_options_sym
                  }
               }
            )
         }|

         ## DGLESScalarEddyViscosity
         ##
         ## Outputs the calculated (isotropic) eddy viscosity term
         ## for a discontinuous velocity field. The subgridscale viscosity
         ## follows that of the Smagorinsky–Lilly model. Requires that both a discontinous
         ## representation of velocity and an les model have been selected.

         element scalar_field {
            attribute rank { "0" },
            attribute name { "DGLESScalarEddyViscosity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         } 

# Insert new diagnostic scalar fields here using the template:
#        element scalar_field {
#            attribute rank { "0" },
#            attribute name { "NewFieldName" },
#            (
#               element diagnostic {
#                  internal_algorithm,
#                  velocity_mesh_choice,
#                  diagnostic_scalar_field
#               }|
#               element aliased {
#                  generic_aliased_field
#               }
#            )
#        }
      )
   )

# This is the choice of additional vector field to be solved for
vector_field_choice =
   (
# The first is a generic field, which may be used for any user-defined field
# that FLUIDITY knows nothing about, or a generic diagnostic
# Prognostic vector fields are not possible (other than velocity and those known fields below).
      (
         ## Generic field variable (vector)
         element vector_field {
            attribute rank { "1" },
            attribute name { xsd:string },
            ## Field type
            (
               element prescribed {
                  mesh_choice,
                  prescribed_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }|
               element diagnostic {
                  vector_diagnostic_algorithms,
                  velocity_mesh_choice,
                  diagnostic_vector_field
               }
            )
         }|
#
# -- List of fields that are primarily prognostic,
#    but can be aliased.
# -- The list is in order of most frequently used.
#
         ## Prescribed vector fields below this
         element ___Prognostic_fields_below___ {
            empty
         }|
         ## As HydrostaticPressure, but solves for the gradient of the pressure
         ## associated with buoyancy. Requires a discontinuous mesh.
         element vector_field {
            attribute rank { "0" },
            attribute name { "HydrostaticPressureGradient" },
            (
               element prognostic {
                  mesh_choice,
                  prognostic_hydrostatic_pressure_gradient_field
               }
            )
         }|
#
# -- List of fields that are primarily prescribed,
#    but can be aliased.
# -- The list is in order of most frequently used.
#
         ## Prescribed vector fields below this
         element ___Prescribed_fields_below___ {
            empty
         }|

         ## MaterialVelocity field.  Used to impose a velocity on a material
         ## in combination with the imposed_material_velocity_absorption_algorithm
         ## for VelocityAbsorption and imposed_material_velocity_source_algorithm 
         ## for VelocitySource.
         element vector_field {
            attribute rank { "1" },
            attribute name { "MaterialVelocity" },
            (
               element prescribed {
                  mesh_choice,
                  prescribed_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
#
# Insert new prescribed vector fields here using the template:
#        element vector_field {
#            attribute rank { "1" },
#            attribute name { "NewFieldName" },
#            (
#               element prescribed {
#                  mesh_choice,
#                  prescribed_vector_field
#               }|
#               element aliased {
#                  generic_aliased_field
#               }
#            )
#        }
#
# -- Last is a list of fields that are primarily diagnostic,
#    but can be aliased. An example is Tidal Range.
# -- The list is in order of most frequently used.
#
         ## Diagnostic vector fields below this
         element ___Diagnostic_Fields_Below___ {
            empty
         }|

         ## Gradient of a scalar field evaluated using the C gradient
         ## matrix constructed using finite elements.
         ## Field must be in this material_phase.
         element vector_field {
            attribute rank { "1" },
            attribute name { "FiniteElementGradient" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name { string },
                  mesh_choice,
                  element integrate_gradient_by_parts {
                     empty
                  }?,
                  diagnostic_gradient_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## Gradient of a scalar field evaluated using the transpose
         ## of the C^T divergence matrix constructed using finite
         ## elements.
         ## Field must be in this material_phase.
         element vector_field {
            attribute rank { "1" },
            attribute name { "FiniteElementDivergenceTransposed" },
            (
               element diagnostic {
                  legacy_internal_algorithm,
                  attribute field_name { string },
                  mesh_choice,
                  element integrate_divergence_by_parts {
                     empty
                  }?,
                  diagnostic_gradient_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## Foam Velocity field. 
         ## Required for simulations of foam flow and liquid drainage in foams. 
         ## If diagnostic, a FoamVelocityPotential field is required, 
         ## and the Foam Velocity is obtained by taking its gradient.
         element vector_field {
            attribute rank { "1" },
            attribute name { "FoamVelocity" },
            (
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_vector_field_no_adapt
               }|
               element diagnostic {
                  internal_algorithm,
                  mesh_choice,
                  element solver {
                     linear_solver_options_sym
                  },
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## Foam Liquid Content times Velocity field.
         ## Its surface integral gives the liquid 
         ## volumetric flow in the flowing foam 
         element vector_field {
            attribute rank { "1" },
            attribute name { "FoamLiquidContentVelocity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_vector_field
               }     
            )
         }|
         
         ## Relative vorticity field - curl of the velocity field
         element vector_field {
            attribute rank { "1" },
            attribute name { "Vorticity" },
            (
               element diagnostic {
                  vorticity_algorithm,
                  velocity_mesh_choice,
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|  

         ## Planetary vorticity
         ##
         ## Limitations:
         ##  - Requires geometry dimension of 3.
         element vector_field {
            attribute rank { "1" },
            attribute name { "PlanetaryVorticity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Absolute vorticity:
         ##
         ##   f + curl u
         ##
         ## Limitations:
         ##  - Requires a geometry dimension of 3.
         element vector_field {
            attribute rank { "1" },
            attribute name { "AbsoluteVorticity" },
            attribute depends { "Velocity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         
         ## Gradient of a scalar field evaluated using the transpose
         ## of the C^T matrix constructed using control volumes.
         ## Field must be in this material_phase.
         element vector_field {
            attribute rank { "1" },
            attribute name { "ControlVolumeDivergenceTransposed" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name { string },
                  velocity_mesh_choice,
                  diagnostic_cv_gradient_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## LinearMomentum field.
         ##  p = \rho*u 
         ## (where p is the linear momentum, \rho the density and u the velocity)
         element vector_field {
             attribute rank { "1" },
             attribute name { "LinearMomentum" },
             (
                element diagnostic {
                   internal_algorithm,
                   velocity_mesh_choice,
                   diagnostic_vector_field
                }|
                element aliased {
                   generic_aliased_field
                }
             )
         }|

         ## Absolute Difference between two vector fields.
         ##
         ## Both fields must be in this material_phase.
         ## Assumes both fields are on the same mesh as the AbsoluteDifference field.
         element vector_field {
            attribute rank { "1" },
            attribute name { "AbsoluteDifference" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name_a { string },
                  attribute field_name_b { string },
                  mesh_choice,
                  diagnostic_vector_field,
                  ## Evaluate the absolute difference once the average difference has been removed?
                  element relative_to_average {
                    empty
                  }?
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## Absolute Difference between two vector fields.
         ##
         ## Both fields must be in this material_phase.
         ## Assumes both fields are on the same mesh as the AbsoluteDifference field.
         element vector_field {
            attribute rank { "1" },
            attribute name { "VectorAbsoluteDifference" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name_a { string },
                  attribute field_name_b { string },
                  mesh_choice,
                  diagnostic_vector_field,
                  ## Evaluate the absolute difference once the average difference has been removed?
                  element relative_to_average {
                    empty
                  }?
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         
         ## Bed Shear Stress
         ## 
         ## This diagnostic vector field is only calculated over surface elements/nodes, 
         ## interior nodes will have zero value.
         element vector_field {
             attribute rank { "1" },
             attribute name { "BedShearStress" },
             (
                element diagnostic {
                   internal_algorithm,
                   velocity_mesh_choice,
                   diagnostic_vector_field_bed_shear_stress
                }|
                element prescribed {
                   velocity_mesh_choice,
                   prescribed_vector_field_no_adapt
                }|
                element aliased {
                   generic_aliased_field
                }
             )
         }|
         ## Max Bed Shear Stress.
         ##
         ## Note that you need BedShearStress turned on for this to work.
         element vector_field {
             attribute rank { "1" },
             attribute name { "MaxBedShearStress" },
             (
                element diagnostic {
                   internal_algorithm,
                   velocity_mesh_choice
                   #diagnostic_vector_field_max_bed_shear_stress
                }|
                element aliased {
                   generic_aliased_field
                }
             ),
            ## This is the time after which the max operator is
            ## applied to the bed shear stress.
            element spin_up_time {
               real
            }
         }|

         ## Coordinate field remapped to the mesh of your choice.
         element vector_field {
             attribute rank { "1" },
             attribute name { "DiagnosticCoordinate" },
             (
                element diagnostic {
                   internal_algorithm,
                   mesh_choice,
                   diagnostic_vector_field
                }|
                element aliased {
                   generic_aliased_field
                }
             )
         }|
         ## Galerkin projection of one field onto another mesh.
         ##
         ## The field must be in this material_phase.
         ## 
         ## NOTE: you need the solver options if the mesh
         ## of this field is continuous.
         element vector_field {
            attribute rank { "1" },
            attribute name { "GalerkinProjection" },
            (
               element diagnostic {
                  legacy_internal_algorithm,
                  element source_field_name { string },
                  mesh_choice,
                  ## Lump the mass matrix of the galerkin projection
                  ## less accurate but faster and might give smoother result.                  
                  element lump_mass {
                     empty
                  }?,                  
                  element solver {
                  linear_solver_options_sym
                  }?,
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Computes the buoyancy term
         element vector_field {
            attribute rank { "1" },
            attribute name { "Buoyancy" },
            (
               element diagnostic {
                  buoyancy_algorithm,
                  velocity_mesh_choice,
                  diagnostic_vector_field
               }
            )
         }|
         ## Projects the Coriolis term onto the mesh of this diagnostic field.
         ## Note that multiple projection methods are available (under the
         ## algorithm option).
         element vector_field {
            attribute rank { "1" },
            attribute name { "Coriolis" },
            (
               element diagnostic {
                  coriolis_algorithm,
                  velocity_mesh_choice,
                  diagnostic_vector_field
               }
            )
         }

# Insert new diagnostic vector field here using the template:
#        element vector_field {
#            attribute rank { "1" },
#            attribute name { "NewFieldName" },
#            (
#               element diagnostic {
#                  internal_algorithm,
#                  mesh_choice,
#                  diagnostic_vector_field
#               }|
#               element aliased {
#                  generic_aliased_field
#               }
#            )
#        }
      )
   )  

# This is the choice of additional tensor fields
tensor_field_choice =
   (
# The first is a generic field, which may be used for any user-defined field
# that FLUIDITY knows nothing about, or a generic diagnostic
# Prognostic tensor fields are not possible.
      (
         ## Generic field variable (tensor)
         element tensor_field {
            attribute rank { "2" },
            attribute name { xsd:string },
            ## Field type
            (
               element prescribed {
                  mesh_choice,
                  prescribed_tensor_field
               }|
               element aliased {
                  generic_aliased_field
               }|
               element diagnostic {
                  tensor_diagnostic_algorithms,
                  velocity_mesh_choice,
                  diagnostic_tensor_field
               }
            )
         }|
#
# -- Second is a list of tensor fields that are primarily prescribed,
#    but can be aliased.
# -- The list is in order of most frequently used.
#
         ## Prescribed scalar fields below this
         element ___Prescribed_fields_below___ {
            empty
         }|

         ## MaterialViscosity field:
         ##
         ## Field for the viscosity of this material.
         ## Required if using a diagnostic bulk viscosity
         ## in a multimaterial simulation.
         element tensor_field {
            attribute rank { "2" },
            attribute name { "MaterialViscosity" },
            (
               element prescribed {
                  mesh_choice,
                  prescribed_tensor_field
               }|
               element diagnostic {
                  mesh_choice,
                  tensor_python_diagnostic_algorithm,
                  diagnostic_tensor_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

#
# Insert new prescribed tensor fields here using the template:
#        element tensor_field {
#            attribute rank { "2" },
#            attribute name { "NewFieldName" },
#            (
#               element prescribed {
#                  mesh_choice,
#                  prescribed_tensor_field
#               }|
#               element aliased {
#                  generic_aliased_field
#               }
#            )
#        }|
#
# -- Last is a list of fields that are primarily diagnostic,
#    but can be aliased.
# -- The list is in order of most frequently used.
#
         ## Diagnostic tensor fields below this
         element ___Diagnostic_Fields_Below___ {
            empty
         }

# Insert new diagnostic tensor field here using the template:
#        element tensor_field {
#            attribute name { "NewFieldName" },
#            (
#               element diagnostic {
#                  internal_algorithm,
#                  mesh_choice,
#                  diagnostic_tensor_field
#               }|
#               element aliased {
#                  generic_aliased_field
#               }
#            )
#        }
      )
   )

cap_option =
   (
      ## Cap the min and max values of this field when using
      ## it as a volume fraction to work out bulk material
      ## properties.
      ## No capping used if not selected.
      element cap_values {
         ## Set the upper bound on the field.
         ## Defaults to huge(0.0)*epsilon(0.0) if not set.
         element upper_cap {
            real
         }?,
         ## Set the lower bound on the field.
         ## Defaults to -huge(0.0)*epsilon(0.0) if not set.
         element lower_cap {
            real
         }?
      }
   )

surface_tension_option =
   (
      element surface_tension {
        ## Surface tension coefficient
        element surface_tension_coefficient {
          real
        },
        ## The equilibrium contact angle (in radians) with the boundaries identified by the surface ids
        element equilibrium_contact_angle {
          real,
          ## Surface ids:
          element surface_ids {
              integer_vector
          }
        }?
      }
   )

limiter_options =
  (
      (
        ## Limit the face value to satisfy a boundedness criterion.
        element limit_face_value{
          (
            sweby_limiter|
            ultimate_limiter
          )
        }|
        ## Do not limit the face value
        element do_not_limit_face_value{
          empty
        }
      )
  )

# cv limiter options that 
# have NO reliance on any CFL field. 
limiter_options_excluding_cfl =
  (
      (
        ## Limit the face value to satisfy a boundedness criterion.
        element limit_face_value{
          (
            sweby_limiter
          )
        }|
        ## Do not limit the face value
        element do_not_limit_face_value{
          empty
        }
      )
  )

sweby_limiter = 
  ## See "High-Resolution Schemes Using Flux Limiters for
  ## Hyperbolic Conservation-Laws", P. K. Sweby, 1984, Siam
  ## Journal on Numerical Analysis, 21, 995-1011
  element limiter {
    attribute name {"Sweby"},
    slope_options?,
    upwind_value_options?
  }

ultimate_limiter =
  ## See "The Ultimate Conservative Difference Scheme Applied
  ## to Unsteady One-Dimensional Advection", B. P. Leonard,
  ## 1991, Computer Methods in Applied Mechanics and
  ## Engineering, 88, 17-74
  element limiter {
    attribute name {"Ultimate"},
    field_based_cfl_number_options,
    upwind_value_options?
  }

slope_options =
   (
      ## Control the upper and lower slopes of the NVD limiter
      element slopes {
         ## Defaults to Sweby, 1984 limiter (= 1.0) if unselected
         element lower {
            real
         }?,
         ## Defaults to Sweby, 1984 limiter (= 2.0) if unselected
         element upper {
            real
         }?
      }
   )

upwind_value_options =
   (
      (
         ## Select the method to be used for calculating the upwind value.
         ## If not selected will default to project_upwind_value_from_point for
         ## simplex element meshes and to a locally_bound_upwind_value for cube
         ## element meshes.
         ##
         ## This method projects the upwind value from a point in the element just
         ## upwind of the node pair straddling the face.  It is otherwise known as 
         ## anisotropic limiting.
         ## This is only available on simplex meshes as it involes a search around
         ## the donor node to find the upwind element.
         element project_upwind_value_from_point {
            ## When the donor node is on a domain boundary reflect the projection
            ## back into the mesh.
            element reflect_off_domain_boundaries {
               empty
            }?,
            ## Constrain the projected value to be between the min and max of the
            ## element values which it was found from.
            element bound_projected_value_locally {
               empty
            }?,
            ## Store the locations of the elements where the upwind values
            ## are projected from for each node pair.
            ## This inserts an integer csr matrix into state so is memory expensive but
            ## saves a significant amount of time (searching around the neighbouring elements).
            ## This is unsafe for moving meshes but should be ok for adaptive meshes.
            element store_upwind_elements {
               ## Store the quadrature locations within the elements
               ## where the upwind values
               ## are projected from for each node pair.
               ## This inserts a real block csr matrix into state so is even more memory
               ## expensive than just storing the upwind elements and
               ## only saves a comparitively
               ## marginal amount of time (as actually searching the
               ## neighbouring elements is the
               ## slowest bit, finding the quadrature is relatively easy).
               element store_upwind_quadrature {
                  empty
               }?
            }?
         }|
         ## Select the method to be used for calculating the upwind value.
         ## If not selected will default to project_upwind_value_from_point for
         ## simplex element meshes and to a locally_bound_upwind_value for cube
         ## element meshes.
         ##
         ## Projects the value of the advected variable from the downwind or donor node
         ## using the interpolated gradient at the donor node in the
         ## direction of the vector
         ## connecting the node pair straddling the face.
         ## This is available on all meshes (except if bounding the values).
         element project_upwind_value_from_gradient {
            (
               ## Select which node to project from:
               ## Project from the downwind node (Jasak et al., 1999) so that:
               ## upwind_value = downwind_value - 2*gradient.vector
               element project_from_downwind_value {
                  comment
               }|
               ## Select which node to project from:
               ## Project from the donor node so that:
               ## upwind_value = donor_value - gradient.vector
               element project_from_donor_value {
                  comment
               }
            ),
            ## When the donor node is on a domain boundary reflect the projection
            ## back into the mesh.
            element reflect_off_domain_boundaries {
               empty
            }?,
            ## Constrain the projected value to be between the min and max of the
            ## element values which surround it.
            ## This is only available on simplex meshes as it involes a search around
            ## the donor node to find the upwind element.
            element bound_projected_value_locally {
               ## Store the locations of the elements closest to the project value.
               ## This inserts an integer csr matrix into state so is
               ## memory expensive but
               ## saves a significant amount of time (searching around
               ## the neighbouring elements).
               ## This is unsafe for moving meshes but should be ok for adaptive meshes.
               element store_upwind_elements {
                  comment
               }?
            }?
         }|
         ## Select the method to be used for calculating the upwind value.
         ## If not selected will default to project_upwind_value_from_point for
         ## simplex element meshes and to a locally_bound_upwind_value for cube
         ## element meshes.
         ##
         ## Chooses an upwind value by selecting the maximum or minimum of the neighbouring
         ## nodes depending on the local slope of the donor and downwind values.
         ## Otherwise known as isotropic limiting.
         ## This is available on all meshes except periodic domains.
         element locally_bound_upwind_value {
            empty
         }|
         ## Select the method to be used for calculating the upwind value.
         ## If not selected will default to project_upwind_value_from_point for
         ## simplex element meshes and to a locally_bound_upwind_value for cube
         ## element meshes.
         ##
         ## Chooses an upwind value by selecting the value at the node most directy
         ## upwind from the vector connecting the donor and downwind nodes.
         ## This is available on all meshes.
         element pseudo_structured_upwind_value {
            empty
         }
      )
   )

field_based_cfl_number_options =
   (
      (
         ## Select the Courant Number definition to be used.
         ##
         ## This uses the control volume definition of the CFL Number.
         element courant_number {
            attribute name { "ControlVolumeCFLNumber" },
            empty
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the Courant Number definition to be used.
         ##
         ## This uses the control volume definition of the CFL Number.
         element courant_number {
            attribute name { "CVMaterialDensityCFLNumber" },
            empty
         }|
         ## Select the Courant Number definition to be used.
         ##
         ## This uses the finite element approximation of the CFL Number.
         element courant_number {
            attribute name { "CFLNumber" },
            empty
         }|
         ## Select the Courant Number definition to be used.
         element courant_number {
            attribute name { string },
            empty
         }
      )
   )

dg_field_based_cfl_number_options =
   (
      (
         ## IF NOT SELECTED THEN THE DEFAULT DG_CourantNumber IS USED.
         ## 
         ## Select the Courant Number definition to be used.
         ## 
         ## This uses the DG finite element approximation of the CFL Number.
         element courant_number {
            attribute name { "DG_CourantNumber" },
            empty
         }|
         ## IF NOT SELECTED THEN THE DEFAULT DG_CourantNumber IS USED.
         ## 
         ## Select the Courant Number definition to be used.
         element courant_number {
            attribute name { string },
            empty
         }
      )?
   )

cv_face_cfl_number_options =
   (
      (
         ## Select the Courant Number definition to be used in the slope of
         ## the NVD diagram upper bound.
         ## This uses the control volume definition of the CFL Number.
         element courant_number {
            attribute name { "ControlVolumeCFLNumber" },
            empty
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the Courant Number definition to be used in the slope of
         ## the NVD diagram upper bound.
         ## This uses a control volume definition of the CFL Number
         ## that incorporates the MaterialDensity.
         ## Requires a MaterialDensity field in this material_phase!
         element courant_number {
            attribute name { "CVMaterialDensityCFLNumber" },
            empty
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the Courant Number definition to be used in the slope of
         ## the NVD diagram upper bound.
         ## This uses the finite element approximation of the CFL Number.
         element courant_number {
            attribute name { "CFLNumber" },
            empty
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the Courant Number definition to be used in the slope of
         ## the NVD diagram upper bound.
         element courant_number {
            attribute name { string },
            empty
         }
      )
   )

timestep_cfl_number_options =
   (
      (
         ## Select the Courant Number definition to be used for adaptive timestepping.
         ## This uses the finite element approximation of the CFL Number.
         element courant_number {
            attribute name { "CFLNumber" },
              ## Select the mesh on which you wish to evaluate the CFLNumber.
              velocity_mesh_choice
         }|
         ## Select the Courant Number definition to be used for adaptive timestepping.
         ## This uses the control volume definition of the CFL Number.
         element courant_number {
            attribute name { "ControlVolumeCFLNumber" },
              ## Select the mesh on which you wish to evaluate the ControlVolumeCFLNumber.
              velocity_mesh_choice
         }
      )
   )

mixing_stats =
   (
      ## Enable to include in the .stat file the fractions of the
      ## scalar field contained in
      ## bins specified by the user. This allows mixing of the field to be quantified.
      ## Replaces and expands upon the old heaviside.dat file
      element include_mixing_stats{
         attribute name { xsd:string },
         (
            ## Select whether to evaluate the volume fraction over the finite element
            ## (continuous galerkin) or within the control volume (control_volumes).
            ##
            ## NOTE: continuous_galerkin only works with linear tets
            ##
            ## NOTE: continuous_galerkin is not fully validated yet
            element continuous_galerkin {
               ## if select normalise the volume fractions will be
               ## divided by the total volume of the domain
               element normalise {
                  empty
               }?
            }|
            ## Select whether to evaluate the volume fraction over the finite element
            ## (continuous galerkin) or within the control volume (control_volumes).
            element control_volumes {
               ## if select normalise the volume fractions will be divided by the total volume of the domain  
               element normalise {
                  empty
               }?
            }
         ),
         ## The values of the bounds of the bins 
         ## e.g. the values -1.5 0.0 1.5 2.0 will return 4 bins 
         ## and the fraction of the field in each bin with,
         ## -1.5<=field<0.0, 0.0<=field<1.5, 1.5<=field<2.0, 2.0<=field, 
         ## will be calculated.  
         element mixing_bin_bounds { 
            (
                 ## list of bin bounds
                 element constant { 
                     real_vector
                 }|
                 ## Python function prescribing bin bounds. Functions should be of the form:
                 ##
                 ##  def val(t):
                 ##     # Function code
                 ##     return # Return value that should be an array of reals
                 ##
                 ## 
                 element python {
                     python_code
                   }
              )
         },
         ## Define the tolerance beneath the specified bins that should be included.
         ## Defaults to zero at machine tolerance (epsilon(0.0)) if not selected.
         element tolerance {
            real
         }?
      }
   )

cv_stats =
   (
      ## Include statistics evaluated on the control volume mesh.
      element include_cv_stats {
         empty
      }
   )

# Options for inclusion of calculations of surface integrals in the .stat file   
surface_integral_stats_base.surface_integral =
   (
      attribute name { xsd:string },
      ## Surface IDs defining the surface over which to integrate. If disabled, integrates over the whole surface.
      element surface_ids {
         integer_vector
      }?,
      ## Enable to normalise the integral by dividing by the surface area
      element normalise {
         comment
      }?
   )
surface_integral_stats_scalar =
   (
      ## Surface integral calculations. The following integral types are available:
      ##  value: Integrates the field
      ##  gradient_normal: Integrates the normal component of the gradient of the field
      element surface_integral {
         surface_integral_stats_scalar.surface_integral
      }
   )
surface_integral_stats_scalar.surface_integral = surface_integral_stats_base.surface_integral
surface_integral_stats_scalar.surface_integral &=
   (
      attribute type { "value" | "gradient_normal" }
   )
surface_integral_stats_vector =
   (
      ## Surface integral calculations. The following integral types are available:
      ##  normal: Integrates the normal component of the field
      element surface_integral {
         surface_integral_stats_vector.surface_integral
      }
   )
surface_integral_stats_vector.surface_integral = surface_integral_stats_base.surface_integral
surface_integral_stats_vector.surface_integral &=
   (
      attribute type { "normal" }
   )


forcing =
  (
    ## Force the boundary conditions of various fields in ocean simulations
    element ocean_forcing{
        element iceshelf_meltrate{
           ## The type of parameterisation 
           element Holland08{
           ## Specific heat capacity of water [Jkg^{-1}]. Default value 3974.
           element c0 {real},  
           ## Specific heat capacity of ice [Jkg^{-1}] Default value 2009.
           element cI {real},
           ## Transfer of heat and salt through the boundary layer [Jkg^{-1}]. Default value 3.35e5
           element L {real},
           ## Temperature of ice [C]. Default value -25, pretty cold 
           element TI {real},
           ## Tb=aSb+b+cB, See equation (4) of Holland08. Default value -0.0573C. 
           element a {real},
           ## Tb=aSb+b+cB, See equation (4) of Holland08. Default value 0.0832C. 
           element b {real},
           ## Drag coefficient, 1.5e-3, see Holland and Jenkins's(1999) Cd. 
           element Cd {real},
           ## Surface ID for the ice
           element melt_surfaceID {integer_vector}?,
           ## Distance for the far field
           element melt_LayerLength {real},
           ## Meltrate, calculated from the three equations.
           element scalar_field {
                     attribute rank { "0" },
                     attribute name { "MeltRate" },
                     (
                       element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field
                        }
                     )
           }?,
           ## Temperature of the ice-ocean interface, calculated from the three equations.
           element scalar_field {
                     attribute rank { "0" },
                     attribute name { "Tb" },
                     (
                       element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field
                        }
                     )
           }?,
           ## Salinity of the ice-ocean interface, calculated from the three equations.
           element scalar_field {
                     attribute rank { "0" },
                     attribute name { "Sb" },
                     (
                       element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field
                        }
                     )
           }?,
           ## Heat flux at the ice-ocean interface.
           element scalar_field {
                     attribute rank { "0" },
                     attribute name { "Heat_flux" },
                     (
                       element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field
                        }
                     )
           }?,
           ## Salt flux at the ice-ocean interface.
           element scalar_field {
                     attribute rank { "0" },
                     attribute name { "Salt_flux" },
                     (
                       element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field
                        }
                     )
           }?,
           ## Temperature of the far field, which is "melt_LayerLength" away from the ice. 
           ## This variable feeds into the three equations. 
           element scalar_field {
                     attribute rank { "0" },
                     attribute name { "Tloc" },
                     (
                       element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field
                        }
                     )
           }?,
           ## Salinity of the far field, which is "melt_LayerLength" away from the ice.
           ## This variable feeds into the three equations. 
           element scalar_field {
                     attribute rank { "0" },
                     attribute name { "Sloc" },
                     (
                       element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field
                        }
                     )
           }?,
           ## Pressure used in the three equations. 
           ## We need to use the pressure at the ice-ocean interface.
           element scalar_field {
                     attribute rank { "0" },
                     attribute name { "Ploc" },
                     (
                       element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field
                        }
                     )
           }?,
           ## Velocity of the far field, which is "melt_LayerLength" away from the ice.
           ## This variable feeds into the three equations. 
           element vector_field {
                     attribute rank { "0" },
                     attribute name { "Vloc" },
                     (
                       element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_vector_field
                        }
                     )
           }?,
           ## Location of your far field, which is "melt_LayerLength" away from the ice.
           ## Location should be perpendicular to the ice-ocean boundary.
           ## melt_LayerLength = ||Location - Location_org||
           element vector_field {
                     attribute rank { "0" },
                     attribute name { "Location" },
                     (
                       element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_vector_field
                        }
                     )
           }?,    
           ## Location of your ice-ocean boundary.
           element vector_field {
                     attribute rank { "0" },
                     attribute name { "Location_org" },
                     (
                       element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_vector_field
                        }
                     )
           }?, 
            ## Boundary stuff   
           element calculate_boundaries {
                     element string_value {
                     "neumann"|"dirichlet"
                     }
                  }?                       
           }?      
        }?,
        element shelf {
           ## Multiplying factor on equilibrium pressure
           element amplitude {
             real
           }?,
           ## Change sign of depth / perturbation from original eta_0
           element y_sign {
             empty
           }?,
           ## Only calculate the field, don't add to pressure (RHS of Momentum)
           element calculate_only {
             empty
           }?,
           ## DeltaS for constant change over ocean,
           ## or the difference between min and max for a linear variation in S
           element salinity_change_constant {
             real
           }?,
           ## If the ice shelf is not included in the original mesh and a free-surface initial condition is used to generate the shelf, the hydrostatic pressure where the shelf now exists needs to be included
           element add_pressure_from_ice {
             empty
           }?,
           ## If using a linear variation in S, include the associated change in density of the ice shelf
           ## NOT arranged to work with add_pressure_from_ice - TODO
           element include_density_change_of_ice {
             empty
           }?
        }?,
        ## Check this option to use an external data set for the t(0)
        ## and in/out boundary conditions. This sets up NEMO-type forcing.
        element external_data_boundary_conditions {
            ## Path of the file containing the external data set
            element input_file {
                attribute file_name {xsd:string} 
            }
        }?,
        ## Bulk formulae allow the boundary conditions of Temperature, Velocity, Salinty
        ## and photosynthetic radiation to be set using data obtained from
        ## the ERA-40 reanalysis (http://data-portal.ecmwf.int/data/d/era40_daily/)
        element bulk_formulae {
            ## The type of bulk formulae to use. Default is NCAR if no selection is made.
            element bulk_formulae {
                ## Bulk formulae from Large and Yeager (2004)
                ## Large, W. G. & Yeager, S. G. Diurnal to decadal global forcing 
                ## for ocean and sea-ice models: The data sets and flux climatologies. NCAR TR.
                element type {
                    attribute name { "NCAR" },
                    empty
                }|
                ## Bulk formulae from Fairall et al (2003). 
                ## Bulk Parameterization of Air–Sea Fluxes: Updates and 
                ## Verification for the COARE Algorithm Journal of Climate, 2003, 16, 571-591
                element type {
                    attribute name { "COARE" },
                    empty
                }|
                ## Bulk formulae from Kara et al (2005).
                ## Stability-Dependent Exchange Coefficients for Air–Sea Fluxes
                ## Journal of Atmospheric and Oceanic Technology, 2005, 22, 1080-1094
                element type {
                    attribute name { "Kara" },
                    empty
                }
            }?,
            ## The netCDF data file downloaded from ERA-40 reanalysis website
            ## (see above)
            element input_file {
               attribute file_name {xsd:string} 
            },
            element input_file_type {
                ## What kind of file is this? Currently only ERA40 files are supported
                element type {
                    attribute name { "ERA40" },
                    ## The data from ERA40 website included accumulated values (ppt, ro, ssrd, strd).
                    ## If these values have been already ammended to instantaneous values, then switch this
                    ## flag on and the accumulation correction will not be applied.
                    element no_accumulation {
                        empty
                    }?
                }|
                element type {
                    attribute name { "NSEP" },
                    empty
                }|
                element type {
                    attribute name { "ICOM" },
                    empty
                }
            }?,
            ## Adding a latitude and longitude here (specified as two real numbers)
            ## will obtain data from the forcing file at that location. The
            ## mesh is not translated nor is the mesh put on the sphere, instead the
            ## specified lat/long is translated into cartesian coordinates and this is simply
            ## added to the surface mesh node coordinates when fluxes are calculated.
            element position {
                real_vector,
                ## Turning on this option will cause all nodes on the surface mesh to
                ## experience the same forcing, regardless of position. Only really
                ## useful for psuedo-1D simulations.
                element single_location {
                    empty
                }?
            }?,
            ## Ouput some extra diagnostic fields for the momentum, temperature and salinity fluxes
            ## These <b>must</b> be on the velocity mesh
            element output_fluxes_diagnostics {
                element vector_field {
                    attribute rank { "1" },
                    attribute name { "MomentumFlux" },
                    (
                        element diagnostic {
                            internal_algorithm,
                            velocity_mesh_choice,
                            diagnostic_scalar_field
                        }|
                        element aliased {
                            generic_aliased_field
                        }
                    )
                }?, 
                element scalar_field {
                    attribute rank { "0" },
                    attribute name { "HeatFlux" },
                    (
                        element diagnostic {
                            internal_algorithm,
                            velocity_mesh_choice,
                            diagnostic_scalar_field
                        }|
                        element aliased {
                            generic_aliased_field
                        }
                    )
                }?,            
                element scalar_field {
                    attribute rank { "0" },
                    attribute name { "SalinityFlux" },
                    (
                        element diagnostic {
                            internal_algorithm,
                            velocity_mesh_choice,
                            diagnostic_scalar_field
                        }|
                        element aliased {
                            generic_aliased_field
                        }
                    )
                }?,
                element scalar_field {
                    attribute rank { "0" },
                    attribute name { "PhotosyntheticRadiationDownward" },
                    (
                        element diagnostic {
                            internal_algorithm,
                            velocity_mesh_choice,
                            diagnostic_scalar_field
                        }|
                        element aliased {
                            generic_aliased_field
                        }
                    )
                }?
            }
        }?,
        ## Tidal forcing options 
        element tidal_forcing {
               ## M2
               element M2 {
                 ## Ancient frequencies:
                 ##
                 ## Ma     Frequency
                 ##
                 ## 0        1.405e-4
                 ##
                 ## 10      1.408e-4
                 ##
                 ## 50      1.416e-4
                 ##
                 ## 100    1.423e-4
                 ##
                 ## 200    1.428e-4
                 ##
                 ## 300    1.431e-4
                 ##
                 ## 350    1.434e-4
                 ##
                 ## 400    1.436e-4
                 ##
                 ## 450    1.441e-4
                 ##
                 ## 500    1.451e-4
                 ##
                 ## 570    1.469e-4
                 ##
                 ## From Poliakow,2005
                 element frequency {
                    real
                  }?,
                 ## Ma     Amplitude
                 ##
                 ## 0        0.2423
                 ##
                 ## 10      0.2428
                 ##
                 ## 50      0.2446
                 ##
                 ## 100    0.2460
                 ##
                 ## 200    0.2471
                 ##
                 ## 300    0.2477
                 ##
                 ## 350    0.2485
                 ##
                 ## 400    0.2489
                 ##
                 ## 450    0.2499
                 ##
                 ## 500    0.2521
                 ##
                 ## 570    0.2560
                 ##
                 ## From Poliakow,2005
                 element amplitude {
                    real
                 }?
               }?,
               ## S2
               element S2 {
                  empty
               }?,
               ## N2
               element N2 {
                  empty
               }?,
               ## K2
               element K2 {
                  empty
               }?,
               ## K1
               element K1 {
                  empty
               }?,
               ## O1
               element O1 {
                  empty
               }?,
               ## P1
               element P1 {
                  empty
               }?,
               ## Q1
               element Q1 {
                  empty
               }?,        
               ## Mf
               element Mf {
                  empty
               }?,  
               ## Mm
               element Mm {
                  empty
               }?,  
               ## Ssa
               element Ssa {
                  empty
               }?,  
               ## Switch on all tidal components
               element all_tidal_components {
                  empty
               }?,
               ## Sets a user defined Love number. If not active a value of 1.0 is used.
               ## Recommended value is 0.693
               element love_number {
                  element value {
                     real
                  }
               }?,
               ## Self attraction and loading term (SAL). This is a simple implementation
               ## that uses a constant beta value
               element sal {
                element beta {
                    real
                }
              }?,
              ## Add the CHI term to the harmonics. You probably only want to do this is you're simulating a specific time
              element chi {
                empty
              }?

        }?
    }
  )

biology =
   (
      ## Model of biological processes in the ocean.
      element ocean_biology{
         ## A simple model of phytoplankton, zooplankton, general nutrient and detritus. 
         element pznd {
            (
               ## Python code specifying the source and sink relationships 
               ## between the biological tracers. This is usually achieved by 
               ## importing fluidity.ocean_biology and calling a scheme from there. 
               element source_and_sink_algorithm {
                  python_code
               }|
               ## Do not calculate sources and sinks. 
               ## This option is generally only useful for testing. 
               element disable_sources_and_sinks {
                  empty
               }
             ),
            ## Photosynthetically Active Radiation (PAR)
            element scalar_field {
               attribute rank { "0" },
               attribute name { "PhotosyntheticRadiation" },
               (
                  element prognostic {
                     velocity_mesh_choice,
                     prognostic_photosynthetic_radiation
                  }|
                  element prescribed {
                     velocity_mesh_choice,
                     prescribed_scalar_field
                  }
               )
            }
         }|
         ## 6 component biology model, which models Nitrates, Ammonium, 
         ## Phytoplankton, Zooplankton, Detritus and Chlorophyll.
         ##
         ## These fields must be enabled in the material phase
         ##
         ## Based on the equations in
         ## Popova, E. E.; Coward, A. C.; Nurser, G. A.; de Cuevas, B.; Fasham, M. J. R. & Anderson, T. R. 
         ## Mechanisms controlling primary and new production in a global ecosystem model - Part I: 
         ## Validation of the biological simulation Ocean Science, 2006, 2, 249-266. 
         ## DOI: 10.5194/os-2-249-2006
         element six_component {
            (
               ## Python code specifying the biology model. This takes
               ## in velocity and light and outputs Phytoplankton, if
               ## those fields exist.
               element source_and_sink_algorithm {
                  python_code
               }|
               ## Do not calculate biology
               ## This option is generally only useful for testing. 
               element disable_sources_and_sinks {
                  empty
               }
             ),
            ## Photosynthetically Active Radiation (PAR)
            element scalar_field {
               attribute rank { "0" },
               attribute name { "PhotosyntheticRadiation" },
               (
                  element prognostic {
                     velocity_mesh_choice,
                     prognostic_photosynthetic_radiation
                  }|
                  element prescribed {
                     velocity_mesh_choice,
                     prescribed_scalar_field
                  }
               )
            }?
         }
      }
   )


prognostic_photosynthetic_radiation =
   (
      ## PAR equation.
      element equation { 
         attribute name { "PhotosyntheticRadiation" }
      },
      ## Spatial discretisation options
      element spatial_discretisation {
         (
            ## Discontinuous galerkin formulation. You can also use this
            ## formulation with a continuous field in which case a simple
            ## galerkin formulation will result. 
            element discontinuous_galerkin {
               empty
            }
         )
      },
      (
         ## Solver
         element solver {
            linear_solver_options_asym_scalar
         }
      ),
      # Alas, no initial_condition either, so we'd better not checkpoint it...
      ## Disables checkpointing of this field
      element exclude_from_checkpointing {
        comment
      },
      ## Coefficients of absorption of photosynthetically active
      ## radiation for water and phytoplankton.
      element absorption_coefficients {
         ## Photosynthetically active radiation absorption coefficient for water.
         element water {
            real
         },
         ## Photosynthetically active radiation absorption coefficient for water.
         element phytoplankton {
            real
         }
      },
      ## Boundary conditions
      element boundary_conditions {
         attribute name { string },
         ## Surface id:
         element surface_ids {
            integer_vector
         },
         ## Type
         (
            element type {
               attribute name { "dirichlet" },
               ## Apply the dirichlet bc weakly.  Only available with
               ## discontinuous_galerkin and control_volume
               ## spatial_discretisations.
               ##
               ## If not selected boundary conditions are applied strongly.
               element apply_weakly {
                 empty
               },
               input_choice_real
            }|
            element type {
               attribute name { "bulk_formulae" },
               empty
            }|
            element type {
               attribute name { "neumann" },
               input_choice_real
            }
         )
      }?,
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options,
      scalar_convergence_options,
      prognostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_scalar_field,
      ## Set the priority of this field
      ## This determines the order in which scalar_fields are solved for:
      ##  - higher numbers have the highest priority
      ##  - lower numbers (including negative) have the lowest priority
      ##  - default if not set is 0
      element priority {
         integer
      }?
   )

recalculation_options =
   (
      ## Prevent this field from being recalculated at every timestep.
      ## This is cheaper especially if you are enforcing discrete properties on the field.
      element do_not_recalculate {
        empty
      }
   )

discrete_properties_algorithm_scalar =
   (
      ## Select discrete properties to enforce on the field
      ## either after being prescribed or interpolated
      element enforce_discrete_properties {
        ## Update this field using the lagrangian multiplier
        ## calculated in the solenoidal projection of a
        ## scalar field.
        ##
        ## Note this field must be specified as the update field
        ## underneath that vector field too.
        ##
        ## Note also this only really makes sense for coupled
        ## fields like velocity and pressure.
        element solenoidal_lagrange_update {
          empty
        }?
      }
   )

discrete_properties_algorithm_vector =
   (
      ## Select discrete properties to enforce on the field
      ## either after being prescribed or interpolated
      element enforce_discrete_properties {
        solenoidal_options?
      }
   )

solenoidal_options =
    ## Constrained divergence-free projection.
    ## This adds an additional constraint that ensures that the field
    ## is solenoidal, i.e. divergence-free.
    ## This is equivalent in cost to a pressure solve.
    ## This is expensive, and thus best left until
    ## needed.
    ##
    ## Note well: this only makes sense for nondivergent
    ## vector fields, such as incompressible velocity!
    element solenoidal {
      ## Options for the mass matrix of the field being interpolated
      element interpolated_field {
        (
          element continuous {
            ## Lump the mass matrix for the assembly of the projection matrix
            ## (not for the initial galerkin projection)
            ##
            ## Required when using interpolating continuous fields
            element lump_mass_matrix {
              ## Lump on the submesh.
              ## This only works for simplex meshes and is only
              ## strictly valid on 2d meshes.
              element use_submesh {
                empty
              }?
            }
          }|
          element discontinuous {
            ## Lump the mass matrix for the assembly of the projection matrix
            ## (not for the initial galerkin projection)
            element lump_mass_matrix {
              empty
            }?
          }
        )
      },
      ## Options for the lagrange multiplier
      ##
      ## Must be on a continuous mesh!
      element lagrange_multiplier {
        pressure_mesh_choice,
        element spatial_discretisation {
          (
            element continuous_galerkin {
              ## Remove the stabilisation term from the projection operator.
              ##
              ## Automatic when not using P1P1.
              element remove_stabilisation_term {
                empty
              }?,
              ## Integrate the divergence operator by parts.
              ##
              ## Automatic when projecting a discontinuous field
              element integrate_divergence_by_parts {
                empty
              }?,
              ## Test the divergence equation with the control volume dual mesh
              ## of the finite element lagrange multiplier mesh.
              ## This will make the lagrange multiplier matrix non symmetric which must be 
              ## considered when selecting the lagrange multiplier  solver options.
              element test_divergence_with_cv_dual {
                 comment
              }?
            }|
            element control_volumes {
              empty
            }
          )
        },
        element reference_node {
          integer
        }?,
        ## **UNDER DEVELOPMENT**
        ## This searches the CMC matrix diagonal looking for nodes that are less than the maximum value time epsilon(0.0) (i.e. nodes that are effectively zero).
        ## It then zeros that row and column and places a one on the diagonal and a zero on the rhs.
        ## At a debug level of 2 it also prints out the value and the sum of the row values.
        ## This is useful as a debugging tool if PETSc complains about zeros on the diagonal (i.e. if you have a stiff node in your mesh) but doesn't necessary produce nice answers at the end.
        element repair_stiff_nodes {
           empty
        }?,
        (
          ## Update a scalar field using the lagrange multiplier from
          ## the divergence free projection of this field.  The selected
          ## scalar field must have solenoidal selected in its interpolation
          ## options too and it must be on the same mesh as used for the
          ## solenoidal projection above.
          ##
          ## Note well: this only really makes sense for scalar fields linked to nondivergent
          ## vector fields, such as pressure to incompressible velocity!                  
          element update_scalar_field {
            attribute name { "Pressure" },
            empty
          }|
          ## Update a scalar field using the lagrange multiplier from
          ## the divergence free projection of this field.  The selected
          ## scalar field must have solenoidal selected in its interpolation
          ## options too and it must be on the same mesh as used for the
          ## solenoidal projection above.
          ##
          ## Note well: this only really makes sense for scalar fields linked to nondivergent
          ## vector fields, such as pressure to incompressible velocity!                  
          element update_scalar_field {
            attribute name { string },
            empty
          }
        )?,
        ## Solver options for the linear solve.
        ## This method requires the inversion of a projection matrix.
        element solver {
          linear_solver_options_sym
        }
      }
    }

represcribe_before_interpolation =
    ## Represcribe the field before interpolation.
    ##
    ## This means the interpolation will not be conservative from the previous mesh so be careful what you're trying to achieve!
    element represcribe_before_interpolation {
      empty
    }