! file src_equations/equation_module_template.f90
!
! Copyright 2009-2018 Dalton Harvie (daltonh@unimelb.edu.au)
! 
! This file is part of arb finite volume solver, referred to as `arb'.
! 
! arb is a software package designed to solve arbitrary partial
! differential equations on unstructured meshes using the finite volume
! method.  Primarily it consists of fortran source code, perl source
! code and shell scripts.  arb replies on certain third party software
! to run, most notably the computer algebra system maxima
! <http://maxima.sourceforge.net/> which is released under the GNU GPL.
! 
! The original copyright of arb is held by Dalton Harvie, however the
! project is now under collaborative development.
! 
! arb is released under the GNU GPL.  arb is free software: you can
! redistribute it and/or modify it under the terms of the GNU General
! Public License (version 3) as published by the Free Software Foundation.
! You should have received a copy of the GNU General Public Licence
! along with arb (see file licence/gpl.txt after unpacking).  If not,
! see <http://www.gnu.org/licences/>.
! 
! arb is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
! FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public Licence
! for more details.
! 
! For full details of arb's licence see the licence directory.
! 
! The current homepage for the arb finite volume solver project is
! <http://people.eng.unimelb.edu.au/daltonh/downloads/arb>.
!
!-------------------------------------------------------------------------
!
! this file has strings replaced by the script setup_equations.pl to form
! the fortran code equation_module.f90
!
! Simulation info generated by setup_equations when equation_module.f90 was created
!<sub_string:simulation_info>

module equation_module

! drag in only these data types from general_module
use general_module, only : vector_list_type, scalar_list_type, integer_list_type

implicit none
private
public allocate_meta_arrays, update_derived_and_equations, update_constants, update_unknowns, update_newtients, &
  update_initial_newtients, update_transients, update_initial_transients, update_outputs, &
  check_variable_constraints, check_condition, find_condition, var_value, varcdivgrad, varcgrad, varcgrad_nodelimited, &
  setup_external_functions, read_initial_outputs, update_region

! include external functions preambles here
!<arb_external_preamble>

contains

!----------------------------------------------------------------------------

subroutine allocate_meta_arrays

! here we allocate array elements of the meta data that corresponds to the input files

use general_module
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine allocate_meta_arrays'

! allocate meta arrays for variables, compounds and region_links
!<sub_string:allocate_meta_arrays>

! allocate meta arrays for regions
!<sub_string:allocate_regions>

! allocate var_lists which are lists of variables and regions used for fast looping
!<sub_string:allocate_var_lists>

! set these commonly used variable-only var_list_numbers
var_list_number_unknown = var_list_number(centring="all",type="unknown",include_regions=.false.)
var_list_number_equation = var_list_number(centring="all",type="equation",include_regions=.false.)
var_list_number_condition = var_list_number(centring="all",type="condition",include_regions=.false.)
! and these variable and region ones
var_list_number_all_region = var_list_number(centring="all",type="all",include_regions=.true.)
var_list_number_derived_region = var_list_number(centring="all",type="derived",include_regions=.true.)
var_list_number_equation_region = var_list_number(centring="all",type="equation",include_regions=.true.)
var_list_number_constant_region = var_list_number(centring="all",type="constant",include_regions=.true.)
var_list_number_unknown_region = var_list_number(centring="all",type="unknown",include_regions=.true.)
var_list_number_newtient_region = var_list_number(centring="all",type="newtient",include_regions=.true.)
var_list_number_transient_region = var_list_number(centring="all",type="transient",include_regions=.true.)
var_list_number_output_region = var_list_number(centring="all",type="output",include_regions=.true.)

! setup msomeloop which is used to allocate the funks within the someloop (thread) containers
!<sub_string:set_msomeloop>

! setup mseparation_list which is used to allocate the separation_lists within the someloop (thread) containers
!<sub_string:set_mseparation_list>

! set transient_simulation logical
!<sub_string:transient_simulation>

! set newtient_simulation logical
!<sub_string:newtient_simulation>

! set kernel availability
!<sub_string:kernel_availability>

if (debug) write(*,'(a/80(1h-))') 'subroutine allocate_meta_arrays'

end subroutine allocate_meta_arrays

!----------------------------------------------------------------------------

subroutine setup_external_functions

! hook subroutine that will include setup for any external functions

use general_module
logical, parameter :: debug = .false.
integer :: n
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine setup_external_functions'

! include external functions setup here
!<arb_external_setup>

if (debug) write(*,'(a/80(1h-))') 'subroutine setup_external_functions'

end subroutine setup_external_functions

!-----------------------------------------------------------------

subroutine update_derived_and_equations(setup)

! here we update all the deriveds and equations

use general_module
use gmesh_module
!$ use omp_lib
logical :: setup ! now needs to know whether we are in the setup phase or not - if so, posibly read in data
integer :: nvar, m, ns, i, j, k
logical :: region_l
integer :: thread = 1
double precision :: derivative_multiplier
character(len=1000) :: formatline
character(len=1000) :: error_string
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine update_derived_and_equations'

if (debug) write(*,'(40(1h+)/a)') 'derived updates'

do nvar = 1, allocatable_size(var_list(var_list_number_derived_region)%list)
  m = var_list(var_list_number_derived_region)%list(nvar)
  region_l = var_list(var_list_number_derived_region)%region(nvar)

  if (region_l) then

    if (debug) write(*,*) 'updating dynamic derived region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
      region(m)%centring

    if (output_region_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! derived region
!    <sub_string:derived_region>

    if (output_region_update_times) call time_variable_update(thread,1,m,region_l=region_l)

  else

    i = 0
    j = 0
    k = 0

    if (debug) write(*,*) 'updating var derived: m = ',m,': name = ',trim(var(m)%name),': centring = ',var(m)%centring
    error_string = "Error occurred while updating derived "//trim(var(m)%name)

    if (output_variable_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! derived
!    <sub_string:derived>

    if (output_variable_update_times) call time_variable_update(thread,1,m,region_l=region_l)

  ! read from data file only (possibly) during variable setup
    if (setup) call read_gmesh(contents='data',var_number=m)

    if (debug) then
      formatline = '(a,'//trim(indexformat)//',a,i3,a,a)'
      do ns = 1, ubound(var(m)%funk,1)
        write(*,fmt=formatline) '(ns,m) = (',ns,',',m,'): var = ',trim(print_funk(var(m)%funk(ns)))
      end do
    end if

  end if

end do

if (debug) write(*,'(a/40(1h-))') 'derived updates'

!-----------------------------------------------------

if (debug) write(*,'(40(1h+)/a)') 'equation updates'

do nvar = 1, allocatable_size(var_list(var_list_number_equation_region)%list)
  m = var_list(var_list_number_equation_region)%list(nvar)
  region_l = var_list(var_list_number_equation_region)%region(nvar)

  if (region_l) then

    if (debug) write(*,*) 'updating dynamic equation region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
      region(m)%centring

    if (output_region_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! equation region
!    <sub_string:equation_region>

    if (output_region_update_times) call time_variable_update(thread,1,m,region_l=region_l)

  else

    i = 0
    j = 0
    k = 0

    if (debug) write(*,*) 'updating var equation: m = ',m,': name = ',trim(var(m)%name),': centring = ',var(m)%centring
    error_string = "Error occurred while updating equation "//trim(var(m)%name)

    if (output_variable_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! equations
!    <sub_string:equation>

    if (output_variable_update_times) call time_variable_update(thread,1,m,region_l=region_l)

  ! read from data file only (possibly) during variable setup
    if (setup) call read_gmesh(contents='data',var_number=m)

    if (debug) then
      formatline = '(a,'//trim(indexformat)//',a,i3,a,a)'
      do ns = 1, ubound(var(m)%funk,1)
        write(*,fmt=formatline) '(ns,m) = (',ns,',',m,'): var = ',trim(print_funk(var(m)%funk(ns)))
      end do
    end if

  end if

end do

if (debug) write(*,'(a/40(1h-))') 'equation updates'

if (debug) write(*,'(a/80(1h-))') 'subroutine update_derived_and_equations'

end subroutine update_derived_and_equations

!-----------------------------------------------------------------

recursive subroutine update_someloop(thread,m,ilast,jlast,klast,error_string)

! here we update each of all the possible someloops
! bug wife again!!!

use general_module
integer :: m, ilast, jlast, klast, i, j, k, ns, region_number, flux_direction, to_ns, from_ns, thread, static_ns, ii, &
  nseparation, i2, ii2, iistart, iiend
double precision :: derivative_multiplier, reflect_multiplier
character(len=1000) :: error_string
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine update_someloop'
!if (.not.present(error_string)) call error_stop("something is calling update_someloop without an error_string")

if (debug) then
  write(*,*) 'm (someloop) = ',m
  write(*,*) 'thread = ',thread
  write(*,*) 'calling error_string = '//trim(error_string)
end if

! ilast and jlast these hold the indices from the calling routine, and cannot be changed since they are used in the
! loop of the calling routine - so save a copy of these that can be altered within this routine
i = ilast 
j = jlast
k = klast

if (debug) then
  if (.not.allocated(someloop)) call error_stop("someloop not allocated")
  if (thread > ubound(someloop,1) .or. thread < lbound(someloop,1)) call error_stop("someloop not in range")
  if (.not.allocated(someloop(thread)%funk)) call error_stop("someloop%funk not allocated")
  if (m > ubound(someloop(thread)%funk,1) .or. m < lbound(someloop(thread)%funk,1)) call error_stop("someloop%funk not in range")
  someloop(thread)%funk(m)%v = 0.d0
  someloop(thread)%funk(m)%ndv = 0
else
  call reset_funk(someloop(thread)%funk(m))
end if

!<sub_string:someloop>

if (debug) write(*,'(a/80(1h-))') 'subroutine update_someloop'

end subroutine update_someloop

!-----------------------------------------------------------------

subroutine update_constants

! here we update all the constants that are calculated using equations

use general_module
use gmesh_module
!$ use omp_lib
integer :: nvar, m, ns, i, j, k
logical :: region_l
integer :: thread = 1
double precision :: derivative_multiplier
character(len=1000) :: formatline
character(len=1000) :: error_string
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine update_constants'

do nvar = 1, allocatable_size(var_list(var_list_number_constant_region)%list)
  m = var_list(var_list_number_constant_region)%list(nvar)
  region_l = var_list(var_list_number_constant_region)%region(nvar)

  if (region_l) then

    if (debug) write(*,*) 'updating dynamic constant region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
      region(m)%centring

    if (output_region_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! constant region
!    <sub_string:constant_region>

    if (output_region_update_times) call time_variable_update(thread,1,m,region_l=region_l)

  else

    i = 0
    j = 0
    k = 0

    if (debug) write(*,*) 'updating var constant: m = ',m,': name = ',trim(var(m)%name),': centring = ',var(m)%centring
    error_string = "Error occurred while updating constant "//trim(var(m)%name)

    if (output_variable_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! constants
!    <sub_string:constant>

    if (output_variable_update_times) call time_variable_update(thread,1,m,region_l=region_l)

! read from data file
    call read_gmesh(contents='data',var_number=m)

    if (debug) then
      formatline = '(a,'//trim(indexformat)//',a,i3,a,a)'
      do ns = 1, ubound(var(m)%funk,1)
        write(*,fmt=formatline) '(ns,m) = (',ns,',',m,'): var = ',trim(print_funk(var(m)%funk(ns)))
      end do
    end if

  end if

end do

if (debug) write(*,'(a/80(1h-))') 'subroutine update_constants'

end subroutine update_constants

!-----------------------------------------------------------------

subroutine update_unknowns(initial,lambda)

! here we update all the unknown initial values

use general_module
use gmesh_module
!$ use omp_lib
double precision, optional :: lambda ! if initial=.false. then lambda is the backstepping parameter and needs to be specified
logical :: initial ! whether this is the first initialisation of the unknowns (true) or the backstepping update of them (false)
integer :: nvar, m, ns, i, j, k, p
logical :: region_l
integer :: thread = 1
double precision :: derivative_multiplier
character(len=1000) :: formatline
character(len=1000) :: error_string
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine update_unknowns'

do nvar = 1, allocatable_size(var_list(var_list_number_unknown_region)%list)
  m = var_list(var_list_number_unknown_region)%list(nvar)
  region_l = var_list(var_list_number_unknown_region)%region(nvar)

  if (region_l) then

    if (debug) write(*,*) 'updating dynamic unknown region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
      region(m)%centring

    if (output_region_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! unknown region
!    <sub_string:unknown_region>

    if (output_region_update_times) call time_variable_update(thread,1,m,region_l=region_l)

  else

    if (initial) then

      i = 0
      j = 0
      k = 0

      if (debug) write(*,*) 'updating var unknown: m = ',m,': name = ',trim(var(m)%name),': centring = ',var(m)%centring
      error_string = "Error occurred while updating unknown "//trim(var(m)%name)

      if (output_variable_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! unknowns, updated via user-entered initial expressions
!      <sub_string:unknown>

      if (output_variable_update_times) call time_variable_update(thread,1,m,region_l=region_l)

! read from data file
      call read_gmesh(contents='data',var_number=m)

    else

! backstepping update
      do ns = 1, ubound(var(m)%funk,1)
        p = var(m)%funk(ns)%pp(1)
        var(m)%funk(ns)%v = phiold(p) + lambda*delphi(p)
      end do

    end if 

    if (debug) then
      formatline = '(a,'//trim(indexformat)//',a,i3,a,a)'
      do ns = 1, ubound(var(m)%funk,1)
        write(*,fmt=formatline) '(ns,m) = (',ns,',',m,'): var = ',trim(print_funk(var(m)%funk(ns)))
      end do
    end if

  end if

end do

if (debug) write(*,'(a/80(1h-))') 'subroutine update_unknowns'

end subroutine update_unknowns

!-----------------------------------------------------------------

subroutine update_newtients

! here we update all the newtients for a new timestep

use general_module
use gmesh_module
!$ use omp_lib
integer :: nvar, m, ns, i, j, k, relstep
logical :: region_l
integer :: thread = 1
double precision :: derivative_multiplier
character(len=1000) :: formatline
character(len=1000) :: error_string
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine update_newtients'

do relstep = newtient_relstepmax, 0, -1 ! look through newtient variables in reverse order
  do nvar = 1, allocatable_size(var_list(var_list_number_newtient_region)%list)
    m = var_list(var_list_number_newtient_region)%list(nvar)
    region_l = var_list(var_list_number_newtient_region)%region(nvar)

    if (region_l) then

      if (debug) then
        if (region(m)%relstep /= relstep) then
          write(*,*) 'skipping dynamic newtient region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
            trim(region(m)%centring),': relstep = ',region(m)%relstep
          cycle ! skip this region if it isn't the correct relstep
        else
          write(*,*) 'updating dynamic newtient region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
            trim(region(m)%centring),': relstep = ',region(m)%relstep
        end if
      else
        if (region(m)%relstep /= relstep) cycle ! skip this region if it isn't the correct relstep
      end if

      if (output_region_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! newtient region
!      <sub_string:newtient_region>

      if (output_region_update_times) call time_variable_update(thread,1,m,region_l=region_l)

    else

      if (var(m)%relstep /= relstep) cycle ! skip this variable if it isn't the correct relstep

      i = 0
      j = 0
      k = 0
      error_string = "Error occurred while updating newtient "//trim(var(m)%name)

      if (debug) write(*,*) 'updating var newtient: m = ',m,': name = ',trim(var(m)%name),': centring = ',var(m)%centring, &
        ": relstep = ",relstep

      if (output_variable_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! newtients
!      <sub_string:newtient>

      if (output_variable_update_times) call time_variable_update(thread,1,m,region_l=region_l)

      if (debug) then
        formatline = '(a,'//trim(indexformat)//',a,i3,a,a)'
        do ns = 1, ubound(var(m)%funk,1)
          write(*,fmt=formatline) '(ns,m) = (',ns,',',m,'): var = ',trim(print_funk(var(m)%funk(ns)))
        end do
      end if

    end if

  end do
end do

if (debug) write(*,'(a/80(1h-))') 'subroutine update_newtients'

end subroutine update_newtients

!-----------------------------------------------------------------

subroutine update_initial_newtients

! here we initialised all the newtients

use general_module
use gmesh_module
!$ use omp_lib
integer :: nvar, m, ns, i, j, k, relstep
logical :: region_l
integer :: thread = 1
double precision :: derivative_multiplier
character(len=1000) :: formatline
character(len=1000) :: error_string
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine update_initial_newtients'

do relstep = 0, newtient_relstepmax ! look through newtient variables in forward order
  do nvar = 1, allocatable_size(var_list(var_list_number_newtient_region)%list)
    m = var_list(var_list_number_newtient_region)%list(nvar)
    region_l = var_list(var_list_number_newtient_region)%region(nvar)

    if (region_l) then

      if (debug) then
        if (region(m)%relstep /= relstep) then
          write(*,*) 'skipping dynamic initial newtient region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
            trim(region(m)%centring),': relstep = ',region(m)%relstep
          cycle ! skip this region if it isn't the correct relstep
        else
          write(*,*) 'updating dynamic initial newtient region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
            trim(region(m)%centring),': relstep = ',region(m)%relstep
        end if
      else
        if (region(m)%relstep /= relstep) cycle ! skip this region if it isn't the correct relstep
      end if

      if (output_region_update_times.and..not.ignore_initial_update_times) call time_variable_update(thread,0,m,region_l=region_l)
 
! newtient region
!      <sub_string:initial_newtient_region>

      if (output_region_update_times.and..not.ignore_initial_update_times) call time_variable_update(thread,1,m,region_l=region_l)

    else

      if (var(m)%relstep /= relstep) cycle ! skip this variable if it isn't the correct relstep

      i = 0
      j = 0
      k = 0

      if (debug) write(*,*) 'initialising var newtient: m = ',m,': name = ',trim(var(m)%name),': centring = ',var(m)%centring,": relstep = ",relstep
      error_string = "Error occurred while updating initial_newtient "//trim(var(m)%name)

      if (output_variable_update_times.and..not.ignore_initial_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! initial_newtients
!      <sub_string:initial_newtient>

      if (output_variable_update_times.and..not.ignore_initial_update_times) call time_variable_update(thread,1,m,region_l=region_l)

! read from data file
      call read_gmesh(contents='data',var_number=m)

      if (debug) then
        formatline = '(a,'//trim(indexformat)//',a,i3,a,a)'
        do ns = 1, ubound(var(m)%funk,1)
          write(*,fmt=formatline) '(ns,m) = (',ns,',',m,'): var = ',trim(print_funk(var(m)%funk(ns)))
        end do
      end if

    end if

  end do
end do

if (debug) write(*,'(a/80(1h-))') 'subroutine update_initial_newtients'

end subroutine update_initial_newtients

!-----------------------------------------------------------------

subroutine update_transients

! here we update all the transients for a new timestep

use general_module
use gmesh_module
!$ use omp_lib
integer :: nvar, m, ns, i, j, k, relstep
logical :: region_l
integer :: thread = 1
double precision :: derivative_multiplier
character(len=1000) :: formatline
character(len=1000) :: error_string
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine update_transients'

do relstep = transient_relstepmax, 0, -1 ! look through transient variables in reverse order
  do nvar = 1, allocatable_size(var_list(var_list_number_transient_region)%list)
    m = var_list(var_list_number_transient_region)%list(nvar)
    region_l = var_list(var_list_number_transient_region)%region(nvar)

    if (region_l) then

      if (debug) then
        if (region(m)%relstep /= relstep) then
          write(*,*) 'skipping dynamic transient region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
            trim(region(m)%centring),': relstep = ',region(m)%relstep
          cycle ! skip this region if it isn't the correct relstep
        else
          write(*,*) 'updating dynamic transient region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
            trim(region(m)%centring),': relstep = ',region(m)%relstep
        end if
      else
        if (region(m)%relstep /= relstep) cycle ! skip this region if it isn't the correct relstep
      end if

      if (output_region_update_times) call time_variable_update(thread,0,m,region_l=region_l)
 
! transient region
!      <sub_string:transient_region>

      if (output_region_update_times) call time_variable_update(thread,1,m,region_l=region_l)
 
    else

      if (var(m)%relstep /= relstep) cycle ! skip this variable if it isn't the correct relstep

      i = 0
      j = 0
      k = 0

      if (debug) write(*,*) 'updating var transient: m = ',m,': name = ',trim(var(m)%name),': centring = ',var(m)%centring, &
        ": relstep = ",relstep
      error_string = "Error occurred while updating transient "//trim(var(m)%name)

      if (output_variable_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! transients
!      <sub_string:transient>

      if (output_variable_update_times) call time_variable_update(thread,1,m,region_l=region_l)

      if (debug) then
        formatline = '(a,'//trim(indexformat)//',a,i3,a,a)'
        do ns = 1, ubound(var(m)%funk,1)
          write(*,fmt=formatline) '(ns,m) = (',ns,',',m,'): var = ',trim(print_funk(var(m)%funk(ns)))
        end do
      end if

    end if

  end do
end do

if (debug) write(*,'(a/80(1h-))') 'subroutine update_transients'

end subroutine update_transients

!-----------------------------------------------------------------

subroutine update_initial_transients

! here we initialised all the transients

use general_module
use gmesh_module
!$ use omp_lib
integer :: nvar, m, ns, i, j, k, relstep
logical :: region_l
integer :: thread = 1
double precision :: derivative_multiplier
character(len=1000) :: formatline
character(len=1000) :: error_string
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine update_initial_transients'

if (debug) write(*,*) 'var_list_number_transient_region = ',var_list_number_transient_region
do relstep = 0, transient_relstepmax ! look through transient variables in forward order
  do nvar = 1, allocatable_size(var_list(var_list_number_transient_region)%list)
    m = var_list(var_list_number_transient_region)%list(nvar)
    region_l = var_list(var_list_number_transient_region)%region(nvar)

    if (debug) write(*,*) 'relstep = ',relstep,': nvar = ',nvar,': m = ',m,': region_l = ',region_l
    if (region_l) then

      if (debug) then
        if (region(m)%relstep /= relstep) then
          write(*,*) 'skipping dynamic initial transient region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
            trim(region(m)%centring),': relstep = ',region(m)%relstep
          cycle ! skip this region if it isn't the correct relstep
        else
          write(*,*) 'updating dynamic initial transient region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
            trim(region(m)%centring),': relstep = ',region(m)%relstep
        end if
      else
        if (region(m)%relstep /= relstep) cycle ! skip this region if it isn't the correct relstep
      end if

      if (output_region_update_times.and..not.ignore_initial_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! transient region
!      <sub_string:initial_transient_region>

      if (output_region_update_times.and..not.ignore_initial_update_times) call time_variable_update(thread,1,m,region_l=region_l)

    else

      if (var(m)%relstep /= relstep) cycle ! skip this variable if it isn't the correct relstep

      i = 0
      j = 0
      k = 0

      if (debug) write(*,*) 'initialising var transient: m = ',m,': name = ',trim(var(m)%name),': centring = ',var(m)%centring,": relstep = ",relstep
      error_string = "Error occurred while updating initial_transient "//trim(var(m)%name)

      if (output_variable_update_times.and..not.ignore_initial_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! initial_transients
!      <sub_string:initial_transient>

      if (output_variable_update_times.and..not.ignore_initial_update_times) call time_variable_update(thread,1,m,region_l=region_l)

  ! read from data file
      call read_gmesh(contents='data',var_number=m)

      if (debug) then
        formatline = '(a,'//trim(indexformat)//',a,i3,a,a)'
        do ns = 1, ubound(var(m)%funk,1)
          write(*,fmt=formatline) '(ns,m) = (',ns,',',m,'): var = ',trim(print_funk(var(m)%funk(ns)))
        end do
      end if

    end if

  end do
end do

if (debug) write(*,'(a/80(1h-))') 'subroutine update_initial_transients'

end subroutine update_initial_transients

!-----------------------------------------------------------------

subroutine update_outputs(stepoutput)

! here we update all the outputs that are calculated using equations

use general_module
!$ use omp_lib
logical, optional, intent(in) :: stepoutput
logical :: stepoutput_local
integer :: nvar, m, ns, i, j, k
logical :: region_l
integer :: thread = 1
double precision :: derivative_multiplier
character(len=1000) :: formatline
character(len=1000) :: error_string
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine update_outputs'

stepoutput_local = .false.
if (present(stepoutput)) then
  stepoutput_local = stepoutput
end if

do nvar = 1, allocatable_size(var_list(var_list_number_output_region)%list)
  m = var_list(var_list_number_output_region)%list(nvar)
  region_l = var_list(var_list_number_output_region)%region(nvar)

  if (region_l) then

    if (debug) write(*,*) 'updating dynamic output region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
      region(m)%centring

! TODO: use options in regions to allow dynamic regions to be updated using stepoutput
! right now output dynamic regions are never updated during setupoutput
    if (stepoutput_local) cycle

    if (output_region_update_times) call time_variable_update(thread,0,m,region_l=region_l)
 
! output region
!    <sub_string:output_region>

    if (output_region_update_times) call time_variable_update(thread,1,m,region_l=region_l)

  else

    i = 0
    j = 0
    k = 0

    if (debug) write(*,*) 'updating var output: m = ',m,': name = ',trim(var(m)%name),': centring = ',var(m)%centring
    error_string = "Error occurred while updating output "//trim(var(m)%name)

  ! if stepoutput is not on for either the component or compound and this is being called from output_step, then skip
  !  updating this variable
    if (stepoutput_local) then
      if (.not.( trim(check_option(compound(var(m)%compound_number)%options,stepoutput_options)) == "stepoutput" .or. &
                 trim(check_option(var(m)%options,stepoutput_options)) == "stepoutput" )) cycle
    end if

    if (output_variable_update_times) call time_variable_update(thread,0,m,region_l=region_l)

! equations
!    <sub_string:output>

    if (output_variable_update_times) call time_variable_update(thread,1,m,region_l=region_l)

    if (debug) then
      formatline = '(a,'//trim(indexformat)//',a,i3,a,a)'
      do ns = 1, ubound(var(m)%funk,1)
        write(*,fmt=formatline) '(ns,m) = (',ns,',',m,'): var = ',trim(print_funk(var(m)%funk(ns)))
      end do
    end if

  end if

end do

if (debug) write(*,'(a/80(1h-))') 'subroutine update_outputs'

end subroutine update_outputs

!-----------------------------------------------------------------

subroutine read_initial_outputs

! here we read any output data from previous files

use general_module
use gmesh_module
integer :: nvar, m, ns
logical :: region_l
character(len=1000) :: formatline
character(len=1000) :: error_string
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine read_initial_outputs'

do nvar = 1, allocatable_size(var_list(var_list_number_output_region)%list)
  m = var_list(var_list_number_output_region)%list(nvar)
  region_l = var_list(var_list_number_output_region)%region(nvar)

  if (region_l) then

    if (debug) write(*,*) 'updating dynamic output region: m = ',m,': name = ',trim(region(m)%name),': centring = ', &
      region(m)%centring

! output region
!    <sub_string:output_region>

  else

    if (debug) write(*,*) 'updating var output: m = ',m,': name = ',trim(var(m)%name),': centring = ',var(m)%centring
    error_string = "Error occurred while updating output "//trim(var(m)%name)

  ! read from data file
    call read_gmesh(contents='data',var_number=m)

    if (debug) then
      formatline = '(a,'//trim(indexformat)//',a,i3,a,a)'
      do ns = 1, ubound(var(m)%funk,1)
        write(*,fmt=formatline) '(ns,m) = (',ns,',',m,'): var = ',trim(print_funk(var(m)%funk(ns)))
      end do
    end if

  end if

end do

if (debug) write(*,'(a/80(1h-))') 'subroutine read_initial_outputs'

end subroutine read_initial_outputs

!-----------------------------------------------------------------

function find_condition(condition_type)

! here we check whether a condition is satisfied for action = output|stop|error|convergence|bell|(time|newt)steprewind
! do this in series as if in parallel all threads would be attempting to alter the one logical
! could be done but not worth it (conditions are usually none centred anyway)
! returns fortran number of condition that is true

use general_module
character(len=*) :: condition_type
integer :: find_condition
integer :: nvar, m, ns, i, j, k
integer :: thread = 1
double precision :: derivative_multiplier
character(len=1000) :: formatline
character(len=1000) :: error_string
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'function find_condition'

find_condition = 0 ! indicating that the condition has not been satisfied

if (var_list_number_condition < 0) var_list_number_condition = var_list_number(centring="all",type="condition",include_regions=.false.)
do nvar = 1, allocatable_size(var_list(var_list_number_condition)%list)
  m = var_list(var_list_number_condition)%list(nvar)

  i = 0
  j = 0
  k = 0

  if (debug) write(*,*) 'checking condition var: m = ',m,': name = ',trim(var(m)%name),': centring = ',var(m)%centring
  error_string = "Error occurred while updating condition "//trim(var(m)%name)

  if (output_variable_update_times) call time_variable_update(thread,0,m,region_l=.false.)

! conditions
!  <sub_string:condition>

  if (output_variable_update_times) call time_variable_update(thread,1,m,region_l=.false.)

  if (debug) then
    formatline = '(a,'//trim(indexformat)//',a,i3,a,a)'
    do ns = 1, ubound(var(m)%funk,1)
      write(*,fmt=formatline) '(ns,m) = (',ns,',',m,'): var = ',trim(print_funk(var(m)%funk(ns)))
    end do
  end if

end do

if (debug) write(*,'(a/80(1h-))') 'function find_condition'

end function find_condition

!-----------------------------------------------------------------

function check_condition(condition_type)

! wrapper function that returns a check on a condition as a logical
use general_module
character(len=*) :: condition_type
logical :: check_condition

if (find_condition(condition_type) == 0) then
  check_condition = .false.
else
  check_condition = .true.
end if

end function check_condition

!-----------------------------------------------------------------

subroutine check_variable_constraints(type,ierror)

! here we check on the magnitudes of any variables against specified options

use general_module
character(len=*) :: type
integer :: ns, ierror
character(len=1000) :: formatline
integer :: thread = 1
logical, parameter :: debug = .false.
                  
if (debug) write(*,'(80(1h+)/a)') 'subroutine check_variable_constraints'

ierror = 0

if (trim(type) == 'unknown') then
  if (debug) write(*,*) 'subroutine check_variable_constraints called with type = '//trim(type)
!  <sub_string:unknown_constraints>
else if (trim(type) == 'derived') then
  if (debug) write(*,*) 'subroutine check_variable_constraints called with type = '//trim(type)
!  <sub_string:derived_constraints>
else if (trim(type) == 'equation') then
  if (debug) write(*,*) 'subroutine check_variable_constraints called with type = '//trim(type)
!  <sub_string:equation_constraints>
else if (trim(type) == 'local') then
  if (debug) write(*,*) 'subroutine check_variable_constraints called with type = '//trim(type)
!  <sub_string:local_constraints>
else
  call error_stop("subroutine check_variable_constraints called with variable type = "//trim(type))
end if

if (debug) write(*,'(a/80(1h-))') 'subroutine check_variable_constraints'

end subroutine check_variable_constraints

!-----------------------------------------------------------------

function var_value(m,ns,noerror)

! little function to find value of var, irrespective of whether it is local (ie not stored) or not
! if noerror is omitted, then the default is not to issue an error message but instead return a 0.d0 result

use general_module
integer, intent(in) :: m ! var number
integer, intent(in) :: ns ! ns reference for region of var
logical, optional, intent(in) :: noerror ! if ns is passed as zero just output result as zero rather than flagging an error
logical :: noerror_l
integer :: i, j, k
integer :: thread = 1
double precision :: var_value
character(len=1000) :: error_string

if (ns == 0) then

  noerror_l = .true.
  if (present(noerror)) noerror_l = noerror
  if (noerror_l) then
    var_value = 0.d0
  else
    call error_stop('ns is equal to 0 in var_value:  Generally this means that an attempt is being made '// &
      'to reference '//trim(var(m)%centring)//' centred variable '//trim(var(m)%name)//' outside of the region '// &
      trim(var(m)%region)//' in which it is defined. '// &
      'Look at the use of this variable in the equations and check that its region context is correct.')
  end if
    
else if (var(m)%type == 'local') then

  i = 0
  j = 0
  k = 0
  if (var(m)%centring == 'cell') then
    i = region(var(m)%region_number)%ijk(ns)
    if (i == 0) call error_stop('i is equal to 0 in var_value: problem with use of '//trim(var(m)%name))
  else if (var(m)%centring == 'face') then
    j = region(var(m)%region_number)%ijk(ns)
    if (j == 0) call error_stop('j is equal to 0 in var_value: problem with use of '//trim(var(m)%name))
  else if (var(m)%centring == 'node') then
    k = region(var(m)%region_number)%ijk(ns)
    if (k == 0) call error_stop('k is equal to 0 in var_value: problem with use of '//trim(var(m)%name))
  end if
  error_string = "Error occurred while updating someloop "//trim(var(m)%name)//" from subroutine var_value.  It is likely that "// &
    "this is being called from an output routine while outputting a local variable."

  call update_someloop(thread,var(m)%someloop,i,j,k,error_string)
  var_value = someloop(thread)%funk(var(m)%someloop)%v

else

  var_value = var(m)%funk(ns)%v

end if

end function var_value

!-----------------------------------------------------------------

function varcdivgrad(i,m)

! same as cellgrad but references a var function (which must be cell centred)
! calculate gradient (vector) of cell-centred data var(m) in cell i

use general_module
integer, intent(in) :: i, m
integer :: ii, j, jj, ik, ns
double precision, dimension(totaldimensions) :: varcdivgrad
logical, parameter :: debug = .false.

if (debug) write(*,'(80(1h+)/a)') 'function varcdivgrad'

varcdivgrad = [ 0.d0, 0.d0, 0.d0 ]

do jj = 1, ubound(cell(i)%jface,1)
  j = cell(i)%jface(jj)
  do ii = 1, ubound(face(j)%kernel(0)%ijk,1)
    ik = face(j)%kernel(0)%ijk(ii) 
    if (ik == 0) cycle
    if (region(var(m)%update_region_number)%ns(ik) == 0) ik = i ! if point doesn't lie in (update_)region then use centre
    ns = region(var(m)%region_number)%ns(ik) ! find real ns which is based on parent region
    varcdivgrad = varcdivgrad + face(j)%norm(:,1)*divop(i,j)*face(j)%kernel(0)%v(ii)*var_value(m,ns)
  end do
end do

if (debug) write(*,'(a/80(1h-))') 'function varcdivgrad'

end function varcdivgrad

!-----------------------------------------------------------------

function varcgrad(i,m)

! same as cellgrad but references a var function (which must be cell centred)
! calculate gradient (vector) of cell-centred data var(m) in cell i

use general_module
integer, intent(in) :: i, m
integer :: l, ii, ik, ns, nsc
double precision, dimension(totaldimensions) :: varcgrad
logical, parameter :: debug = .false.

if (debug) write(*,'(80(1h+)/a)') 'function varcgrad'

varcgrad = [ 0.d0, 0.d0, 0.d0 ]
nsc = region(var(m)%region_number)%ns(i)
if (nsc == 0) call error_stop("trying to evaluate a cell centred variable in varcgrad in a cell in which it is not defined")

do l = 1, totaldimensions
  do ii = 1, allocatable_integer_size(cell(i)%kernel(l)%ijk)
    if (debug) then
      write(90,*) 'i = ',i,': l = ',l,': ii = ',ii,': m = ',m,': centring = ',trim(var(m)%centring)
      flush(90)
    end if
    ik = cell(i)%kernel(l)%ijk(ii) 
    if (ik == 0) cycle
    if (region(var(m)%update_region_number)%ns(ik) == 0) then
      ns = nsc ! if point doesn't lie in region then use centre
    else
      ns = region(var(m)%region_number)%ns(ik)
    end if
    varcgrad(l) = varcgrad(l) + cell(i)%kernel(l)%v(ii)*var_value(m,ns)
  end do
end do

if (debug) write(*,'(a/80(1h-))') 'function varcgrad'

end function varcgrad

!-----------------------------------------------------------------

function varcgrad_nodelimited(i,m)

! same as varcgrad but limits the gradient so that all node values lie within the range of 
!  node-neighbouring surrounding cell values
! calculate limited gradient (vector) of cell-centred data var(m) in cell i

use general_module
integer, intent(in) :: i, m
integer :: kk, k, ii2, i2, ns
double precision, dimension(totaldimensions) :: varcgrad_nodelimited
double precision :: varmax, varmin, varc, varl, alpha
logical, parameter :: debug = .false.

if (debug) write(*,'(80(1h+)/a)') 'function varcgrad_nodelimited'

varcgrad_nodelimited = varcgrad(i,m)

! find var in the central cell and set default limits
ns = region(var(m)%region_number)%ns(i)
if (ns == 0) call error_stop("trying to evaluate a cell centred variable in varcgrad in a cell in which it is not defined")
varc = var_value(m,ns)
varmax = 0.d0
varmin = 0.d0

! cycle through each node of the cell, calculating the (relative) maximum and minimum node values
do kk = 1, allocatable_integer_size(cell(i)%knode)
  k = cell(i)%knode(kk)
  do ii2 = 1, allocatable_integer_size(node(k)%icell)
    i2 = node(k)%icell(ii2)
    if (i2 == i) cycle ! central cell already included
    if (region(var(m)%update_region_number)%ns(i2) == 0) cycle
    ns = region(var(m)%region_number)%ns(i2)
    varl = var_value(m,ns) - varc
    varmax = max(varmax,varl)
    varmin = min(varmin,varl)
  end do
end do
    
! cycle through the nodes again, calculating the non-limited var values at each
alpha = 1.d0 ! this is the multiplier that will eventually limit the gradient
do kk = 1, allocatable_integer_size(cell(i)%knode)
  k = cell(i)%knode(kk)
  varl = dot_product(node(k)%x-cell(i)%x,varcgrad_nodelimited)
  if (varl > varmax) then
    alpha = min(alpha,varmax/max(varl,tinyish))
  else if (varl < varmin) then
    alpha = min(alpha,-varmin/max(-varl,tinyish))
  end if
end do

!alpha = max(alpha,0.d0) 
if (alpha < 1.d0) varcgrad_nodelimited = alpha*varcgrad_nodelimited 

if (debug) write(*,'(a/80(1h-))') 'function varcgrad_nodelimited'

end function varcgrad_nodelimited

!-----------------------------------------------------------------

subroutine update_region(m,initial)

! here we setup the regions by finding the i or j indices for each

use general_module
integer :: m ! region number to be updated
logical :: initial ! whether the the initial or normal location string is to be used
type(region_location_type) :: local_location ! set to either initial or normal location
integer :: i, j, k, n, nregion, ijkregion, nsregion, ns, ii, jj, kk, ijk, l, ijktotal, ns2, donor_region, part_of_region, &
  nsnext, iinext, inext, separation, iimax, glue_ijk
double precision :: tmp, tmpmax
double precision, dimension(totaldimensions) :: x, xmin, xmax, unitnormal ! a single location
character(len=1000) :: formatline
logical :: compoundtype, compoundadd, setijk, setns, faceseparation
logical, dimension(:), allocatable :: elementisin
logical :: debug_sparse = .false.
logical, parameter :: debug = .false.

if (debug) debug_sparse = .true.
                  
if (debug) write(82,'(80(1h+)/a)') 'subroutine update_region'

if (debug) write(82,*) 'Processing region m = ',m,': name = '//trim(region(m)%name)// &
  ': centring = '//trim(region(m)%centring)//': initial = ',initial,': part_of = ', &
  region(m)%part_of,': parent = ',region(m)%parent

! ns array is allocated always, but is zeroed here (possibly again)
region(m)%ns = 0

! specify which of ijk and ns are already (or about to be) set, and from which the other needs to be (possibly) calculated
! for gmsh and system variables, ijk is read in/set, so make this the default
setijk = .true.
setns = .false. ! also need to calculate nsregion, which is the number of elements in the region

!-----------------------------------------------------------------
if (trim(region(m)%type) == 'system') then

  deallocate(region(m)%ijk) ! undo temporary allocation from setup

  if (trim(region(m)%name) == '<allcells>') then
    allocate(region(m)%ijk(itotal))
    do n = 1, itotal
      region(m)%ijk(n) = n
    end do
  else if (trim(region(m)%name) == '<domain>') then
    allocate(region(m)%ijk(idomain))
    n = 0
    do i = 1, itotal
      if (cell(i)%type == 1) then
        n = n + 1
        region(m)%ijk(n) = i
      end if
    end do
  else if (trim(region(m)%name) == '<boundarycells>') then
    allocate(region(m)%ijk(iboundary))
    n = 0
    do i = 1, itotal
      if (cell(i)%type == 2) then
        n = n + 1
        region(m)%ijk(n) = i
      end if
    end do
  else if (trim(region(m)%name) == '<allfaces>') then
    allocate(region(m)%ijk(jtotal))
    do n = 1, jtotal
      region(m)%ijk(n) = n
    end do
  else if (trim(region(m)%name) == '<domainfaces>') then
    allocate(region(m)%ijk(jdomain))
    n = 0
    do j = 1, jtotal
      if (face(j)%type == 1) then
        n = n + 1
        region(m)%ijk(n) = j
      end if
    end do
  else if (trim(region(m)%name) == '<boundaries>') then
    allocate(region(m)%ijk(jboundary))
    n = 0
    do j = 1, jtotal
      if (face(j)%type == 2) then
        n = n + 1
        region(m)%ijk(n) = j
      end if
    end do
  else if (trim(region(m)%name) == '<allnodes>') then
    allocate(region(m)%ijk(ktotal))
    do n = 1, ktotal
      region(m)%ijk(n) = n
    end do
  else if (trim(region(m)%name) == '<domainnodes>') then
    allocate(region(m)%ijk(kdomain))
    n = 0
    do k = 1, ktotal
      if (node(k)%type == 1) then
        n = n + 1
        region(m)%ijk(n) = k
      end if
    end do
  else if (trim(region(m)%name) == '<boundarynodes>') then
    allocate(region(m)%ijk(kboundary))
    n = 0
    do k = 1, ktotal
      if (node(k)%type == 2) then
        n = n + 1
        region(m)%ijk(n) = k
      end if
    end do
  else
    call error_stop("subroutine update_system_region called with incorrect region name "//trim(region(m)%name))
  end if

!-----------------------------------------------------------------
else if (trim(region(m)%type) /= 'gmsh') then
! now deal with all user regions which have locations and/or initial locations to be processed

! note: this may mean that the region is already defined but defining twice won't hurt if the definition is the same

  if (initial) then
    local_location = region(m)%initial_location
    if (.not.local_location%active) call error_stop("region "//trim(region(m)%name)// &
      " is trying to be updated but its initial location it isn't active")
  else
    local_location = region(m)%location
    if (.not.local_location%active) call error_stop("region "//trim(region(m)%name)// &
      " is trying to be updated but its location it isn't active")
  end if

  if (debug) then
    write(82,*) "now processing local_location"
    write(82,*) "type = "//trim(local_location%type)
    write(82,*) "description = "//trim(local_location%description)
    write(82,*) "glue_face = ",local_location%glue_face
    if (allocated(local_location%floats)) write(82,*) "floats = ",local_location%floats
    if (allocated(local_location%integers)) write(82,*) "integers = ",local_location%integers
    if (allocated(local_location%regions)) write(82,*) "regions = ",local_location%regions
  end if

! TODO:
! expand region

! deallocate any allocated regions, starting afresh each time the region is calculated
! will have to think about this more for separation regions etc
! if (allocated(region(m)%ijk)) then
  if (allocatable_integer_size(region(m)%ijk) > 0) then
    if (.not.region(m)%dynamic) write(*,'(a)') "NOTE: a "//trim(local_location%type)// &
      " region operator is acting on region "//trim(region(m)%name)// &
      " that already contains elements: the previous element will be overwritten with the new"
  end if

! think about sizing of ijk for dynamic elements - maybe allocate once based on parent size? - and then deal with zero ijk elements when used?

! set ijktotal as it is used for most of these location types
  if (region(m)%centring == 'cell') then
    ijktotal = itotal
  else if (region(m)%centring == 'face') then
    ijktotal = jtotal
  else
    ijktotal = ktotal
  end if

!---------------------
! ref: at region
! a user defined region from the arb input file that is a single point
! setting ijk but not ns

  if (trim(local_location%type) == "at") then

    call resize_integer_array(array=region(m)%ijk,new_size=1,keep_data=.false.,default_value=0)
    region(m)%ijk(1) = 0
    x = local_location%floats ! the perl has ensured that this has exactly 3 elements

    tmpmax = 1.d+20
    do ns=1,allocatable_integer_size(region(region(m)%part_of)%ijk)
      ijk = region(region(m)%part_of)%ijk(ns)
      if (ijk == 0) cycle ! allow for zero elements in ijk
      if (region(m)%centring == "cell") then
        tmp = distance(x,cell(ijk)%x)
      else if (region(m)%centring == "face") then
        tmp = distance(x,face(ijk)%x)
      else
        tmp = distance(x,node(ijk)%x)
      end if
      if (tmp < tmpmax) then
        region(m)%ijk(1) = ijk
        tmpmax = tmp
      end if
    end do

    if (region(m)%ijk(1) == 0) call resize_integer_array(array=region(m)%ijk,new_size=0,keep_data=.false.)

!---------------------
! ref: withinbox region ref: within region
! a user defined region from the arb input file that is any elements within a box
! TODO: deal with other geometries
! setting ns but not ijk

  else if (trim(local_location%type) == "withinbox") then

    setijk = .false.
    setns = .true.

! check that points are in min and max order and otherwise reorder

    xmin = local_location%floats(1:3) ! the perl has ensured that this has exactly 3 elements
    xmax = local_location%floats(4:6) ! the perl has ensured that this has exactly 3 elements
    do l = 1, 3
      if (xmin(l) > xmax(l)) then
        write(*,'(a,i1,a)') 'WARNING: dimension ',l,' of the points that define the BOX geometry in region '// &
          trim(region(m)%name)//' were incorrectly ordered: they should be in the order of the minimum coordinate values '// &
          'in each dimension, followed by the maximum coordinate values in each dimension.'
        write(*,'(2(a,i1,a,g14.6))') ' xmin(',l,') = ',xmin(l),' xmax(',l,') = ',xmax(l)
        write(*,*) 'These values will be swapped.'
        tmp = xmax(l)
        xmax(l) = xmin(l)
        xmin(l) = tmp
      end if
    end do

    nsregion = 0
    outer: do ns = 1, allocatable_integer_size(region(region(m)%part_of)%ijk) ! just set the elements within the part_of region to be true
      ijk = region(region(m)%part_of)%ijk(ns)
      if (ijk == 0) cycle ! allow for zero elements in ijk
      if (region(m)%centring == "cell") then
        x = cell(ijk)%x
      else if (region(m)%centring == "face") then
        x = face(ijk)%x
      else
        x = node(ijk)%x
      end if
      inner: do l = 1, 3
        if ((x(l)-xmin(l))*(xmax(l)-x(l)) < 0.d0) cycle outer
      end do inner
      nsregion = nsregion + 1
      region(m)%ns(ijk) = nsregion
    end do outer

!---------------------
! ref: normal region, only face centred
! a user defined region from the arb input file that is any faces that have a normal that is within a certain dot-product of the specified unit normal
! the first three components define the normal, and will be normalised
! the fourth component is the maximum deviation from the normal that is allowed, expressed as a dot-product (deviation, ie, valid range is between 0 (exact) and 2 (all faces))
! setting ns but not ijk

  else if (trim(local_location%type) == "normal") then

    setijk = .false.
    setns = .true.

! this should already be checked from the perl, but...
    if (region(m)%centring /= "face") call error_stop("normal region "//trim(region(m)%name)//" must be face centred")

! create unit normal from three components
    tmp = vector_magnitude(local_location%floats(1:3))
    if (tmp < tinyish) call error_stop("the normal for normal region "//trim(region(m)%name)// &
      " is numerically too small:  have you entered the three components as the first three floats within the region description?")
    unitnormal = local_location%floats(1:3)/tmp

    nsregion = 0
    do ns = 1, allocatable_integer_size(region(region(m)%part_of)%ijk) ! just set the elements within the part_of region to be true
      ijk = region(region(m)%part_of)%ijk(ns)
      if (ijk == 0) cycle ! allow for zero elements in ijk
      if (1.d0 - dot_product(face(ijk)%norm(:,1),unitnormal) > local_location%floats(4)) cycle
      nsregion = nsregion + 1
      region(m)%ns(ijk) = nsregion
    end do

!---------------------
! ref: associatedwith ref: boundaryof ref: domainof ref: surrounds
! a new region composed of the boundary to another region, or similar type of related domain
! setting ns but not ijk

  else if (trim(local_location%type) == "boundaryof".or.trim(local_location%type) == "domainof".or.trim(local_location%type) == "associatedwith".or. &
    trim(local_location%type) == "surrounds") then

    setijk = .false.
    setns = .true.

! check centring of requested related region
    nregion = local_location%regions(1) ! perl only allows one region to be defined
    if (region(nregion)%centring == 'none') &
      call error_stop('incorrect centring for '//trim(local_location%type)//' consitutent region '//trim(region(nregion)%name))

    nsregion = 0
    do ns2 = 1,allocatable_integer_size(region(nregion)%ijk)  ! loop through all ijk indices in constituent region
      ijkregion = region(nregion)%ijk(ns2) ! NB, ijkregion and region(m)%ijk may actually be i or j values depending on centring
!     if (ijkregion == 0) cycle ! in case this is empty due to a dynamic region - no, presently all ijk values have to be non-zero
      if (region(nregion)%centring == 'cell') then
        if (region(m)%centring == 'cell') then
! create cell region from cell region
          do ii = 1, ubound(cell(ijkregion)%jface,1)+1 ! loop around cells that border cells, and itself
            i = cell(ijkregion)%icell(ii)
            if (region(region(m)%part_of)%ns(i) == 0) cycle ! check that element is a member of the part_of region
            if (trim(local_location%type) == "boundaryof".and.cell(i)%type /= 2) cycle
            if (trim(local_location%type) == "domainof".and.cell(i)%type /= 1) cycle
            if (trim(local_location%type) == "surrounds".and.region(nregion)%ns(i) /= 0) cycle
            if (region(m)%ns(i) == 0) then
              nsregion = nsregion + 1
              region(m)%ns(i) = nsregion
            end if
          end do
        else if (region(m)%centring == 'face') then
! create face region from cell region
          do jj = 1, ubound(cell(ijkregion)%jface,1)
            j = cell(ijkregion)%jface(jj)
            if (region(region(m)%part_of)%ns(j) == 0) cycle ! check that element is a member of the part_of region
            if (trim(local_location%type) == "boundaryof".and.face(j)%type /= 2) cycle
            if (trim(local_location%type) == "domainof".and.face(j)%type /= 1) cycle
            if (trim(local_location%type) == "surrounds") then
              if (ijkregion == face(j)%icell(2).and.region(nregion)%ns(face(j)%icell(1)) /= 0) cycle
              if (ijkregion == face(j)%icell(1).and.region(nregion)%ns(face(j)%icell(2)) /= 0) cycle
            end if
            if (region(m)%ns(j) == 0) then
              nsregion = nsregion + 1
              region(m)%ns(j) = nsregion
            end if
          end do
        else
! create node region from cell region
! surrounds not implemented
          do kk = 1, ubound(cell(ijkregion)%knode,1)
            k = cell(ijkregion)%knode(kk)
            if (region(region(m)%part_of)%ns(k) == 0) cycle ! check that element is a member of the part_of region
            if (trim(local_location%type) == "boundaryof".and.node(k)%type /= 2) cycle
            if (trim(local_location%type) == "domainof".and.node(k)%type /= 1) cycle
            if (trim(local_location%type) == "surrounds") &
              call error_stop("surrounds not implemented for constructing node from cell region "//trim(region(nregion)%name))
            if (region(m)%ns(k) == 0) then
              nsregion = nsregion + 1
              region(m)%ns(k) = nsregion
            end if
          end do
        end if
      else if (region(nregion)%centring == 'face') then
! create cell region from face region
! boundaryof will pick out boundary cells coincident with boundary faces
! surrounds not implemented
        if (region(m)%centring == 'cell') then
          do ii = 1, 2
            i = face(ijkregion)%icell(ii)
            if (region(region(m)%part_of)%ns(i) == 0) cycle ! check that element is a member of the part_of region
            if (trim(local_location%type) == "boundaryof".and.cell(i)%type /= 2) cycle
            if (trim(local_location%type) == "domainof".and.cell(i)%type /= 1) cycle
            if (trim(local_location%type) == "surrounds") &
              call error_stop("surrounds not implemented for constructing node from cell region "//trim(region(nregion)%name))
            if (region(m)%ns(i) == 0) then
              nsregion = nsregion + 1
              region(m)%ns(i) = nsregion
            end if
          end do
        else if (region(m)%centring == 'face') then
! create face region from face region
! boundaryof will pick out faces that are on the boundary
! associatedwith is nonsense - will just copy region
          j = ijkregion
          if (region(region(m)%part_of)%ns(j) == 0) cycle ! check that element is a member of the part_of region
          if (trim(local_location%type) == "boundaryof".and.face(j)%type /= 2) cycle
          if (trim(local_location%type) == "domainof".and.face(j)%type /= 1) cycle
          if (region(m)%ns(j) == 0) then
            nsregion = nsregion + 1
            region(m)%ns(j) = nsregion
          end if
        else
! create node region from face region
          do kk = 1, ubound(face(ijkregion)%knode,1)
            k = face(ijkregion)%knode(kk)
            if (region(region(m)%part_of)%ns(k) == 0) cycle ! check that element is a member of the part_of region
            if (trim(local_location%type) == "boundaryof".and.node(k)%type /= 2) cycle
            if (trim(local_location%type) == "domainof".and.node(k)%type /= 1) cycle
            if (trim(local_location%type) == "surrounds") &
              call error_stop("surrounds not implemented for constructing node from face region "//trim(region(nregion)%name))
            if (region(m)%ns(k) == 0) then
              nsregion = nsregion + 1
              region(m)%ns(k) = nsregion
            end if
          end do
        end if
      else
! TODO: fix this node region stuff
        call error_stop(trim(local_location%type)//" not implemented for constructing node region "//trim(region(nregion)%name))
      end if
    end do

!---------------------
! ref: compound ref: common ref: intersection ref: union
! a region composed of a compound list of other regions, or a region that has all of the listed regions in common
! we use elementisin here as the ns index of the elements is unknown when they are being added
! ns and ijk will be calculated below from elementisin
! NB, alternative but equivalent names supported by setup_equations: intersection = common and union = compound

  else if (trim(local_location%type) == "compound" .or. trim(local_location%type) == "common") then

    setijk = .false.
    setns = .false.

    allocate(elementisin(ijktotal)) 
    elementisin = .false.
    if (trim(local_location%type) == "compound") then
      compoundtype = .true. ! use this for fast lookup later
    else
      compoundtype = .false. 
      do ns = 1, allocatable_integer_size(region(region(m)%part_of)%ijk) ! just set the elements within the part_of region to be true
        ijk = region(region(m)%part_of)%ijk(ns)
        if (ijk == 0) cycle
        elementisin(ijk) = .true.
      end do
    end if

! loop through all regions identified in the list
    do n = 1, allocatable_integer_size(local_location%regions)
      nregion = local_location%regions(n)
      if (compoundtype) then
        compoundadd = .false.
        if (local_location%integers(n) == 1) compoundadd = .true.
      end if

      if (debug) write(82,*) '++ in compound loop processing region = '//trim(region(nregion)%name)//' with integer = ', &
        local_location%integers(n)

! check centring of requested related region is consistent with the parent
      if (region(nregion)%centring /= region(m)%centring) &
        call error_stop('incorrect centring for '//trim(local_location%type)//' consitutent region '//trim(region(nregion)%name)// &
        ' that is being used in region '//trim(region(m)%name))

! now loop through all elements within part_of region

      do ns = 1, allocatable_integer_size(region(region(m)%part_of)%ijk) ! just set the elements within the part_of region to be true
        ijk = region(region(m)%part_of)%ijk(ns)
        if (debug) write(82,*) 'before: ns = ',ns,': ijk = ',ijk,': elementisin(ijk) = ',elementisin(ijk)
        if (ijk == 0) cycle
        if (compoundtype) then
          if (region(nregion)%ns(ijk) == 0) cycle ! if this element is not in the location region then don't do anything - cycle
          if (compoundadd) then
            if (.not.elementisin(ijk)) elementisin(ijk) = .true. ! adding element
          else
            if (elementisin(ijk)) elementisin(ijk) = .false. ! subtracting element
          end if
        else
          if (elementisin(ijk)) then
            if (region(nregion)%ns(ijk) == 0) elementisin(ijk) = .false.
          end if
        end if
        if (debug) write(82,*) 'after: ns = ',ns,': ijk = ',ijk,': elementisin(ijk) = ',elementisin(ijk)
      end do

    end do

! deal with possible glued faces (and nodes), making sure that they are also in the region if also in the part_of region
    if (local_location%glue_face) then
! sanity check on centring
      if (region(m)%centring /= 'face'.and.region(m)%centring /= 'node') &
        call error_stop("subroutine update_system_region has been called with glueface on, but region "//trim(region(m)%name)// &
        " has incorrect "//trim(region(m)%centring)//" centring (must be either face or node centred to use glueface option)")
      do ijk = 1, ijktotal
        if (.not.elementisin(ijk)) cycle ! face/node isn't in region
! there is only one possible glued face, but an array of possible glued nodes, so deal with separately
        if (region(m)%centring == 'face') then
          glue_ijk = face(ijk)%glue_jface
          if (glue_ijk == 0) cycle ! face/node isn't glued to another
          if (elementisin(glue_ijk)) cycle ! face/node is already in region
          if (region(region(m)%part_of)%ns(glue_ijk) == 0) cycle ! glued face/node isn't in part_of region
! otherwise add face
          nsregion = nsregion + 1
          elementisin(glue_ijk) = .true.
        else
          do kk = 1, allocatable_integer_size(node(ijk)%glue_knode)
            glue_ijk = node(ijk)%glue_knode(kk)
!           if (glue_ijk == 0) cycle ! face/node isn't glued to another - assume that this array doesn't have zero elements
            if (elementisin(glue_ijk)) cycle ! face/node is already in region
            if (region(region(m)%part_of)%ns(glue_ijk) == 0) cycle ! glued face/node isn't in part_of region
! otherwise add node
            nsregion = nsregion + 1
            elementisin(glue_ijk) = .true.
          end do
        end if
      end do
    end if

!---------------------
! ref: variable
! a region which is only true where a variable is greater than zero
! setting ns but not ijk

  else if (trim(local_location%type) == "variable") then

    setijk = .false.
    setns = .true.

    nsregion = 0
    do ns = 1, allocatable_integer_size(region(region(m)%part_of)%ijk) ! just set the elements within the part_of region to be true
      ijk = region(region(m)%part_of)%ijk(ns)
      if (ijk == 0) cycle ! allow for zero elements in ijk
      ns2 = region(var(local_location%variables(1))%region_number)%ns(ijk)
      if (ns2 == 0) cycle
!     if (var(local_location%variables(1))%funk(ns2)%v > 0.d0) then
      if (var_value(local_location%variables(1),ns2,noerror=.false.) > 0.d0) then
        nsregion = nsregion + 1
        region(m)%ns(ijk) = nsregion
      end if
    end do

!---------------------
! ref: expand
! expand a region outwards through a set number of separation levels
! setting both ns and ijk

  else if (trim(local_location%type) == "expand") then

    setijk = .true.
    setns = .true.

    donor_region = local_location%regions(1)
    part_of_region = region(m)%part_of
! need to copy ijk and ns indicies from donor region into new region while respecting part_of size of new region
! first increase the new region ijk array size to the maximum size possible, which is that of its part_of region
    call resize_integer_array(array=region(m)%ijk,new_size=allocatable_integer_size(region(part_of_region)%ns),keep_data=.false.)
    region(m)%ijk = 0
    region(m)%ns = 0
    region(m)%nslast = 0 ! this is the last index of the cells one separation less than the maximum in region(m)
    nsregion = 0
    do ns = 1, allocatable_integer_size(region(donor_region)%ijk)
      ijk = region(donor_region)%ijk(ns)
      if (region(part_of_region)%ns(ijk) /= 0) then
        nsregion = nsregion + 1
        region(m)%ijk(nsregion) = ijk
        region(m)%ns(ijk) = nsregion
        if (ns <= region(donor_region)%nslast) region(m)%nslast = nsregion
      end if
    end do

! at this stage we have all elements from donor region in the possibly resized new array

    if (local_location%integers(1) < 0) then
      faceseparation = .true.
    else
      faceseparation = .false.
    end if

! loop through separation levels
    do separation = 1, abs(local_location%integers(1))
      nsnext = nsregion
! loop through indicies in outer layer of region
      do ns2 = region(m)%nslast+1, nsnext
        i = region(m)%ijk(ns2)
! loop through cells that surround this cell
        if (faceseparation) then
! for a faceseparation loop break after the surround face cells have been done
! a faceseparation loop is indicated by a negative maximumseparation integer
          iimax = ubound(cell(i)%jface,1) + 1
        else
          iimax = ubound(cell(i)%icell,1) ! for cellseparation
        end if
        do iinext = 1, iimax
          inext = cell(i)%icell(iinext) ! we need to check whether this cell is in the part_of_region and already in the new region
          if (region(m)%ns(inext) == 0) then
            if (region(part_of_region)%ns(inext) /= 0) then
              nsregion = nsregion + 1
              region(m)%ijk(nsregion) = inext
              region(m)%ns(inext) = nsregion
            end if
          end if
        end do
      end do
! update outer region indicies
      region(m)%nslast=nsnext
    end do

! finally resize ijk array so that it doesn't include any zero elements
    call resize_integer_array(array=region(m)%ijk,new_size=nsregion,keep_data=.true.)

!---------------------
! ref: all
! a region composed of all elements - actually a straight copy from region(m)%part_of
! setting both ns and ijk

  else if (trim(local_location%type) == "all") then

    setijk = .true.
    setns = .true.

    nsregion = allocatable_integer_size(region(region(m)%part_of)%ijk)
    call resize_integer_array(array=region(m)%ijk,new_size=nsregion,keep_data=.false.)
    if (nsregion > 0) region(m)%ijk = region(region(m)%part_of)%ijk
    region(m)%ns = region(region(m)%part_of)%ns

!---------------------
  else if (trim(local_location%type) /= "none") then ! anything but none is an error

    call error_stop('location type for region '//trim(region(m)%name)//' is not understood: type = '//trim(local_location%type))

  end if

!-----------------------------------------------------------------
end if

if (allocated(elementisin)) then
  setijk = .true.
  setns = .true.
! now loop through the elementisin list creating a new ijk list from it
  nsregion = 0 ! to avoid multiple calls to push_array allocate the array using needed size
  do ijk = 1, ijktotal
    if (elementisin(ijk)) nsregion = nsregion + 1
  end do
  call resize_integer_array(array=region(m)%ijk,new_size=nsregion,keep_data=.false.)
  nsregion = 0
  do ijk = 1, ijktotal
    if (elementisin(ijk)) then
      nsregion = nsregion + 1
      region(m)%ns(ijk) = nsregion
      region(m)%ijk(nsregion) = ijk
    end if
  end do
  deallocate(elementisin)
end if


! if ns isn't already set, find ns indicies which give the data number corresponding to location i, j or k
! alternatively, if ijk isn't already set, set it using reverse lookup
! also make sure that nsregion is set to be the number of elements in the region
if (.not.setns) then
  nsregion = allocatable_size(region(m)%ijk)
  if (setijk) then
    do ns = 1, allocatable_size(region(m)%ijk)
      region(m)%ns(region(m)%ijk(ns)) = ns
    end do
  end if
else if (.not.setijk.and.setns) then

! deal with glued faces only for regions that use setijk off and setns on (compound/common regions are handled separately in their section)
! when glueface is on, check that faces/nodes glued to the faces/nodes within the region are also included, if they are within the part_of region
  if (local_location%glue_face) then
! sanity check on centring
    if (region(m)%centring /= 'face'.and.region(m)%centring /= 'node') &
      call error_stop("subroutine update_system_region has been called with glueface on, but region "//trim(region(m)%name)// &
      " has incorrect "//trim(region(m)%centring)//" centring (must be either face or node centred to use glueface option)")
! now add any missing faces/nodes
    do ijk = 1, ijktotal
      if (region(m)%ns(ijk) == 0) cycle ! face/node isn't in region
! there is only one possible glued face, but an array of possible glued nodes, so deal with separately
      if (region(m)%centring == 'face') then
        glue_ijk = face(ijk)%glue_jface
        if (glue_ijk == 0) cycle ! face/node isn't glued to another
        if (region(m)%ns(glue_ijk) /= 0) cycle ! glued face/node is already included
        if (region(region(m)%part_of)%ns(glue_ijk) == 0) cycle ! glued face/node isn't in part_of region
! otherwise add face
        nsregion = nsregion + 1
        region(m)%ns(glue_ijk) = nsregion
      else
        do kk = 1, allocatable_integer_size(node(ijk)%glue_knode)
          glue_ijk = node(ijk)%glue_knode(kk)
!         if (glue_ijk == 0) cycle ! face/node isn't glued to another - assume that this array doesn't have zero elements
          if (region(m)%ns(glue_ijk) /= 0) cycle ! glued face/node is already included
          if (region(region(m)%part_of)%ns(glue_ijk) == 0) cycle ! glued face/node isn't in part_of region
! otherwise add node
          nsregion = nsregion + 1
          region(m)%ns(glue_ijk) = nsregion
        end do
      end if
    end do
  end if

! if setns is true, then nsregion must hold number of elements within region
  call resize_integer_array(array=region(m)%ijk,new_size=nsregion,keep_data=.false.)
  do ijk = 1, ijktotal
    if (region(m)%ns(ijk) /= 0) region(m)%ijk(region(m)%ns(ijk)) = ijk
  end do
end if

if (debug_sparse) then
  if (region(m)%dynamic) then
    formatline = "(a,"//trim(dindexformat(nsregion))//",a)"
    write(*,fmt=formatline) "INFO: updated dynamic "//trim(region(m)%centring)//" "//trim(region(m)%type)//" region "// &
      trim(region(m)%name)//" which now has ",nsregion," elements"
  end if
end if

! check that each region contains some elements, and that it is allocated (even if zero length)
! ie, ensure that all region%ijk arrays are allocated

if (nsregion == 0) then
  if (.not.region(m)%dynamic) &
    write(*,'(a)') "WARNING: the "//trim(region(m)%centring)//" "//trim(region(m)%type)//" region "//trim(region(m)%name)// &
    " contains no elements (none allocated)"
  if (.not.allocated(region(m)%ijk)) allocate(region(m)%ijk(0)) ! don't think that this is required in the code anymore, but it is valid fortran now
end if

    
if (debug) then
  write(82,*) '# region = '//trim(region(m)%name)//': centring = '//region(m)%centring
  write(82,*) '# '//ijkstring(region(m)%centring)//', ns'
  do ijk = 1, ijktotal
    write(82,*) ijk, region(m)%ns(ijk)
  end do
end if

!---------------------

if (debug) write(82,'(a/80(1h-))') 'subroutine update_region'

end subroutine update_region

!-----------------------------------------------------------------

! include external functions contents here
!<arb_external_contents>

!-----------------------------------------------------------------

end module equation_module

!----------------------------------------------------------------------------
