text:
module LakeModule

using HierarchicalStateMachineModule: State, Event
using UserExceptionModule: UserError
using LakeModelParametersModule: LakeModelParameters
using MergeTypesModule
using HDModule: RiverParameters, PrognosticFields, RiverPrognosticFields, RiverDiagnosticFields
using HDModule: RiverDiagnosticOutputFields, PrintResults, PrintSection
using HDModule: print_river_results, get_river_parameters, get_river_fields
using CoordsModule: Coords,LatLonCoords,LatLonSectionCoords,is_lake
using CoordsModule: Generic1DCoords
using GridModule: Grid, LatLonGrid, for_all, for_all_fine_cells_in_coarse_cell, for_all_with_line_breaks
using GridModule: for_section_with_line_breaks,find_coarse_cell_containing_fine_cell
using FieldModule: Field, set!,elementwise_divide
import HDModule: water_to_lakes,water_from_lakes
import HierarchicalStateMachineModule: handle_event

abstract type Lake <: State end

struct RunLakes <: Event end

struct AddWater <: Event
 inflow::Float64
end

struct RemoveWater <: Event
  outflow::Float64
end

struct StoreWater <: Event
 inflow::Float64
end

struct AcceptMerge <: Event
  redirect_coords::Coords
end

struct AcceptSplit <: Event end

struct DrainExcessWater <: Event end

struct ReleaseNegativeWater <: Event end

struct WriteLakeNumbers <: Event
  timestep::Int64
end

struct WriteDiagnosticLakeVolumes <: Event
  timestep::Int64
end

struct WriteLakeVolumes <: Event end

struct SetupLakes <: Event
  initial_water_to_lake_centers::Field{Float64}
end

struct DistributeSpillover <: Event
  initial_spillover_to_rivers::Field{Float64}
end

struct CheckWaterBudget <: Event
  total_initial_water_volume::Float64
end

CheckWaterBudget() = CheckWaterBudget(0.0)

abstract type GridSpecificLakeParameters end

struct LakeParameters
  lake_centers::Field{Bool}
  connection_volume_thresholds::Field{Float64}
  flood_volume_thresholds::Field{Float64}
  flood_only::Field{Bool}
  flood_local_redirect::Field{Bool}
  connect_local_redirect::Field{Bool}
  additional_flood_local_redirect::Field{Bool}
  additional_connect_local_redirect::Field{Bool}
  merge_points::Field{MergeTypes}
  basin_numbers::Field{Int64}
  number_fine_grid_cells::Field{Int64}
  basins::Vector{Vector{Coords}}
  grid::Grid
  hd_grid::Grid
  surface_model_grid::Grid
  grid_specific_lake_parameters::GridSpecificLakeParameters
  lake_model_parameters::LakeModelParameters
  function LakeParameters(lake_centers::Field{Bool},
                          connection_volume_thresholds::Field{Float64},
                          flood_volume_thresholds::Field{Float64},
                          flood_local_redirect::Field{Bool},
                          connect_local_redirect::Field{Bool},
                          additional_flood_local_redirect::Field{Bool},
                          additional_connect_local_redirect::Field{Bool},
                          merge_points::Field{MergeTypes},
                          grid::Grid,
                          hd_grid::Grid,
                          surface_model_grid::Grid,
                          grid_specific_lake_parameters::GridSpecificLakeParameters)
    flood_only::Field{Bool} =  Field{Bool}(grid,false)
    number_fine_grid_cells = Field{Int64}(surface_model_grid,0)
    for_all(grid) do coords::Coords
      if connection_volume_thresholds(coords) == -1.0
        set!(flood_only,coords,true)
      end
      surface_model_coords::Coords =
        get_corresponding_surface_model_grid_cell(coords,grid_specific_lake_parameters)
      set!(number_fine_grid_cells,
           surface_model_coords,
           number_fine_grid_cells(surface_model_coords) + 1)
    end
    basins::Vector{Vector{Coords}} = Vector{Vector{Coords}}[]
    basin_numbers::Field{Int64} = Field{Int64}(hd_grid)
    basin_number::Int64 = 1
    for_all(hd_grid) do coords::Coords
      basins_in_coarse_cell::Vector{Coords} = Vector{Coords}[]
      basins_found::Bool = false
      for_all_fine_cells_in_coarse_cell(grid,hd_grid,coords) do fine_coords::Coords
        if lake_centers(fine_coords)
          push!(basins_in_coarse_cell,fine_coords)
          basins_found = true
        end
      end
      if basins_found
        push!(basins,basins_in_coarse_cell)
        set!(basin_numbers,coords,basin_number)
        basin_number += 1
      end
    end
    return new(lake_centers,
               connection_volume_thresholds,
               flood_volume_thresholds,
               flood_only,
               flood_local_redirect,
               connect_local_redirect,
               additional_flood_local_redirect,
               additional_connect_local_redirect,
               merge_points,
               basin_numbers,
               number_fine_grid_cells,
               basins,grid,hd_grid,
               surface_model_grid,
               grid_specific_lake_parameters,
               LakeModelParameters())
  end
end

struct LatLonLakeParameters <: GridSpecificLakeParameters
  flood_next_cell_lat_index::Field{Int64}
  flood_next_cell_lon_index::Field{Int64}
  connect_next_cell_lat_index::Field{Int64}
  connect_next_cell_lon_index::Field{Int64}
  flood_force_merge_lat_index::Field{Int64}
  flood_force_merge_lon_index::Field{Int64}
  connect_force_merge_lat_index::Field{Int64}
  connect_force_merge_lon_index::Field{Int64}
  flood_redirect_lat_index::Field{Int64}
  flood_redirect_lon_index::Field{Int64}
  connect_redirect_lat_index::Field{Int64}
  connect_redirect_lon_index::Field{Int64}
  additional_flood_redirect_lat_index::Field{Int64}
  additional_flood_redirect_lon_index::Field{Int64}
  additional_connect_redirect_lat_index::Field{Int64}
  additional_connect_redirect_lon_index::Field{Int64}
  corresponding_surface_cell_lat_index::Field{Int64}
  corresponding_surface_cell_lon_index::Field{Int64}
end

struct UnstructuredLakeParameters <: GridSpecificLakeParameters
  flood_next_cell_index::Field{Int64}
  connect_next_cell_index::Field{Int64}
  flood_force_merge_index::Field{Int64}
  connect_force_merge_index::Field{Int64}
  flood_redirect_index::Field{Int64}
  connect_redirect_index::Field{Int64}
  additional_flood_redirect_index::Field{Int64}
  additional_connect_redirect_index::Field{Int64}
  corresponding_surface_cell_index::Field{Int64}
end

struct LakeFields
  connected_lake_cells::Field{Bool}
  flooded_lake_cells::Field{Bool}
  lake_numbers::Field{Int64}
  buried_lake_numbers::Field{Int64}
  water_to_lakes::Field{Float64}
  water_to_hd::Field{Float64}
  lake_water_from_ocean::Field{Float64}
  number_lake_cells::Field{Int64}
  cells_with_lakes::Vector{Coords}
  function LakeFields(lake_parameters::LakeParameters)
    connected_lake_cells = Field{Bool}(lake_parameters.grid,false)
    flooded_lake_cells = Field{Bool}(lake_parameters.grid,false)
    lake_numbers = Field{Int64}(lake_parameters.grid,0)
    buried_lake_numbers = Field{Int64}(lake_parameters.grid,0)
    water_to_lakes = Field{Float64}(lake_parameters.hd_grid,0.0)
    water_to_hd = Field{Float64}(lake_parameters.hd_grid,0.0)
    lake_water_from_ocean = Field{Float64}(lake_parameters.hd_grid,0.0)
    number_lake_cells = Field{Int64}(lake_parameters.surface_model_grid,0)
    cells_with_lakes::Vector{Coords} = Vector{Coords}[]
    for_all(lake_parameters.hd_grid) do coords::Coords
      contains_lake::Bool = false
      for_all_fine_cells_in_coarse_cell(lake_parameters.grid,
                                        lake_parameters.hd_grid,coords) do fine_coords::Coords
        if lake_parameters.lake_centers(fine_coords) contains_lake = true end
      end
      if contains_lake
        push!(cells_with_lakes,coords)
      end
    end
    new(connected_lake_cells,flooded_lake_cells,
        lake_numbers,buried_lake_numbers,water_to_lakes,water_to_hd,
        lake_water_from_ocean,number_lake_cells,cells_with_lakes)
  end
end

struct LakePrognostics
  lakes::Vector{Lake}
  function LakePrognostics(lake_parameters::LakeParameters,
                           lake_fields::LakeFields)
    lakes::Vector{Lake} = Vector{Lake}[]
    lake_number::Int64 = 1
    for_all(lake_parameters.grid) do coords::Coords
      if lake_parameters.lake_centers(coords)
        coarse_coords = find_coarse_cell_containing_fine_cell(lake_parameters.grid,
                                                              lake_parameters.hd_grid,
                                                              coords)
        push!(lakes,FillingLake(lake_parameters,
                                LakeVariables(coords,lakes,
                                              lake_number,
                                              coarse_coords),
                                lake_fields))
        set!(lake_fields.lake_numbers,coords,lake_number)
        lake_number += 1
      end
    end
    return new(lakes)
  end
end

mutable struct LakeVariables
  center_cell::Coords
  lake_volume::Float64
  secondary_lake_volume::Float64
  unprocessed_water::Float64
  current_cell_to_fill::Coords
  previous_cell_to_fill::Coords
  other_lakes::Vector{Lake}
  lake_number::Int64
  filled_lake_cells::Vector{Coords}
  center_cell_coarse_coords::Coords
  LakeVariables(center_cell::Coords,
                other_lakes::Vector{Lake},
                lake_number::Int64,
                center_cell_coarse_coords::Coords) =
    new(center_cell,0.0,0.0,0.0,deepcopy(center_cell),deepcopy(center_cell),
        other_lakes,lake_number,Vector{Coords}[],center_cell_coarse_coords)
end

get_lake_parameters(obj::T) where {T <: Lake} =
  obj.lake_parameters::LakeParameters

get_lake_variables(obj::T) where {T <: Lake} =
  obj.lake_variables::LakeVariables

get_lake_fields(obj::T) where {T <: Lake} =
  obj.lake_fields::LakeFields

function handle_event(lake::Lake,inflow::Float64)
  throw(UserError())
end

mutable struct FillingLakeVariables
  primary_merge_completed::Bool
  use_additional_fields::Bool
end

struct FillingLake <: Lake
  lake_parameters::LakeParameters
  lake_variables::LakeVariables
  lake_fields::LakeFields
  filling_lake_variables::FillingLakeVariables
end

FillingLake(lake_parameters::LakeParameters,
            lake_variables::LakeVariables,
            lake_fields::LakeFields) =
  FillingLake(lake_parameters,
              lake_variables,
              lake_fields,
              FillingLakeVariables(false,false))

mutable struct OverflowingLakeVariables
  excess_water::Float64
end

OverflowingLakeVariables() = OverflowingLakeVariables(0.0)

struct OverflowingLake <: Lake
  lake_parameters::LakeParameters
  lake_variables::LakeVariables
  lake_fields::LakeFields
  outflow_redirect_coords::Coords
  local_redirect::Bool
  lake_retention_coefficient::Float64
  next_merge_target_coords::Coords
  overflowing_lake_variables::OverflowingLakeVariables
end

struct SubsumedLake <: Lake
  lake_parameters::LakeParameters
  lake_variables::LakeVariables
  lake_fields::LakeFields
  primary_lake_number::Int64
end

mutable struct LakeDiagnosticVariables
  total_lake_volume::Float64
end

LakeDiagnosticVariables() = LakeDiagnosticVariables(0.0)

struct RiverAndLakePrognosticFields <: PrognosticFields
  river_parameters::RiverParameters
  river_fields::RiverPrognosticFields
  river_diagnostic_fields::RiverDiagnosticFields
  river_diagnostic_output_fields::RiverDiagnosticOutputFields
  lake_parameters::LakeParameters
  lake_prognostics::LakePrognostics
  lake_fields::LakeFields
  lake_diagnostics_variables::LakeDiagnosticVariables
  using_lakes::Bool
  function RiverAndLakePrognosticFields(river_parameters::RiverParameters,
                               river_fields::RiverPrognosticFields,
                               lake_parameters::LakeParameters,
                               lake_prognostics::LakePrognostics,
                               lake_fields::LakeFields)
    for_all(river_parameters.grid) do coords::Coords
      if is_lake(river_parameters.flow_directions(coords))
        set!(river_parameters.cascade_flag,coords,false)
      end
    end
    new(river_parameters,river_fields,RiverDiagnosticFields(river_parameters),
        RiverDiagnosticOutputFields(river_parameters),lake_parameters,
        lake_prognostics,lake_fields,LakeDiagnosticVariables(),true)
  end
end

get_lake_prognostics(river_and_lake_prognostics::RiverAndLakePrognosticFields) =
                     river_and_lake_prognostics.lake_prognostics::LakePrognostics

get_lake_fields(river_and_lake_prognostics::RiverAndLakePrognosticFields) =
                     river_and_lake_prognostics.lake_fields::LakeFields

get_lake_parameters(river_and_lake_prognostics::RiverAndLakePrognosticFields) =
                     river_and_lake_prognostics.lake_parameters::LakeParameters

get_lake_diagnostic_variables(river_and_lake_prognostics::RiverAndLakePrognosticFields) =
  river_and_lake_prognostics.lake_diagnostics_variables::LakeDiagnosticVariables

function water_to_lakes(prognostic_fields::RiverAndLakePrognosticFields,coords::Coords,
                        inflow::Float64,step_length::Float64)
  lake_fields = get_lake_fields(prognostic_fields)
  set!(lake_fields.water_to_lakes,coords,inflow*step_length)
end

function water_from_lakes(prognostic_fields::RiverAndLakePrognosticFields,
                          step_length::Float64)
  lake_fields = get_lake_fields(prognostic_fields)
  return lake_fields.water_to_hd/step_length,
         lake_fields.lake_water_from_ocean/step_length
end

function handle_event(prognostic_fields::RiverAndLakePrognosticFields,setup_lakes::SetupLakes)
  lake_fields::LakeFields = get_lake_fields(prognostic_fields)
  lake_prognostics::LakePrognostics = get_lake_prognostics(prognostic_fields)
  lake_parameters::LakeParameters = get_lake_parameters(prognostic_fields)
  for_all(lake_parameters.grid) do coords::Coords
    if lake_parameters.lake_centers(coords)
      initial_water_to_lake_center::Float64 = setup_lakes.initial_water_to_lake_centers(coords)
      if  initial_water_to_lake_center > 0.0
        add_water::AddWater = AddWater(initial_water_to_lake_center)
        lake_index::Int64 = lake_fields.lake_numbers(coords)
        lake::Lake = lake_prognostics.lakes[lake_index]
        lake_prognostics.lakes[lake_index] = handle_event(lake,add_water)
      end
    end
  end
  return prognostic_fields
end

function handle_event(prognostic_fields::RiverAndLakePrognosticFields,
                      distribute_spillover::DistributeSpillover)
  lake_fields::LakeFields = get_lake_fields(prognostic_fields)
  lake_parameters::LakeParameters = get_lake_parameters(prognostic_fields)
  for_all(lake_parameters.hd_grid) do coords::Coords
    initial_spillover::Float64 = distribute_spillover.initial_spillover_to_rivers(coords)
    if initial_spillover > 0.0
      set!(lake_fields.water_to_hd,coords,lake_fields.water_to_hd(coords) + initial_spillover)
    end
  end
  return prognostic_fields
end

function handle_event(prognostic_fields::RiverAndLakePrognosticFields,run_lake::RunLakes)
  river_parameters::RiverParameters = get_river_parameters(prognostic_fields)
  river_fields::RiverPrognosticFields = get_river_fields(prognostic_fields)
  lake_fields::LakeFields = get_lake_fields(prognostic_fields)
  lake_parameters::LakeParameters = get_lake_parameters(prognostic_fields)
  lake_prognostics::LakePrognostics = get_lake_prognostics(prognostic_fields)
  drain_excess_water::DrainExcessWater = DrainExcessWater()
  fill!(lake_fields.water_to_hd,0.0)
  fill!(lake_fields.lake_water_from_ocean,0.0)
  local basins_in_cell::Array{Coords,1}
  for coords::Coords in lake_fields.cells_with_lakes
    if lake_fields.water_to_lakes(coords) > 0.0
      basins_in_cell =
        lake_parameters.basins[lake_parameters.basin_numbers(coords)]
      share_to_each_lake = lake_fields.water_to_lakes(coords)/length(basins_in_cell)
      add_water::AddWater = AddWater(share_to_each_lake)
      for basin_center::Coords in basins_in_cell
        lake_index::Int64 = lake_fields.lake_numbers(basin_center)
        lake::Lake = lake_prognostics.lakes[lake_index]
        lake_prognostics.lakes[lake_index] = handle_event(lake,add_water)
      end
    elseif lake_fields.water_to_lakes(coords) < 0.0
      basins_in_cell =
        lake_parameters.basins[lake_parameters.basin_numbers(coords)]
      share_to_each_lake = -1.0*lake_fields.water_to_lakes(coords)/length(basins_in_cell)
      remove_water::RemoveWater = RemoveWater(share_to_each_lake)
      for basin_center::Coords in basins_in_cell
        lake_index::Int64 = lake_fields.lake_numbers(basin_center)
        lake::Lake = lake_prognostics.lakes[lake_index]
        lake_prognostics.lakes[lake_index] = handle_event(lake,remove_water)
      end
    end
  end
  for lake::Lake in lake_prognostics.lakes
    lake_variables::LakeVariables = get_lake_variables(lake)
    lake_index::Int64 = lake_variables.lake_number
    if lake_variables.unprocessed_water > 0.0
      lake_prognostics.lakes[lake_index] = handle_event(lake,AddWater(0.0))
    end
    if lake_variables.lake_volume < 0.0
      lake_prognostics.lakes[lake_index] = handle_event(lake,ReleaseNegativeWater())
    end
    if isa(lake,OverflowingLake)
      handle_event(lake,drain_excess_water)
    end
  end
  return prognostic_fields
end

function handle_event(lake::FillingLake,add_water::AddWater)
  inflow::Float64 = add_water.inflow + lake.lake_variables.unprocessed_water
  lake.lake_variables.unprocessed_water = 0.0
  while inflow > 0.0
    inflow,filled::Bool = fill_current_cell(lake,inflow)
    if filled
      merge_type::SimpleMergeTypes = get_merge_type(lake)
      if merge_type != no_merge
        while true
          if ! (merge_type == primary_merge &&
                lake.filling_lake_variables.primary_merge_completed)
            if (merge_type == double_merge &&
                lake.filling_lake_variables.primary_merge_completed)
              merge_type = secondary_merge
              lake.filling_lake_variables.use_additional_fields = true
            end
            local merge_possible::Bool
            local already_merged::Bool
            merge_possible,already_merged = check_if_merge_is_possible(lake,merge_type)
            local subsumed_lake::SubsumedLake
            if merge_possible
              if merge_type == secondary_merge
                subsumed_lake = perform_secondary_merge(lake)
                subsumed_lake = handle_event(subsumed_lake,
                                             StoreWater(inflow))
                return subsumed_lake
              else
                perform_primary_merge(lake)
              end
            elseif ! already_merged
              overflowing_lake::Lake = change_to_overflowing_lake(lake,merge_type)
              overflowing_lake = handle_event(overflowing_lake,StoreWater(inflow))
              return overflowing_lake
            elseif merge_type == secondary_merge
              target_cell::Coords = get_secondary_merge_coords(lake)
              other_lake_number::Integer = lake.lake_fields.lake_numbers(target_cell)
              other_lake::SubsumedLake = lake.lake_variables.other_lakes[other_lake_number]
              other_lake_as_filling_lake::FillingLake = change_to_filling_lake(other_lake)
              lake.lake_variables.other_lakes[other_lake_number] =
                change_to_overflowing_lake(other_lake_as_filling_lake,primary_merge)
              subsumed_lake = perform_secondary_merge(lake)
              subsumed_lake = handle_event(subsumed_lake,
                                             StoreWater(inflow))
              return subsumed_lake
            elseif merge_type == double_merge
              lake.filling_lake_variables.primary_merge_completed = true
            end
          end
          if merge_type != double_merge
            break
          end
        end
      end
      update_filling_cell(lake)
    end
  end
  return lake
end

function handle_event(lake::FillingLake,remove_water::RemoveWater)
  if remove_water.outflow <= lake.lake_variables.unprocessed_water
    lake.lake_variables.unprocessed_water -= remove_water.outflow
    return lake
  end
  outflow::Float64 = remove_water.outflow - lake.lake_variables.unprocessed_water
  lake.lake_variables.unprocessed_water = 0.0
  while outflow > 0.0
    outflow,drained::Bool = drain_current_cell(lake,outflow)
    if drained
      rollback_filling_cell(lake)
      merge_type::SimpleMergeTypes = get_merge_type(lake)
      if merge_type == primary_merge || merge_type == double_merge
        new_outflow::Float64 = outflow/2.0
        rollback_primary_merge(lake,new_outflow)
        outflow -= new_outflow
        lake.filling_lake_variables.primary_merge_completed = false
        lake.filling_lake_variables.use_additional_fields = false
      end
    end
  end
  return lake
end

function handle_event(lake::OverflowingLake,add_water::AddWater)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  inflow::Float64 = add_water.inflow +
                    lake_variables.unprocessed_water
  lake_variables.unprocessed_water = 0.0
  if lake.local_redirect
    other_lake_number::Int64 = lake_fields.lake_numbers(lake.outflow_redirect_coords)
    other_lake::Lake =
      lake_variables.other_lakes[other_lake_number]
    if lake_parameters.lake_model_parameters.instant_throughflow
      lake_variables.other_lakes[other_lake_number] = handle_event(other_lake,
                                                             AddWater(inflow))
    else
      lake_variables.other_lakes[other_lake_number] = handle_event(other_lake,
                                                             StoreWater(inflow))
    end
  else
    lake.overflowing_lake_variables.excess_water += inflow
  end
  return lake_variables.other_lakes[lake_variables.lake_number]
end

function handle_event(lake::OverflowingLake,remove_water::RemoveWater)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  if remove_water.outflow <= lake_variables.unprocessed_water
    lake_variables.unprocessed_water -= remove_water.outflow
    return lake
  end
  outflow::Float64 = remove_water.outflow - lake_variables.unprocessed_water
  lake_variables.unprocessed_water = 0.0
  if outflow <= lake.overflowing_lake_variables.excess_water
    lake.overflowing_lake_variables.excess_water -= outflow
    return lake
  end
  outflow -= lake.overflowing_lake_variables.excess_water
  lake.overflowing_lake_variables.excess_water = 0
  if lake_fields.flooded_lake_cells(lake_variables.current_cell_to_fill)
    set!(lake_fields.flooded_lake_cells,lake_variables.current_cell_to_fill,false)
    surface_model_coords::Coords =
      get_corresponding_surface_model_grid_cell(lake_variables.current_cell_to_fill,
                                                lake_parameters.grid_specific_lake_parameters)
    set!(lake_fields.number_lake_cells,surface_model_coords,
         lake_fields.number_lake_cells(surface_model_coords) - 1)
  else
    set!(lake_fields.connected_lake_cells,lake_variables.current_cell_to_fill,false)
  end
  lake_as_filling_lake::FillingLake = change_to_filling_lake(lake)
  merge_type::SimpleMergeTypes = get_merge_type(lake)
  if merge_type == double_merge
    if lake_as_filling_lake.filling_lake_variables.primary_merge_completed
        new_outflow::Float64 = outflow/2.0
        rollback_primary_merge(lake,new_outflow)
        outflow -= new_outflow
    end
  end
  lake_as_filling_lake.filling_lake_variables.primary_merge_completed = false
  lake_as_filling_lake.filling_lake_variables.use_additional_fields = false
  new_remove_water::RemoveWater = RemoveWater(outflow)
  return handle_event(lake_as_filling_lake,new_remove_water)
end

function handle_event(lake::Lake,store_water::StoreWater)
  lake.lake_variables.unprocessed_water += store_water.inflow
  return lake
end

function handle_event(lake::SubsumedLake,add_water::AddWater)
  lake_variables::LakeVariables = get_lake_variables(lake)
  inflow::Float64 = add_water.inflow + lake_variables.unprocessed_water
  lake_variables.unprocessed_water = 0.0
  lake_variables.other_lakes[lake.primary_lake_number] =
    handle_event(lake_variables.other_lakes[lake.primary_lake_number],
                 AddWater(inflow))
  return lake
end

function handle_event(lake::SubsumedLake,remove_water::RemoveWater)
  lake_variables::LakeVariables = get_lake_variables(lake)
  if remove_water.outflow <= lake_variables.unprocessed_water
    lake_variables.unprocessed_water -= remove_water.outflow
    return lake
  end
  outflow::Float64 = remove_water.outflow - lake_variables.unprocessed_water
  lake_variables.unprocessed_water = 0.0
  lake_variables.other_lakes[lake.primary_lake_number] =
    handle_event(lake_variables.other_lakes[lake.primary_lake_number],
                 RemoveWater(outflow))
  return lake_variables.other_lakes[lake_variables.lake_number]::Lake
end

function handle_event(lake::Lake,release_negative_water::ReleaseNegativeWater)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  set!(lake_fields.lake_water_from_ocean,
       lake_variables.center_cell_coarse_coords,
       lake_fields.lake_water_from_ocean(lake_variables.center_cell_coarse_coords) -
       lake_variables.lake_volume)
  lake_variables.lake_volume = 0.0
  return lake
end

function handle_event(lake::Lake,accept_merge::AcceptMerge)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  if ! (lake_fields.connected_lake_cells(lake_variables.current_cell_to_fill) ||
        lake_parameters.flood_only(lake_variables.current_cell_to_fill))
    set!(lake_fields.connected_lake_cells,lake_variables.current_cell_to_fill,true)
  elseif ! lake_fields.flooded_lake_cells(lake_variables.current_cell_to_fill)
    set!(lake_fields.flooded_lake_cells,lake_variables.current_cell_to_fill,true)
    surface_model_coords::Coords =
      get_corresponding_surface_model_grid_cell(lake_variables.current_cell_to_fill,
                                                lake_parameters.grid_specific_lake_parameters)
    set!(lake_fields.number_lake_cells,surface_model_coords,
         lake_fields.number_lake_cells(surface_model_coords) + 1)
  end
  subsumed_lake::SubsumedLake = SubsumedLake(lake_parameters,
                                             lake_variables,
                                             lake_fields,
                                             lake_fields.
                                             lake_numbers(accept_merge.
                                                          redirect_coords))
  if isa(lake,OverflowingLake)
    subsumed_lake = handle_event(subsumed_lake,StoreWater(lake.
                                                          overflowing_lake_variables.
                                                          excess_water))
  end
  return subsumed_lake
end

function handle_event(lake::SubsumedLake,accept_split::AcceptSplit)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  if lake_fields.flooded_lake_cells(lake_variables.current_cell_to_fill)
      set!(lake_fields.flooded_lake_cells,lake_variables.current_cell_to_fill,false)
      surface_model_coords::Coords =
        get_corresponding_surface_model_grid_cell(lake_variables.current_cell_to_fill,
                                                  lake_parameters.grid_specific_lake_parameters)
      set!(lake_fields.number_lake_cells,surface_model_coords,
           lake_fields.number_lake_cells(surface_model_coords) - 1)
  else
    set!(lake_fields.connected_lake_cells,lake_variables.current_cell_to_fill,false)
  end
  filling_lake::FillingLake = FillingLake(lake_parameters,
                                          lake_variables,
                                          lake.lake_fields,
                                          FillingLakeVariables(false,
                                                               false))
  return filling_lake
end

function change_to_overflowing_lake(lake::FillingLake,merge_type::SimpleMergeTypes)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  if ! (lake_fields.connected_lake_cells(lake_variables.current_cell_to_fill) ||
        lake_parameters.flood_only(lake_variables.current_cell_to_fill))
    set!(lake_fields.connected_lake_cells,lake_variables.current_cell_to_fill,true)
  elseif ! lake_fields.flooded_lake_cells(lake_variables.current_cell_to_fill)
    set!(lake_fields.flooded_lake_cells,lake_variables.current_cell_to_fill,true)
    surface_model_coords::Coords =
      get_corresponding_surface_model_grid_cell(lake_variables.current_cell_to_fill,
                                                  lake_parameters.grid_specific_lake_parameters)
    set!(lake_fields.number_lake_cells,surface_model_coords,
         lake_fields.number_lake_cells(surface_model_coords) + 1)
  end
  outflow_redirect_coords,local_redirect =
    get_outflow_redirect_coords(lake)
  local next_merge_target_coords::Coords
  if merge_type == secondary_merge
    next_merge_target_coords = get_secondary_merge_coords(lake)
  else
    next_merge_target_coords = get_primary_merge_coords(lake)
  end
  return OverflowingLake(lake.lake_parameters,
                         lake_variables,
                         lake_fields,
                         outflow_redirect_coords,
                         local_redirect,
                         lake.lake_parameters.lake_model_parameters.lake_retention_constant,
                         next_merge_target_coords,
                         OverflowingLakeVariables())
end

function handle_event(lake::OverflowingLake,drain_excess_water::DrainExcessWater)
  lake_fields::LakeFields = get_lake_fields(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  if lake.overflowing_lake_variables.excess_water > 0.0
    lake.overflowing_lake_variables.excess_water +=
      lake_variables.unprocessed_water
    lake_variables.unprocessed_water = 0.0
    flow = (lake.overflowing_lake_variables.excess_water+
            lake_variables.lake_volume+
            lake_variables.secondary_lake_volume)/
           (lake.lake_retention_coefficient + 1.0)
    flow = min(flow,lake.overflowing_lake_variables.excess_water)
    set!(lake_fields.water_to_hd,lake.outflow_redirect_coords,
         lake_fields.water_to_hd(lake.outflow_redirect_coords)+flow)
    lake.overflowing_lake_variables.excess_water -= flow
  end
  return lake
end

function fill_current_cell(lake::FillingLake,inflow::Float64)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  current_cell_to_fill::Coords = lake_variables.current_cell_to_fill
  new_lake_volume::Float64 = inflow + lake_variables.lake_volume
  maximum_new_lake_volume::Float64 =
    (lake_fields.connected_lake_cells(current_cell_to_fill) ||
    lake_parameters.flood_only(current_cell_to_fill)) ?
    lake_parameters.flood_volume_thresholds(current_cell_to_fill) :
    lake_parameters.connection_volume_thresholds(current_cell_to_fill)
  if new_lake_volume <= maximum_new_lake_volume
    lake_variables.lake_volume = new_lake_volume
    return 0.0,false
  else
    inflow = new_lake_volume - maximum_new_lake_volume
    lake_variables.lake_volume = maximum_new_lake_volume
    return inflow,true
  end
end

function drain_current_cell(lake::FillingLake,outflow::Float64)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  if size(lake_variables.filled_lake_cells,1) == 0
    lake_variables.lake_volume -= outflow
    return 0.0,false
  end
  current_cell_to_fill::Coords = lake_variables.current_cell_to_fill
  previous_cell_to_fill::Coords = lake_variables.previous_cell_to_fill
  new_lake_volume::Float64 = lake_variables.lake_volume - outflow
  minimum_new_lake_volume::Float64 =
    lake_fields.flooded_lake_cells(previous_cell_to_fill)  ?
    lake_parameters.flood_volume_thresholds(previous_cell_to_fill) :
    lake_parameters.connection_volume_thresholds(previous_cell_to_fill)
  if new_lake_volume >= minimum_new_lake_volume
    lake_variables.lake_volume = new_lake_volume
    return 0.0,false
  else
    outflow = minimum_new_lake_volume - new_lake_volume
    lake_variables.lake_volume = minimum_new_lake_volume
    return outflow,true
  end
end

function get_merge_type(lake::Lake)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  current_cell::Coords = lake_variables.current_cell_to_fill
  extended_merge_type::MergeTypes = lake_parameters.merge_points(current_cell)
  if lake_fields.connected_lake_cells(current_cell) ||
     lake_parameters.flood_only(current_cell)
    return convert_to_simple_merge_type_flood[Int(extended_merge_type)+1]
  else
    return convert_to_simple_merge_type_connect[Int(extended_merge_type)+1]
  end
end

function check_if_merge_is_possible(lake::Lake,merge_type::SimpleMergeTypes)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  local target_cell::Coords
  if merge_type == secondary_merge
    target_cell = get_secondary_merge_coords(lake)
  else
    target_cell = get_primary_merge_coords(lake)
  end
  if (! (lake_fields.connected_lake_cells(target_cell) ||
         lake_parameters.flood_only(target_cell))) ||
     ( lake_parameters.flood_only(target_cell) &&
      ! lake_fields.flooded_lake_cells(target_cell))
    return false,false
  end
  other_lake::Lake = lake_variables.other_lakes[lake_fields.lake_numbers(target_cell)]
  other_lake = find_true_primary_lake(other_lake)
  if (get_lake_variables(other_lake)::LakeVariables).lake_number == lake_variables.lake_number
    return false,true
  elseif isa(other_lake,OverflowingLake)
    other_lake_target_lake_number::Int64 = lake_fields.lake_numbers(other_lake.next_merge_target_coords)
    while true
      if other_lake_target_lake_number == 0
        return false,false
      end
      other_lake_target_lake::Lake = lake_variables.other_lakes[other_lake_target_lake_number]
      if (get_lake_variables(other_lake_target_lake)::LakeVariables).lake_number == lake_variables.lake_number
        return true,false
      else
        if isa(other_lake_target_lake,OverflowingLake)
          other_lake_target_lake_number =
            lake_fields.lake_numbers(other_lake_target_lake.next_merge_target_coords)
        elseif isa(other_lake_target_lake,SubsumedLake)
          other_lake_target_lake_number =
            (get_lake_variables(find_true_primary_lake(other_lake_target_lake)::Lake)
                                ::LakeVariables).lake_number
        else
          return false,false
        end
      end
    end
  else
    return false,false
  end
end

function perform_primary_merge(lake::FillingLake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  target_cell::Coords = get_primary_merge_coords(lake)
  other_lake_number::Integer = lake_fields.lake_numbers(target_cell)
  other_lake::Lake = lake_variables.other_lakes[other_lake_number]
  other_lake = find_true_primary_lake(other_lake)
  lake_variables.secondary_lake_volume += (other_lake.lake_variables.lake_volume +
    other_lake.lake_variables.secondary_lake_volume)
  other_lake_number = other_lake.lake_variables.lake_number
  lake.filling_lake_variables.primary_merge_completed = true
  accept_merge::AcceptMerge = AcceptMerge(lake_variables.center_cell)
  lake_variables.other_lakes[other_lake_number] =
    handle_event(other_lake,accept_merge)
end

function perform_secondary_merge(lake::FillingLake)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  target_cell::Coords = get_secondary_merge_coords(lake)
  local surface_model_coords::Coords
  if ! (lake_fields.connected_lake_cells(lake_variables.current_cell_to_fill) ||
        lake_parameters.flood_only(lake_variables.current_cell_to_fill))
    set!(lake_fields.connected_lake_cells,
         lake_variables.current_cell_to_fill,true)
  elseif ! (lake_fields.flooded_lake_cells(lake_variables.current_cell_to_fill))
    set!(lake_fields.flooded_lake_cells,
         lake_variables.current_cell_to_fill,true)
    surface_model_coords =
      get_corresponding_surface_model_grid_cell(lake_variables.current_cell_to_fill,
                                                lake_parameters.grid_specific_lake_parameters)
    set!(lake_fields.number_lake_cells,surface_model_coords,
         lake_fields.number_lake_cells(surface_model_coords) + 1)
  end
  other_lake_number::Integer = lake_fields.lake_numbers(target_cell)
  other_lake::Lake = lake_variables.other_lakes[other_lake_number]
  other_lake = find_true_primary_lake(other_lake)
  other_lake.lake_variables.secondary_lake_volume += (lake_variables.lake_volume +
    lake_variables.secondary_lake_volume)
  other_lake_number = other_lake.lake_variables.lake_number
  other_lake_current_cell_to_fill::Coords = other_lake.lake_variables.current_cell_to_fill
  if ! (lake_fields.flooded_lake_cells(other_lake_current_cell_to_fill) ||
        lake_parameters.flood_only(other_lake_current_cell_to_fill))
    set!(lake_fields.connected_lake_cells,other_lake_current_cell_to_fill,false)
  else
    set!(lake_fields.flooded_lake_cells,other_lake_current_cell_to_fill,false)
    surface_model_coords =
      get_corresponding_surface_model_grid_cell(other_lake_current_cell_to_fill,
                                                lake_parameters.grid_specific_lake_parameters)
    set!(lake_fields.number_lake_cells,surface_model_coords,
         lake_fields.number_lake_cells(surface_model_coords) - 1)
  end
  other_lake_as_filling_lake::FillingLake = change_to_filling_lake(other_lake)
  lake_variables.other_lakes[other_lake_number] = other_lake_as_filling_lake
  accept_merge::AcceptMerge =
    AcceptMerge(get_lake_variables(other_lake_as_filling_lake).center_cell)
  return handle_event(lake,accept_merge)
end

function rollback_primary_merge(lake::FillingLake,other_lake_outflow::Float64)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  target_cell::Coords = get_primary_merge_coords(lake)
  other_lake_number::Integer = lake_fields.lake_numbers(target_cell)
  other_lake::Lake = lake_variables.other_lakes[other_lake_number]
  other_lake = find_true_rolledback_primary_lake(other_lake)
  lake_variables.secondary_lake_volume -=
    (other_lake.lake_variables.lake_volume +
     other_lake.lake_variables.secondary_lake_volume)
  other_lake_number = other_lake.lake_variables.lake_number
  accept_split::AcceptSplit = AcceptSplit()
  lake_variables.other_lakes[other_lake_number] =
    handle_event(other_lake,accept_split)
  other_lake = lake_variables.other_lakes[other_lake_number]
  remove_water::RemoveWater = RemoveWater(other_lake_outflow)
  lake_variables.other_lakes[other_lake_number] =
    handle_event(other_lake,remove_water)
end

function change_to_filling_lake(lake::OverflowingLake)
  lake.lake_variables.unprocessed_water +=
    lake.overflowing_lake_variables.excess_water
  use_additional_fields = get_merge_type(lake) == double_merge
  filling_lake::FillingLake = FillingLake(lake.lake_parameters,
                                          lake.lake_variables,
                                          lake.lake_fields,
                                          FillingLakeVariables(true,
                                                               use_additional_fields))
  return filling_lake
end

function change_to_filling_lake(lake::SubsumedLake)
  filling_lake::FillingLake = FillingLake(lake.lake_parameters,
                                          lake.lake_variables,
                                          lake.lake_fields,
                                          FillingLakeVariables(false,
                                                               false))
end

function update_filling_cell(lake::FillingLake)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  lake.filling_lake_variables.primary_merge_completed = false
  coords::Coords = lake_variables.current_cell_to_fill
  local new_coords::Coords
  if lake_fields.connected_lake_cells(coords) || lake_parameters.flood_only(coords)
    new_coords = get_flood_next_cell_coords(lake_parameters,coords)
  else
    new_coords = get_connect_next_cell_coords(lake_parameters,coords)
  end
  old_lake_number::Int64 = lake_fields.lake_numbers(new_coords)
  if old_lake_number != 0 &&
     old_lake_number != lake_variables.lake_number
     set!(lake_fields.buried_lake_numbers,new_coords,old_lake_number)
  end
  set!(lake_fields.lake_numbers,new_coords,lake_variables.lake_number)
  if ! (lake_fields.connected_lake_cells(coords) ||
        lake_parameters.flood_only(coords))
    set!(lake_fields.connected_lake_cells,coords,true)
  elseif ! lake_fields.flooded_lake_cells(coords)
    set!(lake_fields.flooded_lake_cells,coords,true)
    surface_model_coords::Coords =
      get_corresponding_surface_model_grid_cell(coords,
                                                lake_parameters.grid_specific_lake_parameters)
    set!(lake_fields.number_lake_cells,surface_model_coords,
         lake_fields.number_lake_cells(surface_model_coords) + 1)
  end
  push!(lake_variables.filled_lake_cells,coords)
  lake_variables.previous_cell_to_fill = coords
  lake_variables.current_cell_to_fill = new_coords
end

function rollback_filling_cell(lake::FillingLake)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  lake.filling_lake_variables.primary_merge_completed = false
  coords::Coords = lake_variables.current_cell_to_fill
  new_coords::Coords = lake_variables.previous_cell_to_fill
  if lake_parameters.flood_only(coords) ||
     ! lake_fields.connected_lake_cells(coords)
    set!(lake_fields.lake_numbers,coords,0)
  elseif lake_fields.buried_lake_numbers(coords) != 0
    set!(lake_fields.lake_numbers,coords,lake_fields.buried_lake_numbers(coords))
    set!(lake_fields.buried_lake_numbers,coords,0)
  end
  if ! (lake_fields.flooded_lake_cells(new_coords) ||
        lake_parameters.flood_only(new_coords))
    set!(lake_fields.connected_lake_cells,new_coords,false)
  else
    set!(lake_fields.flooded_lake_cells,new_coords,false)
    surface_model_coords::Coords =
      get_corresponding_surface_model_grid_cell(new_coords,
                                                lake_parameters.grid_specific_lake_parameters)
    set!(lake_fields.number_lake_cells,surface_model_coords,
         lake_fields.number_lake_cells(surface_model_coords) - 1)
  end
  lake_variables.current_cell_to_fill = new_coords
  pop!(lake_variables.filled_lake_cells)
  if size(lake_variables.filled_lake_cells,1) != 0
    lake_variables.previous_cell_to_fill = lake_variables.filled_lake_cells[end]
  else
    lake_variables.previous_cell_to_fill = new_coords
  end
end

function get_primary_merge_coords(lake::Lake)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  current_cell::Coords = lake_variables.current_cell_to_fill
  return get_primary_merge_coords(lake_parameters,
                                  current_cell,
                                  lake_fields.connected_lake_cells(current_cell) ||
                                  lake_parameters.flood_only(current_cell))
end

function get_secondary_merge_coords(lake::Lake)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  current_cell::Coords = lake_variables.current_cell_to_fill
  return get_secondary_merge_coords(lake_parameters,
                                    current_cell,
                                    lake_fields.connected_lake_cells(current_cell) ||
                                    lake_parameters.flood_only(current_cell))
end

function get_outflow_redirect_coords(lake::FillingLake)
  lake_parameters::LakeParameters = get_lake_parameters(lake)
  lake_variables::LakeVariables = get_lake_variables(lake)
  lake_fields::LakeFields = get_lake_fields(lake)
  current_cell::Coords = lake_variables.current_cell_to_fill
  connected_lake_cell::Bool = lake_fields.connected_lake_cells(current_cell) ||
                              lake_parameters.flood_only(current_cell)
  use_additional_fields = lake.filling_lake_variables.use_additional_fields
  local local_redirect::Bool
  if use_additional_fields
    local_redirect = connected_lake_cell ?
      lake_parameters.additional_flood_local_redirect(current_cell) :
      lake_parameters.additional_connect_local_redirect(current_cell)
  else
    local_redirect = connected_lake_cell ?
      lake_parameters.flood_local_redirect(current_cell) :
      lake_parameters.connect_local_redirect(current_cell)
  end
  return get_outflow_redirect_coords(lake_parameters,
                                     current_cell,
                                     connected_lake_cell,
                                     use_additional_fields),local_redirect
end

function find_true_primary_lake(lake::Lake)
  if isa(lake,SubsumedLake)
    return find_true_primary_lake(lake.lake_variables.other_lakes[lake.primary_lake_number])
  else
    return lake
  end
end

function find_true_rolledback_primary_lake(lake::Lake)
  next_lake::Lake = lake.lake_variables.other_lakes[lake.primary_lake_number]
  if isa(next_lake,SubsumedLake)
    return find_true_primary_lake(next_lake)
  else
    return lake
  end
end

function get_flood_next_cell_coords(lake_parameters::LakeParameters,initial_coords::Coords)
  throw(UserError())
end

function get_connect_next_cell_coords(lake_parameters::LakeParameters,initial_coords::Coords)
  throw(UserError())
end

function get_primary_merge_coords(lake_parameters::LakeParameters,initial_coords::Coords,
                                  connected_lake_cell::Bool)
  throw(UserError())
end

function get_secondary_merge_coords(lake_parameters::LakeParameters,initial_coords::Coords,
                                    connected_lake_cell::Bool)
  throw(UserError())
end

function get_outflow_redirect_coords(lake_parameters::LakeParameters,
                                     initial_coords::Coords,
                                     use_flood_redirect::Bool,
                                     use_additional_fields::Bool)
  throw(UserError())
end

function get_corresponding_surface_model_grid_cell(coords::Coords,
                                                   grid_specific_lake_parameters::GridSpecificLakeParameters)
 throw(UsersError())
end

function get_primary_merge_coords(lake_parameters::LakeParameters,initial_coords::LatLonCoords,
                                  connected_lake_cell::Bool)
  if connected_lake_cell
    return LatLonCoords(
           lake_parameters.grid_specific_lake_parameters.flood_force_merge_lat_index(initial_coords),
           lake_parameters.grid_specific_lake_parameters.flood_force_merge_lon_index(initial_coords))
  else
    return LatLonCoords(
           lake_parameters.grid_specific_lake_parameters.connect_force_merge_lat_index(initial_coords),
           lake_parameters.grid_specific_lake_parameters.connect_force_merge_lon_index(initial_coords))
  end
end

function get_primary_merge_coords(lake_parameters::LakeParameters,initial_coords::Generic1DCoords,
                                  connected_lake_cell::Bool)
  if connected_lake_cell
    return Generic1DCoords(
           lake_parameters.grid_specific_lake_parameters.flood_force_merge_index(initial_coords))
  else
    return Generic1DCoords(
           lake_parameters.grid_specific_lake_parameters.connect_force_merge_index(initial_coords))
  end
end

function get_secondary_merge_coords(lake_parameters::LakeParameters,initial_coords::LatLonCoords,
                                    connected_lake_cell::Bool)
  if connected_lake_cell
    return LatLonCoords(
           lake_parameters.grid_specific_lake_parameters.flood_next_cell_lat_index(initial_coords),
           lake_parameters.grid_specific_lake_parameters.flood_next_cell_lon_index(initial_coords))
  else
    return LatLonCoords(
           lake_parameters.grid_specific_lake_parameters.connect_next_cell_lat_index(initial_coords),
           lake_parameters.grid_specific_lake_parameters.connect_next_cell_lon_index(initial_coords))
  end
end

function get_secondary_merge_coords(lake_parameters::LakeParameters,initial_coords::Generic1DCoords,
                                    connected_lake_cell::Bool)
  if connected_lake_cell
    return Generic1DCoords(
           lake_parameters.grid_specific_lake_parameters.flood_next_cell_index(initial_coords))
  else
    return Generic1DCoords(
           lake_parameters.grid_specific_lake_parameters.connect_next_cell_index(initial_coords))
  end
end

function get_flood_next_cell_coords(lake_parameters::LakeParameters,initial_coords::LatLonCoords)
  return LatLonCoords(
    lake_parameters.grid_specific_lake_parameters.flood_next_cell_lat_index(initial_coords),
    lake_parameters.grid_specific_lake_parameters.flood_next_cell_lon_index(initial_coords))
end

function get_flood_next_cell_coords(lake_parameters::LakeParameters,initial_coords::Generic1DCoords)
  return Generic1DCoords(
    lake_parameters.grid_specific_lake_parameters.flood_next_cell_index(initial_coords))
end

function get_connect_next_cell_coords(lake_parameters::LakeParameters,initial_coords::LatLonCoords)
  return LatLonCoords(
    lake_parameters.grid_specific_lake_parameters.connect_next_cell_lat_index(initial_coords),
    lake_parameters.grid_specific_lake_parameters.connect_next_cell_lon_index(initial_coords))
end

function get_connect_next_cell_coords(lake_parameters::LakeParameters,initial_coords::Generic1DCoords)
  return Generic1DCoords(
    lake_parameters.grid_specific_lake_parameters.connect_next_cell_index(initial_coords))
end

function get_outflow_redirect_coords(lake_parameters::LakeParameters,
                                     initial_coords::LatLonCoords,
                                     use_flood_redirect::Bool,
                                     use_additional_fields::Bool)
  local lat::Int64
  local lon::Int64
  if use_additional_fields
    lat = use_flood_redirect ?
      lake_parameters.grid_specific_lake_parameters.additional_flood_redirect_lat_index(initial_coords) :
      lake_parameters.grid_specific_lake_parameters.additional_connect_redirect_lat_index(initial_coords)

    lon = use_flood_redirect ?
      lake_parameters.grid_specific_lake_parameters.additional_flood_redirect_lon_index(initial_coords) :
      lake_parameters.grid_specific_lake_parameters.additional_connect_redirect_lon_index(initial_coords)
  else
    lat = use_flood_redirect ?
      lake_parameters.grid_specific_lake_parameters.flood_redirect_lat_index(initial_coords) :
      lake_parameters.grid_specific_lake_parameters.connect_redirect_lat_index(initial_coords)

    lon = use_flood_redirect ?
      lake_parameters.grid_specific_lake_parameters.flood_redirect_lon_index(initial_coords) :
      lake_parameters.grid_specific_lake_parameters.connect_redirect_lon_index(initial_coords)
  end
  return LatLonCoords(lat,lon)
end

function get_outflow_redirect_coords(lake_parameters::LakeParameters,
                                     initial_coords::Generic1DCoords,
                                     use_flood_redirect::Bool,
                                     use_additional_fields::Bool)
  local index::Int64
  if use_additional_fields
    index = use_flood_redirect ?
      lake_parameters.grid_specific_lake_parameters.additional_flood_redirect_index(initial_coords) :
      lake_parameters.grid_specific_lake_parameters.additional_connect_redirect_index(initial_coords)
  else
    index = use_flood_redirect ?
      lake_parameters.grid_specific_lake_parameters.flood_redirect_index(initial_coords) :
      lake_parameters.grid_specific_lake_parameters.connect_redirect_index(initial_coords)
  end
  return Generic1DCoords(index)
end

function get_corresponding_surface_model_grid_cell(coords::LatLonCoords,
                                                   grid_specific_lake_parameters::LatLonLakeParameters)
  return LatLonCoords(grid_specific_lake_parameters.corresponding_surface_cell_lat_index(coords),
                      grid_specific_lake_parameters.corresponding_surface_cell_lon_index(coords))
end

function get_corresponding_surface_model_grid_cell(coords::Generic1DCoords,
                                                   grid_specific_lake_parameters::UnstructuredLakeParameters)
  return Generic1DCoords(grid_specific_lake_parameters.corresponding_surface_cell_index(coords))
end

function calculate_lake_fraction_on_surface_grid(lake_parameters::LakeParameters,
                                                 lake_fields::LakeFields)
  return elementwise_divide(lake_fields.number_lake_cells,
                            lake_parameters.number_fine_grid_cells)
end

function handle_event(prognostic_fields::RiverAndLakePrognosticFields,
                      print_results::PrintResults)
  lake_parameters::LakeParameters = get_lake_parameters(prognostic_fields)
  lake_prognostics::LakePrognostics = get_lake_prognostics(prognostic_fields)
  lake_fields::LakeFields = get_lake_fields(prognostic_fields)
  println("Timestep: $(print_results.timestep)")
  print_river_results(prognostic_fields)
  println("")
  println("Water to HD")
  println(lake_fields.water_to_hd)
  for lake in lake_prognostics.lakes
    println("Lake Center: $(lake.lake_variables.center_cell) "*
            "Lake Volume: $(lake.lake_variables.lake_volume)")
  end
  println("")
  print_lake_types(lake_parameters.grid,lake_prognostics,lake_fields)
  return prognostic_fields
end

function print_lake_types(grid::LatLonGrid,lake_prognostics::LakePrognostics,
                          lake_fields::LakeFields)
  for_all_with_line_breaks(grid) do coords::Coords
    print_lake_type(coords,lake_prognostics,lake_fields)
  end
end

function print_lake_type(coords::Coords,lake_prognostics::LakePrognostics,
                         lake_fields::LakeFields)
  lake_number::Int64 = lake_fields.lake_numbers(coords)
  if lake_number == 0
    print("- ")
  else
    lake::Lake = lake_prognostics.lakes[lake_number]
    if isa(lake,FillingLake)
      print("F ")
    elseif isa(lake,OverflowingLake)
      print("O ")
    elseif isa(lake,SubsumedLake)
      print("S ")
    else
      print("U ")
    end
  end
end

function handle_event(prognostic_fields::RiverAndLakePrognosticFields,
                      print_results::PrintSection)
  lake_parameters::LakeParameters = get_lake_parameters(prognostic_fields)
  lake_prognostics::LakePrognostics = get_lake_prognostics(prognostic_fields)
  lake_fields::LakeFields = get_lake_fields(prognostic_fields)
  println("")
  print_lake_types_section(lake_parameters.grid,lake_prognostics,lake_fields)
  return prognostic_fields
end

function print_lake_types_section(grid::LatLonGrid,lake_prognostics::LakePrognostics,
                                  lake_fields::LakeFields)
  section_coords::LatLonSectionCoords = LatLonSectionCoords(225,295,530,580)
  for_section_with_line_breaks(grid,section_coords) do coords::Coords
    print_lake_type(coords,lake_prognostics,lake_fields)
  end
end

function handle_event(prognostic_fields::RiverAndLakePrognosticFields,
                      write_lake_numbers::WriteLakeNumbers)
  lake_parameters::LakeParameters = get_lake_parameters(prognostic_fields)
  lake_fields::LakeFields = get_lake_fields(prognostic_fields)
  write_lake_numbers_field(lake_parameters,lake_fields,timestep=write_lake_numbers.timestep)
  return prognostic_fields
end

function write_lake_numbers_field(lake_parameters::LakeParameters,lake_fields::LakeFields;
                                  timestep::Int64=-1)
  throw(UserError())
end

function handle_event(prognostic_fields::RiverAndLakePrognosticFields,
                      write_lake_volumes::WriteLakeVolumes)
  lake_parameters::LakeParameters = get_lake_parameters(prognostic_fields)
  lake_prognostics::LakePrognostics = get_lake_prognostics(prognostic_fields)
  lake_volumes::Field{Float64} = Field{Float64}(lake_parameters.grid,0.0)
  for lake::Lake in lake_prognostics.lakes
    lake_variables::LakeVariables = get_lake_variables(lake)
    lake_center_cell::Coords = lake_variables.center_cell
    set!(lake_volumes,lake_center_cell,lake_variables.lake_volume +
                                      lake_variables.unprocessed_water)
  end
  write_lake_volumes_field(lake_parameters,lake_volumes)
  return prognostic_fields
end

function handle_event(prognostic_fields::RiverAndLakePrognosticFields,
                      write_diagnostic_lake_volumes::WriteDiagnosticLakeVolumes)
  lake_parameters::LakeParameters = get_lake_parameters(prognostic_fields)
  lake_fields::LakeFields = get_lake_fields(prognostic_fields)
  lake_prognostics::LakePrognostics = get_lake_prognostics(prognostic_fields)
  diagnostic_lake_volumes::Field{Float64} = calculate_diagnostic_lake_volumes_field(lake_parameters,
                                                                                lake_fields,
                                                                                lake_prognostics)
  write_diagnostic_lake_volumes_field(lake_parameters,diagnostic_lake_volumes,
                                      timestep=write_diagnostic_lake_volumes.timestep)
  return prognostic_fields
end

function calculate_diagnostic_lake_volumes_field(lake_parameters::LakeParameters,
                                                 lake_fields::LakeFields,
                                                 lake_prognostics::LakePrognostics)
  lake_volumes_by_lake_number::Vector{Float64} = zeros(Float64,size(lake_prognostics.lakes,1))
  diagnostic_lake_volumes::Field{Float64} = Field{Float64}(lake_parameters.grid,0.0)
  for lake::Lake in lake_prognostics.lakes
    lake_variables::LakeVariables = get_lake_variables(lake)
    original_lake_index::Int64 = lake_variables.lake_number
    if isa(lake,SubsumedLake)
      lake = find_true_primary_lake(lake)
    end
    lake_volumes_by_lake_number[original_lake_index] =
      lake.lake_variables.lake_volume + lake.lake_variables.secondary_lake_volume
  end
  for_all(lake_parameters.grid) do coords::Coords
    lake_number = lake_fields.lake_numbers(coords)
    if (lake_number > 0)
      set!(diagnostic_lake_volumes,coords,lake_volumes_by_lake_number[lake_number])
    end
  end
  return diagnostic_lake_volumes
end

function write_lake_volumes_field(lake_parameters::LakeParameters,lake_volumes::Field{Float64})
  throw(UserError())
end

function write_diagnostic_lake_volumes_field(lake_parameters::LakeParameters,
                                             diagnostic_lake_volumes::Field{Float64};
                                             timestep::Int64=-1)
  throw(UserError())
end

function handle_event(prognostic_fields::RiverAndLakePrognosticFields,
                      check_water_budget::CheckWaterBudget)
  lake_prognostics::LakePrognostics = get_lake_prognostics(prognostic_fields)
  lake_fields::LakeFields = get_lake_fields(prognostic_fields)
  lake_diagnostic_variables = get_lake_diagnostic_variables(prognostic_fields)
  new_total_lake_volume::Float64 = 0.0
  for lake::Lake in lake_prognostics.lakes
    lake_variables::LakeVariables = get_lake_variables(lake)
    new_total_lake_volume += lake_variables.unprocessed_water + lake_variables.lake_volume
    if isa(lake,OverflowingLake)
      new_total_lake_volume += lake.overflowing_lake_variables.excess_water
    end
  end
  change_in_total_lake_volume::Float64 = new_total_lake_volume -
                                         lake_diagnostic_variables.total_lake_volume -
                                         check_water_budget.
                                         total_initial_water_volume
  total_water_to_lakes::Float64 = sum(lake_fields.water_to_lakes)
  total_inflow_minus_outflow::Float64 = total_water_to_lakes +
                                        sum(lake_fields.lake_water_from_ocean) -
                                        sum(lake_fields.water_to_hd)
  difference::Float64 = change_in_total_lake_volume - total_inflow_minus_outflow
  tolerance::Float64 = 10e-16*max(new_total_lake_volume, total_water_to_lakes)
  if ! isapprox(difference,0,atol=tolerance)
    println("*** Lake Water Budget ***")
    println("Total lake volume: $(new_total_lake_volume)")
    println("Total inflow - outflow: $(total_inflow_minus_outflow)")
    println("Change in lake volume: $(change_in_total_lake_volume)")
    println("Difference: $(difference)")
    println(total_water_to_lakes)
    println(sum(lake_fields.lake_water_from_ocean))
    println(sum(lake_fields.water_to_hd))
    println(lake_diagnostic_variables.total_lake_volume)
  end
  lake_diagnostic_variables.total_lake_volume = new_total_lake_volume
  return prognostic_fields
end

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
struct JeffreysPrior{T<:Real,TD<:AbstractDesign} <: Prior
    pivots::Vector{T}
    weights::Vector{T}
    design::TD
    pdf::Vector{T}
    pdf_function::Function
    low::T
    high::T
end

function JeffreysPrior(design::TD; low::Real = 0, high::Real = 1)::JeffreysPrior{Float64,TD} where
    {TD<:AbstractDesign}

    eps          = 1e-3
    low, high    = max(low, eps), min(high, 1 - eps)
    p, ω         = gauss_legendre_25(low, high)
    z            = quadgk( p -> fisher_information(p, design), low, high; atol = 1e-5 )[1]
    pdf_function = p -> ( (p < low) | (p > high) ) ? 0.0 : fisher_information(p, design) / z
    pdf          = pdf_function.(p)
    zz           = (high - low)/2 * sum(pdf .* ω)
    pdf          = pdf ./ zz # normalize
    return JeffreysPrior{Float64,TD}(
        convert(Vector{Float64}, p),
        convert(Vector{Float64}, ω),
        design,
        convert(Vector{Float64}, pdf),
        pdf_function,
        convert(Float64, low),
        convert(Float64, high)
    )
end

function pdf(p::TR, prior::JeffreysPrior{TR,TD})::TR where
    {TR<:Real,TD<:AbstractDesign}

    return prior.pdf_function(p)
end

function cdf(p::TR, prior::JeffreysPrior{TR,TD})::TR where
    {TR<:Real,TD<:AbstractDesign}

    return min(1.0, max(0.0, quadgk(prior.pdf_function, 0, p; atol = 1e-5)[1]))
end

function expectation(f::Function, prior::JeffreysPrior{TR,TD})::TR where
    {TR<:Real,TD<:AbstractDesign}

    return (prior.high - prior.low)/2 * sum( prior.pdf .* f.(prior.pivots) .* prior.weights )
end

function mean(prior::JeffreysPrior{TR,TD})::TR where
    {TR<:Real,TD<:AbstractDesign}

    return expectation(p -> p, prior)
end

function condition(prior::JeffreysPrior{TR,TD}; low::TR = prior.low, high::TR = prior.high)::JeffreysPrior{TR,TD} where
    {TR<:Real,TD<:AbstractDesign}

    return JeffreysPrior{TR,TD}(TD; low = max(low, prior.low), high = min(high, prior.high))
end

function update(prior::JeffreysPrior{TR,TD}, x::TI, n::TI)::GenericDistribution{TR} where
    {TR<:Real,TI<:Integer,TD<:AbstractDesign}

    return GenericDistribution( p -> pdf(p, prior) * pmf(x, n, p) )
end

function string(prior::JeffreysPrior)
    eps          = 1e-3
    low, high    = max(prior.low, eps), min(prior.high, 1 - eps)
    if (prior.low > low) | (prior.high < high)
        @sprintf "JeffreysPrior|[%.2f,%.2f](%s)" prior.low prior.high prior.design
    else
        @sprintf "JeffreysPrior(%s)" prior.design
    end
end

ext: jl


text:
__precompile__()

module IntervalConstraintProgramming

using   IntervalArithmetic,
        IntervalRootFinding,
        IntervalContractors

using MacroTools

import Base:
    show, ∩, ∪, !, ⊆, setdiff

import IntervalArithmetic: sqr, setindex

export
    @contractor,
    Separator, separator, @separator, @constraint,
    @function,
    SubPaving, Paving,
    pave, refine!,
    Vol,
    show_code

const reverse_operations = IntervalContractors.reverse_operations

include("ast.jl")
include("code_generation.jl")
include("contractor.jl")
include("separator.jl")
include("paving.jl")
include("setinversion.jl")
include("volume.jl")
include("functions.jl")

end # module


text:
# Belief Propagation                #
# --------------------------------- #
@rule Wishart(:out, Marginalisation) (m_ν::PointMass, m_S::PointMass) = Wishart(mean(m_ν), mean(m_S))

# Variational                       #
# --------------------------------- #
@rule Wishart(:out, Marginalisation) (q_ν::Any, m_S::PointMass) = Wishart(mean(q_ν), mean(m_S))
@rule Wishart(:out, Marginalisation) (m_ν::PointMass, q_S::Any) = Wishart(mean(m_ν), mean(q_S))
@rule Wishart(:out, Marginalisation) (q_ν::Any, q_S::Any)       = Wishart(mean(q_ν), mean(q_S))


text:
using InterpNd
using PyPlot

function add_partial_waves(data, orbitals)
    ml = maximum(map(d -> size(d,2), data))
    mr = maximum(map(d -> maximum(d[1,:]), data))
    rsq = linspace(0,mr,ml)
    core = zeros(ml)
    for i in eachindex(orbitals)
        w = degeneracy(orbitals[i])
        ii = interp1d(vec(data[i][2,:]))
        core += w*abs(ii(vec(data[i][1,:]),rsq)).^2
    end

    rsq,core
end

function plot_wfn(config::Config, term::Term, ncorr)
    filename = active_file(config, term)
    conf_orbs = map(config) do o
        Orbital((o[1],o[2],1,"*"))
    end
    dir = "$filename/$ncorr"
    sep = "    0.0000       0.000"
    cp = r"([0-9]+)([a-z])\(([ 0-9]+?)\)"
    orbitals,data = dir_run(dir) do
        clist = map(filter(l -> ismatch(cp, l), eachline("$(filename).c"))) do l
            map(matchall(cp,l)) do m
                ref_set_list("$(m[1])$(m[2])")[1]
            end
        end
        orbitals = sort(unique(vcat(conf_orbs...,clist...)))
        rd,wr = redirect_stdin()
        for o in orbitals
            write(wr, " \n")
        end
        run(pipeline(rd, `$(atsp)/plotw $(filename).w`))
        orbitals,map(split(readall("plot.dat"), sep)[2:end]) do pw
            v = map([sep; split(pw, "\n")[2:end-4]]) do l
                map(f -> parse(Float64, f), split(l))
            end
            hcat(v...)
        end
    end
    subplot(211)
    for i in eachindex(orbitals)
        PyPlot.plot(data[i][1,:]',data[i][2,:]'.^2,label=string(orbitals[i]))
    end
    legend(framealpha=0.75)
    margins(0,0.1)
    subplot(212)
    PyPlot.plot(add_partial_waves(data,orbitals)...)
    margins(0,0.1)
end

export plot_wfn


text:
import MPI

include("../types/types.jl")
include("../setting/constants.jl")
include("../utils/readfiles.jl")
include("../utils/kdtree.jl")
include("../utils/parse_commandline.jl")

# ! this program uses Float64, which may cause the out of memory problem when In values_from_mesh.h
# ! approximate total number of points in entire mesh > 151837623.000000.
function run_interp(myrank::Int64, nrank::Int64, command_args::Dict{String,Any})
    # * init some variables
    isroot = (myrank == 0)
    nproc_old = command_args["nproc_old"]
    old_mesh_dir = command_args["old_mesh_dir"]
    old_model_dir = command_args["old_model_dir"]
    nproc_new = command_args["nproc_new"]
    new_mesh_dir = command_args["new_mesh_dir"]
    new_model_dir = command_args["new_model_dir"]
    model_tags = command_args["model_tags"]
    output_dir = command_args["output_dir"]

    # * init some variables
    mesh_new = sem_mesh_data()
    mesh_old = sem_mesh_data()

    # * parse model tags
    model_names = String.(split(model_tags, ","))
    nmodel = length(model_names)
    if isroot
        @info "[$(myrank)]# nmodel=$(nmodel)"
        @info "[$(myrank)]# model names=$(model_names)"
    end

    # * interpolate old mesh/model onto new mesh

    # typical element size at surface for old mesh
    # read from constants.jl for old mesh
    typical_size = deg2rad(max(ANGULAR_WIDTH_XI_IN_DEGREES_VAL / NEX_XI_VAL, ANGULAR_WIDTH_ETA_IN_DEGREES_VAL / NEX_ETA_VAL)) * R_UNIT_SPHERE
    max_search_dist = 10.0 * typical_size
    max_misloc = typical_size / 4.0

    # * loop each new mesh slice (where MPI parallels)
    for iproc_new in myrank:nrank:(nproc_new - 1)
        # if iproc_new != 386
        #     continue
        # end
        @info "[$(myrank)]# iproc_new=$(iproc_new)"

        # * read new mesh slice
        mesh_new = sem_mesh_read(new_mesh_dir, iproc_new)
        nspec_new = mesh_new.nspec
        ngll_new = NGLLX * NGLLY * NGLLZ * nspec_new

        xyz_new = zeros(Float64, 3, ngll_new)
        idoubling_new = zeros(Int64, ngll_new)
        xyz_center_new = zeros(Float64, 3, nspec_new)

        for ispec = 1:nspec_new
            iglob = mesh_new.ibool[MIDX,MIDY,MIDZ,ispec]
            xyz_center_new[:,ispec] = mesh_new.xyz_glob[:,iglob]
            for igllz = 1:NGLLZ
                for iglly = 1:NGLLY
                    for igllx = 1:NGLLX
                        igll = igllx + NGLLX * ( (iglly - 1) + NGLLY * ( (igllz - 1) + NGLLZ * ( (ispec - 1))))
                        iglob = mesh_new.ibool[igllx,iglly,igllz,ispec]
                        xyz_new[:,igll] = mesh_new.xyz_glob[:,iglob]
                        idoubling_new[igll] = mesh_new.idoubling[ispec]
                    end
                end
            end
        end

        # * initialize variables for interpolation
        location_1slice = Array{sem_mesh_location}(undef, ngll_new)
        stat_final = zeros(Int64, ngll_new)
        misloc_final = zeros(Float64, ngll_new)

        stat_final .= -1
        misloc_final .= HUGEVAL

        model_interp = zeros(Float64, nmodel, ngll_new)

        # * read in the new model as the background model
        model_gll_new = zeros(Float64, nmodel, NGLLX, NGLLY, NGLLZ, nspec_new)
        sem_io_read_gll_file_n!(new_model_dir, iproc_new, model_names, nmodel, model_gll_new)

        # * loop each slices of the old mesh
        flag = true
        for iproc_old = 0:nproc_old - 1
            @info "[$(myrank)]# iproc_old=$(iproc_old)"
            # read old mesh slice
            mesh_old = sem_mesh_read(old_mesh_dir, iproc_old)
            nspec_old = mesh_old.nspec
            # test if the new and old mesh slices are separated apart
            min_dist = HUGEVAL
            for ispec = 1:nspec_old
                iglob = mesh_old.ibool[MIDX,MIDY,MIDZ,ispec]
                min_dist = min(min_dist, sqrt(minimum((xyz_center_new[1,:] .- mesh_old.xyz_glob[1,iglob]).^2 + (xyz_center_new[2,:] .- mesh_old.xyz_glob[2,iglob]).^2 + (xyz_center_new[3,:] .- mesh_old.xyz_glob[3,iglob]).^2)))
            end
            if min_dist > max_search_dist
                @info "[$(myrank)]# $(iproc_old)/$(iproc_new) slices too far away, skip"
                continue
            end

            # read old model
            model_gll_old = zeros(Float64, nmodel, NGLLX, NGLLY, NGLLZ, nspec_old)
            sem_io_read_gll_file_n!(old_model_dir, iproc_old, model_names, nmodel, model_gll_old)

            # locate points in this mesh slice
            nnearest = 10
        #     location_1slice = sem_mesh_locate_kdtree2!(mesh_old, ngll_new, xyz_new, idoubling_new, nnearest, max_search_dist, max_misloc, iproc_old)
        #     for igll = 1:ngll_new
        #         if (stat_final[igll] == 1 && location_1slice[igll].stat == 1)
        #             @info "[$(myrank)]# multi-located, $(xyz_new[:,igll])"
        #             continue
        #         end
        #         # for point located inside one element in the first time or closer to one element than located before
        #         if location_1slice[igll].stat == 1 || (location_1slice[igll].stat == 0 && location_1slice[igll].misloc < misloc_final[igll])
        #             for imodel = 1:nmodel
        #                 model_interp[imodel,igll] = sum(location_1slice[igll].lagrange .* model_gll_old[imodel,:,:,:,location_1slice[igll].eid])
        #             end
        #             stat_final[igll] = location_1slice[igll].stat
        #             misloc_final[igll] = location_1slice[igll].misloc
        #         end
        #         # if location_1slice[igll].stat == 1 || (location_1slice[igll].stat == 0 && location_1slice[igll].misloc < misloc_final[igll])
        #             # flag = false
        #         # end
        #     end
        # end

        # * saved memory version
            kdtree, xyz_elem, xyz_anchor, xigll, yigll, zigll, hlagx, hlagy, hlagz = prepare_kdtree(mesh_old)
            for igll = 1:ngll_new
                # * loop for all the points
                sem_mesh_locate_kdtree_series!(igll,nnearest, max_search_dist, max_misloc, nmodel, kdtree, mesh_old, xyz_new, idoubling_new, xyz_elem, xyz_anchor, xigll, yigll, zigll, model_gll_old,hlagx,
                # modified values
                hlagy, hlagz, stat_final, misloc_final, model_interp)
            end
        end

        # * write out gll files for this new mesh slice

        # reshape model_interp to model_gll
        for ispec = 1:nspec_new
            for igllz = 1:NGLLZ
                for iglly = 1:NGLLY
                    for igllx = 1:NGLLX
                        igll = igllx + NGLLX * ( (iglly - 1) + NGLLY * ( (igllz - 1) + NGLLZ * ( (ispec - 1))))
                        if stat_final[igll] != -1
                            model_gll_new[:,igllx,iglly,igllz,ispec] = model_interp[:,igll]
                        end
                    end
                end
            end
        end
        sem_io_write_gll_file_n(output_dir, iproc_new, model_names, nmodel, model_gll_new)
    end
end

function main()
    MPI.Init()

    comm = MPI.COMM_WORLD
    myrank = MPI.Comm_rank(comm)
    nrank = MPI.Comm_size(comm)
    command_args = parse_commandline()

    run_interp(myrank, nrank, command_args)

    MPI.Barrier(comm)
    MPI.Finalize()
end

main()

max_stars_repo_name: ziyixi/seisflow

max_issues_repo_name: ziyixi/seisflow

max_forks_repo_name: ziyixi/seisflow


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct BgpPeerStatusListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{BgpPeerStatus} } # spec name: value

    function BgpPeerStatusListResult(;value=nothing)
        o = new()
        validate_property(BgpPeerStatusListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type BgpPeerStatusListResult

const _property_map_BgpPeerStatusListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"))
const _property_types_BgpPeerStatusListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{BgpPeerStatus}")
Base.propertynames(::Type{ BgpPeerStatusListResult }) = collect(keys(_property_map_BgpPeerStatusListResult))
Swagger.property_type(::Type{ BgpPeerStatusListResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_BgpPeerStatusListResult[name]))}
Swagger.field_name(::Type{ BgpPeerStatusListResult }, property_name::Symbol) =  _property_map_BgpPeerStatusListResult[property_name]

function check_required(o::BgpPeerStatusListResult)
    true
end

function validate_property(::Type{ BgpPeerStatusListResult }, name::Symbol, val)
end


text:
"""
    from_dict(::Type{T}, d::AbstractDict{String}; kw...) where T

Convert dictionary `d` to an option type `T`, the value of valid fields of `T`
in this dictionary `d` can be override by keyword arguments.

!!! compat "Configurations 0.17"

    `convert_to_option` interface is deprecated for conversion, please overload
    the 3-arg or 4-arg `from_dict` instead. See also [Type Conversion](@ref type-conversion)
    section.

# Example

```julia-repl
julia> @option struct OptionA
           name::String = "Sam"
           age::Int = 25
       end

julia> d = Dict{String, Any}(
           "name" => "Roger",
           "age" => 10,
       );

julia> from_dict(OptionA, d; age=25)
OptionA(;
    name = "Roger",
    age = 25,
)
```
"""
function from_dict(::Type{OptionType}, d::AbstractDict{String}; kw...) where {OptionType}
    if !isconcretetype(OptionType)
        throw(ArgumentError("expect a concrete type, got $OptionType"))
    end

    # deepcopy is quite expensive
    # error earlier before that
    assert_field_match_exactly(OptionType, d)

    if !isempty(kw)
        d = from_underscore_kwargs!(deepcopy(d), OptionType; kw...)
    end

    return from_dict_specialize(OptionType, d)
end

# An already-parsed option should be kept in nested parsing.
from_dict(::Type{T}, t::T) where {T} = t

"""
    ignore_extra(option_type) -> Bool

Return `true` if the option type ignores extra fields when
read from a dict-like object.

!!! note
    Normally, we require the dict-like object to have exactly the same
    number of fields with the option type. However, it could be useful
    to have ignore extra fields when wrapping network work services to
    ignore some irrelavent optional fields.

!!! note
    Unlike [pydantic](https://pydantic-docs.helpmanual.io/usage/model_config/),
    we do not allow dynamically adding fields to a given type. One should manually
    define fields you would like to include in a struct type and let it to have
    type `Dict{String, Any}`.

# Example

```julia
julia> Configurations.ignore_extra(::Type{MyOption}) = true
```
"""
function ignore_extra(::Type{OptionType}) where {OptionType}
    is_option(OptionType) || error("expect an option type")
    return false
end

function assert_field_match_exactly(::Type{OptionType}, d::AbstractDict{String}) where {OptionType}
    ignore_extra(OptionType) && return

    nf = fieldcount(OptionType)
    option_keys = [string(fieldname(OptionType, idx)) for idx in 1:nf]
    for key in keys(d)
        key == "#metadata#" && continue
        key in option_keys || throw(InvalidKeyError(key, option_keys))
    end
    return
end

"""
    OptionField{name}
    OptionField(name::Symbol)

Trait type that denotes a field of an option type. Mainly used
for dispatch purpose. `name` should be a `Symbol`.
"""
struct OptionField{name} end

OptionField(name::Symbol) = OptionField{name}()

"""
    from_dict(::Type{OptionType}, ::OptionField{f_name}, ::Type{T}, x) where {OptionType, f_name, T}

For option type `OptionType`, convert the object `x` to the field type `T` and assign it to the field
`f_name`.
Raise `FieldTypeConversionError`s errors if `Base.convert` raises exception
```
ERROR: MethodError: Cannot `convert` an object of type ...
```
"""
function from_dict(::Type{OptionType}, optionfield::OptionField{f_name}, ::Type{T}, x
      ) where {OptionType,f_name,T}
    try
        return from_dict(OptionType, T, x)
    catch err
        if err isa MethodError && err.f === convert
            throw(FieldTypeConversionError(typeof(x), f_name, T, OptionType))
        else
            throw(err)
        end
    end
end

function from_dict(
    ::Type{OptionType}, of::OptionField, ::Type{T}, x
) where {OptionType,T<:AbstractVector}
    if eltype(T) isa Union
        return map(x) do each
            from_dict_union_type(OptionType, of, eltype(T), each)
        end
    else
        return map(x) do each
            from_dict(OptionType, of, eltype(T), each)
        end
    end
end

# default conversions
"""
    from_dict(::Type{OptionType}, ::Type{T}, x) where {OptionType, T}

For option type `OptionType`, convert the object `x` to type `T`. This is
similar to `Base.convert(::Type{T}, x)` and will fallback to `Base.convert`
if not defined.
"""
function from_dict(::Type{OptionType}, ::Type{T}, x) where {OptionType,T}
    is_option(T) && return from_dict(T, x)
    # TODO: deprecate convert_to_option
    # then just use the following
    # return convert(T, x)
    return deprecated_conversion(OptionType, T, x)
    # return convert(T, x)
end

from_dict(::Type, ::Type{VersionNumber}, x) = VersionNumber(x)

function deprecated_conversion(::Type{OptionType}, ::Type{T}, x) where {OptionType,T}
    ret = convert_to_option(OptionType, T, x)
    if ret === ConvertNotFound()
        return convert(T, x)
    else
        return ret
    end
end

@generated function from_dict_union_type(
    ::Type{OptionType}, ::OptionField{f_name}, ::Type{FieldType}, value
) where {OptionType,f_name,FieldType}
    types = Base.uniontypes(FieldType)
    return from_dict_union_type_generated(OptionType, OptionField(f_name), types, :value)
end

function from_dict_union_type_dynamic(
    ::Type{OptionType}, of::OptionField{f_name}, ::Type{FieldType}, value
) where {OptionType,f_name,FieldType}
    FieldType isa Union || return from_dict(OptionType, of, FieldType, value)
    assert_duplicated_alias_union(FieldType)

    types = Base.uniontypes(FieldType)
    if Nothing in types
        value === nothing && return nothing # happy path
        types = filter(x -> x !== Nothing, types)
        if length(types) == 1
            return from_dict(OptionType, of, types[1], value)
        end
    end

    for T in types
        if is_option(T)
            alias = type_alias(T)
            reflect_idx = find_reflect_field(T)
            # type can be determined by alias
            if haskey(value, alias)
                return from_dict(T, value[alias])
            elseif reflect_idx !== nothing # type can be determined by reflect field
                reflect_f_name = fieldname(T, reflect_idx)
                reflect_key = string(reflect_f_name)
                if haskey(value, reflect_key)
                    if alias == value[reflect_key] # find alias in reflect
                        return from_dict(T, value)
                    else
                        type = tryparse_jltype(value[reflect_key], get_type_alias_map(OptionType))
                        type === nothing && continue
                        # NOTE: type is always more specialized
                        type <: T && return from_dict(type, value)
                    end
                end
                continue
            end

            if alias === nothing && reflect_idx === nothing
                # not much information, try parse
                try
                    return from_dict(T, value)
                catch
                    continue
                end
            end
        else
            try
                return from_dict(OptionType, of, T, value)
            catch
                continue
            end
        end
    end
    return error("cannot parse field $f_name, expect $FieldType, got $(typeof(value))")
end

function find_reflect_field(::Type{OptionType}) where {OptionType}
    nf = fieldcount(OptionType)
    for f_idx in 1:nf
        f_type = fieldtype(OptionType, f_idx)
        if f_type <: Reflect
            return f_idx
        end
    end
    return nothing
end

function assert_duplicated_alias_union(::Type{UnionType}) where {UnionType}
    set = String[]
    for T in Base.uniontypes(UnionType)
        if is_option(T)
            alias = type_alias(T)
            if alias !== nothing
                alias in set && throw(DuplicatedAliasError(alias))
                push!(set, alias)
            end
        end
    end
    return nothing
end

"""
    from_dict_specialize(::Type{OptionType}, x) where {OptionType}

A specialized [`from_dict`](@ref) method for option type `OptionType`.
This is usually generated by [`@option`](@ref), but one may also specialized
this manually to acheive maximal performance.
"""
function from_dict_specialize(::Type{OptionType}, x) where {OptionType}
    error(
        "specialized from_dict is not generated for `$OptionType` " *
        "consider define your own `from_dict` " *
        "otherwise use the `@option` macro to create the option type."
    )
end

"""
    from_dict_generated(::Type{OptionType}, value::Symbol) where {OptionType}

Generate a specialized Julia expression to convert an `AbstractDict{String}` to
our `OptionType`.
"""
function from_dict_generated(::Type{OptionType}, value::Symbol) where {OptionType}
    nf = fieldcount(OptionType)
    ret = Expr(:block)
    construct = Expr(:call, OptionType)
    for f_idx in 1:nf
        f_name = fieldname(OptionType, f_idx)
        f_type = fieldtype(OptionType, f_idx)
        f_default = field_default(OptionType, f_name)

        var = f_name
        key = string(f_name)
        field_value = gensym(:field_value)
        push!(construct.args, var)

        jl = JLIfElse()

        if f_default === no_default
            err_msg = "expect key: $key"
            jl[:(!haskey($value, $key))] = :(error($err_msg))
        elseif f_default isa PartialDefault
            jl[:(!haskey($value, $key))] = :($var = $(f_default.lambda)($(f_default.vars...)))
        else
            jl[:(!haskey($value, $key))] = :($var = $f_default)
        end

        body = from_dict_generated(OptionType, OptionField(f_name), f_type, field_value)
        # Maybe{option} type wants to treat empty dict
        # as all default value
        types = Base.uniontypes(f_type)
        types = filter(x -> x !== Nothing, types)

        if f_default === nothing && length(types) == 1 && is_option(types[1])
            jl.otherwise = quote
                $field_value = $value[$key]
                $var = $body
            end
        elseif f_default === nothing
            jl.otherwise = quote
                $field_value = $value[$key]
                if $field_value isa AbstractDict && isempty($field_value)
                    $var = nothing
                else
                    $var = $body
                end
            end
        else
            jl.otherwise = quote
                $field_value = $value[$key]
                $var = $body
            end
        end

        push!(ret.args, codegen_ast(jl))
    end
    push!(ret.args, construct)
    return ret
end

function from_dict_generated(
    option_type, of::OptionField, f_type::Type, field_value::Symbol
)
    if is_option(f_type)
        quote
            $field_value isa AbstractDict || $field_value isa $f_type ||
                error("expect an AbstractDict or $($f_type), got $(typeof($field_value))")
            # NOTE: we want to allow user overloaded from_dict here
            # thus we don't use $(from_dict_generated(f_type, value))
            $Configurations.from_dict($f_type, $field_value)
        end
    elseif f_type isa Union
        types = Base.uniontypes(f_type)
        from_dict_union_type_generated(option_type, of, types, field_value)
    elseif f_type <: Reflect # check if the reflect value match current type
        msg = "type mismatch, expect $option_type got"
        alias = type_alias(option_type)
        alias_map = get_type_alias_map(option_type)
        if alias === nothing
            quote
                $Configurations.parse_jltype($field_value, $alias_map) <: $option_type ||
                    throw(ArgumentError($msg * " $($field_value)"))
                Reflect()
            end
        else
            quote
                $field_value == $alias ||
                    $Configurations.parse_jltype($field_value, $alias_map) <: $option_type ||
                    throw(ArgumentError($msg * " $($field_value)"))
                Reflect()
            end
        end
    else
        quote
            $Configurations.from_dict($option_type, $of, $f_type, $field_value)
        end
    end
end

function from_dict_union_type_generated(
    option_type, of::OptionField, types::Vector{Any}, value::Symbol
)
    if Nothing in types
        from_dict_maybe_type_generated(option_type, of, types, value)
    elseif has_same_reflect_field(types)
        _from_dict_union_type_similar_reflect_field(types, value)
    else # fallback to dynamic
        FieldType = Union{types...}
        return quote
            $Configurations.from_dict_union_type_dynamic(
                $option_type, $of, $FieldType, $value
            )
        end
    end
end

function from_dict_maybe_type_generated(
    option_type, of::OptionField, types::Vector{Any}, value::Symbol
)
    types = filter(x -> x !== Nothing, types)
    if length(types) == 1 # Maybe{T}
        return quote
            if $value === nothing
                nothing
            else
                $(from_dict_generated(option_type, of, types[1], value))
            end
        end
    else # Maybe{Union{A, B, C...}}
        return quote
            if $value === nothing
                nothing
            else
                $(from_dict_union_type_generated(option_type, of, types, value))
            end
        end
    end
end

function _from_dict_union_type_similar_reflect_field(types::Vector{Any}, value::Symbol)
    T = first(types)
    f_alias_map = alias_map(types)
    f_alias_map = isempty(f_alias_map) ? nothing : f_alias_map
    idx = find_reflect_field(T)
    reflect_key = string(fieldname(T, idx))
    msg = "expect key: $reflect_key"

    type_err_msg = "expect one of the following option type: " *
        join(map(x->string("`", x, "`"), types), ", ")
    @gensym type
    return quote
        haskey($value, $reflect_key) || error($msg)
        $type = $Configurations.parse_jltype($value[$reflect_key], $f_alias_map)
        $type <: $(Union{types...}) || error($type_err_msg * " got `$($type)`")
        $Configurations.from_dict_specialize($type, $value)
    end
end

"""
    has_same_reflect_field(types::Vector{Any})

Check if all types has the same reflect field name.
"""
function has_same_reflect_field(types::Vector{Any})
    idx = find_reflect_field(first(types))
    idx === nothing && return false

    f_reflect = fieldname(first(types), idx)
    for t_idx in 2:length(types)
        T = types[t_idx]
        t_reflect_idx = find_reflect_field(T)
        t_reflect_idx === nothing && return false

        t_f_reflect = fieldname(T, t_reflect_idx)
        if f_reflect !== t_f_reflect
            return false
        end
    end
    return true
end

"""
    alias_map(types::Vector{Any})

Create a `Dict` mapping type alias to a type.
"""
function alias_map(types::Vector{Any})
    d = Dict{String,Any}()
    for t in types
        if is_option(t)
            alias_map = get_type_alias_map(t)::Dict{String, Any}
            merge!(d, alias_map)
        end
    end
    return d
end


text:
"""
    CoulombReactionField(; dist_cutoff, solvent_dielectric, nl_only, weight_14,
                            coulomb_const, force_unit, energy_unit)

The Coulomb electrostatic interaction modified using the reaction field approximation.
"""
struct CoulombReactionField{D, S, W, T, F, E} <: GeneralInteraction
    dist_cutoff::D
    solvent_dielectric::S
    nl_only::Bool
    weight_14::W
    coulomb_const::T
    force_unit::F
    energy_unit::E
end

const solventdielectric = 78.3

function CoulombReactionField(;
                    dist_cutoff,
                    solvent_dielectric=solventdielectric,
                    nl_only=false,
                    weight_14=1.0,
                    coulomb_const=coulombconst,
                    force_unit=u"kJ * mol^-1 * nm^-1",
                    energy_unit=u"kJ * mol^-1")
    return CoulombReactionField{typeof(dist_cutoff), typeof(solvent_dielectric), typeof(weight_14),
                                typeof(coulomb_const), typeof(force_unit), typeof(energy_unit)}(
        dist_cutoff, solvent_dielectric, nl_only, weight_14,
        coulomb_const, force_unit, energy_unit)
end

@inline @inbounds function force(inter::CoulombReactionField,
                                    coord_i,
                                    coord_j,
                                    atom_i,
                                    atom_j,
                                    box_size,
                                    weight_14::Bool=false)
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)

    if r2 > (inter.dist_cutoff ^ 2)
        return ustrip.(zero(coord_i)) * inter.force_unit
    end

    coulomb_const = inter.coulomb_const
    qi, qj = atom_i.charge, atom_j.charge
    r = √r2
    i, j = atom_i.index, atom_j.index
    if weight_14
        # 1-4 interactions do not use the reaction field approximation
        krf = (1 / (inter.dist_cutoff ^ 3)) * 0
    else
        # These values could be pre-computed but this way is easier for AD
        krf = (1 / (inter.dist_cutoff ^ 3)) * ((inter.solvent_dielectric - 1) / (2 * inter.solvent_dielectric + 1))
    end

    f = (coulomb_const * qi * qj) * (inv(r) - 2 * krf * r2) * inv(r2)

    if weight_14
        return f * dr * inter.weight_14
    else
        return f * dr
    end
end

@inline @inbounds function potential_energy(inter::CoulombReactionField,
                                    s::System,
                                    i::Integer,
                                    j::Integer,
                                    weight_14::Bool=false)
    dr = vector(s.coords[i], s.coords[j], s.box_size)
    r2 = sum(abs2, dr)

    if r2 > (inter.dist_cutoff ^ 2)
        return ustrip(zero(s.box_size[1])) * inter.energy_unit
    end

    coulomb_const = inter.coulomb_const
    qi, qj = s.atoms[i].charge, s.atoms[j].charge
    r = √r2
    if weight_14
        # 1-4 interactions do not use the reaction field approximation
        krf = (1 / (inter.dist_cutoff ^ 3)) * 0
        crf = (1 /  inter.dist_cutoff     ) * 0
    else
        krf = (1 / (inter.dist_cutoff ^ 3)) * ((inter.solvent_dielectric - 1) / (2 * inter.solvent_dielectric + 1))
        crf = (1 /  inter.dist_cutoff     ) * ((3 * inter.solvent_dielectric) / (2 * inter.solvent_dielectric + 1))
    end

    pe = (coulomb_const * qi * qj) * (inv(r) + krf * r2 - crf)

    if weight_14
        return pe * inter.weight_14
    else
        return pe
    end
end