text:
!##############################################################################
!# ****************************************************************************
!# <name> cc2dmedium_callback </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains callback functions for the cc2d-mini problem that are
!# used during the matrix/vector assembly for specifying analytical data.
!# There are three callback functions involved, which may be called depending
!# on the situation. All of them correspond to a specific interface for
!# callback functions, defined in 'intf_xxxx.inc' files.
!#
!# 1.) coeff_Stokes
!#     -> Returns the coefficients for the Laplace matrix. This routine is
!#        only used if the problem to calculate has nonconstant coefficients!
!#        Otherwise the routine is dead.
!#     -> Corresponds to the interface defined in the file
!#        'intf_coefficientMatrixSc.inc'
!#
!# 1.) coeff_Pressure
!#     -> Returns the coefficients for the pressure matrix. This routine is
!#        only used if the problem to calculate has nonconstant coefficients!
!#        Otherwise the routine is dead.
!#     -> Corresponds to the interface defined in the file
!#        'intf_coefficientMatrixSc.inc'
!#
!# 2.) coeff_RHS_x
!#     -> Returns analytical values for the right hand side of the X-velocity
!#        equation.
!#     -> Corresponds to the interface defined in the file
!#        'intf_coefficientVectorSc.inc'
!#
!# 3.) coeff_RHS_y
!#     -> Returns analytical values for the right hand side of the Y-velocity
!#        equation.
!#     -> Corresponds to the interface defined in the file
!#        'intf_coefficientVectorSc.inc'
!#
!# 4.) ffunction_TargetX
!#     -> Returns analytical values for the desired flow field in X-direction.
!#     -> Is used for error analysis during the postprocessing.
!#     -> Corresponds to the interface defined in the file
!#        'intf_coefficientVectorSc.inc'
!#
!# 5.) ffunction_TargetY
!#     -> Returns analytical values for the desired flow field in Y-direction.
!#     -> Is used for error analysis during the postprocessing.
!#     -> Corresponds to the interface defined in the file
!#        'intf_coefficientVectorSc.inc'
!#
!# 6.) ffunction_TargetP
!#     -> Returns analytical values for the desired pressure.
!#     -> Is used for error analysis during the postprocessing.
!#     -> Corresponds to the interface defined in the file
!#        'intf_coefficientVectorSc.inc'
!#
!# 7.) getBoundaryValues
!#     -> Returns analitical values on the boundary of the
!#        problem to solve.
!#     -> Corresponds to the interface defined in the file
!#        'intf_bcassembly.inc'
!#
!# 8.) getBoundaryValuesFBC
!#     -> Returns analytical values on the fictitious boundary components
!#     -> Corresponds to the interface defined in the file
!#        'intf_fbcassembly.inc'
!#
!# 9.) c2d2_initCollectForAssembly
!#     -> Is called prior to the assembly process.
!#     -> Stores some information from the problem structure to a collection
!#        such that it can be accessed in callback routines
!#
!# 10.) c2d2_doneCollectForAssembly
!#      -> Is called after the assembly process.
!#      -> Releases information stored in the collection by
!#         c2d2_initCollectForAssembly.
!#
!# For nonstationary simulation, it might be neccessary in these routines
!# to access the current simulation time. Before the assembly process, the cc2d
!# framework calls c2d2_initCollectForAssembly to stores the current point
!# in time (and probably other necessary information) to the quickaccess-array
!# in the collection which is passed to the callback routines. The callback
!# routines can access this as follows:
!#
!# -> p_rcollection%IquickAccess(1)   = 0: stationary,
!#                                      1: nonstationary with explicit time stepping
!# -> p_rcollection%DquickAccess(1)   = current simulation time
!# -> p_rcollection%DquickAccess(2)   = minimum simulation time
!# -> p_rcollection%DquickAccess(3)   = maximum simulation time
!#
!# After the assembly, c2d2_doneCollectForAssembly is called to clean up.
!# Note: Information stored in the quick-access array are of temporary
!# nature and does not have to be cleaned up.
!#
!# </purpose>
!##############################################################################

module cc2dmedium_callback

  use fsystem
  use storage
  use linearsolver
  use boundary
  use bilinearformevaluation
  use linearformevaluation
  use cubature
  use matrixfilters
  use vectorfilters
  use bcassembly
  use mprimitives

  use cc2dmediumm2basic
  use cc2dmediumm2boundarydef

  implicit none

contains

! ***************************************************************************

!<subroutine>

  subroutine c2d2_initCollectForAssembly (rproblem,rcollection)

!<description>
  ! This subroutine is an auxiliary subroutine called by the CC2D framework
  ! and has usually not to be changed by the user.
  !
  ! The subroutine prepares the collection rcollection to be passed to callback
  ! routines for assembling boundary conditions or RHS vectors. It is
  ! called directly prior to the assembly to store problem-specific information
  ! in the quick-access arrays of the collection.
  ! Basically speaking, the routine stores information about whether thw problem
  ! is stationary, nonstationary or about the current simulation time.
!</description>

!<input>
  ! Problem structure with all problem relevant data.
  type(t_problem), intent(IN) :: rproblem
!</input>

!<inputoutput>
  ! Collection structure to be initialised.
  type(t_collection), intent(INOUT) :: rcollection
!</inputoutput>

!</subroutine>

    ! In a nonstationary simulation, save the simulation time as well as the
    ! minimum and maximum time to the quick-access array of the collection,
    ! so it can be accessed in the callback routines!
    rcollection%Iquickaccess(1) = rproblem%itimedependence
    select case (rproblem%itimedependence)
    case (0)
      ! Stationary simulation
      rcollection%Dquickaccess(1) = 0.0_DP
      rcollection%Dquickaccess(2) = 0.0_DP
      rcollection%Dquickaccess(3) = 0.0_DP
    case (1)
      rcollection%Dquickaccess(1) = rproblem%rtimedependence%dtime
      rcollection%Dquickaccess(2) = rproblem%rtimedependence%dtimeInit
      rcollection%Dquickaccess(3) = rproblem%rtimedependence%dtimeMax
    end select

  end subroutine

! ***************************************************************************

!<subroutine>

  subroutine c2d2_doneCollectForAssembly (rproblem,rcollection)

!<description>
  ! This subroutine is an auxiliary subroutine called by the CC2D framework
  ! and has usually not to be changed by the user.
  !
  ! After the assembly process, this subroutine is called to release temporary
  ! information from the collection which was stored there by
  ! c2d2_initCollectForAssembly.
!</description>

!<input>
  ! Problem structure with all problem relevant data.
  type(t_problem), intent(IN) :: rproblem
!</input>

!<inputoutput>
  ! Collection structure to be cleaned up.
  type(t_collection), intent(INOUT) :: rcollection
!</inputoutput>

!</subroutine>

    ! Currently, this subroutine is empty as all information stored in
    ! the collection in c2d2_initCollectForAssembly is put to the quick-access
    ! arrays -- which do not have to be cleaned up.
    ! This might change in future...

  end subroutine

! ***************************************************************************
  !<subroutine>

  subroutine coeff_Stokes (rdiscretisation,rform, &
                  nelements,npointsPerElement,Dpoints, &
                  IdofsTrial,IdofsTest,rdomainIntSubset, p_rcollection,&
                  Dcoefficients)

    use basicgeometry
    use triangulation
    use collection
    use scalarpde
    use domainintegration
    use dofmapping

  !<description>
    ! This subroutine is called during the matrix assembly. It has to compute
    ! the coefficients in front of the terms of the bilinear form.
    !
    ! The routine accepts a set of elements and a set of points on these
    ! elements (cubature points) in real coordinates.
    ! According to the terms in the bilinear form, the routine has to compute
    ! simultaneously for all these points and all the terms in the bilinear form
    ! the corresponding coefficients in front of the terms.
  !</description>

  !<input>
    ! The discretisation structure that defines the basic shape of the
    ! triangulation with references to the underlying triangulation,
    ! analytic boundary boundary description etc.
    type(t_spatialDiscretisation), intent(IN)                   :: rdiscretisation

    ! The bilinear form which is currently being evaluated:
    type(t_bilinearForm), intent(IN)                            :: rform

    ! Number of elements, where the coefficients must be computed.
    integer, intent(IN)                                         :: nelements

    ! Number of points per element, where the coefficients must be computed
    integer, intent(IN)                                         :: npointsPerElement

    ! This is an array of all points on all the elements where coefficients
    ! are needed.
    ! Remark: This usually coincides with rdomainSubset%p_DcubPtsReal.
    ! DIMENSION(dimension,npointsPerElement,nelements)
    real(DP), dimension(:,:,:), intent(IN)  :: Dpoints

    ! An array accepting the DOF's on all elements trial in the trial space.
    ! DIMENSION(\#local DOF's in trial space,nelements)
    integer(PREC_DOFIDX), dimension(:,:), intent(IN) :: IdofsTrial

    ! An array accepting the DOF's on all elements trial in the trial space.
    ! DIMENSION(\#local DOF's in test space,nelements)
    integer(PREC_DOFIDX), dimension(:,:), intent(IN) :: IdofsTest

    ! This is a t_domainIntSubset structure specifying more detailed information
    ! about the element set that is currently being integrated.
    ! It's usually used in more complex situations (e.g. nonlinear matrices).
    type(t_domainIntSubset), intent(IN)              :: rdomainIntSubset

    ! A pointer to a collection structure to provide additional
    ! information to the coefficient routine. May point to NULL() if not defined.
    type(t_collection), pointer                      :: p_rcollection

  !</input>

  !<output>
    ! A list of all coefficients in front of all terms in the bilinear form -
    ! for all given points on all given elements.
    !   DIMENSION(itermCount,npointsPerElement,nelements)
    ! with itermCount the number of terms in the bilinear form.
    real(DP), dimension(:,:,:), intent(OUT)                      :: Dcoefficients
  !</output>

  !</subroutine>

    ! Routine is not called anyway in case of constant coefficients!
    Dcoefficients = 1.0_DP

  end subroutine

! ***************************************************************************
  !<subroutine>

  subroutine coeff_Pressure (rdiscretisation,rform, &
                  nelements,npointsPerElement,Dpoints, &
                  IdofsTrial,IdofsTest,rdomainIntSubset, p_rcollection,&
                  Dcoefficients)

    use basicgeometry
    use triangulation
    use collection
    use scalarpde
    use domainintegration
    use dofmapping

  !<description>
    ! This subroutine is called during the matrix assembly. It has to compute
    ! the coefficients in front of the terms of the bilinear form.
    !
    ! The routine accepts a set of elements and a set of points on these
    ! elements (cubature points) in real coordinates.
    ! According to the terms in the bilinear form, the routine has to compute
    ! simultaneously for all these points and all the terms in the bilinear form
    ! the corresponding coefficients in front of the terms.
  !</description>

  !<input>
    ! The discretisation structure that defines the basic shape of the
    ! triangulation with references to the underlying triangulation,
    ! analytic boundary boundary description etc.
    type(t_spatialDiscretisation), intent(IN)                   :: rdiscretisation

    ! The bilinear form which is currently being evaluated:
    type(t_bilinearForm), intent(IN)                            :: rform

    ! Number of elements, where the coefficients must be computed.
    integer, intent(IN)                                         :: nelements

    ! Number of points per element, where the coefficients must be computed
    integer, intent(IN)                                         :: npointsPerElement

    ! This is an array of all points on all the elements where coefficients
    ! are needed.
    ! Remark: This usually coincides with rdomainSubset%p_DcubPtsReal.
    ! DIMENSION(dimension,npointsPerElement,nelements)
    real(DP), dimension(:,:,:), intent(IN)  :: Dpoints

    ! An array accepting the DOF's on all elements trial in the trial space.
    ! DIMENSION(\#local DOF's in trial space,nelements)
    integer(PREC_DOFIDX), dimension(:,:), intent(IN) :: IdofsTrial

    ! An array accepting the DOF's on all elements trial in the trial space.
    ! DIMENSION(\#local DOF's in test space,nelements)
    integer(PREC_DOFIDX), dimension(:,:), intent(IN) :: IdofsTest

    ! This is a t_domainIntSubset structure specifying more detailed information
    ! about the element set that is currently being integrated.
    ! It's usually used in more complex situations (e.g. nonlinear matrices).
    type(t_domainIntSubset), intent(IN)              :: rdomainIntSubset

    ! A pointer to a collection structure to provide additional
    ! information to the coefficient routine. May point to NULL() if not defined.
    type(t_collection), pointer                      :: p_rcollection

  !</input>

  !<output>
    ! A list of all coefficients in front of all terms in the bilinear form -
    ! for all given points on all given elements.
    !   DIMENSION(itermCount,npointsPerElement,nelements)
    ! with itermCount the number of terms in the bilinear form.
    real(DP), dimension(:,:,:), intent(OUT)                      :: Dcoefficients
  !</output>

  !</subroutine>

    ! Routine is not called anyway in case of constant coefficients!
    Dcoefficients = 1.0_DP

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine coeff_RHS_x (rdiscretisation,rform, &
                  nelements,npointsPerElement,Dpoints, &
                  IdofsTest,rdomainIntSubset,p_rcollection, &
                  Dcoefficients)

    use basicgeometry
    use triangulation
    use collection
    use scalarpde
    use domainintegration
    use dofmapping

  !<description>
    ! This subroutine is called during the vector assembly. It has to compute
    ! the coefficients in front of the terms of the linear form of the
    ! X-velocity part of the right hand side vector.
    !
    ! The routine accepts a set of elements and a set of points on these
    ! elements (cubature points) in real coordinates.
    ! According to the terms in the linear form, the routine has to compute
    ! simultaneously for all these points and all the terms in the linear form
    ! the corresponding coefficients in front of the terms.
  !</description>

  !<input>
    ! The discretisation structure that defines the basic shape of the
    ! triangulation with references to the underlying triangulation,
    ! analytic boundary boundary description etc.
    type(t_spatialDiscretisation), intent(IN)                   :: rdiscretisation

    ! The linear form which is currently to be evaluated:
    type(t_linearForm), intent(IN)                              :: rform

    ! Number of elements, where the coefficients must be computed.
    integer, intent(IN)                                         :: nelements

    ! Number of points per element, where the coefficients must be computed
    integer, intent(IN)                                         :: npointsPerElement

    ! This is an array of all points on all the elements where coefficients
    ! are needed.
    ! Remark: This usually coincides with rdomainSubset%p_DcubPtsReal.
    ! DIMENSION(dimension,npointsPerElement,nelements)
    real(DP), dimension(:,:,:), intent(IN)  :: Dpoints

    ! An array accepting the DOF's on all elements trial in the trial space.
    ! DIMENSION(\#local DOF's in test space,nelements)
    integer(PREC_DOFIDX), dimension(:,:), intent(IN) :: IdofsTest

    ! This is a t_domainIntSubset structure specifying more detailed information
    ! about the element set that is currently being integrated.
    ! It's usually used in more complex situations (e.g. nonlinear matrices).
    type(t_domainIntSubset), intent(IN)              :: rdomainIntSubset

    ! A pointer to a collection structure to provide additional
    ! information to the coefficient routine. May point to NULL() if not defined.
    type(t_collection), pointer                      :: p_rcollection

  !</input>

  !<output>
    ! A list of all coefficients in front of all terms in the linear form -
    ! for all given points on all given elements.
    !   DIMENSION(itermCount,npointsPerElement,nelements)
    ! with itermCount the number of terms in the linear form.
    real(DP), dimension(:,:,:), intent(OUT)                      :: Dcoefficients
  !</output>

  !</subroutine>

    real(DP) :: dtime

    ! In a nonstationary simulation, one can get the simulation time
    ! with the quick-access array of the collection.
    if (associated(p_rcollection)) then
      dtime = p_rcollection%Dquickaccess(1)
    else
      dtime = 0.0_DP
    end if

    Dcoefficients(:,:,:) = 0.0_DP

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine coeff_RHS_y (rdiscretisation,rform, &
                  nelements,npointsPerElement,Dpoints, &
                  IdofsTest,rdomainIntSubset,p_rcollection, &
                  Dcoefficients)

    use basicgeometry
    use triangulation
    use collection
    use scalarpde
    use domainintegration
    use dofmapping

  !<description>
    ! This subroutine is called during the vector assembly. It has to compute
    ! the coefficients in front of the terms of the linear form of the
    ! Y-velocity part of the right hand side vector.
    !
    ! The routine accepts a set of elements and a set of points on these
    ! elements (cubature points) in real coordinates.
    ! According to the terms in the linear form, the routine has to compute
    ! simultaneously for all these points and all the terms in the linear form
    ! the corresponding coefficients in front of the terms.
  !</description>

  !<input>
    ! The discretisation structure that defines the basic shape of the
    ! triangulation with references to the underlying triangulation,
    ! analytic boundary boundary description etc.
    type(t_spatialDiscretisation), intent(IN)                   :: rdiscretisation

    ! The linear form which is currently to be evaluated:
    type(t_linearForm), intent(IN)                              :: rform

    ! Number of elements, where the coefficients must be computed.
    integer, intent(IN)                                         :: nelements

    ! Number of points per element, where the coefficients must be computed
    integer, intent(IN)                                         :: npointsPerElement

    ! This is an array of all points on all the elements where coefficients
    ! are needed.
    ! Remark: This usually coincides with rdomainSubset%p_DcubPtsReal.
    ! DIMENSION(dimension,npointsPerElement,nelements)
    real(DP), dimension(:,:,:), intent(IN)  :: Dpoints

    ! An array accepting the DOF's on all elements trial in the trial space.
    ! DIMENSION(\#local DOF's in test space,nelements)
    integer(PREC_DOFIDX), dimension(:,:), intent(IN) :: IdofsTest

    ! This is a t_domainIntSubset structure specifying more detailed information
    ! about the element set that is currently being integrated.
    ! It's usually used in more complex situations (e.g. nonlinear matrices).
    type(t_domainIntSubset), intent(IN)              :: rdomainIntSubset

    ! A pointer to a collection structure to provide additional
    ! information to the coefficient routine. May point to NULL() if not defined.
    type(t_collection), pointer                      :: p_rcollection

  !</input>

  !<output>
    ! A list of all coefficients in front of all terms in the linear form -
    ! for all given points on all given elements.
    !   DIMENSION(itermCount,npointsPerElement,nelements)
    ! with itermCount the number of terms in the linear form.
    real(DP), dimension(:,:,:), intent(OUT)                      :: Dcoefficients
  !</output>

  !</subroutine>

    real(DP) :: dtime

    ! In a nonstationary simulation, one can get the simulation time
    ! with the quick-access array of the collection.
    if (associated(p_rcollection)) then
      dtime = p_rcollection%Dquickaccess(1)
    else
      dtime = 0.0_DP
    end if

    Dcoefficients(:,:,:) = 0.0_DP

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine ffunction_TargetX (cderivative,rdiscretisation, &
                nelements,npointsPerElement,Dpoints, &
                IdofsTest,rdomainIntSubset,p_rcollection, &
                Dvalues)

  use basicgeometry
  use triangulation
  use collection
  use scalarpde
  use domainintegration
  use derivatives
  use dofmapping

!<description>
  ! This subroutine is called during the postprocessing.
  ! It should return values of the analytical solution (if it is known).
  ! These are compared with the calculated solution to calculate the
  ! error in the X-velocity.
  !
  ! If the analytical solution is unknown, this routine doesn't make sense.
  ! In this case, error analysis should be deactivated in the .DAT files!
!</description>

!<input>
  ! This is a DER_xxxx derivative identifier (from derivative.f90) that
  ! specifies what to compute: DER_FUNC=function value, DER_DERIV_X=x-derivative,...
  ! The result must be written to the Dvalue-array below.
  integer, intent(IN)                                         :: cderivative

  ! The discretisation structure that defines the basic shape of the
  ! triangulation with references to the underlying triangulation,
  ! analytic boundary boundary description etc.
  type(t_spatialDiscretisation), intent(IN)                   :: rdiscretisation

  ! Number of elements, where the coefficients must be computed.
  integer, intent(IN)                                         :: nelements

  ! Number of points per element, where the coefficients must be computed
  integer, intent(IN)                                         :: npointsPerElement

  ! This is an array of all points on all the elements where coefficients
  ! are needed.
  ! DIMENSION(NDIM2D,npointsPerElement,nelements)
  ! Remark: This usually coincides with rdomainSubset%p_DcubPtsReal.
  real(DP), dimension(:,:,:), intent(IN)  :: Dpoints

  ! An array accepting the DOF's on all elements trial in the trial space.
  ! DIMENSION(\#local DOF's in trial space,Number of elements)
  integer(PREC_DOFIDX), dimension(:,:), intent(IN) :: IdofsTest

  ! This is a t_domainIntSubset structure specifying more detailed information
  ! about the element set that is currently being integrated.
  ! It's usually used in more complex situations (e.g. nonlinear matrices).
  type(t_domainIntSubset), intent(IN)              :: rdomainIntSubset

  ! A pointer to a collection structure to provide additional
  ! information to the coefficient routine. May point to NULL() if not defined.
  type(t_collection), pointer                      :: p_rcollection

!</input>

!<output>
  ! This array has to receive the values of the (analytical) function
  ! in all the points specified in Dpoints, or the appropriate derivative
  ! of the function, respectively, according to cderivative.
  !   DIMENSION(npointsPerElement,nelements)
  real(DP), dimension(:,:), intent(OUT)                      :: Dvalues
!</output>

!</subroutine>

    real(DP) :: dtime,dtimeMax
    integer :: itimedependence

    ! In a nonstationary simulation, one can get the simulation time
    ! with the quick-access array of the collection.
    if (associated(p_rcollection)) then
      dtime = p_rcollection%Dquickaccess(1)
      dtimeMax = p_rcollection%Dquickaccess(3)
      itimedependence = p_rcollection%Iquickaccess(1)
    else
      itimedependence = 0
      dtime = 0.0_DP
      dtimeMax = 0.0_DP
    end if

    Dvalues(:,:) = 0.0_DP

    if (cderivative .eq. DER_FUNC) then
      Dvalues(:,:) = (-dtime**2/100.+dtime/5.)*(Dpoints(1,:,:))
    end if

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine ffunction_TargetY (cderivative,rdiscretisation, &
                nelements,npointsPerElement,Dpoints, &
                IdofsTest,rdomainIntSubset,p_rcollection, &
                Dvalues)

  use basicgeometry
  use triangulation
  use collection
  use scalarpde
  use domainintegration
  use derivatives
  use dofmapping

!<description>
  ! This subroutine is called during the postprocessing.
  ! It should return values of the analytical solution (if it is known).
  ! These are compared with the calculated solution to calculate the
  ! error in the Y-velocity.
  !
  ! If the analytical solution is unknown, this routine doesn't make sense.
  ! In this case, error analysis should be deactivated in the .DAT files!
!</description>

!<input>
  ! This is a DER_xxxx derivative identifier (from derivative.f90) that
  ! specifies what to compute: DER_FUNC=function value, DER_DERIV_X=x-derivative,...
  ! The result must be written to the Dvalue-array below.
  integer, intent(IN)                                         :: cderivative

  ! The discretisation structure that defines the basic shape of the
  ! triangulation with references to the underlying triangulation,
  ! analytic boundary boundary description etc.
  type(t_spatialDiscretisation), intent(IN)                   :: rdiscretisation

  ! Number of elements, where the coefficients must be computed.
  integer, intent(IN)                                         :: nelements

  ! Number of points per element, where the coefficients must be computed
  integer, intent(IN)                                         :: npointsPerElement

  ! This is an array of all points on all the elements where coefficients
  ! are needed.
  ! DIMENSION(NDIM2D,npointsPerElement,nelements)
  ! Remark: This usually coincides with rdomainSubset%p_DcubPtsReal.
  real(DP), dimension(:,:,:), intent(IN)  :: Dpoints

  ! An array accepting the DOF's on all elements trial in the trial space.
  ! DIMENSION(\#local DOF's in trial space,Number of elements)
  integer(PREC_DOFIDX), dimension(:,:), intent(IN) :: IdofsTest

  ! This is a t_domainIntSubset structure specifying more detailed information
  ! about the element set that is currently being integrated.
  ! It's usually used in more complex situations (e.g. nonlinear matrices).
  type(t_domainIntSubset), intent(IN)              :: rdomainIntSubset

  ! A pointer to a collection structure to provide additional
  ! information to the coefficient routine. May point to NULL() if not defined.
  type(t_collection), pointer                      :: p_rcollection

!</input>

!<output>
  ! This array has to receive the values of the (analytical) function
  ! in all the points specified in Dpoints, or the appropriate derivative
  ! of the function, respectively, according to cderivative.
  !   DIMENSION(npointsPerElement,nelements)
  real(DP), dimension(:,:), intent(OUT)                      :: Dvalues
!</output>

!</subroutine>

    real(DP) :: dtime,dtimeMax
    integer :: itimedependence

    ! In a nonstationary simulation, one can get the simulation time
    ! with the quick-access array of the collection.
    if (associated(p_rcollection)) then
      dtime = p_rcollection%Dquickaccess(1)
      dtimeMax = p_rcollection%Dquickaccess(3)
      itimedependence = p_rcollection%Iquickaccess(1)
    else
      itimedependence = 0
      dtime = 0.0_DP
      dtimeMax = 0.0_DP
    end if

    Dvalues(:,:) = 0.0_DP

    if (cderivative .eq. DER_FUNC) then
      Dvalues(:,:) = (-dtime**2/100.+dtime/5.)*(-Dpoints(2,:,:))
    end if

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine ffunction_TargetP (cderivative,rdiscretisation, &
                nelements,npointsPerElement,Dpoints, &
                IdofsTest,rdomainIntSubset,p_rcollection, &
                Dvalues)

  use basicgeometry
  use triangulation
  use collection
  use scalarpde
  use domainintegration
  use derivatives
  use dofmapping

!<description>
  ! This subroutine is called during the postprocessing.
  ! It should return values of the analytical solution (if it is known).
  ! These are compared with the calculated solution to calculate the
  ! error in the pressure
  !
  ! If the analytical solution is unknown, this routine doesn't make sense.
  ! In this case, error analysis should be deactivated in the .DAT files!
!</description>

!<input>
  ! This is a DER_xxxx derivative identifier (from derivative.f90) that
  ! specifies what to compute: DER_FUNC=function value, DER_DERIV_X=x-derivative,...
  ! The result must be written to the Dvalue-array below.
  integer, intent(IN)                                         :: cderivative

  ! The discretisation structure that defines the basic shape of the
  ! triangulation with references to the underlying triangulation,
  ! analytic boundary boundary description etc.
  type(t_spatialDiscretisation), intent(IN)                   :: rdiscretisation

  ! Number of elements, where the coefficients must be computed.
  integer, intent(IN)                                         :: nelements

  ! Number of points per element, where the coefficients must be computed
  integer, intent(IN)                                         :: npointsPerElement

  ! This is an array of all points on all the elements where coefficients
  ! are needed.
  ! DIMENSION(NDIM2D,npointsPerElement,nelements)
  ! Remark: This usually coincides with rdomainSubset%p_DcubPtsReal.
  real(DP), dimension(:,:,:), intent(IN)  :: Dpoints

  ! An array accepting the DOF's on all elements trial in the trial space.
  ! DIMENSION(\#local DOF's in trial space,Number of elements)
  integer(PREC_DOFIDX), dimension(:,:), intent(IN) :: IdofsTest

  ! This is a t_domainIntSubset structure specifying more detailed information
  ! about the element set that is currently being integrated.
  ! It's usually used in more complex situations (e.g. nonlinear matrices).
  type(t_domainIntSubset), intent(IN)              :: rdomainIntSubset

  ! A pointer to a collection structure to provide additional
  ! information to the coefficient routine. May point to NULL() if not defined.
  type(t_collection), pointer                      :: p_rcollection

!</input>

!<output>
  ! This array has to receive the values of the (analytical) function
  ! in all the points specified in Dpoints, or the appropriate derivative
  ! of the function, respectively, according to cderivative.
  !   DIMENSION(npointsPerElement,nelements)
  real(DP), dimension(:,:), intent(OUT)                      :: Dvalues
!</output>

!</subroutine>

    real(DP) :: dtime,dtimeMax
    integer :: itimedependence

    ! In a nonstationary simulation, one can get the simulation time
    ! with the quick-access array of the collection.
    if (associated(p_rcollection)) then
      dtime = p_rcollection%Dquickaccess(1)
      dtimeMax = p_rcollection%Dquickaccess(3)
      itimedependence = p_rcollection%Iquickaccess(1)
    else
      itimedependence = 0
      dtime = 0.0_DP
      dtimeMax = 0.0_DP
    end if

    Dvalues(:,:) = 0.0_DP

    if (cderivative .eq. DER_FUNC) then
      ! ...
    end if

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine getBoundaryValues (Icomponents,rdiscretisation,rbcRegion,ielement, &
                                cinfoNeeded,iwhere,dwhere, p_rcollection, Dvalues)

  use collection
  use spatialdiscretisation
  use discretebc

!<description>
  ! This subroutine is called during the discretisation of boundary
  ! conditions. It calculates a special quantity on the boundary, which is
  ! then used by the discretisation routines to generate a discrete
  ! 'snapshot' of the (actually analytic) boundary conditions.
!</description>

!<input>
  ! Component specifier.
  ! For Dirichlet boundary:
  !   Icomponents(1) defines the number of the boundary component, the value
  !   should be calculated for (e.g. 1=1st solution component, e.g. X-velocitry,
  !   2=2nd solution component, e.g. Y-velocity,...)
  integer, dimension(:), intent(IN)                           :: Icomponents

  ! The discretisation structure that defines the basic shape of the
  ! triangulation with references to the underlying triangulation,
  ! analytic boundary boundary description etc.
  type(t_spatialDiscretisation), intent(IN)                   :: rdiscretisation

  ! Boundary condition region that is currently being processed.
  ! (This e.g. defines the type of boundary conditions that are
  !  currently being calculated, as well as information about the current
  !  boundary segment 'where we are at the moment'.)
  type(t_bcRegion), intent(IN)                                :: rbcRegion

  ! The element number on the boundary which is currently being processed
  integer(I32), intent(IN)                                    :: ielement

  ! The type of information, the routine should calculate. One of the
  ! DISCBC_NEEDxxxx constants. Depending on the constant, the routine has
  ! to return one or multiple information value in the result array.
  integer, intent(IN)                                         :: cinfoNeeded

  ! A reference to a geometric object where information should be computed.
  ! cinfoNeeded=DISCBC_NEEDFUNC :
  !   iwhere = number of the point in the triangulation or
  !          = 0, if only the parameter value of the point is known; this
  !               can be found in dwhere,
  ! cinfoNeeded=DISCBC_NEEDDERIV :
  !   iwhere = number of the point in the triangulation or
  !          = 0, if only the parameter value of the point is known; this
  !               can be found in dwhere,
  ! cinfoNeeded=DISCBC_NEEDINTMEAN :
  !   iwhere = number of the edge where the value integral mean value
  !            should be computed
  integer, intent(IN)                                         :: iwhere

  ! A reference to a geometric object where information should be computed.
  ! cinfoNeeded=DISCBC_NEEDFUNC :
  !   dwhere = parameter value of the point where the value should be computed,
  ! cinfoNeeded=DISCBC_NEEDDERIV :
  !   dwhere = parameter value of the point where the value should be computed,
  ! cinfoNeeded=DISCBC_NEEDINTMEAN :
  !   dwhere = 0 (not used)
  real(DP), intent(IN)                                        :: dwhere

  ! A pointer to a collection structure to provide additional
  ! information to the coefficient routine. May point to NULL() if not defined.
  type(t_collection), pointer                  :: p_rcollection

!</input>

!<output>
  ! This array receives the calculated information. If the caller
  ! only needs one value, the computed quantity is put into Dvalues(1).
  ! If multiple values are needed, they are collected here (e.g. for
  ! DISCBC_NEEDDERIV: Dvalues(1)=x-derivative, Dvalues(2)=y-derivative,...)
  !
  ! The function may return SYS_INFINITY_DP as a value. This indicates the
  ! framework to ignore the node and treat it as 'natural boundary condition'
  ! node.
  real(DP), dimension(:), intent(OUT)                         :: Dvalues
!</output>

!</subroutine>

    integer :: icomponent,iexprtyp

    real(DP) :: dtime

    ! In a nonstationary simulation, one can get the simulation time
    ! with the quick-access array of the collection.
    if (associated(p_rcollection)) then
      dtime = p_rcollection%Dquickaccess(1)
    else
      dtime = 0.0_DP
    end if

    ! Use boundary conditions from DAT files.
    select case (cinfoNeeded)

    case (DISCBC_NEEDFUNC,DISCBC_NEEDFUNCMID,DISCBC_NEEDDERIV, &
          DISCBC_NEEDINTMEAN,DISCBC_NEEDNORMALSTRESS)

      ! Dirichlet boundary conditions

      ! Get from the current component of the PDE we are discretising:
      icomponent = Icomponents(1)

      ! -> 1=X-velocity, 2=Y-velocity.

      ! Return zero Dirichlet boundary values for all situations by default.
      Dvalues(1) = 0.0_DP

      ! Now, depending on the problem, calculate the return value.

      ! Get the type of the expression to evaluate from the
      ! integer tag of the BC-region - if there is an expression to evaluate
      ! at all.
      iexprtyp = rbcRegion%ibdrexprtype

      ! Now, which boundary condition do we have here?
      select case (rbcRegion%ctype)
      case (BC_DIRICHLET)
        ! Simple Dirichlet BC's. Evaluate the expression iexprtyp.
        Dvalues(1) = evalBoundary (rdiscretisation, rbcRegion%rboundaryRegion, &
                                    iexprtyp, rbcRegion%itag, rbcRegion%dtag, &
                                    dwhere, rbcRegion%stag,dtime,&
                                    p_rcollection)

      case (BC_PRESSUREDROP)
        ! Normal stress / pressure drop. Evaluate Evaluate the
        ! expression iexprtyp.
        Dvalues(1) = evalBoundary (rdiscretisation, rbcRegion%rboundaryRegion, &
                                    iexprtyp, rbcRegion%itag, rbcRegion%dtag, &
                                    dwhere, rbcRegion%stag,dtime,&
                                    p_rcollection)
      end select

    end select

  contains

    ! Auxiliary function: Evaluate a scalar expression on the boundary.

    real(DP) function evalBoundary (rdiscretisation, rboundaryRegion, &
                                    ityp, ivalue, dvalue, dpar, stag, dtime, p_rcollection)

    ! Discretisation structure of the underlying discretisation
    type(t_spatialDiscretisation), intent(IN) :: rdiscretisation

    ! Current boundary region
    type(t_boundaryRegion), intent(IN) :: rboundaryRegion

    ! Type of expression to evaluate.
    ! One of the BDC_xxxx constants from cc2dmediumm2boundarydef.f90.
    integer, intent(IN) :: ityp

    ! Integer tag. If ityp=BDC_EXPRESSION, this must specify the number of
    ! the expression in the expression object to evaluate.
    ! Otherwise unused.
    integer, intent(IN) :: ivalue

    ! Double precision parameter for simple expressions
    real(DP), intent(IN) :: dvalue

    ! Current parameter value of the point on the boundary.
    ! 0-1-parametrisation.
    real(DP), intent(IN) :: dpar

    ! String tag that defines more complicated BC's.
    character(LEN=*), intent(IN) :: stag

    ! For nonstationary simulation: Simulation time.
    ! =0 for stationary simulations.
    real(DP), intent(IN) :: dtime

    ! A compiled expression for evaluation at runtime
    type(t_fparser), pointer :: p_rparser

    ! A pointer to a collection structure to provide additional
    ! information to the coefficient routine. May point to NULL() if not defined.
    type(t_collection), pointer                  :: p_rcollection

      ! local variables
      real(DP) :: d,dx,dy
      character(LEN=PARLST_MLDATA) :: sexpr
      real(DP), dimension(size(SEC_EXPRVARIABLES)) :: Rval

      select case (ityp)
      case (BDC_USERDEF)
        ! This is a hardcoded, user-defined identifier.
        ! In stag, the name of the identifier is noted.
        ! Get the identifier itself from the collection.
        call collct_getvalue_string (p_rcollection, stag, sexpr, &
                                     0, SEC_SBDEXPRESSIONS)

        ! Now we can decide on 'sexpr' how to evaluate.
        ! By default, we return 0.0. A user defined calculation can be added here!
        evalBoundary = 0.0_DP

      case (BDC_VALDOUBLE)
        ! A simple constant, given by dvalue
        evalBoundary = dvalue

      case (BDC_EXPRESSION)
        ! A complex expression.
        ! Get the expression object from the collection.

        p_rparser => collct_getvalue_pars (p_rcollection, BDC_BDPARSER, &
                                   0, SEC_SBDEXPRESSIONS)

        ! Set up an array with variables for evaluating the expression.
        ! Give the values in exactly the same order as specified
        ! by SEC_EXPRVARIABLES!
        Rval = 0.0_DP

        call boundary_getCoords(rdiscretisation%p_rboundary, &
                                rboundaryRegion%iboundCompIdx, &
                                dpar, dx, dy)

        ! Get the local parameter value 0 <= d <= 1.
        ! Note that if dpar < rboundaryRegion%dminParam, we have to add the maximum
        ! parameter value on the boundary to dpar as normally 0 <= dpar < max.par.
        ! although 0 <= dminpar <= max.par
        !      and 0 <= dmaxpar <= max.par!
        d = dpar
        if (d .lt. rboundaryRegion%dminParam) &
          d = d + boundary_dgetMaxParVal(rdiscretisation%p_rboundary,&
                                         rboundaryRegion%iboundCompIdx)
        d = d - rboundaryRegion%dminParam

        Rval(1) = dx
        Rval(2) = dy
        ! Rval(3) = .
        Rval(4) = d
        Rval(5) = dpar
        Rval(6) = boundary_convertParameter(rdiscretisation%p_rboundary, &
                                            rboundaryRegion%iboundCompIdx, dpar, &
                                            BDR_PAR_01, BDR_PAR_LENGTH)
        Rval(7) = dtime

        ! Evaluate the expression. ivalue is the number of
        ! the expression to evaluate.
        call fparser_evalFunction (p_rparser, ivalue, Rval, evalBoundary)

      case (BDC_VALPARPROFILE)
        ! A parabolic profile. dvalue expresses the
        ! maximum value of the profile.
        !
        ! Get the local parameter value 0 <= d <= 1.
        ! Note that if dpar < rboundaryRegion%dminParam, we have to add the maximum
        ! parameter value on the boundary to dpar as normally 0 <= dpar < max.par.
        ! although 0 <= dminpar <= max.par
        !      and 0 <= dmaxpar <= max.par!
        d = dpar
        if (d .lt. rboundaryRegion%dminParam) &
          d = d + boundary_dgetMaxParVal(rdiscretisation%p_rboundary,&
                                         rboundaryRegion%iboundCompIdx)
        d = d - rboundaryRegion%dminParam

        evalBoundary = mprim_getParabolicProfile (d,1.0_DP,dvalue)
      end select

    end function

  end subroutine

  ! ***************************************************************************
  ! Values in a fictitious boundary component:

!<subroutine>

  subroutine getBoundaryValuesFBC (Icomponents,rdiscretisation,rbcRegion, &
                                   Revaluation, p_rcollection)

  use collection
  use spatialdiscretisation
  use discretefbc

!<description>
  ! This subroutine is called during the discretisation of boundary
  ! conditions on fictitious boundary components. It calculates a special quantity
  ! on the boundary, which is then used by the discretisation routines to
  ! generate a discrete 'snapshot' of the (actually analytic) boundary conditions.
  !
  ! The routine must calculate the values on all elements of the element
  ! list Ielements simultaneously. Iwhere is a list with vertex or edge numbers
  ! where information is to be retrieved. Dvalues is filled with function values
  ! while Binside is set to TRUE for every vertex/edge that is inside of the
  ! corresponding fictitious boundary region (identified by rbcRegion).
!</description>

!<input>
  ! Component specifier.
  ! For Dirichlet boundary:
  !   Icomponents(1..SIZE(Icomponents)) defines the number of the solution component,
  !   the value should be calculated for
  !   (e.g. 1=1st solution component, e.g. X-velocity,
  !         2=2nd solution component, e.g. Y-velocity,...,
  !         3=3rd solution component, e.g. pressure)
  !   Example: Icomponents(:) = [1,2] -> Compute values for X- and Y-velocity
  !     (1=x, 2=y component)
  integer, dimension(:), intent(IN)                           :: Icomponents

  ! The discretisation structure that defines the basic shape of the
  ! triangulation with references to the underlying triangulation,
  ! analytic boundary boundary description etc.
  type(t_blockDiscretisation), intent(IN)                     :: rdiscretisation

  ! Boundary condition region that is currently being processed.
  ! (This e.g. defines the type of boundary conditions that are
  !  currently being calculated (Dirichlet, Robin,...), as well as information
  !  about the current boundary region 'what is discretised at the moment'.)
  type(t_bcRegion), intent(IN)                                :: rbcRegion

  ! A pointer to a collection structure to provide additional
  ! information to the coefficient routine. May point to NULL() if not defined.
  type(t_collection), pointer                                 :: p_rcollection

!</input>

!<inputoutput>
  ! A t_discreteFBCevaluation structure array that defines what to evaluate,
  ! where to evaluate and which accepts the return values.
  ! This callback routine must check out the cinfoNeeded-entry in this structure
  ! to find out what to evaluate.
  ! The other entries in this structure describe where to evaluate.
  ! The result of the evaluation must be written into the p_Dvalues array entry
  ! in this structure.
  !
  ! The number of structures in this array depend on what to evaluate:
  !
  ! For Dirichlet boundary:
  !   revaluation contains as many entries as Icomponents; every entry in
  !   Icomponent corresponds to one entry in revaluation
  !   (so Icomponent(1)=1 defines to evaluate the X-velocity while the
  !    values for the X-velocity are written to revaluation(1)\%p_Dvalues;
  !    Icomponent(2)=2 defines to evaluate the Y-velocity while the values
  !    for the Y-velocity are written to revaluation(2)\%p_Dvalues, etc).
  !
  type(t_discreteFBCevaluation), dimension(:), intent(INOUT) :: Revaluation
!</inputoutput>

!</subroutine>

    ! local variables
    real(DP) :: ddistance, dxcenter, dycenter, dradius, dx, dy
    real(DP), dimension(:,:), pointer :: p_DvertexCoordinates
    integer(PREC_POINTIDX), dimension(:,:), pointer :: p_IverticesAtElement
    integer(PREC_POINTIDX), dimension(:,:), pointer :: p_IverticesAtEdge
    type(t_triangulation), pointer :: p_rtriangulation
    integer :: ipoint,idx

    ! Note: the definition of (analytic) fictitious boundary components
    ! is performed in 'c2d2_parseFBDconditions'.

    ! Are we evaluating our fictitious boundary component?
    if (rbcRegion%rfictBoundaryRegion%sname .eq. 'CIRCLE') then

      ! Get the triangulation array for the point coordinates
      p_rtriangulation => rdiscretisation%RspatialDiscr(1)%p_rtriangulation
      call storage_getbase_double2d (p_rtriangulation%h_DvertexCoords,&
                                     p_DvertexCoordinates)
      call storage_getbase_int2d (p_rtriangulation%h_IverticesAtElement,&
                                  p_IverticesAtElement)
      call storage_getbase_int2d (p_rtriangulation%h_IverticesAtEdge,&
                                  p_IverticesAtEdge)

      ! Definition of the circle
      dxcenter = 0.6
      dycenter = 0.2
      dradius  = 0.05

      ! Loop through the points where to evaluate:
      do idx = 1,Revaluation(1)%nvalues

        ! Get the number of the point to process; may also be number of an
        ! edge or element...
        ipoint = Revaluation(1)%p_Iwhere(idx)

        ! Get x- and y-coordinate
        call getXYcoord (Revaluation(1)%cinfoNeeded,ipoint,&
                         p_DvertexCoordinates,&
                         p_IverticesAtElement,p_IverticesAtEdge,&
                         p_rtriangulation%NVT,&
                         dx,dy)

        ! Get the distance to the center
        ddistance = sqrt( (dx-dxcenter)**2 + (dy-dycenter)**2 )

        ! Point inside?
        if (ddistance .le. dradius) then

          ! Denote in the p_Iinside array that we prescribe a value here:
          Revaluation(1)%p_Iinside (idx) = 1
          Revaluation(2)%p_Iinside (idx) = 1

          ! We prescribe 0.0 as Dirichlet value here - for x- and y-velocity
          Revaluation(1)%p_Dvalues (idx,1) = 0.0_DP
          Revaluation(2)%p_Dvalues (idx,1) = 0.0_DP

        end if

      end do

!      ! Definition of a 2nd circle
!      dxcenter = 1.1
!      dycenter = 0.31
!      dradius  = 0.05
!
!      ! Loop through the points where to evaluate:
!      DO idx = 1,Revaluation(1)%nvalues
!
!        ! Get the number of the point to process; may also be number of an
!        ! edge or element...
!        ipoint = Revaluation(1)%p_Iwhere(idx)
!
!        ! Get x- and y-coordinate
!        CALL getXYcoord (Revaluation(1)%cinfoNeeded,ipoint,&
!                         p_DvertexCoordinates,&
!                         p_IverticesAtElement,p_IverticesAtEdge,&
!                         p_rtriangulation%NVT,&
!                         dx,dy)
!
!        ! Get the distance to the center
!        ddistance = SQRT( (dx-dxcenter)**2 + (dy-dycenter)**2 )
!
!        ! Point inside?
!        IF (ddistance .LE. dradius) THEN
!
!          ! Denote in the p_Iinside array that we prescribe a value here:
!          Revaluation(1)%p_Iinside (idx) = 1
!          Revaluation(2)%p_Iinside (idx) = 1
!
!          ! We prescribe 0.0 as Dirichlet value here - vor X- and Y-velocity
!          Revaluation(1)%p_Dvalues (idx,1) = 0.0_DP
!          Revaluation(2)%p_Dvalues (idx,1) = 0.0_DP
!
!        END IF
!
!      END DO

    end if

  contains

    ! ---------------------------------------------------------------
    ! Auxiliary routine: Get X/Y coordinate of a point.
    ! This is either the X/Y coordinate of a corner point or
    ! of an edge-midpoint, depending of cinfoNeeded.

    subroutine getXYcoord (cinfoNeeded,iwhere,DvertexCoords,&
                           IverticesAtElement,IverticesAtEdge,NVT,&
                           dx,dy)

    ! One of the DISCFBC_NEEDxxxx constanr. DISCFBC_NEEDFUNC interprets
    ! iwhere as ivt (vertex number) and returns the coordinates of that
    ! vertex. DISCFBC_NEEDFUNCMID/DISCFBC_NEEDINTMEAN interprets iwhere
    ! as imid (edge number) and returns the coordinates of the edge
    ! midpoint. DISCFBC_NEEDFUNCELMID interprets iwhere as element number
    ! and returns the element midpoint.
    integer, intent(IN) :: cinfoNeeded

    ! Identifier. Either vertex number, edge number or element number,
    ! depending on cinfoNeeded.
    integer(I32), intent(IN) :: iwhere

    ! Array with coordinates of all corner vertices (DCORVG)
    real(DP), dimension(:,:), intent(IN) :: DvertexCoords

    ! Array with numbers of corner coordinates for all elements (KVERT)
    integer(PREC_POINTIDX), dimension(:,:), intent(IN) :: IverticesAtElement

    ! Array with numbers of points adjacent to all edges
    integer(PREC_POINTIDX), dimension(:,:), intent(IN) :: IverticesAtEdge

    ! Number of vertices in the triangulation
    integer(PREC_POINTIDX), intent(IN) :: NVT

    ! Output: X-coordinate
    real(DP), intent(OUT) :: dx

    ! Output: Y-coordinate
    real(DP), intent(OUT) :: dy

      ! local variables
      real(DP) :: dm1,dm2
      integer :: i,j
      integer(PREC_POINTIDX) :: iv1,iv2

      ! Let's see, what do we have...
      select case (cinfoNeeded)
      case (DISCFBC_NEEDFUNC)
        ! That's easy
        dx = DvertexCoords(1,iwhere)
        dy = DvertexCoords(2,iwhere)

      case (DISCFBC_NEEDFUNCMID,DISCFBC_NEEDINTMEAN)
        ! Not much harder; get the two points on the edge and calculate
        ! the mean coordinates.
        iv1 = IverticesAtEdge(1,iwhere)
        iv2 = IverticesAtEdge(2,iwhere)
        dx = 0.5_DP*(DvertexCoords(1,iv1)+DvertexCoords(1,iv2))
        dy = 0.5_DP*(DvertexCoords(2,iv1)+DvertexCoords(2,iv2))

      case (DISCFBC_NEEDFUNCELMID)
        ! A little bit more to do; we have three or four corners
        dx = 0.0_DP
        dy = 0.0_DP
        do i=1,size(IverticesAtElement,1)
          iv1 = IverticesAtElement(i,iwhere)
          if (iv1 .ne. 0) then
            dx = dx + DvertexCoords(1,iv1)
            dy = dy + DvertexCoords(2,iv1)
            j = i
          end if
        end do

        dx = dx / real(j,DP)
        dy = dy / real(j,DP)

      case DEFAULT
        print *,'getBoundaryValuesFBC: Insupported coordinate type to return.'
        stop
      end select

    end subroutine

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
      subroutine sub()
      parameter (n = 3)
      parameter (m = 4)
      integer ka(n, m)
      data ka/11,21,31,12, 22, 32/

      do 10 k = 1,n
         do 20 j = 1,m
            print 9000, k, j, ka(k, j)
  020    continue
   10 continue
 9000 format(1x, 'KA(', i1, ', ', i1, ') = ', i5)

      return
      end

      program init
      integer iarray(3,2)
      data iarray/11,21,31,12,22,32/

      do 10 k = 1,3
         do 20 j = 1,2
            print 9000, k, j, iarray(k, j)
  020    continue
   10 continue
 9000 format(1x, 'X(', i1, ', ', i1, ') = ', i5)

      call sub
      end

ext: f

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
Then describe the Wiki Animals animal here!

Some ideas: What do you like best about your pet? What do you not like? How old is it in human years? What is the funniest/strangest/most remarkable thing it has done? Would you recommend this breed to others? Whats involved in caring for this type of animal (does it have to be walked, milked, sheared)? What costs (food, vaccinations, etc.)? Are there restrictions on where/who can own it? Who should not get this type of pet? Is it childfriendly? Furniturefriendly?

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! PWR001: Declare global variables as function parameters
! https://www.appentra.com/knowledge/checks/pwr001

PROGRAM example
    INTEGER :: global = 10
    CALL s
CONTAINS
    SUBROUTINE s
        IMPLICIT NONE
        global = 20
    END SUBROUTINE
END PROGRAM

max_stars_repo_name: teamappentra/check-examples

max_issues_repo_name: teamappentra/check-examples

max_forks_repo_name: teamappentra/check-examples


text:
module diag_axis_mod

! <CONTACT EMAIL="Giang.Nong@noaa.gov">
!   Giang Nong
! </CONTACT>

! <OVERVIEW>
!   <TT>diag_axis_mod</TT> is an integral part of diag_manager_mod. It helps to create
! axis IDs that are used in register_diag_field.
! </OVERVIEW>

! <DESCRIPTION>
! Users first create axis ID by calling diag_axis_init, then use this axis ID in
! register_diag_field.
!
!</DESCRIPTION>

use mpp_domains_mod, only: domain1d, domain2d, mpp_get_compute_domain, &
                           mpp_get_domain_components, null_domain1d, &
                           null_domain2d, operator(.ne.)
use         fms_mod, only: error_mesg, write_version_number, lowercase, FATAL
use diag_data_mod, only  : diag_axis_type, max_subaxes, max_axes
implicit none

private
public  diag_axis_init, get_diag_axis, get_domain1d, get_domain2d, &
        get_axis_length, get_axis_global_length, diag_subaxes_init, &
        get_diag_axis_cart, get_diag_axis_data, max_axes, get_axis_aux

!counter of number of axes defined
integer,allocatable :: num_subaxes(:)
integer             :: num_def_axes = 0

!-----------------------------------------------------------------------
!
!  name        = short name for axis
!  units       = units for axis
!  long_name   = long name for axis
!  cart_name   = cartesian axis ('x','y','z', or 't')
!  length      = number of coordinates for axis
!  direction   = if +1, data are in a positive direction (default:+1)
!                if -1, data are in a negative direction
!  data        = array of coordinate values for this axis
!
!-----------------------------------------------------------------------

! ---- storage for axis set names ----

integer, parameter :: max_num_axis_sets = 25
integer            :: num_axis_sets = 0
character(len=128) :: Axis_sets(max_num_axis_sets)
! ---- global storage for all defined axes ----
type (diag_axis_type), allocatable, save :: Axes(:)
logical            :: module_is_initialized = .FALSE.
character(len=128) :: &
     version='$Id: diag_axis.f90,v 13.0 2006/03/28 21:37:49 fms Exp $'
character(len=128) :: tagname='$Name: latest $'

contains
!#######################################################################

function diag_axis_init (name, data, units, cart_name, long_name,     &
       direction, set_name, edges, Domain, Domain2, aux) &
       result (indexx)

! increment axis counter and fill in axes
!-----------------------------------------------------------------------
!  name                = short name for axis
!  data                = array of coordinate values for this axis
!  units               = units for axis
!  cart_name           = cartesian axis ("x",'y','z','t')
!  direction(optional) = if +1, data are in a up   direction
!                      = if -1, data are in a down direction
!                      = if  0, neither up or down (default)
!  long_name(optional) = long name for axis (default: name)
!  edges    (optional) = axis id for the previously defined "edges axis"
!  aux      (optional) = auxiliary name, can be only either geolon_t or geolat_t
!-----------------------------------------------------------------------

  character(len=*), intent(in)           :: name
  real            , intent(in)           :: data(:)
  character(len=*), intent(in)           :: units
  character(len=*), intent(in)           :: cart_name
  character(len=*), intent(in), optional :: long_name, set_name
  integer         , intent(in), optional :: direction, edges
  type(domain1d)  , intent(in), optional :: Domain
  type(domain2d)  , intent(in), optional :: Domain2
  character(len=*), intent(in), optional :: aux
  type(domain1d)                         :: domain_x, domain_y
  integer                                :: indexx, ierr, axlen
  integer                                :: i, set

  if ( .not.module_is_initialized ) then
     call write_version_number( version, tagname )
  endif
  if (.not. allocated(Axes)) allocate(Axes(max_axes))
  if (.not. allocated(num_subaxes)) then
     allocate(num_subaxes(max_axes))
     num_subaxes = 0
  endif
!---- is there an axis set? ----
  if ( present(set_name) ) then
     set = get_axis_set_num (set_name)
!    ---- add new set name ----
     if (set == 0) then
        num_axis_sets = num_axis_sets + 1
        if (num_axis_sets > max_num_axis_sets)   &
             call error_mesg('diag_axis_init in diag_axis_mod',  &
             'exceeded max_num_axis_sets', FATAL)
        set = num_axis_sets
        Axis_sets(set) = set_name
     endif
  else
     set = 0
  endif

!---- see if axis already exists --
! if this is time axis, return the ID of a previously defined
! if this is spatial axis, FATAL error
  do i = 1, num_def_axes
     if (trim(name) == Axes(i)%name .and. trim(name) == 'Stations') then
        indexx = i
        return
     endif
     if (trim(name) == Axes(i)%name .and. trim(name) == 'Levels') then
        indexx = i
        return
     endif
     if (name == Axes(i)%name .and. set == Axes(i)%set) then
        if(trim(lowercase(name))=='time'.or.trim(lowercase(cart_name))=='t') then
           indexx = i
           return
        else
           call error_mesg ('diag_axis_init in diag_axis_mod', &
                'axis_name '//trim(name)//' and axis_set already exists', FATAL)
        endif
     endif
  enddo

!---- register axis ----
  num_def_axes = num_def_axes + 1
  if (num_def_axes > max_axes) call error_mesg ('diag_axis_init in diag_axis_mod', &
       'max_axes exceeded, increase it via diag_manager_nml', FATAL)
  indexx = num_def_axes

!---- check and then save cart_name name ----
  if (cart_name == 'x' .or.cart_name == 'X' ) then
     Axes(indexx)%cart_name = 'X'
  else if (cart_name == 'y' .or.cart_name == 'Y' ) then
     Axes(indexx)%cart_name = 'Y'
  else if (cart_name == 'z' .or.cart_name == 'Z' ) then
     Axes(indexx)%cart_name = 'Z'
  else if (cart_name == 't' .or.cart_name == 'T' ) then
      Axes(indexx)%cart_name = 'T'
  else if (cart_name == 'n' .or.cart_name == 'N' ) then
      Axes(indexx)%cart_name = 'N'
  else
      call error_mesg('diag_axis_init in diag_axis_mod', 'Invalid cart_name name.', FATAL)
  endif

!---- allocate storage for coordinate values of axis ----
  axlen = size(data(:))
  if ( Axes(indexx)%cart_name == 'T' ) axlen = 0
  allocate ( Axes(indexx)%data(1:axlen) )
  Axes(indexx)%name   = trim(name)
  Axes(indexx)%data   = data(1:axlen)
  Axes(indexx)%units  = units
  Axes(indexx)%length = axlen
  Axes(indexx)%set    = set
! start and end are used in subaxes information only
  Axes(indexx)%start = -1
  Axes(indexx)%end = -1
  Axes(indexx)%subaxis_name = ""

  if (present(long_name))then
     Axes(indexx)%long_name = long_name
  else
     Axes(indexx)%long_name = name
  endif
  if (present(aux))then
     Axes(indexx)%aux = trim(aux)
  else
     Axes(indexx)%aux = 'none'
  endif

!---- axis direction (-1, 0, or +1) ----
  if (present(direction))then
     if(abs(direction) /= 1) call error_mesg('diag_axis_init in diag_axis_mod', &
          'direction must be either +1 or -1',FATAL)
     Axes(indexx)%direction = direction
  else
     Axes(indexx)%direction = 0
  endif
!---- domain2d type ----
  if ( present(Domain2) .or. present(Domain)) then
    if ( Axes(indexx)%cart_name /= 'X' .and. Axes(indexx)%cart_name /= 'Y') then
      call error_mesg('diag_axis_init in diag_axis_mod', &
           'Domain must not be present for an axis which is not in the X or Y direction', &
            FATAL)
    endif
  endif
  if ( present(Domain2) .and. present(Domain)) then
    call error_mesg('diag_axis_init in diag_axis_mod', &
           'Presence of both Domain and Domain2 at the same time is prohibited', &
            FATAL)
  endif
  if ( present(Domain2) ) then
     Axes(indexx)%Domain2 = Domain2
     call mpp_get_domain_components(Domain2, domain_x, domain_y)
     if ( Axes(indexx)%cart_name == 'X' ) Axes(indexx)%Domain = domain_x
     if ( Axes(indexx)%cart_name == 'Y' ) Axes(indexx)%Domain = domain_y
  else
     Axes(indexx)%Domain2 = null_domain2d
!---- domain1d type ----
     if ( present(Domain)) then
        Axes(indexx)%Domain = Domain
     else
        Axes(indexx)%Domain = null_domain1d
     endif
  endif
!---- have axis edges been defined ? ----
  Axes(indexx)%edges = 0
  if (present(edges))then
     if ( edges > 0 .and. edges < num_def_axes ) then
        ierr=0
        if ( Axes(edges)%cart_name /= Axes(indexx)%cart_name) ierr=1
        if ( Axes(edges)%length    /= Axes(indexx)%length+1 ) ierr=ierr+2
        if ( Axes(edges)%set       /= Axes(indexx)%set      ) ierr=ierr+4
        if ( ierr > 0 )   call error_mesg ('diag_axis_init in diag_axis_mod', &
             'Edges axis does not match axis', FATAL)
        Axes(indexx)%edges = edges
     else
       call error_mesg ('diag_axis_init in diag_axis_mod', &
                        'Edges axis is not defined', FATAL)
     endif
  endif
  module_is_initialized = .TRUE.
!-----------------------------------------------------------------------
end function diag_axis_init

function diag_subaxes_init(axis,subdata,start_indx,end_indx,domain_1d,domain_2d)  result(index)
! Given ID of parent axis and data of subaxis, this function returns ID of corresponding subaxis
! subaxis is defined on parent axis from start_indx to end_indx
  integer, intent(in)            :: axis ! ID of parent axis
  real, dimension(:), intent(in) :: subdata ! data of subaxis
  integer, intent(in)            :: start_indx ! start index of subaxis
  integer, intent(in)            :: end_indx ! end index of subaxis
  type(domain1d), intent(in), optional  :: domain_1d
  type(domain2d), intent(in), optional  :: domain_2d

  integer                        :: index
  integer                        :: i,nsub_axis
  character(len=128)             :: name, nsub_name
  character(len=128)             :: units
  character(len=128)             :: cart_name
  character(len=128)             :: long_name
  logical                        :: subaxis_set

! there may be more than 1 subaxis on a parent axis, check for redundancy
  nsub_axis = 0; subaxis_set = .false.
  do i = 1, num_subaxes(axis)
     if(start_indx == Axes(axis)%start(i) .and. end_indx == Axes(axis)%end(i)) then
        nsub_axis = i
        subaxis_set = .true.    !subaxis already exists
        name = trim(Axes(axis)%subaxis_name(nsub_axis))
        exit
     endif
  enddo
  if (nsub_axis == 0) then  ! create new subaxis
     num_subaxes(axis) = num_subaxes(axis) + 1
     if (num_subaxes(axis) > max_subaxes) &
          call error_mesg ('diag_subaxes_init in diag_axis_mod',' max_subaxes too small, increase max_subaxes', FATAL)
     nsub_axis = num_subaxes(axis)
     Axes(axis)%start(nsub_axis) = start_indx
     Axes(axis)%end(nsub_axis)   = end_indx
  endif

  ! Create new name for the subaxis from name of parent axis

! If subaxis already exists, get the index and return
  if(subaxis_set) then
     index = get_axis_num(name)
  else
! get a new index for subaxis
     write(nsub_name,'(I1)') nsub_axis
     name = trim(Axes(axis)%name)//'_sub'//trim(nsub_name)
     Axes(axis)%subaxis_name(nsub_axis) = name
     long_name = trim(Axes(axis)%long_name)
     units = trim(Axes(axis)%units)
     cart_name = trim(Axes(axis)%cart_name)
     index =  diag_axis_init (trim(name), subdata, trim(units), trim(cart_name), trim(long_name), &
          Domain2=domain_2d)
  endif
end function diag_subaxes_init
!#######################################################################

subroutine get_diag_axis (id, name, units, long_name, cart_name, &
                           direction, edges, Domain, data)

! Return information about the axis with index id

!-----------------------------------------------------------------------
!  id         =  axis number
!  name       = short name for axis
!  units      = units for axis
!  long_name  = long name for axis
!  cart_name  = cartesian axis ("x",'y','z', or 't')
!  direction  = if +1, data are in a positive direction (default:+1)
!               if -1, data are in a negative direction
!  data       = array of coordinate values for this axis
!-----------------------------------------------------------------------

  character(len=*), intent(out) :: name, units, long_name, cart_name
  integer, intent(in)           :: id
  type(domain1d), intent(out)   :: Domain
  integer, intent(out)          :: direction, edges
  real, intent(out)             :: data(:)
  character(len=128)            :: error_msg

  if (id < 1 .or. id > num_def_axes) then
     write(error_msg,'(i2)')id
     call error_mesg('get_diag_axis in diag_axis_mod', &
          trim(error_msg)//' is illegal value for axis_id', FATAL)
  endif
  name      = Axes(id)%name
  units     = Axes(id)%units
  long_name = Axes(id)%long_name
  cart_name = Axes(id)%cart_name
  direction = Axes(id)%direction
  edges     = Axes(id)%edges
  Domain    = Axes(id)%Domain
  if (Axes(id)%length > size(data(:))) call error_mesg ('get_diag_axis in diag_axis_mod', &
       'array data is too small', FATAL)
  data(1:Axes(id)%length) = Axes(id)%data
!#####################################################################
end subroutine get_diag_axis
subroutine get_diag_axis_cart(id, cart_name)
! Return the axis cartesian from  index id
!  id         =  axis number
!  cart_name  = cartesian axis ("x",'y','z', or 't')
  character(len=*), intent(out) :: cart_name
  integer, intent(in)           :: id
  cart_name = Axes(id)%cart_name
end subroutine get_diag_axis_cart

!######################################################################
subroutine get_diag_axis_data(id,data)
! Return the axis data from  index id
!  id         =  axis number
  integer, intent(in) :: id
  real, intent(out)   :: data(:)
  if (Axes(id)%length > size(data(:))) call error_mesg ('get_diag_axis_data in diag_axis_mod', &
       'array data is too small', FATAL)
  data(1:Axes(id)%length) = Axes(id)%data
end subroutine get_diag_axis_data

!######################################################################

function get_axis_length (id) result (length)
  integer, intent(in) :: id
  integer             :: length
  if ( Axes(id)%Domain /= null_domain1d ) then
     call mpp_get_compute_domain(Axes(id)%Domain,size=length)
  else
     length = Axes(id) % length
  endif
end function get_axis_length
!######################################################################
function get_axis_aux (id) result (aux)

  integer, intent(in)   :: id
  character(len=128)    :: aux

  aux =  Axes(id)%aux
end function get_axis_aux

!#######################################################################

function get_axis_global_length (id) result (length)
  integer, intent(in) :: id
  integer             :: length
  length = Axes(id) % length
end function get_axis_global_length
!#######################################################################

function get_domain1d (id) result (Domain1)
  integer, intent(in) :: id
  type(domain1d)      :: Domain1
  if (Axes(id)%Domain .NE. NULL_DOMAIN1D) then
     Domain1 = Axes(id)%Domain
  else
     Domain1 = NULL_DOMAIN1D
  endif
end function get_domain1d
!#######################################################################

function get_domain2d (ids) result (Domain2)
  integer, intent(in) :: ids(:)
  type(domain2d)      :: Domain2
  integer             :: i, id, flag

  if ( size(ids(:)) < 1 .or. size(ids(:)) > 4 ) call error_mesg  &
       ('get_domain2d in diag_axis_mod', &
       'input argument has incorrect size', FATAL)
  Domain2 = null_domain2d
  flag = 0
  do i = 1, size(ids(:))
     id = ids(i)
     if ( Axes(id)%cart_name == 'X' .or.  &
          Axes(id)%cart_name == 'Y' ) flag = flag + 1
!     --- both x/y axes found ---
     if ( flag == 2 ) then
        if (Axes(id)%Domain2 .NE. NULL_DOMAIN2D) then
           Domain2 = Axes(id)%Domain2
        endif
        exit
     endif
  enddo
end function get_domain2d

!#######################################################################
function get_axis_num (axis_name, set_name) result (num)

! Returns index into axis table corresponding to a given axis name
  character(len=*), intent(in)           :: axis_name
  character(len=*), intent(in), optional ::  set_name
  integer                                :: num, set, n

  if (present(set_name)) then
     set = get_axis_set_num (trim(set_name))
  else
     set = 0
  endif
  num = 0
  do n = 1, num_def_axes
     if ( trim(axis_name) == trim(Axes(n)%name) .and. Axes(n)%set == set ) then
        num = n
        return
     endif
  enddo
end function get_axis_num

!#######################################################################
function get_axis_set_num (set_name) result (num)

! Returns index in axis set table corresponding to a given axis set name

  character(len=*), intent(in) :: set_name
  integer                      :: num, iset

  num = 0
  do iset = 1, num_axis_sets
     if (set_name == Axis_sets(iset))then
        num = iset
        return
     endif
  enddo
end function get_axis_set_num

!#######################################################################
end module diag_axis_mod

max_stars_repo_name: JoyMonteiro/FMSDynamicalCore


!
! Distributed under the OSI-approved Apache License, Version 2.0.  See
!  accompanying file Copyright.txt for details.
!
!  adios2_variable_mod.f90 : ADIOS2 Fortran bindings for Variable class
!
!   Created on: Mar 13, 2017
!       Author: William F Godoy godoywf@ornl.gov
!

module adios2_variable_mod
    use adios2_functions_mod
    implicit none

contains

    subroutine adios2_variable_name(variable, name, ierr)
        type(adios2_variable), intent(in) :: variable
        character(len=:), allocatable, intent(out) :: name
        integer, intent(out) :: ierr

        character(len=1024) :: c_name
        integer :: length

        call adios2_variable_name_f2c(variable%f2c, c_name, length, ierr)
        call adios2_StringC2F(c_name, length, name)

    end subroutine

    subroutine adios2_variable_type(variable, type, ierr)
        type(adios2_variable), intent(in) :: variable
        integer, intent(out) :: type
        integer, intent(out) :: ierr

        integer :: c_type

        call adios2_variable_type_f2c(variable%f2c, c_type, ierr)
        call adios2_TypeC2F(c_type, type)

    end subroutine

    subroutine adios2_variable_ndims(variable, ndims, ierr)
        type(adios2_variable), intent(in) :: variable
        integer, intent(out) :: ndims
        integer, intent(out) :: ierr

        call adios2_variable_ndims_f2c(variable%f2c, ndims, ierr)
    end subroutine

    subroutine adios2_variable_shape(variable, ndims, shape_dims, ierr)
        type(adios2_variable), intent(in) :: variable
        integer, intent(out) :: ndims
        integer(kind=8), dimension(:), allocatable, intent(out) :: shape_dims
        integer, intent(out) :: ierr

        call adios2_variable_ndims_f2c(variable%f2c, ndims, ierr)
        allocate (shape_dims(ndims))
        call adios2_variable_shape_f2c(variable%f2c, shape_dims, ierr)

    end subroutine

    subroutine adios2_variable_steps_start(variable, steps_start, ierr)
        type(adios2_variable), intent(in) :: variable
        integer(kind=8), intent(out) :: steps_start
        integer, intent(out) :: ierr

        call adios2_variable_steps_start_f2c(variable%f2c, steps_start, ierr)

    end subroutine

    subroutine adios2_variable_steps(variable, steps_count, ierr)
        type(adios2_variable), intent(in) :: variable
        integer(kind=8), intent(out) :: steps_count
        integer, intent(out) :: ierr

        call adios2_variable_steps_f2c(variable%f2c, steps_count, ierr)

    end subroutine

    subroutine adios2_set_shape(variable, ndims, shape_dims, ierr)
        type(adios2_variable), intent(in) :: variable
        integer, intent(in) :: ndims
        integer(kind=8), dimension(:), intent(in) :: shape_dims
        integer, intent(out) :: ierr

        call adios2_set_shape_f2c(variable%f2c, ndims, shape_dims, ierr)
    end subroutine

    subroutine adios2_set_selection(variable, ndims, start_dims, count_dims, &
                                    ierr)
        type(adios2_variable), intent(in) :: variable
        integer, intent(in) :: ndims
        integer(kind=8), dimension(:), intent(in) :: start_dims
        integer(kind=8), dimension(:), intent(in) :: count_dims
        integer, intent(out) :: ierr

        call adios2_set_selection_f2c(variable%f2c, ndims, start_dims, &
                                      count_dims, ierr)
    end subroutine

    subroutine adios2_set_step_selection(variable, step_start, step_count, ierr)
        type(adios2_variable), intent(in) :: variable
        integer(kind=8), intent(in) :: step_start
        integer(kind=8), intent(in) :: step_count
        integer, intent(out) :: ierr

        call adios2_set_step_selection_f2c(variable%f2c, step_start, &
                                           step_count, ierr)
    end subroutine

    subroutine adios2_variable_check_type(variable, adios2_type, hint, ierr)
        type(adios2_variable), intent(in):: variable
        integer, intent(in):: adios2_type
        character*(*), intent(in):: hint
        integer, intent(out):: ierr
        !local
        integer :: variable_type
        character(len=:), allocatable :: variable_name

        call adios2_variable_type(variable, variable_type, ierr)

        if( variable_type /= adios2_type ) then
            call adios2_variable_name(variable, variable_name, ierr)
            write(0,*) 'ERROR: adios2 variable ', variable_name, &
                       ' type mismatch, in call to adios2_', hint
            ierr = -1
        end if

    end subroutine

end module


text:
      SUBROUTINE mpmul(w,u,v,n,m)
      INTEGER m,n,NMAX
      CHARACTER*1 w(n+m),u(n),v(m)
      DOUBLE PRECISION RX
      PARAMETER (NMAX=8192,RX=256.D0)
CU    USES drealft
      INTEGER j,mn,nn
      DOUBLE PRECISION cy,t,a(NMAX),b(NMAX)
      mn=max(m,n)
      nn=1
1     if(nn.lt.mn) then
        nn=nn+nn
      goto 1
      endif
      nn=nn+nn
      if(nn.gt.NMAX)pause 'NMAX too small in fftmul'
      do 11 j=1,n
        a(j)=ichar(u(j))
11    continue
      do 12 j=n+1,nn
        a(j)=0.D0
12    continue
      do 13 j=1,m
        b(j)=ichar(v(j))
13    continue
      do 14 j=m+1,nn
        b(j)=0.D0
14    continue
      call drealft(a,nn,1)
      call drealft(b,nn,1)
      b(1)=b(1)*a(1)
      b(2)=b(2)*a(2)
      do 15 j=3,nn,2
        t=b(j)
        b(j)=t*a(j)-b(j+1)*a(j+1)
        b(j+1)=t*a(j+1)+b(j+1)*a(j)
15    continue
      call drealft(b,nn,-1)
      cy=0.
      do 16 j=nn,1,-1
        t=b(j)/(nn/2)+cy+0.5D0
        b(j)=mod(t,RX)
        cy=int(t/RX)
16    continue
      if (cy.ge.RX) pause 'cannot happen in fftmul'
      w(1)=char(int(cy))
      do 17 j=2,n+m
        w(j)=char(int(b(j-1)))
17    continue
      return
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/mpmul.for


text: The Charlotte Pehlke Benches Bench Memorials honors the memory of Charlotte Pehlke.  It can be found just off the Bike Paths path running along the east side of West Pond.  It is in a small garden by the parking lot for the UC Davis Medical Group.


text:
 All About Lucky
  The best thing about Lucky is his ability to make you smile no matter how bad youre feeling. He has catlike intuition and loves to cuddle and sleep all day. While he refuses to play with toys or other dogs, he loves the company of doting women...men not so much.

  Lucky was adopted from an abusive home at about 1 year old, and has come a long way from being a cowardly, fearful animal, to the loving, loyal pet we know today.

  Luckys favourite hiding places include: under the bed, behind the table, inside boxes, in his crate, kitchen cabinets, behind the toilet, on closet shelves...

  If Lucky could, he would likely enjoy spending his whole day behind a couch somewhere with no distractions or interruptionsnot even a break for a snack could tear him away.

  Lucky is a light eater, but when he does sit down for a meal, he requires privacy. He will not eat if someone is watching him, and he prefers to take his food out of his bowl in order to eat in a secluded, unreachable place.

  Lucky loves Users/EricWilliamWong Eric Wong, but just doesnt know it yet. :)

  On 1/31/14, Lucky had surgery to remove a bleeding tumor from his small intestine! Since then, he has made a full recovery (yay, go Lucky!)
 > Luckys Photo Gallery!<

  Image(lick.jpg, left, thumbnail, 500, Tasty?)

  Image(steamy.jpg, thumbnail, 1000, Napping. Photo credit:  Users/JenniferGiang Jennifer Giang)

  Image(burried.jpg, thumbnail, 500, Buried in the Sand!)

  Image(lucky vet.jpg, left, thumbnail, 500, PostSurgery Vet Visit!)

  Image(worker dog.jpg, thumbnail, 500)

  Image(sleepy dog.jpg, thumbnail, 500, Exhausted from hiking!)


text:
     subroutine testConversion (messageUnit, status)
!
!
!      USE IFPORT
!      USE IFCORE
!
      implicit none
!

      integer messageUnit, status, ierr
      integer highFunction, lowFunction, errorCode, istatus, severity
       character c*5

      common /lchk/ lcheck
      logical lcheck
!
!
      integer nrec
      integer ifltabSize
      integer(8) ifltab1(600), ifltab2(600)
      ifltabSize = 600
!
!
     write(messageUnit, 10)
 10  format('*** Enter testConversion ***')
     lcheck = .true.
 !   Set the output message level to moderate
    ! call zset('mlvl', ' ', 4)
   !  call zopen(ifltab1, 'C:/temp/sample6.dss', status)
   !  ierr = DELFILESQQ ('C:/temp/sample7convert.dss')
   !  call zopen7(ifltab2, 'C:/temp/sample7convert.dss', status)

     !call zset('mlvl', ' ', 2)
     call zgetfileversion('sample.dss', status)
     if (status.ne.6) then
            write(messageUnit,*)' '
            write(messageUnit,*)'sample.dss is required for conversion test '
            write(messageUnit,*)'It must be the sample.dss file from DSSVue '
            write(messageUnit,*)'And must be a version 6 file'
            if (status.eq.0) then
                write(messageUnit,*)'That file does not exist'
            else if (status.eq.7) then
                write(messageUnit,*)'That file is a version 7 file'
            endif
            write(messageUnit,*)' '
            return
     endif

     call zopen6(ifltab1, 'sample.dss', status)
     if (status.ne.0) return
     c = ' '
     call zinqir(ifltab1, 'NREC', c, nrec)
     if (nrec.eq.0) then
        write(messageUnit,*)' '
        write(messageUnit,*)' First file for converion test is not correct.'
        write(messageUnit,*)' It should be the sample dss 6 file.'
        status = -1
        return
     endif
     call cdelet ('sample7.dss',ierr)
     call zopen7(ifltab2, 'sample7.dss', status)
     call zcopyfile(ifltab1, ifltab2, status)
     if (status.eq.123) then
        call zclose(ifltab1)
        call CDELET ('sample6convert.dss',ierr)
        call zopen6(ifltab1, 'sample6convert.dss', status)
        call zcopyfile(ifltab2, ifltab1, status)
     endif

!
!
 800    continue
        if (status .ne. 0) go to 900
        call zclose(ifltab1)
        call zclose(ifltab2)
        Write(messageUnit,*)'test conversion  passed successfully!'

        return
!
!       Error return
 900    Continue
        Write(messageUnit, 820) status
        !Write(messageUnit, 811) filename1
        !Write(messageUnit, 811) filename2
 811    Format( 'File tested: ', a)
 820    Format ('test conversion FAILED; error: ', I8)
        if (status .gt. 1000) then
            call zerrordecode(status, highFunction, lowFunction, errorCode, istatus, severity)
            write (*,*) 'highFunction   ',highFunction
            write (*,*) 'lowFunction   ',lowFunction
            write (*,*) 'errorCode   ',errorCode
            write (*,*) 'istatus   ',istatus
            write (*,*) 'severity   ',severity
        endif
      return
      end

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


text:
  function checkIndex( dataObj, i ) result(closestBound)

! ==============================================================================
!
! This function returns an index for a valid entry in the target object arrays.
! These arrays include radii and normalized radii, nucleon densities, and nucleon
! Fermi momenta.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: int32

    implicit none
    class(StandardDCMData), intent(inout) :: dataObj
    integer(int32), intent(in   ) :: i
    integer(int32)                :: closestBound

! ==============================================================================

    ! Assume closest bound is what was passed in
    closestBound = i

    ! Check if the element can exist, if not, move to the closest existing bound
    if ( i < minZonesAllowed ) then
       write(dataObj%io%message, 1000) i, minZonesAllowed
       call dataObj%io%print(0, 3, dataObj%io%message)
       closestBound = minZonesAllowed
    else if ( i > maxZonesAllowed ) then
       write(dataObj%io%message, 1000) i, maxZonesAllowed
       call dataObj%io%print(0, 3, dataObj%io%message)
       closestBound = maxZonesAllowed
    end if

    return
! ==============================================================================
1000 format("Standard DCM data array exceeded - out of bounds (element ", i3, &
          & "). Using closest valid index (", i3, ").")
! ==============================================================================
  end function checkIndex

  function zoneBoundR( dataObj, i ) result(boundR)

! ==============================================================================
!
! This function returns to the client the value of the 'zoneBoundR' variable.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64

    implicit none
    class(StandardDCMData), intent(inout) :: dataObj
    integer(int32),         intent(in   ) :: i

    integer(int32) :: indx
    real(real64)   :: boundR

! ==============================================================================

    indx = dataObj%checkIndex( i )

    boundR = dataObj%target%zoneBoundR( indx )

    return
! ==============================================================================
  end function zoneBoundR

  function zoneBRFrac( dataObj, i ) result(boundRF)

! ==============================================================================
!
! This function returns to the client the value of the 'zoneBRFrac' variable.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64

    implicit none
    class(StandardDCMData), intent(inout) :: dataObj
    integer(int32),         intent(in   ) :: i

    integer(int32) :: indx
    real(real64)   :: boundRF

! ==============================================================================

    indx = dataObj%checkIndex( i )

    boundRF = dataObj%target%zoneBRFrac( indx )

    return
! ==============================================================================
  end function zoneBRFrac

  function protonDensity( dataObj, i ) result(density)

! ==============================================================================
!
! This function returns to the client the value of the 'protonDensity' variable.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64

    implicit none
    class(StandardDCMData), intent(inout) :: dataObj
    integer(int32),         intent(in   ) :: i

    integer(int32) :: indx
    real(real64)   :: density

! ==============================================================================

    indx = dataObj%checkIndex( i )

    density = dataObj%target%protonDensity( indx )

    return
! ==============================================================================
  end function protonDensity

  function neutronDensity( dataObj, i ) result(density)

! ==============================================================================
!
! This function returns to the client the value of the 'neutronDensity' variable.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64

    implicit none
    class(StandardDCMData), intent(inout) :: dataObj
    integer(int32),         intent(in   ) :: i

    integer(int32) :: indx
    real(real64)   :: density

! ==============================================================================

    indx = dataObj%checkIndex( i )

    density = dataObj%target%neutronDensity( indx )

    return
! ==============================================================================
  end function neutronDensity

  function coulombPote( dataObj, i ) result(coulV)

! ==============================================================================
!
! This function returns to the client the value of the 'coulombPote' variable.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64

    implicit none
    class(StandardDCMData), intent(inout) :: dataObj
    integer(int32),         intent(in   ) :: i

    integer(int32) :: indx
    real(real64)   :: coulV

! ==============================================================================

    indx = dataObj%checkIndex( i )

    coulV = dataObj%target%coulombPote( indx )

    return
! ==============================================================================
  end function coulombPote

  function protFermiMom( dataObj, i ) result(fProtMom)

! ==============================================================================
!
! This function returns to the client the value of the 'protFermiMom' variable.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64

    implicit none
    class(StandardDCMData), intent(inout) :: dataObj
    integer(int32),         intent(in   ) :: i

    integer(int32) :: indx
    real(real64)   :: fProtMom

! ==============================================================================

    indx = dataObj%checkIndex( i )

    fProtMom = dataObj%target%protFermiMom( indx )

    return
! ==============================================================================
  end function protFermiMom

  function neutFermiMom( dataObj, i ) result(fNeutMom)

! ==============================================================================
!
! This function returns to the client the value of the 'neutFermiMom' variable.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64

    implicit none
    class(StandardDCMData), intent(inout) :: dataObj
    integer(int32),         intent(in   ) :: i

    integer(int32) :: indx
    real(real64)   :: fNeutMom

! ==============================================================================

    indx = dataObj%checkIndex( i )

    fNeutMom = dataObj%target%neutFermiMom( indx )

    return
! ==============================================================================
  end function neutFermiMom

  function numBaryons ( dataObj ) result(aNumber)

! ==============================================================================
!
! This function returns the number of baryons contained by the target that the
! Standard DCM Data object is established for.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: real64

    implicit none
    class(StandardDCMData), intent(in   ) :: dataObj
    real(real64) :: aNumber

! ==============================================================================

    aNumber = dataObj%target%numBaryons
    return
! ==============================================================================
  end function numBaryons

  function numProtons ( dataObj ) result(zNumber)

! ==============================================================================
!
! This function returns the number of protons contained by the target that the
! Standard DCM Data object is established for.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: real64

    implicit none
    class(StandardDCMData), intent(in   ) :: dataObj
    real(real64) :: zNumber

! ==============================================================================

    zNumber = dataObj%target%numProtons
    return
! ==============================================================================
  end function numProtons

  function numZones ( dataObj ) result(nZones)

! ==============================================================================
!
! This function returns the number of zones contained by the target that the
! Standard DCM Data object is established for.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: int32

    implicit none
    class(StandardDCMData), intent(in   ) :: dataObj
    integer(int32) :: nZones

! ==============================================================================

    nZones = dataObj%options%numZones
    return
! ==============================================================================
  end function numZones

  function aTargThrd ( dataObj ) result(athrd)

! ==============================================================================
!
! This function returns the aTarg**(1/3) for the Standard DCM Data object.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: real64

    implicit none
    class(StandardDCMData), intent(in   ) :: dataObj
    real(real64)                          :: athrd

! ==============================================================================

    athrd = dataObj%target%aTargThrd
    return
! ==============================================================================
  end function aTargThrd

  function geomCrossSection ( dataObj ) result(geomSigma)

! ==============================================================================
!
! This function returns the geometric cross section of the target that the
! Standard DCM Data object is established for.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: real64

    implicit none
    class(StandardDCMData), intent(in   ) :: dataObj
    real(real64)                          :: geomSigma

! ==============================================================================

    geomSigma = dataObj%target%geomCrossSection
    return
! ==============================================================================
  end function geomCrossSection

  function pionPote ( dataObj ) result(pionPotential)

! ==============================================================================
!
! This function returns the geometric cross section of the target that the
! Standard DCM Data object is established for.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: real64

    implicit none
    class(StandardDCMData), intent(in   ) :: dataObj
    real(real64)                          :: pionPotential

! ==============================================================================

    pionPotential = dataObj%target%pionPote
    return
! ==============================================================================
  end function pionPote

  function getSepEnergy ( dataObj ) result(sepEnergy)

! ==============================================================================
!
! This function returns the geometric cross section of the target that the
! Standard DCM Data object is established for.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: real64

    implicit none
    class(StandardDCMData), intent(in   ) :: dataObj
    real(real64)                          :: sepEnergy

! ==============================================================================

    sepEnergy = dataObj%target%sepEnergy
    return
! ==============================================================================
  end function getSepEnergy

  subroutine setSepEnergy( dataObj, newSepEnergy )

! ==============================================================================
!
! This function returns the geometric cross section of the target that the
! Standard DCM Data object is established for.
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use, intrinsic:: iso_fortran_env, only: real64

    implicit none
    class(StandardDCMData), intent(inout) :: dataObj
    real(real64),           intent(in   ) :: newSepEnergy

! ==============================================================================

    dataObj%target%sepEnergy = newSepEnergy
    return
! ==============================================================================
  end subroutine setSepEnergy

  function properlyConstructedObject(dataObj) result(constructed)

! ==============================================================================
!
! Returns flag stating if the data object is constructed or not
!
! ==============================================================================

    implicit none
    class(StandardDCMData), intent(in   ) :: dataObj
    logical :: constructed

! ==============================================================================

    constructed = dataObj%objectConstructed

    return
! ==============================================================================
  end function properlyConstructedObject

  function queryOptions( dataObj ) result(options)

! ==============================================================================
!
! Returns the options contained by the sDCM Data object
!
! ==============================================================================

    implicit none
    class(StandardDCMData), intent(in   ) :: dataObj
    type(sDCMDataOptions) :: options

! ==============================================================================

    options = dataObj%options
    return
! ==============================================================================
  end function queryOptions

max_stars_repo_name: lanl/generalized-spallation-model

max_issues_repo_name: lanl/generalized-spallation-model

max_forks_repo_name: lanl/generalized-spallation-model