constructor Function

  function constructor(source_field, direction, extra_derivative, &
                       boundary_locs, boundary_types,             &
                       boundary_operations, coef) result(this)
    !* Author: Chris MacMackin
    !  Date: December 2016
    !
    ! Build a block in a Jacobian matrix, with the form
    ! $$\frac{\partial F}{\partial x_i} + F\Delta_i,$$ where \(F\) is
    ! a scalar field and \(\Delta_i\) is the differentiation operator
    ! in the \(i\)-direction. Additionally, a further differentiation
    ! operator may be added to the right hand side of this matrix
    ! block.  Optional arguments allow for handling of boundary
    ! conditions. See the end of the documentation of the
    ! [[jacobian_block(type)]] type for a description of how boundary
    ! conditions are treated.
    !
    class(scalar_field), intent(in)             :: source_field
      !! A scalar field (\(F\)) making up this block of the Jacobian
    integer, intent(in)                         :: direction
      !! The direction in which field derivatives are taken.
    integer, intent(in), optional               :: extra_derivative
      !! If present, specifies the direction of a differentiation
      !! operator to be added to the right hand side of this matrix
      !! block.
    integer, dimension(:), optional, intent(in) :: boundary_locs
      !! The locations in the raw representation of `rhs` for which
      !! boundary conditions are specified. Defaults to there being
      !! none.
    integer, dimension(:), optional, intent(in) :: boundary_types
      !! Integers specifying the type of boundary condition. The type
      !! of boundary condition corresponding to a given integer is
      !! specified in [[boundary_types_mod]]. Only Dirichlet and
      !! Neumann conditions are supported. Defaults to Dirichlet. The
      !! order in which they are stored must match that of
      !! `boundary_locs`.
    procedure(jacobian_block_bounds),  optional :: boundary_operations
      !! A function specifying the values to place at the boundaries
      !! of the result when using the Jacobian block for
      !! multiplication. By default, all boundaries are set to 0. The
      !! order in which the resulting values are stored should match
      !! that of `boundary_locs`.
    real(r8), optional, intent(in)              :: coef
      !! An optional coefficient by which the the \(\partial
      !! F/\partial x \) term in the operator will be
      !! multipled. Default value is 1.
    type(jacobian_block)                        :: this
      !! A new Jacobian block
    call source_field%guard_temp()
    allocate(this%contents, mold=source_field)
    allocate(this%derivative, mold=source_field)
    this%contents = source_field
    this%direction = direction
    this%derivative = this%contents%d_dx(this%direction)
    if (present(extra_derivative)) this%extra_derivative = extra_derivative
    if (present(coef)) this%coef = coef
    if (present(boundary_locs)) then
      this%boundary_locs = boundary_locs
    else
      allocate(this%boundary_locs(0))
    end if
    if (present(boundary_types)) then
      this%boundary_types = boundary_types
    else
      allocate(this%boundary_types(size(this%boundary_locs)))
      this%boundary_types = dirichlet
    end if
    if (present(boundary_operations)) then
      this%get_boundaries => boundary_operations
    else
      this%get_boundaries => jacobian_block_bounds
    end if
    call source_field%clean_temp()
#ifdef DEBUG
    call logger%debug('jacobian_block','Instantiated a Jacobian block object.')
#endif
  end function constructor