text:
!##############################################################################
!# ****************************************************************************
!# <name> stokes2d_method3_block </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module is a demonstation program how to solve a Stokes
!# problem on a simple domain.
!#
!# The routine uses the BiCGStab solver with a simple-VANKA preconditioner
!# for 2D saddle point problems, Jacobi-Type.
!# The discretisation uses the block assembly method to evaluate the
!# mattrix all-in-one.
!# </purpose>
!##############################################################################

module stokes2d_method3_block

  use fsystem
  use storage
  use genoutput
  use boundary
  use cubature
  use derivatives
  use matrixfilters
  use vectorfilters
  use linearalgebra
  use discretebc
  use bcassembly
  use triangulation
  use element
  use spatialdiscretisation
  use linearsystemscalar
  use linearsystemblock
  use spdiscprojection
  use scalarpde
  use bilinearformevaluation
  use linearformevaluation
  use discretebc
  use filtersupport
  use coarsegridcorrection
  use linearsolver
  use ucd

  use stokes2d_callback

  use blockmatassemblybase
  use blockmatassembly
  use feevaluation2
  use collection

  implicit none

contains

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

!<subroutine>

  subroutine st2d3_fcalc_Stokes(RmatrixData,rassemblyData,rmatrixAssembly,&
      npointsPerElement,nelements,revalVectors,rcollection)

!<description>
    ! Calculates the Stokes operator.
!</description>

!<inputoutput>
    ! Matrix data of all matrices. The arrays p_Dentry of all submatrices
    ! have to be filled with data.
    type(t_bmaMatrixData), dimension(:,:), intent(inout), target :: RmatrixData
!</inputoutput>

!<input>
    ! Data necessary for the assembly. Contains determinants and
    ! cubature weights for the cubature,...
    type(t_bmaMatrixAssemblyData), intent(in) :: rassemblyData

    ! Structure with all data about the assembly
    type(t_bmaMatrixAssembly), intent(in) :: rmatrixAssembly

    ! Number of points per element
    integer, intent(in) :: npointsPerElement

    ! Number of elements
    integer, intent(in) :: nelements

    ! Values of FEM functions automatically evaluated in the
    ! cubature points.
    type(t_fev2Vectors), intent(in) :: revalVectors

    ! User defined collection structure
    type(t_collection), intent(inout), target, optional :: rcollection
!</input>

!<subroutine>

    ! Local variables
    real(DP) :: dbasJ, dval, dbasIx, dbasIy, dbasJx, dbasJy, dnu
    integer :: iel, icubp, idofe, jdofe
    real(DP), dimension(:,:,:), pointer :: p_DlocalMatrixA11,p_DlocalMatrixA22
    real(DP), dimension(:,:,:), pointer :: p_DlocalMatrixA13,p_DlocalMatrixA23
    real(DP), dimension(:,:,:), pointer :: p_DlocalMatrixA31,p_DlocalMatrixA32
    real(DP), dimension(:,:,:,:), pointer :: p_DbasTrialA11,p_DbasTestA11
    real(DP), dimension(:,:,:,:), pointer :: p_DbasTrialA13,p_DbasTestA13
    real(DP), dimension(:,:), pointer :: p_DcubWeight
    type(t_bmaMatrixData), pointer :: p_rmatrixDataA11,p_rmatrixDataA22,p_rmatrixDataA13

    ! Viscosity
    dnu = rcollection%DquickAccess(1)

    ! Get cubature weights data
    p_DcubWeight => rassemblyData%p_DcubWeight
    p_rmatrixDataA11 => RmatrixData(1,1)
    p_rmatrixDataA22 => RmatrixData(2,2)
    p_rmatrixDataA13 => RmatrixData(1,3)
    p_DlocalMatrixA11 => RmatrixData(1,1)%p_Dentry
    p_DlocalMatrixA22 => RmatrixData(2,2)%p_Dentry
    p_DlocalMatrixA13 => RmatrixData(1,3)%p_Dentry
    p_DlocalMatrixA23 => RmatrixData(2,3)%p_Dentry
    p_DlocalMatrixA31 => RmatrixData(3,1)%p_Dentry
    p_DlocalMatrixA32 => RmatrixData(3,2)%p_Dentry
    p_DbasTrialA11 => RmatrixData(1,1)%p_DbasTrial
    p_DbasTestA11 => RmatrixData(1,1)%p_DbasTest
    p_DbasTrialA13 => RmatrixData(1,3)%p_DbasTrial
    p_DbasTestA13 => RmatrixData(1,3)%p_DbasTest

    ! Calculate Laplace in the diagonal blocks

    ! Loop over the elements in the current set.
    do iel = 1,nelements

      ! Loop over all cubature points on the current element
      do icubp = 1,npointsPerElement

        ! Outer loop over the DOFs i=1..ndof on our current element,
        ! which corresponds to the (test) basis functions Phi_i:
        do idofe=1,p_rmatrixDataA11%ndofTest

          ! Fetch the contributions of the (test) basis functions Phi_i
          ! into dbasI
          dbasIx = p_DbasTestA11(idofe,DER_DERIV2D_X,icubp,iel)
          dbasIy = p_DbasTestA11(idofe,DER_DERIV2D_Y,icubp,iel)

          ! Inner loop over the DOFs j=1..ndof, which corresponds to
          ! the basis function Phi_j:
          do jdofe=1,p_rmatrixDataA11%ndofTrial

            ! Fetch the contributions of the (trial) basis function Phi_j
            ! into dbasJ
            dbasJx = p_DbasTrialA11(jdofe,DER_DERIV2D_X,icubp,iel)
            dbasJy = p_DbasTrialA11(jdofe,DER_DERIV2D_Y,icubp,iel)

            ! Multiply the values of the basis functions
            ! (1st derivatives) by the cubature weight and sum up
            ! into the local matrices.
            dval = dnu * p_DcubWeight(icubp,iel) * ( dbasJx*dbasIx + dbasJy*dbasIy )
            p_DlocalMatrixA11(jdofe,idofe,iel) = p_DlocalMatrixA11(jdofe,idofe,iel) + dval
            p_DlocalMatrixA22(jdofe,idofe,iel) = p_DlocalMatrixA22(jdofe,idofe,iel) + dval

          end do ! idofe

        end do ! jdofe

      end do ! icubp

    end do ! iel

    ! Calculate the B1/B2/B1^T/B2^T blocks
    ! Loop over the elements in the current set.
    do iel = 1,nelements

      ! Loop over all cubature points on the current element
      do icubp = 1,npointsPerElement

        ! Outer loop over the DOFs i=1..ndof on our current element,
        ! which corresponds to the (test) basis functions Phi_i:
        do idofe=1,p_rmatrixDataA13%ndofTest

          ! Fetch the contributions of the (test) basis functions Phi_i
          ! into dbasI
          dbasIx = p_DbasTestA13(idofe,DER_DERIV2D_X,icubp,iel)
          dbasIy = p_DbasTestA13(idofe,DER_DERIV2D_Y,icubp,iel)

          ! Inner loop over the DOFs j=1..ndof, which corresponds to
          ! the basis function Phi_j:
          do jdofe=1,p_rmatrixDataA13%ndofTrial

            ! Fetch the contributions of the (trial) basis function Phi_j
            ! into dbasJ
            dbasJ = p_DbasTrialA13(jdofe,DER_FUNC,icubp,iel)

            ! Multiply the values of the basis functions
            ! (1st derivatives) by the cubature weight and sum up
            ! into the local matrices.

            ! B-matrices
            p_DlocalMatrixA13(jdofe,idofe,iel) = p_DlocalMatrixA13(jdofe,idofe,iel) + &
                p_DcubWeight(icubp,iel) * ( -dbasJ*dbasIx )

            p_DlocalMatrixA23(jdofe,idofe,iel) = p_DlocalMatrixA23(jdofe,idofe,iel) + &
                p_DcubWeight(icubp,iel) * ( -dbasJ*dbasIy )

            ! B^T-matrices
            p_DlocalMatrixA31(idofe,jdofe,iel) = p_DlocalMatrixA31(idofe,jdofe,iel) + &
                p_DcubWeight(icubp,iel) * ( -dbasJ*dbasIx )

            p_DlocalMatrixA32(idofe,jdofe,iel) = p_DlocalMatrixA32(idofe,jdofe,iel) + &
                p_DcubWeight(icubp,iel) * ( -dbasJ*dbasIy )

          end do ! idofe

        end do ! jdofe

      end do ! icubp

    end do ! iel

  end subroutine

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

!<subroutine>

  subroutine st2d3_fcalc_rhs(rvectorData,rassemblyData,rvectorAssembly,&
      npointsPerElement,nelements,revalVectors,rcollection)

!<description>
    ! Calculates the Mass operator in all diagonal matrices.
!</description>

!<inputoutput>
    ! Vector data of all subvectors. The arrays p_Dentry of all subvectors
    ! have to be filled with data.
    type(t_bmaVectorData), dimension(:), intent(inout), target :: rvectorData
!</inputoutput>

!<input>
    ! Data necessary for the assembly. Contains determinants and
    ! cubature weights for the cubature,...
    type(t_bmaVectorAssemblyData), intent(in) :: rassemblyData

    ! Structure with all data about the assembly
    type(t_bmaVectorAssembly), intent(in) :: rvectorAssembly

    ! Number of points per element
    integer, intent(in) :: npointsPerElement

    ! Number of elements
    integer, intent(in) :: nelements

    ! Values of FEM functions automatically evaluated in the
    ! cubature points.
    type(t_fev2Vectors), intent(in) :: revalVectors

    ! User defined collection structure
    type(t_collection), intent(inout), target, optional :: rcollection
!</input>

!<subroutine>

    ! Local variables
    real(DP) :: dbasI, dval1, dval2, dval3
    integer :: iel, icubp, idofe
    real(DP), dimension(:,:), pointer :: p_DlocalVector1,p_DlocalVector2,p_DlocalVector3
    real(DP), dimension(:,:,:,:), pointer :: p_DbasTest1,p_DbasTest3
    real(DP), dimension(:,:), pointer :: p_DcubWeight
    type(t_bmaVectorData), pointer :: p_rvectorData1,p_rvectorData2,p_rvectorData3

    ! Get cubature weights data
    p_DcubWeight => rassemblyData%p_DcubWeight
    p_rvectorData1 => RvectorData(1)
    p_rvectorData2 => RvectorData(2)
    p_rvectorData3 => RvectorData(3)
    p_DlocalVector1 => RvectorData(1)%p_Dentry
    p_DlocalVector2 => RvectorData(2)%p_Dentry
    p_DlocalVector3 => RvectorData(3)%p_Dentry
    p_DbasTest1 => RvectorData(1)%p_DbasTest
    p_DbasTest3 => RvectorData(3)%p_DbasTest

    ! Calculate the RHS of the momentum equation

    ! Loop over the elements in the current set.
    do iel = 1,nelements

      ! Loop over all cubature points on the current element
      do icubp = 1,npointsPerElement

        ! Outer loop over the DOFs i=1..ndof on our current element,
        ! which corresponds to the (test) basis functions Phi_i:
        do idofe=1,p_rvectorData1%ndofTest

          ! Fetch the contributions of the (test) basis functions Phi_i
          ! into dbasI
          dbasI = p_DbasTest1(idofe,DER_FUNC,icubp,iel)

          ! Values of the velocity RHS for the X1 and X2 component
          dval1 = 0.0_DP
          dval2 = 0.0_DP

          ! Multiply the values of the basis functions
          ! (1st derivatives) by the cubature weight and sum up
          ! into the local vectors.
          p_DlocalVector1(idofe,iel) = p_DlocalVector1(idofe,iel) + &
              p_DcubWeight(icubp,iel) * dval1 * dbasI
          p_DlocalVector2(idofe,iel) = p_DlocalVector2(idofe,iel) + &
              p_DcubWeight(icubp,iel) * dval2 * dbasI

        end do ! jdofe

      end do ! icubp

    end do ! iel

    ! Calculate the RHS of the continuity equation

    ! Loop over the elements in the current set.
    do iel = 1,nelements

      ! Loop over all cubature points on the current element
      do icubp = 1,npointsPerElement

        ! Outer loop over the DOFs i=1..ndof on our current element,
        ! which corresponds to the (test) basis functions Phi_i:
        do idofe=1,p_rvectorData3%ndofTest

          ! Fetch the contributions of the (test) basis functions Phi_i
          ! into dbasI
          dbasI = p_DbasTest3(idofe,DER_FUNC,icubp,iel)

          ! Value of the pressure RHS
          dval3 = 0.0_DP

          ! Multiply the values of the basis functions
          ! by the cubature weight and sum up into the local vectors.
          p_DlocalVector3(idofe,iel) = p_DlocalVector3(idofe,iel) + &
              p_DcubWeight(icubp,iel) * dval3 * dbasI

        end do ! jdofe

      end do ! icubp

    end do ! iel

  end subroutine

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

!<subroutine>

  subroutine st2d3_fcalc_L2error(Dintvalue,rassemblyData,rintegralAssembly,&
      npointsPerElement,nelements,revalVectors,rcollection)

!<description>
    ! Calculates the (squared) L2 error.
    ! The component to calculate the error for has to be provided
    ! in rcollection%IquickAccess(1). The viscosity has to be stored
    ! in rcollection%DquickAccess(1).
!</description>

!<input>
    ! Data necessary for the assembly. Contains determinants and
    ! cubature weights for the cubature,...
    type(t_bmaIntegralAssemblyData), intent(in) :: rassemblyData

    ! Structure with all data about the assembly
    type(t_bmaIntegralAssembly), intent(in) :: rintegralAssembly

    ! Number of points per element
    integer, intent(in) :: npointsPerElement

    ! Number of elements
    integer, intent(in) :: nelements

    ! Values of FEM functions automatically evaluated in the
    ! cubature points.
    type(t_fev2Vectors), intent(in) :: revalVectors

    ! User defined collection structure
    type(t_collection), intent(inout), target, optional :: rcollection
!</input>

!<output>
    ! Returns the value of the integral
    real(DP), dimension(:), intent(out) :: Dintvalue
!</output>

!<subroutine>

    ! Local variables
    real(DP) :: dval1,dval2,dx,dy,dnu
    integer :: iel, icubp, icomp
    real(DP), dimension(:,:), pointer :: p_DcubWeight
    real(DP), dimension(:,:), pointer :: p_Dfunc
    real(DP), dimension(:,:,:), pointer :: p_Dpoints

    ! Calcel if no FEM function is given.
    if (revalVectors%ncount .eq. 0) return

    ! Skip interleaved vectors.
    if (revalVectors%p_RvectorData(1)%bisInterleaved) return

    ! Get cubature weights data
    p_DcubWeight => rassemblyData%p_DcubWeight

    ! Get the coordinates of the cubature points
    p_Dpoints => rassemblyData%revalElementSet%p_DpointsReal

    ! Get the component from the collection
    icomp = rcollection%IquickAccess(1)

    ! Get the viscosity parameter
    dnu = rcollection%DquickAccess(1)

    ! Get the data array with the values of the FEM function
    ! in the cubature points
    p_Dfunc => revalVectors%p_RvectorData(icomp)%p_Ddata(:,:,DER_FUNC2D)

    Dintvalue = 0.0_DP

    ! Which component to calculate?
    select case (icomp)

    ! ----------
    ! X-velocity
    ! ----------
    case (1)

      ! Loop over the elements in the current set.
      do iel = 1,nelements

        ! Loop over all cubature points on the current element
        do icubp = 1,npointsPerElement

          ! Get the value of the bubble function
          dx = p_Dpoints(1,icubp,iel)
          dy = p_Dpoints(2,icubp,iel)

          dval1 = dy * (1.0_DP-dy)

          ! Get the error of the FEM function to the bubble function
          dval2 = p_Dfunc(icubp,iel)

          ! Multiply the values by the cubature weight and sum up
          ! into the (squared) L2 error:
          Dintvalue = Dintvalue + &
              p_DcubWeight(icubp,iel) * (dval1 - dval2)**2

        end do ! icubp

      end do ! iel

    ! ----------
    ! Y-velocity
    ! ----------
    case (2)

      ! Loop over the elements in the current set.
      do iel = 1,nelements

        ! Loop over all cubature points on the current element
        do icubp = 1,npointsPerElement

          ! Get the value of the bubble function
          dx = p_Dpoints(1,icubp,iel)
          dy = p_Dpoints(2,icubp,iel)

          dval1 = 0.0_DP

          ! Get the error of the FEM function to the bubble function
          dval2 = p_Dfunc(icubp,iel)

          ! Multiply the values by the cubature weight and sum up
          ! into the (squared) L2 error:
          Dintvalue = Dintvalue + &
              p_DcubWeight(icubp,iel) * (dval1 - dval2)**2

        end do ! icubp

      end do ! iel

    ! ----------
    ! Pressure
    ! ----------
    case (3)
      ! Loop over the elements in the current set.
      do iel = 1,nelements

        ! Loop over all cubature points on the current element
        do icubp = 1,npointsPerElement

          ! Get the value of the bubble function
          dx = p_Dpoints(1,icubp,iel)
          dy = p_Dpoints(2,icubp,iel)

          dval1 = 2.0_DP * dnu * (1.0_DP-dx)

          ! Get the error of the FEM function to the bubble function
          dval2 = p_Dfunc(icubp,iel)

          ! Multiply the values by the cubature weight and sum up
          ! into the (squared) L2 error:
          Dintvalue = Dintvalue + &
              p_DcubWeight(icubp,iel) * (dval1 - dval2)**2

        end do ! icubp

      end do ! iel

    end select

  end subroutine

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

!<subroutine>

  subroutine st2d3_fcalc_H1error(Dintvalue,rassemblyData,rintegralAssembly,&
      npointsPerElement,nelements,revalVectors,rcollection)

!<description>
    ! Calculates the (squared) H1 error.
    ! The component to calculate the error for has to be provided
    ! in rcollection%IquickAccess(1). The viscosity has to be stored
    ! in rcollection%DquickAccess(1).
!</description>

!<input>
    ! Data necessary for the assembly. Contains determinants and
    ! cubature weights for the cubature,...
    type(t_bmaIntegralAssemblyData), intent(in) :: rassemblyData

    ! Structure with all data about the assembly
    type(t_bmaIntegralAssembly), intent(in) :: rintegralAssembly

    ! Number of points per element
    integer, intent(in) :: npointsPerElement

    ! Number of elements
    integer, intent(in) :: nelements

    ! Values of FEM functions automatically evaluated in the
    ! cubature points.
    type(t_fev2Vectors), intent(in) :: revalVectors

    ! User defined collection structure
    type(t_collection), intent(inout), target, optional :: rcollection
!</input>

!<output>
    ! Returns the value of the integral
    real(DP), dimension(:), intent(out) :: Dintvalue
!</output>

!<subroutine>

    ! Local variables
    real(DP) :: dderivX1,dderivY1,dderivX2,dderivY2,dx,dy
    integer :: iel, icubp, icomp
    real(DP), dimension(:,:), pointer :: p_DcubWeight
    real(DP), dimension(:,:), pointer :: p_DderivX,p_DderivY
    real(DP), dimension(:,:,:), pointer :: p_Dpoints

    ! Calcel if no FEM function is given.
    if (revalVectors%ncount .eq. 0) return

    ! Skip interleaved vectors.
    if (revalVectors%p_RvectorData(1)%bisInterleaved) return

    ! Get cubature weights data
    p_DcubWeight => rassemblyData%p_DcubWeight

    ! Get the coordinates of the cubature points
    p_Dpoints => rassemblyData%revalElementSet%p_DpointsReal

    ! Get the component from the collection
    icomp = rcollection%IquickAccess(1)

    ! Get the data array with the values of the FEM function
    ! in the cubature points
    p_DderivX => revalVectors%p_RvectorData(icomp)%p_Ddata(:,:,DER_DERIV2D_X)
    p_DderivY => revalVectors%p_RvectorData(icomp)%p_Ddata(:,:,DER_DERIV2D_Y)

    Dintvalue = 0.0_DP

    ! Which component to calculate?
    select case (icomp)

    ! ----------
    ! X-velocity
    ! ----------
    case (1)

      ! Loop over the elements in the current set.
      do iel = 1,nelements

        ! Loop over all cubature points on the current element
        do icubp = 1,npointsPerElement

          ! Get the derivatives of the bubble function in the cubature point
          dx = p_Dpoints(1,icubp,iel)
          dy = p_Dpoints(2,icubp,iel)

          dderivX1 = 0.0_DP
          dderivY1 = 1.0_DP-2.0_DP*dy

          ! Get the error of the FEM function derivatives of the bubble function
          ! in the cubature point
          dderivX2 = p_DderivX(icubp,iel)
          dderivY2 = p_DderivY(icubp,iel)

          ! Multiply the values by the cubature weight and sum up
          ! into the (squared) H1 error:
          Dintvalue = Dintvalue + p_DcubWeight(icubp,iel) * &
              ( (dderivX1 - dderivX2)**2 + (dderivY1 - dderivY2)**2 )

        end do ! icubp

      end do ! iel

    ! ----------
    ! Y-velocity
    ! ----------
    case (2)

      ! Loop over the elements in the current set.
      do iel = 1,nelements

        ! Loop over all cubature points on the current element
        do icubp = 1,npointsPerElement

          ! Get the derivatives of the bubble function in the cubature point
          dx = p_Dpoints(1,icubp,iel)
          dy = p_Dpoints(2,icubp,iel)

          dderivX1 = 0.0_DP
          dderivY1 = 0.0_DP

          ! Get the error of the FEM function derivatives of the bubble function
          ! in the cubature point
          dderivX2 = p_DderivX(icubp,iel)
          dderivY2 = p_DderivY(icubp,iel)

          ! Multiply the values by the cubature weight and sum up
          ! into the (squared) H1 error:
          Dintvalue = Dintvalue + p_DcubWeight(icubp,iel) * &
              ( (dderivX1 - dderivX2)**2 + (dderivY1 - dderivY2)**2 )

        end do ! icubp

      end do ! iel

    ! ----------
    ! Pressure
    ! ----------
    case (3)

      ! not used here.

    end select

  end subroutine

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

!<subroutine>

  subroutine stokes2d_3_block

!<description>
  ! This is an all-in-one stokes solver for directly solving a stokes
  ! problem without making use of special features like collections
  ! and so on. The routine performs the following tasks:
  !
  ! 1.) Read in parametrisation
  ! 2.) Read in triangulation
  ! 3.) Set up RHS
  ! 4.) Set up matrix
  ! 5.) Create solver structure
  ! 6.) Solve the problem
  ! 7.) Write solution to GMV file
  ! 8.) Release all variables, finish
!</description>

!</subroutine>

    ! Definitions of variables.
    !
    ! We need a couple of variables for this problem. Let us see...
    !
    ! An object for saving the domain:
    type(t_boundary) :: rboundary

    ! An object for saving the triangulation on the domain
    type(t_triangulation) :: rtriangulation

    ! An object specifying the discretisation.
    ! This contains also information about trial/test functions,...
    type(t_blockDiscretisation) :: rdiscretisation,rprjDiscretisation

    ! A block matrix and a couple of block vectors. These will be filled
    ! with data for the linear solver.
    type(t_matrixBlock) :: rmatrix
    type(t_vectorBlock) :: rvector,rrhs,rtempBlock,rprjVector

    ! Cubature information structure which defines the cubature formula.
    type(t_scalarCubatureInfo) :: rcubatureInfo

    ! A set of variables describing the analytic and discrete boundary
    ! conditions.
    type(t_boundaryRegion) :: rboundaryRegion
    type(t_discreteBC), target :: rdiscreteBC, rprjDiscreteBC

    ! A solver node that accepts parameters for the linear solver
    type(t_linsolNode), pointer :: p_rsolverNode,p_rpreconditioner

    ! A filter chain that describes how to filter the matrix/vector
    ! before/during the solution process. The filters usually implement
    ! boundary conditions.
    type(t_filterChain), dimension(1), target :: RfilterChain

    ! Number of filters in the filter chain
    integer :: nfilters

    ! NLMAX receives the level where we want to solve.
    integer :: NLMAX

    ! Viscosity parameter nu = 1/Re
    real(DP) :: dnu

    ! Error indicator during initialisation of the solver
    integer :: ierror

    ! Collection structure for callback routines
    type(t_collection) :: rcollection

    ! Output block for UCD output to GMV file
    type(t_ucdExport) :: rexport
    character(len=SYS_STRLEN) :: sucddir
    real(DP), dimension(:), pointer :: p_Ddata,p_Ddata2

    ! Path to the mesh
    character(len=SYS_STRLEN) :: spredir

    ! Vector evaluation structure for the calculation of errors.
    type(t_fev2Vectors) :: revalVectors

    ! Error arrays for post-processing
    real(DP), dimension(2) :: DerrorUL2, DerrorUH1
    real(DP), dimension(1) :: DerrorPL2

    ! Ok, let us start.
    !
    ! We want to solve our Poisson problem on level...
    ! As we do not use a multigrid solver here, we will set the level
    ! to 5 instead of 7.
    NLMAX = 5

    ! Viscosity parameter:
    dnu = 1.0_DP

    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Read the domain, read the mesh, refine
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    ! Get the path $PREDIR from the environment, where to read .prm/.tri files
    ! from. If that does not exist, write to the directory "./pre".
    if (.not. sys_getenv_string("PREDIR", spredir)) spredir = "./pre"

    ! At first, read in the parametrisation of the boundary and save
    ! it to rboundary.
    call boundary_read_prm(rboundary, trim(spredir)//"/QUAD.prm")

    ! Now read in the basic triangulation.
    call tria_readTriFile2D (rtriangulation, trim(spredir)//"/QUAD.tri", rboundary)

    ! Refine the mesh up to the minimum level
    call tria_quickRefine2LevelOrdering(NLMAX-1,rtriangulation,rboundary)

    ! Create information about adjacencies and everything one needs from
    ! a triangulation. Afterwards, we have the coarse mesh.
    call tria_initStandardMeshFromRaw (rtriangulation,rboundary)

    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Set up a discretisation structure which tells the code which
    ! finite element to use
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    ! Now we can start to initialise the discretisation. At first, set up
    ! a block discretisation structure that specifies 3 blocks in the
    ! solution vector.
    call spdiscr_initBlockDiscr (rdiscretisation,3,&
                                 rtriangulation, rboundary)

    ! rdiscretisation%RspatialDiscr is a list of scalar
    ! discretisation structures for every component of the solution vector.
    ! We have a solution vector with three components:
    !  Component 1 = X-velocity
    !  Component 2 = Y-velocity
    !  Component 3 = Pressure
    ! For simplicity, we set up one discretisation structure for the
    ! velocity...
    call spdiscr_initDiscr_simple (rdiscretisation%RspatialDiscr(1),&
                EL_EM30, rtriangulation, rboundary)

    ! ...and copy this structure also to the discretisation structure
    ! of the 2nd component (Y-velocity). This needs no additional memory,
    ! as both structures will share the same dynamic information afterwards.
    call spdiscr_duplicateDiscrSc(rdiscretisation%RspatialDiscr(1),&
        rdiscretisation%RspatialDiscr(2))

    ! For the pressure (3rd component), we set up a separate discretisation
    ! structure, as this uses different finite elements for trial and test
    ! functions.
    call spdiscr_deriveSimpleDiscrSc (rdiscretisation%RspatialDiscr(1), &
        EL_Q0, rdiscretisation%RspatialDiscr(3))

    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Set up an cubature info structure to tell the code which cubature
    ! formula to use
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    ! Create an assembly information structure which tells the code
    ! the cubature formula to use. Standard: Gauss 3x3.
    call spdiscr_createDefCubStructure(&
        rdiscretisation%RspatialDiscr(1),rcubatureInfo,CUB_GEN_AUTO_G3)

    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Create a block matrix structure for the operator
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    ! Initialise the block matrix with default values based on
    ! the discretisation.
    call lsysbl_createMatBlockByDiscr (rdiscretisation,rmatrix)

    ! Inform the matrix that we build a saddle-point problem.
    ! Normally, imatrixSpec has the value LSYSBS_MSPEC_GENERAL,
    ! but probably some solvers can use the special structure later.
    rmatrix%imatrixSpec = LSYSBS_MSPEC_SADDLEPOINT

    ! Now as the discretisation is set up, we can start to generate
    ! the structure of the system matrix which is to solve.
    ! We create that directly in the block (1,1) of the block matrix
    ! using the discretisation structure of the first block.
    !
    ! In the global system, there are two coupling matrices B1 and B2.
    ! Both have the same structure.
    !
    !    ( A         B1 )
    !    (      A    B2 )
    !    ( B1^T B2^T    )
    !
    ! Create the matrix structure of the X-velocity.
    call bilf_createMatrixStructure (rdiscretisation%RspatialDiscr(1),&
                                     LSYSSC_MATRIX9, rmatrix%RmatrixBlock(1,1))

    ! Use it for the Y-velocity
    call lsyssc_duplicateMatrix (rmatrix%RmatrixBlock(1,1),&
        rmatrix%RmatrixBlock(2,2),LSYSSC_DUP_SHARE,LSYSSC_DUP_IGNORE)

    ! Create the structure for the B-matrices using the 3rd
    ! spatial discretisation structure in p_rdiscretisation%RspatialDiscr.
    call bilf_createMatrixStructure (rdiscretisation%RspatialDiscr(3),&
                                     LSYSSC_MATRIX9, rmatrix%RmatrixBlock(1,3),&
                                     rdiscretisation%RspatialDiscr(1))

    call lsyssc_duplicateMatrix (rmatrix%RmatrixBlock(1,3),&
        rmatrix%RmatrixBlock(2,3),LSYSSC_DUP_SHARE,LSYSSC_DUP_IGNORE)

    ! Create the structure for the B^T matrices by transposing the structre
    ! of the B-matrices. The structure of B1^T and B2^T is again the same.
    call lsyssc_transposeMatrix (&
        rmatrix%RmatrixBlock(1,3),rmatrix%RmatrixBlock(3,1),LSYSSC_TR_STRUCTURE)

    call lsyssc_duplicateMatrix (rmatrix%RmatrixBlock(3,1),&
        rmatrix%RmatrixBlock(3,2),LSYSSC_DUP_SHARE,LSYSSC_DUP_IGNORE)

    ! Now re-assign the block discretisation structure to all matrices
    call lsysbl_assignDiscrDirectMat (rmatrix,rdiscretisation)

    ! Allocate memory for the matrix
    call lsysbl_allocEmptyMatrix (rmatrix,LSYSSC_SETM_ZERO)

    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Create a block matrix entries
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    ! Pass dnu via rcollection
    rcollection%DquickAccess(1) = dnu

    call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,&
          st2d3_fcalc_Stokes,rcubatureInfo=rcubatureInfo,rcollection=rcollection)

    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Create RHS and solution vectors
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    ! Create a RHS and a solution vector based on the discretisation.
    ! Fill with zero.
    call lsysbl_createVectorBlock (rdiscretisation,rrhs,.true.)
    call lsysbl_createVectorBlock (rdiscretisation,rvector,.true.)

    ! Assemble the right hand side.
    call bma_buildVector (rrhs,BMA_CALC_STANDARD,&
          st2d3_fcalc_rhs,rcubatureInfo=rcubatureInfo,rcollection=rcollection)

    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Assembly of matrices/vectors finished
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Discretise the boundary conditions and apply them to the matrix/RHS/sol.
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    ! For implementing boundary conditions, we use a `filter technique with
    ! discretised boundary conditions`. This means, we first have to calculate
    ! a discrete version of the analytic BC, which we can implement into the
    ! solution/RHS vectors using the corresponding filter.
    !
    ! Create a t_discreteBC structure where we store all discretised boundary
    ! conditions.
    call bcasm_initDiscreteBC(rdiscreteBC)

    ! We first set up the boundary conditions for the X-velocity, then those
    ! of the Y-velocity.
    !
    ! We "know" already (from the problem definition) that we have four boundary
    ! segments in the domain. Each of these, we want to use for enforcing
    ! some kind of boundary condition.
    !
    ! We ask the boundary routines to create a "boundary region" - which is
    ! simply a part of the boundary corresponding to a boundary segment.
    ! A boundary region roughly contains the type, the min/max parameter value
    ! and whether the endpoints are inside the region or not.
    call boundary_createRegion(rboundary,1,1,rboundaryRegion)

    ! The endpoint of this segment should also be Dirichlet. We set this by
    ! changing the region properties in rboundaryRegion.
    rboundaryRegion%iproperties = BDR_PROP_WITHSTART + BDR_PROP_WITHEND

    ! We use this boundary region and specify that we want to have Dirichlet
    ! boundary there. The following call does the following:
    ! - Create Dirichlet boundary conditions on the region rboundaryRegion.
    !   We specify icomponent="1" to indicate that we set up the
    !   Dirichlet BC`s for the first (here: one and only) component in the
    !   solution vector.
    ! - Discretise the boundary condition so that the BC`s can be applied
    !   to matrices and vectors
    ! - Add the calculated discrete BC`s to rdiscreteBC for later use.
    call bcasm_newDirichletBConRealBD (rdiscretisation,1,&
                                       rboundaryRegion,rdiscreteBC,&
                                       getBoundaryValues_2D)

    ! Edge 2 is Neumann boundary, so it is commented out.
    ! CALL boundary_createRegion(rboundary,1,2,rboundaryRegion)
    ! CALL bcasm_newDirichletBConRealBD (rdiscretisation,1,&
    !                                    rboundaryRegion,rdiscreteBC,&
    !                                    getBoundaryValues_2D)

    ! Edge 3 of boundary component 1.
    call boundary_createRegion(rboundary,1,3,rboundaryRegion)
    call bcasm_newDirichletBConRealBD (rdiscretisation,1,&
                                       rboundaryRegion,rdiscreteBC,&
                                       getBoundaryValues_2D)

    ! Edge 4 of boundary component 1. That is it.
    call boundary_createRegion(rboundary,1,4,rboundaryRegion)
    call bcasm_newDirichletBConRealBD (rdiscretisation,1,&
                                       rboundaryRegion,rdiscreteBC,&
                                       getBoundaryValues_2D)

    ! Now continue with defining the boundary conditions of the Y-velocity:
    !
    ! Define edge 1.
    call boundary_createRegion(rboundary,1,1,rboundaryRegion)

    ! Edge with start- and endpoint.
    rboundaryRegion%iproperties = BDR_PROP_WITHSTART + BDR_PROP_WITHEND

    ! As we define the Y-velocity, we now set icomponent=2 in the following call.
    call bcasm_newDirichletBConRealBD (rdiscretisation,2,&
                                       rboundaryRegion,rdiscreteBC,&
                                       getBoundaryValues_2D)

    ! Edge 2 is Neumann boundary, so it is commented out.
    ! CALL boundary_createRegion(rboundary,1,2,rboundaryRegion)
    ! CALL bcasm_newDirichletBConRealBD (rdiscretisation,2,&
    !                                    rboundaryRegion,rdiscreteBC,&
    !                                    getBoundaryValues_2D)

    ! Edge 3 of boundary component 1.
    call boundary_createRegion(rboundary,1,3,rboundaryRegion)
    call bcasm_newDirichletBConRealBD (rdiscretisation,2,&
                                       rboundaryRegion,rdiscreteBC,&
                                       getBoundaryValues_2D)

    ! Edge 4 of boundary component 1. That is it.
    call boundary_createRegion(rboundary,1,4,rboundaryRegion)
    call bcasm_newDirichletBConRealBD (rdiscretisation,2,&
                                       rboundaryRegion,rdiscreteBC,&
                                       getBoundaryValues_2D)

    ! The pressure does not need boundary conditions.

    ! Next step is to implement boundary conditions into the RHS,
    ! solution and matrix. This is done using a vector/matrix filter
    ! for discrete boundary conditions.
    call vecfil_discreteBCrhs (rrhs,rdiscreteBC)
    call vecfil_discreteBCsol (rvector,rdiscreteBC)
    call matfil_discreteBC (rmatrix,rdiscreteBC)

    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Set up a linear solver
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    ! During the linear solver, the boundary conditions are also
    ! frequently imposed to the vectors. But as the linear solver
    ! does not work with the actual solution vectors but with
    ! defect vectors instead.
    ! So, set up a filter chain that filters the defect vector
    ! during the solution process to implement discrete boundary conditions.
    call filter_initFilterChain (RfilterChain,nfilters)
    call filter_newFilterDiscBCDef (RfilterChain,nfilters,rdiscreteBC)

    ! Create a BiCGStab-solver with VANCA preconditioner.
    ! Attach the above filter chain to the solver, so that the solver
    ! automatically filters the vector during the solution process.
    nullify(p_rpreconditioner)
    call linsol_initVANKA (p_rpreconditioner,1.0_DP,LINSOL_VANKA_GENERAL)
    call linsol_initBiCGStab (p_rsolverNode,p_rpreconditioner,RfilterChain)

    ! Set the output level of the solver to 2 for some output
    p_rsolverNode%ioutputLevel = 2

    ! We will allow the solver to perform 200 iterations
    p_rsolverNode%nmaxIterations = 200

    ! Attach the system matrix to the solver.
    call linsol_setMatrix(p_rsolverNode,rmatrix)

    ! Initialise structure/data of the solver. This allows the
    ! solver to allocate memory / perform some precalculation
    ! to the problem.
    call linsol_initStructure (p_rsolverNode, ierror)

    if (ierror .ne. LINSOL_ERR_NOERROR) then
      call output_line("Matrix structure invalid!",OU_CLASS_ERROR)
      call sys_halt()
    end if

    call linsol_initData (p_rsolverNode, ierror)

    if (ierror .ne. LINSOL_ERR_NOERROR) then
      call output_line("Matrix singular!",OU_CLASS_ERROR)
      call sys_halt()
    end if

    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Solve the system
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    ! Finally solve the system. As we want to solve Ax=b with
    ! b being the real RHS and x being the real solution vector,
    ! we use linsol_solveAdaptively. If b is a defect
    ! RHS and x a defect update to be added to a solution vector,
    ! we would have to use linsol_precondDefect instead.
    call linsol_solveAdaptively (p_rsolverNode,rvector,rrhs,rtempBlock)

    ! -------------------------------------------------------------------------
    ! Solution vector projection for projecting a solution vector of
    ! an arbitrary element to Q1 such that it can be written into a VTK file
    ! -------------------------------------------------------------------------

    ! The solution vector is probably not in the way GMV likes it!
    ! GMV for example does not understand Q1~ vectors!
    ! Therefore, we first have to convert the vector to a form that
    ! GMV understands.
    ! GMV understands only Q1 solutions! So the task is now to create
    ! a Q1 solution from p_rvector and write that out.
    !
    ! For this purpose, first create a "derived" simple discretisation
    ! structure based on Q1 by copying the main guiding block discretisation
    ! structure and modifying the discretisation structures of the
    ! two velocity subvectors:

    call spdiscr_duplicateBlockDiscr (rdiscretisation,rprjDiscretisation)

    call spdiscr_deriveSimpleDiscrSc (rdiscretisation%RspatialDiscr(1), &
                 EL_Q1, rprjDiscretisation%RspatialDiscr(1))

    call spdiscr_deriveSimpleDiscrSc (rdiscretisation%RspatialDiscr(2), &
                 EL_Q1, rprjDiscretisation%RspatialDiscr(2))

    ! The pressure discretisation substructure stays the old.
    !
    ! Now set up a new solution vector based on this discretisation,
    ! allocate memory.
    call lsysbl_createVectorBlock (rprjDiscretisation,rprjVector,.true.)

    ! Then take our original solution vector and convert it according to the
    ! new discretisation:
    call spdp_projectSolution (rvector,rprjVector)

    ! Discretise the boundary conditions according to the Q1/Q1/Q0
    ! discretisation.
    !
    ! Create a t_discreteBC structure where we store all discretised boundary
    ! conditions.
    call bcasm_initDiscreteBC(rprjDiscreteBC)
    !
    ! Edge 1 of boundary component 1, X-velocity.
    call boundary_createRegion(rboundary,1,1,rboundaryRegion)

    ! Edge with start- and endpoint.
    rboundaryRegion%iproperties = BDR_PROP_WITHSTART + BDR_PROP_WITHEND

    call bcasm_newDirichletBConRealBD (rprjDiscretisation,1,&
                                       rboundaryRegion,rprjDiscreteBC,&
                                       getBoundaryValues_2D)

    ! Edge 2 is Neumann boundary, so it is commented out.
    ! CALL boundary_createRegion(rboundary,1,2,rboundaryRegion)
    ! CALL bcasm_newDirichletBConRealBD (rprjDiscretisation,1,&
    !                                    rboundaryRegion,rprjDiscreteBC,&
    !                                    getBoundaryValues_2D)

    ! Edge 3 of boundary component 1.
    call boundary_createRegion(rboundary,1,3,rboundaryRegion)
    call bcasm_newDirichletBConRealBD (rprjDiscretisation,1,&
                                       rboundaryRegion,rprjDiscreteBC,&
                                       getBoundaryValues_2D)

    ! Edge 4 of boundary component 1. That is it.
    call boundary_createRegion(rboundary,1,4,rboundaryRegion)
    call bcasm_newDirichletBConRealBD (rprjDiscretisation,1,&
                                       rboundaryRegion,rprjDiscreteBC,&
                                       getBoundaryValues_2D)

    ! Edge 1 of boundary component 1, Y-velocity.
    call boundary_createRegion(rboundary,1,1,rboundaryRegion)

    ! Edge with start- and endpoint.
    rboundaryRegion%iproperties = BDR_PROP_WITHSTART + BDR_PROP_WITHEND

    ! As we define the Y-velocity, we now set icomponent=2 in the following call.
    call bcasm_newDirichletBConRealBD (rprjDiscretisation,2,&
                                       rboundaryRegion,rprjDiscreteBC,&
                                       getBoundaryValues_2D)

    ! Edge 2 is Neumann boundary, so it is commented out.
    ! CALL boundary_createRegion(rboundary,1,2,rboundaryRegion)
    ! CALL bcasm_newDirichletBConRealBD (rprjDiscretisation,2,&
    !                                    rboundaryRegion,rprjDiscreteBC,&
    !                                    getBoundaryValues_2D)

    ! Edge 3 of boundary component 1.
    call boundary_createRegion(rboundary,1,3,rboundaryRegion)
    call bcasm_newDirichletBConRealBD (rprjDiscretisation,2,&
                                       rboundaryRegion,rprjDiscreteBC,&
                                       getBoundaryValues_2D)

    ! Edge 4 of boundary component 1. That is it.
    call boundary_createRegion(rboundary,1,4,rboundaryRegion)
    call bcasm_newDirichletBConRealBD (rprjDiscretisation,2,&
                                       rboundaryRegion,rprjDiscreteBC,&
                                       getBoundaryValues_2D)

    ! Send the vector to the boundary-condition implementation filter.
    ! This modifies the vector according to the discrete boundary
    ! conditions.
    call vecfil_discreteBCsol (rprjVector,rprjDiscreteBC)

    ! Get the path for writing postprocessing files from the environment variable
    ! $UCDDIR. If that does not exist, write to the directory "./gmv".
    if (.not. sys_getenv_string("UCDDIR", sucddir)) sucddir = "./gmv"

    ! Now we have a Q1/Q1/Q0 solution in rprjVector.
    ! We can now start the postprocessing.
    ! Start UCD export to GMV file:
    call ucd_startVTK (rexport,UCD_FLAG_STANDARD,rtriangulation,&
        trim(sucddir)//"/u2d_3_block.vtk")

    ! Write velocity field
    call lsyssc_getbase_double (rprjVector%RvectorBlock(1),p_Ddata)
    call lsyssc_getbase_double (rprjVector%RvectorBlock(2),p_Ddata2)

    ! In case we use the VTK exporter, which supports vector output, we will
    ! pass the X- and Y-velocity at once to the ucd module.
    call ucd_addVarVertBasedVec(rexport,"velocity",p_Ddata,p_Ddata2)

    ! If we use the GMV exporter, we might replace the line above by the
    ! following two lines:
    !CALL ucd_addVariableVertexBased (rexport,"X-vel",UCD_VAR_XVELOCITY, p_Ddata)
    !CALL ucd_addVariableVertexBased (rexport,"Y-vel",UCD_VAR_YVELOCITY, p_Ddata2)

    ! Write pressure
    call lsyssc_getbase_double (rprjVector%RvectorBlock(3),p_Ddata)
    call ucd_addVariableElementBased (rexport,"pressure",UCD_VAR_STANDARD, p_Ddata)

    ! Write the file to disc, that is it.
    call ucd_write (rexport)
    call ucd_release (rexport)

    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Projection and VTK export finished.
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Error calculation
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    ! Store the viscosity parameter nu in the collection"s quick access array
    rcollection%DquickAccess(1) = dnu

    ! Set up revalVectors with the velocity/pressure vectors.
    call fev2_addVectorToEvalList(revalVectors,rvector%RvectorBlock(1),1)
    call fev2_addVectorToEvalList(revalVectors,rvector%RvectorBlock(2),1)
    call fev2_addVectorToEvalList(revalVectors,rvector%RvectorBlock(3),0)

    ! Calculate errors of velocity and pressure against analytic solutions.
    rcollection%IquickAccess(1) = 1  ! X-component, L2-error
    call bma_buildIntegral(DerrorUL2(1),BMA_CALC_STANDARD,st2d3_fcalc_L2error,&
        rcollection=rcollection,revalVectors=revalVectors,&
        rcubatureInfo=rcubatureInfo)

    rcollection%IquickAccess(1) = 2  ! Y-component, L2-error
    call bma_buildIntegral(DerrorUL2(2),BMA_CALC_STANDARD,st2d3_fcalc_L2error,&
        rcollection=rcollection,revalVectors=revalVectors,&
        rcubatureInfo=rcubatureInfo)

    rcollection%IquickAccess(1) = 3  ! Pressure, L2-error
    call bma_buildIntegral(derrorPL2(1),BMA_CALC_STANDARD,st2d3_fcalc_L2error,&
        rcollection=rcollection,revalVectors=revalVectors,&
        rcubatureInfo=rcubatureInfo)

    rcollection%IquickAccess(1) = 1  ! X-component, H1-error
    call bma_buildIntegral(DerrorUH1(1),BMA_CALC_STANDARD,st2d3_fcalc_H1error,&
        rcollection=rcollection,revalVectors=revalVectors,&
        rcubatureInfo=rcubatureInfo)

    rcollection%IquickAccess(1) = 2  ! Y-component, H1-error
    call bma_buildIntegral(DerrorUH1(2),BMA_CALC_STANDARD,st2d3_fcalc_H1error,&
        rcollection=rcollection,revalVectors=revalVectors,&
        rcubatureInfo=rcubatureInfo)

    ! Print the errors.
    call output_lbrk()
    call output_line("|u - u_h|_L2 = " // trim(sys_sdEL(sqrt(DerrorUL2(1)), 10)) &
                                // " " // trim(sys_sdEL(sqrt(DerrorUL2(2)), 10)))
    call output_line("|u - u_h|_H1 = " // trim(sys_sdEL(sqrt(DerrorUH1(1)), 10)) &
                                // " " // trim(sys_sdEL(sqrt(DerrorUH1(2)), 10)))
    call output_line("|p - p_h|_L2 = " // trim(sys_sdEL(sqrt(derrorPL2(1)), 10)))

    ! Cleanup
    call fev2_releaseVectorList(revalVectors)

    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! Clean up
    ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    ! We are finished - but not completely!
    ! Now, clean up so that all the memory is available again.
    !
    ! Release solver data and structure
    call linsol_doneData (p_rsolverNode)
    call linsol_doneStructure (p_rsolverNode)

    ! Release the solver node and all subnodes attached to it (if at all):
    call linsol_releaseSolver (p_rsolverNode)

    ! Release the filter chain
    call filter_doneFilterChain (RfilterChain,nfilters)

    ! Release the block matrix/vectors
    call lsysbl_releaseVector (rprjVector)
    call lsysbl_releaseVector (rtempBlock)
    call lsysbl_releaseVector (rvector)
    call lsysbl_releaseVector (rrhs)
    call lsysbl_releaseMatrix (rmatrix)

    ! Release our discrete version of the boundary conditions
    call bcasm_releaseDiscreteBC (rprjDiscreteBC)
    call bcasm_releaseDiscreteBC (rdiscreteBC)

    ! Release the discretisation structure and all spatial discretisation
    ! structures in it.
    call spdiscr_releaseBlockDiscr(rprjDiscretisation)
    call spdiscr_releaseBlockDiscr(rdiscretisation)

    ! Release the triangulation.
    call tria_done (rtriangulation)

    ! Finally release the domain, that is it.
    call boundary_release (rboundary)

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


!***********************************************************************************************************************************
! UNIT:
!
!  (MODULE) MR_DEF_RANKS
!
! PURPOSE:
!
!   TO DEFINE THE RANK PARAMETERS DETERMINING THE SIZE OF ARRAYS.
!
! DEFINITION OF VARIABLES:
!
!                             NND    :    NUMBER OF NODES
!                             NEM    :    NUMBER OF ELEMENTS
!                          NI, NJ    :    NUMBER OF COLUMNS AND ROWS, RESPECTIVELY
!                              NK    :    NUMBER OF LAYERS
!                             NKS    :    NUMBER OF SEDIMENT SIZES
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  MODULE MR_DEF_RANKS

    USE MR_KINDS

    IMPLICIT NONE

    INTEGER(NDID_KIND) :: NND
    INTEGER(EMID_KIND) :: NEM
    INTEGER(IJID_KIND) :: NI , NJ
    INTEGER(KKID_KIND) :: NK = 10
    INTEGER(KSID_KIND) :: NKS = 0

  END MODULE MR_DEF_RANKS


text:
Address(Los Robles Street) is a residential Culdesacs culdesac in South Davis.

 Intersecting Streets

  Lillard Drive
  Monte Vista Avenue and Monte Vista Place
  Monterey Avenue
  Laguna Avenue

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! Copyright (C) 2015, 2016, Simone Marsili
! All rights reserved.
! License: BSD 3 clause

module scrs
  use nrtype
  implicit none
  private

  public :: compute_scores
  public :: print_scores

  real(DP), allocatable, save :: scores(:,:)

contains

  subroutine compute_scores(skip_gaps)
    use data, only: nv
    use model, only: couplings
    logical, intent(in) :: skip_gaps
    integer :: iv,jv,k
    integer(I4B) :: err

    ! at the very end of the run
    allocate(scores(nv,nv),stat=err)
    scores = 0.0_DP

    ! compute scores
    k = 0
    do jv = 1,nv-1
       do iv = jv+1,nv
          k = k + 1
          if (skip_gaps) then
             scores(iv,jv) = frobenius(couplings(2:,2:,k))
          else
             scores(iv,jv) = frobenius(couplings(:,:,k))
          end if
          scores(jv,iv) = scores(iv,jv)
       end do
    end do

    call apc_correction(scores)

  contains

    real(DP) function frobenius(a)
      real(DP), intent(in) :: a(:,:)

      frobenius = sqrt(sum(a**2))

    end function frobenius

    subroutine apc_correction(a)
      real(DP), intent(inout) :: a(:,:)
      real(DP), allocatable :: sums(:)
      real(DP) :: totsum
      integer(I4B) :: err
      integer(I4B) :: i,j
      integer(I4B) :: n

      n = size(a,1)
      allocate(sums(n),stat=err)

      do i = 1,n
         sums(i) = sum(a(:,i))
      end do
      totsum = sum(sums)

      do j = 1,n
         do i = j,n
            a(i,j) = a(i,j) - sums(i)*sums(j)/totsum
            a(j,i) = a(i,j)
         end do
      end do

      deallocate(sums)

    end subroutine apc_correction

  end subroutine compute_scores

  subroutine print_scores(uscrs)
    use units
    use data, only: nv
    integer(I4B), intent(in) :: uscrs
    integer(I4B) :: iv,jv
    real(8) :: sij

    do iv = 1,nv-1
       do jv = iv+1,nv
          sij = scores(iv,jv)
          write(uscrs,'(i5,1x,i5,1x,f8.5)') iv,jv,sij
       end do
    end do

  end subroutine print_scores

end module scrs

max_stars_repo_name: simomarsili/mpl

max_issues_repo_name: simomarsili/mpl

max_forks_repo_name: simomarsili/mpl


!
! Copyright (C) by Argonne National Laboratory
!     See COPYRIGHT in top-level directory
!

subroutine MPI_Free_mem_f08(base, ierror)
    use, intrinsic :: iso_c_binding, only : c_int
    use :: mpi_c_interface, only : MPIR_Free_mem_c

    implicit none

    type(*), dimension(..), intent(in), asynchronous :: base
    integer, optional, intent(out) :: ierror

    integer(c_int) :: ierror_c

    ierror_c = MPIR_Free_mem_c(base)

    if (present(ierror)) ierror = ierror_c

end subroutine MPI_Free_mem_f08


text:
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
! test for comparison operation take quad precision

program test
  integer, parameter :: n = 2
  logical :: r1 = 1.1_16 > 1.2_16
  logical :: r2
  integer :: result(n), expect(n)
  expect = 1
  result = 0

  r2 = 621424.2346523454241479464_16 <= 621424.234652345424147954_16
  if (r1 .eq. .false.) result(1) = 1
  if (r2 .eq. .true.) result(2) = 1
  call check(result, expect, n)

end program

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


subroutine psb_cdall(ictxt, desc, info,mg,ng,parts,vg,vl,flag,nl,repl, globalcheck,lidx)
  use psb_desc_mod
  use psb_serial_mod
  use psb_const_mod
  use psb_error_mod
  use psb_penv_mod
  use psb_cd_tools_mod, psb_protect_name => psb_cdall
  use psi_mod
  implicit None
  procedure(psb_parts)               :: parts
  integer(psb_ipk_), intent(in)               :: mg,ng,ictxt, vg(:), vl(:),nl,lidx(:)
  integer(psb_ipk_), intent(in)               :: flag
  logical, intent(in)               :: repl, globalcheck
  integer(psb_ipk_), intent(out)              :: info
  type(psb_desc_type), intent(out)  :: desc

  optional :: mg,ng,parts,vg,vl,flag,nl,repl, globalcheck,lidx

  interface
    subroutine psb_cdals(m, n, parts, ictxt, desc, info)
      use psb_desc_mod
      procedure(psb_parts)               :: parts
      integer(psb_ipk_), intent(in)                 :: m,n,ictxt
      Type(psb_desc_type), intent(out)    :: desc
      integer(psb_ipk_), intent(out)                :: info
    end subroutine psb_cdals
    subroutine psb_cdalv(v, ictxt, desc, info, flag)
      use psb_desc_mod
      integer(psb_ipk_), intent(in)               :: ictxt, v(:)
      integer(psb_ipk_), intent(in), optional     :: flag
      integer(psb_ipk_), intent(out)              :: info
      Type(psb_desc_type), intent(out)  :: desc
    end subroutine psb_cdalv
    subroutine psb_cd_inloc(v, ictxt, desc, info, globalcheck,idx)
      use psb_desc_mod
      implicit None
      integer(psb_ipk_), intent(in)               :: ictxt, v(:)
      integer(psb_ipk_), intent(out)              :: info
      type(psb_desc_type), intent(out)  :: desc
      logical, intent(in), optional     :: globalcheck
      integer(psb_ipk_), intent(in), optional     :: idx(:)
    end subroutine psb_cd_inloc
    subroutine psb_cdrep(m, ictxt, desc,info)
      use psb_desc_mod
      integer(psb_ipk_), intent(in)               :: m,ictxt
      Type(psb_desc_type), intent(out)  :: desc
      integer(psb_ipk_), intent(out)              :: info
    end subroutine psb_cdrep
  end interface
  character(len=20)   :: name
  integer(psb_ipk_) :: err_act, n_, flag_, i, me, np, nlp, nnv, lr
  integer(psb_ipk_), allocatable :: itmpsz(:)
  integer(psb_mpik_) :: iictxt

  if (psb_get_errstatus() /= 0) return
  info=psb_success_
  name = 'psb_cdall'
  call psb_erractionsave(err_act)

  call psb_info(ictxt, me, np)
  iictxt  = ictxt
  if (count((/ present(vg),present(vl),&
       &  present(parts),present(nl), present(repl) /)) /= 1) then
    info=psb_err_no_optional_arg_
    call psb_errpush(info,name,a_err=" vg, vl, parts, nl, repl")
    goto 9999
  endif

  desc%base_desc => null()
  if (allocated(desc%indxmap)) then
    write(0,*) 'Allocated on an intent(OUT) var?'
  end if

  if (present(parts)) then

    if (.not.present(mg)) then
      info=psb_err_no_optional_arg_
      call psb_errpush(info,name)
      goto 9999
    end if
    if (present(ng)) then
      n_ = ng
    else
      n_ = mg
    endif
    call  psb_cdals(mg, n_, parts, ictxt, desc, info)

  else if (present(repl)) then

    if (.not.present(mg)) then
      info=psb_err_no_optional_arg_
      call psb_errpush(info,name)
      goto 9999
    end if
    if (.not.repl) then
      info=psb_err_no_optional_arg_
      call psb_errpush(info,name)
      goto 9999
    end if

    call  psb_cdrep(mg, ictxt, desc, info)

  else if (present(vg)) then

    if (present(flag)) then
      flag_=flag
    else
      flag_=0
    endif
    if (present(mg)) then
      nnv = min(mg,size(vg))
    else
      nnv = size(vg)
    end if

    call psb_cdalv(vg(1:nnv), ictxt, desc, info, flag=flag_)

  else if (present(vl)) then

    if (present(nl)) then
      nnv = min(nl,size(vl))
    else
      nnv = size(vl)
    end if

    call psb_cd_inloc(vl(1:nnv),ictxt,desc,info, globalcheck=globalcheck,idx=lidx)

  else if (present(nl)) then

    if (np == 1) then
      allocate(psb_repl_map      :: desc%indxmap, stat=info)
    else
      allocate(psb_gen_block_map :: desc%indxmap, stat=info)
    end if
    if (info == psb_success_) then
      select type(aa => desc%indxmap)
      type is (psb_repl_map)
        call aa%repl_map_init(iictxt,nl,info)
      type is (psb_gen_block_map)
        call aa%gen_block_map_init(iictxt,nl,info)
      class default
        ! This cannot happen
        info = psb_err_internal_error_
        goto 9999
      end select
    end if

    call psb_realloc(1,itmpsz, info)
    if (info /= 0) then
      write(0,*) 'Error reallocating itmspz'
      goto 9999
    end if
    itmpsz(:) = -1
    call psi_bld_tmpovrl(itmpsz,desc,info)

  endif

  if (info /= psb_success_) goto 9999

  ! Finish off
  lr = desc%indxmap%get_lr()
  call psb_realloc(max(1,lr/2),desc%halo_index, info)
  if (info == psb_success_) call psb_realloc(1,desc%ext_index, info)
  if (info /= psb_success_) then
    info=psb_err_from_subroutine_
    call psb_errpush(info,name,a_err='psb_realloc')
    Goto 9999
  end if
  desc%halo_index(:)           = -1
  desc%ext_index(:)            = -1
  call psb_cd_set_bld(desc,info)
  if (info /= psb_success_) goto 9999

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(ictxt,err_act)

  return

end subroutine psb_cdall


! ====================================================================
!     subroutine ETp
!
!     purpose: 1. calculate the ET0 and Ep and Tp in the ??.et0 file.
!              2. calculate the Crop growth in the simulation period
!                 in the ??.crp file.
!              3. calculate the Epi[mm] and Tpi[mm] in every 1d layers
!                 in the ??.eti file.
! ====================================================================
! =========================Incoming variables=========================
!       num            The number of the 1D soil columns[-]
!       sp(:,:)        The discreted function for soil evaporation[-]
!       dx(:,:)        The z coordinates of  the discreted layers[m]
! =========================Outcoming variables========================
!       ??.ETi         file stored Epi[mm] and Tpi[mm] in every layer
! =========================related files==============================
!       ??.et0         ET0 file with Ep[mm] and Tp[mm]
!       ??.crpcrop     file with root depth[m]
!       cropdat.dat    crop datebase with the discreted function of
!                      crop transpiration
!       ??.ETi         file stored Epi[mm] and Tpi[mm] in every layer
! =========================related functions==========================
!       cdatebase      The datebase to find discreted function for
!                      soil evaporation and crop transpiration
!       crop           Find root depth[m]
!       Fdense         Find the contribution coefficients of the
!                      discreted layers
!       Ep_Tp          Call for Ep[mm] and Tp[mm]
! ====================================================================
    SUBROUTINE ETp(num,datapath)
    USE parm
    IMPLICIT NONE

    CHARACTER (LEN=6) :: cha1, cha2, cha3
    CHARACTER (LEN=8) :: date1
    CHARACTER (LEN=8), DIMENSION(5) :: da
    CHARACTER (LEN=100) :: datapath
    INTEGER (KIND=4) :: i,j,k,ii,jj,num,lenpath
    REAL (KIND=KR), DIMENSION(Nlayer) :: dz1,sed,dense
    REAL (KIND=KR), DIMENSION(21) :: crpdat
    REAL (KIND=KR), DIMENSION(Nlayer) :: work
    REAL (KIND=KR) :: ET0, EP, TP, Nouse, rd, dtp

    lenpath = Len_Trim(datapath)
    cha1="??.et0"
    cha2="??.crp"
    cha3="??.eti"

    date1 = date

    DO i=1,Nlayer
        dz1(i)=zx(i+1)
    ENDDO

    CALL Caldays(datapath, MaxAL)

    IF (lCrop) THEN
        CALL Crop(datapath,num,numc,date1,MaxAL)
        IF (Terr.ne.0) RETURN
    ENDIF
    CALL ep_tp(datapath,num,date1,MaxAL,IfPM)
    IF (Terr.ne.0) RETURN

    WRITE(cha1(1:2),"(i2.2)")num
    WRITE(cha3(1:2),"(i2.2)")num
    OPEN(4,file=datapath(1:lenpath)//'/'//cha1,status="unknown",err=901)
    OPEN(6,file=datapath(1:lenpath)//'/'//cha3,status="unknown",err=901)
    WRITE(6,*,err=901)"date ordinal Es(1:Nlayer)/mm    TS(1:Nlayer)/mm"
    READ(4,*,err=901)cha1
    IF (lCrop) THEN
        WRITE(cha2(1:2),"(i2.2)")num
        OPEN(5,file=datapath(1:lenpath)//'/'//cha2,status="unknown",err=901)
        READ(5,*,err=901)cha1
    ENDIF

    DO j=1,MaxAL
        READ(4,*,err=901)date1,Nouse,ET0,Ep,Tp

        dtp=max(0.05_KR*ET0, 0.05_KR)
        dtp=dtp/4.0_KR
        DO ii=1,NMat
            DO jj=1,4
                work(jj)=sp(jj,ii)
            ENDDO
        ENDDO
        CALL Fdense(Nlayer,dz1,dtp,work,sed,4)

        IF (lCrop) THEN
            READ(5,*,err=901)Nouse,Nouse,rd

            CALL cdatebase(datapath,numc,da,crpdat)
            IF (Terr.ne.0) RETURN
            DO jj=1,10
                work(jj)=crpdat(jj+10)
            ENDDO
            ptab = crpdat(21)
            CALL Fdense(Nlayer,dz1,0.1_KR*rd,work,dense,10)
        ENDIF
        DO k=1,Nlayer
            sed(k)=sed(k)*ep
            IF (lCrop) THEN
                dense(k)=dense(k)*tp
            ELSE
                dense(k)=0.0_KR
            ENDIF
        ENDDO
        WRITE(6,"(A9,I6,260f6.3)",err=901)date1,i,sed(1:Nlayer),dense(1:Nlayer)
    ENDDO
    CLOSE(4)
    IF (lCrop) CLOSE(5)
    CLOSE(6)
    RETURN

901 Terr=3
    RETURN
    END SUBROUTINE ETp

    SUBROUTINE Caldays(datapath, MaxAL)
    USE parm, ONLY : KI
    IMPLICIT NONE
    INTEGER (KIND=KI) :: MaxAL
    CHARACTER (LEN=100) :: datapath, cha1
    INTEGER (KIND=4) ::lenpath, Terr

    lenpath = Len_Trim(datapath)
    cha1 = '01.wea'
    MaxAL = 0
    OPEN(101,file=datapath(1:lenpath)//'/'//cha1,status="old",err=901)
    DO WHILE (.not. eof(101))
        READ(101,*)
        MaxAL = MaxAL + 1
    ENDDO
    MaxAL = MaxAL - 1
    CLOSE(101)
    RETURN

901 Terr=1
    RETURN
    END SUBROUTINE Caldays

! ====================================================================
!     Subroutine Steady ET
!
!     Purpose: Calculate the ET under steady condition.
! ====================================================================
    SUBROUTINE Steady_ET(num,datapath,res)
    USE parm
    IMPLICIT NONE
    CHARACTER (LEN=100) :: datapath
    CHARACTER (LEN=8), DIMENSION(5) :: da
    REAL (KIND=KR) :: dtp, ckc, rd, clai, claif, ETc, Ep, Tp
    REAL (KIND=KR) :: f,ckc1,ckc2,ckc3,rd1,rd2,clai1,clai2,clai3,clai4
    REAL (KIND=KR), DIMENSION(21) :: crpdat
    REAL (KIND=KR), DIMENSION(Nlayer) :: work,dz1,sed,dense,res
    INTEGER (KIND=4) :: ii, jj, num, jd1, Jdate
    INTEGER (KIND=4), DIMENSION(5) :: jd

    IF (lCrop) THEN
        CALL cdatebase(datapath,num,da,crpdat)
        IF (Terr.NE.0) RETURN
        ! Calculate the potential Ep and Tp.
        f=crpdat(1)
        ckc1=crpdat(2)
        ckc2=crpdat(3)
        ckc3=crpdat(4)
        rd1=crpdat(5)
        rd2=crpdat(6)
        clai1=crpdat(7)
        clai2=crpdat(8)
        clai3=crpdat(9)
        clai4=crpdat(10)

        jd1=Jdate(date)
        DO ii=1,5
            jd(ii)=Jdate(da(ii))
        ENDDO
        IF (jd1<jd(1).or.jd1>jd(5))THEN
            ckc=0.9
            rd=0
            clai=0
        ELSEIF (jd1.ge.jd(1).and.jd1<jd(2)) THEN
            ckc=ckc1
            rd=rd1+(rd2-rd1)*(jd1-jd(1))/(jd(3)-jd(1))
            clai=CLAI1+(CLAI2-CLAI1)*(jd1-jd(1))/(jd(2)-jd(1))
        ELSEIF (jd1.ge.jd(2).and.jd1<jd(3)) THEN
            ckc=ckc1+(ckc2-ckc1)*(jd1-jd(2))/(jd(3)-jd(2))
            rd=rd1+(rd2-rd1)*(jd1-jd(1))/(jd(3)-jd(1))
            clai=CLAI2+(CLAI3-CLAI2)*(jd1-jd(2))/(jd(3)-jd(2))
        ELSEIF (jd1.ge.jd(3).and.jd1<jd(4)) THEN
            ckc=ckc2
            rd=rd2
            clai=clai3
        ELSEIF (jd1.ge.jd(4).and.jd1.le.jd(5)) THEN
	        ckc=ckc2-(ckc2-ckc3)*(jd1-jd(4))/(jd(5)-jd(4))
            rd=rd2
            clai=CLAI3-(CLAI3-CLAI4)*(jd1-jd(4))/(jd(5)-jd(4))
        ENDIF
        Claif=exp(-f*clai)
        ETc=ckc*rET
        Ep = Claif*ETc
        Tp = ETc-Ep
    ELSE
        Ep = rET
        Tp = 0.0_KR
    ENDIF

    ! Redistribution to each layer.
    dtp = max(0.05_KR*rET, 0.05_KR)
    dtp = dtp/4.0_KR
    DO ii=1,NMat
        DO jj=1,4
            work(jj)=sp(jj,ii)
        ENDDO
    ENDDO
    DO ii=1,Nlayer
        dz1(ii)=zx(ii+1)
    ENDDO
    CALL Fdense(Nlayer,dz1,dtp,work,sed,4)

    IF (lCrop) THEN
        DO jj=1,10
            work(jj)=crpdat(jj+10)
        ENDDO
        ptab = crpdat(21)
        CALL Fdense(Nlayer,dz1,0.1_KR*rd,work,dense,10)
    ENDIF
    DO ii=1,Nlayer
        sed(ii)=sed(ii)*ep
        IF (lCrop) THEN
            dense(ii)=dense(ii)*tp
        ELSE
            dense(ii)=0.0_KR
        ENDIF
    ENDDO
    res = sed+dense

    RETURN
    END SUBROUTINE Steady_ET

! ====================================================================
!   Subroutine ep_tp
!
!   Purpose: calculatin the potential Ep and Tp
! ====================================================================
! =========================Incoming variables=========================
!   num            The number of the crop[-]
!   dateini        The beginning day of simulation[yyyymmdd]
!   interval       The duration or length of the simulaiton day[-]
! =========================Outcoming variables========================
!   ??.ET0         A file stored ET0[m],Ep[m],Tp[m],PE[m]
! =========================related files==============================
!   ??.crp         The file stored information about Kc and LAIF is
!                  needed.
!   ??.ET0         A file stored ET0[m],Ep[m],Tp[m],PE[m]
! =========================related functions==========================
!   None.
! ====================================================================
    SUBROUTINE ep_tp(datapath,num,dateini,interval,IfPM)
    USE parm, ONLY : KR, KI, Terr, lCrop
    IMPLICIT NONE

    LOGICAL (KIND=4) :: IfPM
    CHARACTER (LEN=6) :: cha1, cha2, cha3
    CHARACTER (LEN=8) :: dateini
    CHARACTER (LEN=100) :: datapath
    INTEGER (KIND=KI) ::interval
    !INTEGER (kind=KI) :: yearini, monthini, dayini
    INTEGER (KIND=4) :: num, ierr, i, year, month, day, julie, DM, lenpath
    REAL (KIND=KR) :: tmean,tmax,tmin,nact,rhmean,uh,ea,deta,dr,si,tansi,&
    tmean1,ws,ra,np,rns,eamax,eamin,ed,rnl,rn,g,lamta,u2,gama,ato_p,et0,dt,PE
    REAL (KIND=KR) :: fai,hhh,hu
    REAL (KIND=KR) :: aa, ckc, claif, etc, ep, tp

    lenpath = Len_Trim(datapath)

    IF (lCrop) THEN
        cha1="??.wea"
        cha2="??.et0"
        cha3="??.crp"
        WRITE(cha1(1:2),"(i2.2)")num
        WRITE(cha2(1:2),"(i2.2)")num
        WRITE(cha3(1:2),"(i2.2)")num
        OPEN(1,file=datapath(1:lenpath)//'/'//cha1,status="old",err=901)
        OPEN(3,file=datapath(1:lenpath)//'/'//cha3,status="old",err=902)
        OPEN(2,file=datapath(1:lenpath)//'/'//cha2,status="unknown",err=902)
        WRITE(2,*)"year month day ordinal et0 ep tp	pe"
        read(3,*,err=902) !/lai and kc/

        IF (IfPM) THEN
            read(1,*,err=901)
            read(1,*,err=901)fai,hhh,hu
            read(1,*,err=901)
        ELSE
            read(1,*,err=901)
        ENDIF

        DO i=1,interval
            IF (IfPM) THEN
                READ(1,*,err=901)year,month,day,tmean,tmax,tmin,nact,rhmean,uh,ato_p,PE
                CALL Refet(year,month,day,Tmean,tmax,tmin,nact,RHmean,uh,hu,ato_p,hhh,fai,ET0,tmean1)
            ELSEIF (.not.IfPM) THEN
                READ(1,*,err=901)year,month,day,et0,PE
            ENDIF

            READ(3,*,err=902)aa,aa,ckc,aa,aa,CLAif
            ETc=ckc*ET0

!            CLAIF = 0.25    !2017-06-02

            Ep=cLAIf*ETc    !CLAI exp[-f*LAI]
            Tp=ETc-Ep
            WRITE(2,"(i4,i2.2,i2.2,i6,4f8.3)",err=902)year,month,day,i,et0,ep,tp,PE
        ENDDO

        CLOSE(1)
        CLOSE(2)
        CLOSE(3)
    ELSE
        cha1="??.wea"
        cha2="??.et0"
        WRITE(cha1(1:2),"(i2.2)")num
        WRITE(cha2(1:2),"(i2.2)")num
        OPEN(1,file=datapath(1:lenpath)//'/'//cha1,status="old",err=901)
        OPEN(2,file=datapath(1:lenpath)//'/'//cha2,status="unknown",err=902)
        WRITE(2,*)"year month day ordinal et0 ep tp	pe"

        IF (IfPM) THEN
            read(1,*,err=901)
            read(1,*,err=901)fai,hhh,hu
            read(1,*,err=901)
        ELSE
            read(1,*,err=901)
        ENDIF

        DO i=1,interval
            IF (IfPM) THEN
                READ(1,*,err=901)year,month,day,tmean,tmax,tmin,nact,rhmean,uh,ato_p,PE
                CALL Refet(year,month,day,Tmean,tmax,tmin,nact,RHmean,uh,hu,ato_p,hhh,fai,ET0,tmean1)
            ELSEIF (.not.IfPM) THEN
                READ(1,*,err=901)year,month,day,et0,PE
            ENDIF

            Ep=ET0
            Tp=0.0_KR
            WRITE(2,"(i4,i2.2,i2.2,i6,4f8.3)",err=902)year,month,day,i,et0,ep,tp,PE
        ENDDO

        CLOSE(1)
        CLOSE(2)
    ENDIF
    RETURN

901 Terr=1
    RETURN
902 Terr=3
    RETURN

    END SUBROUTINE ep_tp

! ====================================================================
!   Subroutine CROP
!
!   Purpose: Calculation Kc, RootDepth, LAI, CLAI for each day, by an
!            interpolation method.
! ====================================================================
! =========================Incoming variables=========================
!   num        The number of the 1D soil columns[-]
!   numc       The number of the kind of crop[-]
!   dateini    The beginning day of simulation[yyyymmdd]
!   interval   The duration or length of the simulaiton day[-]
! =========================Outcoming variables========================
!   ??.crp     A file stored Kc[-],RootDepth[m],LAI[-],CLAI[-].
! =========================related files==============================
!   ??.crp     cropdat.dat    crop datebase.
! =========================related functions==========================
!   None.
! ====================================================================
    SUBROUTINE Crop(datapath,num,numc,dateini,interval)
    USE parm, ONLY : KR, KI, Terr
    IMPLICIT NONE

    CHARACTER (len=8) :: Dateini, Date, Date1, Db
    CHARACTER (len=8), DIMENSION(5) :: da
    CHARACTER (len=6) :: nameF
    CHARACTER (len=10) :: cha1
    CHARACTER (LEN=100) :: datapath
    INTEGER (kind=KI) interval	! the simulation period
    INTEGER (KIND=4), DIMENSION(5) :: jd
    INTEGER (KIND=4) :: numc, nyear, nyear1, ndyear, JDATE, jd1, jmax, nyrr, dd, ifrun
    REAL (KIND=KR), DIMENSION(21) :: crpdat
    REAL (KIND=KR) :: f, ckc1, ckc2, ckc3, rd1, rd2, clai1, clai2, clai3, clai4
    REAL (KIND=KR) :: ckc, rd, clai, claif
    INTEGER (KIND=4) :: lenpath, num, i, j, k

    lenpath = Len_Trim(datapath)
    namef="??.crp"
    WRITE(namef(1:2),"(i2.2)")num
    OPEN(11,file=datapath(1:lenpath)//'/'//nameF,status="unknown",err=901)
    WRITE(11,*)"date    ordinal Kc  rootdepth   LAI LAIF    rp(1:Nlayer)"
    CALL cdatebase(datapath,numc,da,crpdat)
    IF (Terr.ne.0) RETURN

    f=crpdat(1)
    ckc1=crpdat(2)
    ckc2=crpdat(3)
    ckc3=crpdat(4)
    rd1=crpdat(5)
    rd2=crpdat(6)
    clai1=crpdat(7)
    clai2=crpdat(8)
    clai3=crpdat(9)
    clai4=crpdat(10)

    ! Reset the year of the dataset to the calculation year.
    READ(da(1)(1:4),"(i4)")nyear
    READ(Dateini(1:4),"(i4)")nyear1
    ndyear=nyear1-nyear
    DO i=1,5
        READ(da(i)(1:4),"(i4)")nyear
        WRITE(da(i)(1:4),"(i4)")nyear+ndyear
    ENDDO

    DO i=1,5
        jd(i)=Jdate(da(i))
    ENDDO

    date=dateini
    jd1=Jdate(date)

    DO i=1,interval
        IF (jd1<jd(1).or.jd1>jd(5))THEN
            ckc=0.9
            rd=0
            clai=0
        ELSEIF (jd1.ge.jd(1).and.jd1<jd(2)) THEN
            ckc=ckc1
            rd=rd1+(rd2-rd1)*(jd1-jd(1))/(jd(3)-jd(1))
            clai=CLAI1+(CLAI2-CLAI1)*(jd1-jd(1))/(jd(2)-jd(1))
        ELSEIF (jd1.ge.jd(2).and.jd1<jd(3)) THEN
            ckc=ckc1+(ckc2-ckc1)*(jd1-jd(2))/(jd(3)-jd(2))
            rd=rd1+(rd2-rd1)*(jd1-jd(1))/(jd(3)-jd(1))
            clai=CLAI2+(CLAI3-CLAI2)*(jd1-jd(2))/(jd(3)-jd(2))
        ELSEIF (jd1.ge.jd(3).and.jd1<jd(4)) THEN
            ckc=ckc2
            rd=rd2
            clai=clai3
        ELSEIF (jd1.ge.jd(4).and.jd1.le.jd(5)) THEN
	        ckc=ckc2-(ckc2-ckc3)*(jd1-jd(4))/(jd(5)-jd(4))
            rd=rd2
            clai=CLAI3-(CLAI3-CLAI4)*(jd1-jd(4))/(jd(5)-jd(4))
        ENDIF
        Claif=exp(-f*Clai)

        WRITE(11,"(A9,I6,14f6.3)")date,i,ckc,rd,clai,claif !
        CALL dateadd(date,1,date1)
        date=date1
        jd1=jd1+1

        ! If the calculation date is out of the database range !
        jmax=jd(5) ! the maximum day number
        IF (jd1>jmax) THEN
            READ(da(1)(1:4),"(i4)")nyear
            READ(da(5)(1:4),"(i4)")nyear1
            IF (nyear == nyear1) THEN
                IF (ifrun(nyear) == 0) THEN
                    nyrr = 365
                ELSE
                    Db = da(1)
                    WRITE(Db(5:8),"(i4)")0229
                    dd = Jdate(Db)
                    IF (jd(1)>dd) THEN
                        nyrr = 365
                    ELSE
                        nyrr = 366
                    ENDIF
                ENDIF
            ELSE
                IF ((ifrun(nyear)==0) .and. (ifrun(nyear+1)==0)) THEN
                    nyrr = 365
                ELSEIF (ifrun(nyear)==1) THEN
                    Db = da(1)
                    WRITE(Db(5:8),"(i4)")0229
                    dd = Jdate(Db)
                    IF (jd(1)>dd) THEN
                        nyrr = 365
                    ELSE
                        nyrr = 366
                    ENDIF
                ELSE
                    Db = da(5)
                    WRITE(Db(5:8),"(i4)")0229
                    dd = Jdate(Db)
                    IF (dd>jd(5)) THEN
                        nyrr = 365
                    ELSE
                        nyrr = 366
                    ENDIF
                ENDIF
            ENDIF
            jd1=jd1-Nyrr
        ENDIF

    ENDDO
    CLOSE(11)
    RETURN

901 Terr=3
    RETURN

    END SUBROUTINE Crop

! ====================================================================
!   Subroutine Fdense
!
!   Purpose: Calculation the root density function and evaporation
!            cumulative distribution function.
! ====================================================================
! =========================Incoming variables=========================
!   n1         the number for interpolation,to T is 10, to E is 4.[-]
!   n          the layer number.
!   dz1(:)     the depth of the lower boundary of the layer in its
!              controled area.[m]
!   re         the actual length of root or depth of evaporation.[m]
!   work(:)    the discreted function for interpolation.[-]
! =========================Outcoming variables========================
!   dense(:)   the interpolation results.[-]
! =========================related files==============================
!   None.
! =========================related functions==========================
!   None.
! ====================================================================
    SUBROUTINE Fdense(mlayer,dz1,re,work,dense,N1)
    USE parm
    IMPLICIT NONE

    REAL (KIND=KR), DIMENSION(Nlayer) :: dz1
    REAL (KIND=KR), DIMENSION(Nlayer) :: dense
    REAL (KIND=KR), DIMENSION(10) :: work
    REAL (KIND=KR) :: r2, re, r1
    INTEGER (KIND=KI) :: mlayer
    INTEGER (KIND=4) ::  i, j, L, N1

    DO i=1,mlayer
        dense(i)=0.0_KR
    ENDDO

    r2=0.0_KR

    DO j=1,mlayer
        L=int(dz1(j)/re)

        IF (l-0.0_KR<Tol) THEN
            r1=dz1(j)/re*work(L+1)
        ELSEIF (L>0.0_KR.and.L<N1) THEN
            r1=work(L)+(dz1(j)-re*L)/re*(work(L+1)-work(L))
        ELSEIF (L.GE.N1) THEN
            r1=1.0_KR
        ENDIF

        dense(j)=r1
        dense(j)=r1-r2
        r2=r1
    ENDDO

    IF (re-0.0_KR<Tol) THEN
        dense = 0.0_KR
    ENDIF

    END SUBROUTINE Fdense

! ====================================================================
!   Subroutine cdatebase
!
!   Purpose: read the information in crop datebase for crop NumC
! ====================================================================
! =========================Incoming variables=========================
!   Numc            number of the crop[-]
! =========================Outcoming variables========================
!   da(:)           date of the 5 stage of crop[yyyymmdd]
!   crpdat(:)       other numerical information such as Kc,LAI,f,
!                   root depth[m].
! =========================related files==============================
!   cropdat.dat     crop datebase with the discreted function of
!                   transpiration
! =========================related functions==========================
!   None.
! ====================================================================
    SUBROUTINE cdatebase(datapath,Numc,da,crpdat)
    USE parm, ONLY : KR, KI, Terr
    IMPLICIT NONE

    REAL (kind=KR), DIMENSION(21) :: crpdat
    CHARACTER (len=8), DIMENSION(5) :: da
    CHARACTER (len=10) :: nameC
    CHARACTER (LEN=100) :: datapath
    INTEGER (KIND=4) :: lenpath, i, numc

    lenpath = Len_Trim(datapath)
    OPEN(1,file=datapath(1:lenpath)//'/cropdat.dat',status="old",err=901)

    DO i=1,numc-1
        READ(1,*,err=901)
    ENDDO
    READ(1,*,err=901)nameC,da(1:5),crpdat(1:21)

    CLOSE(1)
    RETURN

901 Terr=2
    RETURN
    END SUBROUTINE cdatebase

! ====================================================================
!     subroutine refet
!
!     purpose: this is a subrutine to cacalate the reference et.the
!              equation comes from the PM equation(1992) recommended
!              by FAO(1998)
! ====================================================================
! =========================Incoming variables=========================
!   year        year[yyyy]
!   month       month[mm]
!   day         day[dd]
!   Tmean       the average temperature[C]
!   Tmax        the maximum temperature[C]
!   Tmin        the minimum temperature[C]
!   Nact        the number of hours of actual sunshine[h]
!   RHmean      the average humidty[%]
!   Uh          the speed of wind[m/s]
!   ato_p       the atmospheric pressure[Kpa]
!   hhh         the sea level elevation[m]
!   fai         the latitude[radian]
!   Hu          the Wind speed measuring height[m].
! =========================Outcoming variables========================
!   Et0         the reference evapotranspiration[m/d]
! =========================related files==============================
!   ??.wea      weather information for ET.
! =========================related functions==========================
!   None.
! ====================================================================
    SUBROUTINE Refet(year,month,day,Tmean,tmax,tmin,nact,RHmean,uh,hu,ato_p,hhh,fai,ET0,tmean1)
    USE parm, ONLY : KR, KI
    IMPLICIT NONE

    INTEGER (KIND=4), DIMENSION(12) :: Dm ! Calculation ordinal number
    INTEGER (KIND=4) :: i,year,month,day,julie
    REAL (KIND=KR) :: tmean,tmax,tmin,nact,rhmean,uh,ea,deta,dr,si,tansi, &
         tmean1,ws,ra,np,rns,eamax,eamin,ed,rnl,rn,g,lamta,u2,gama,ato_p,et0
    REAL (KIND=KR) :: fai, hhh, hu, amta

    DATA DM/0,31,59,90,120,151,181,212,243,273,304,334/

    julie=day+DM(MONTH)

    ea=0.611*EXP(17.27*tmean/(tmean+237.3))
    deta=4098*ea/(tmean+237.3)**2
    dr=1+0.033*COS(0.0172*julie)
    si=0.409*SIN(0.0172*julie-1.39)
    tansi=tan(si)
    ws=acos(-TAN(fai)*tansi)
    ra=37.6*dr*(ws*SIN(fai)*SIN(si)+COS(fai)*COS(si)*SIN(ws))
    Np=7.64*ws
    rns=0.77*(0.25+0.5*nact/Np)*ra
    eamax=0.611*EXP(17.27*tmax/(tmax+237.3))
    eamin=0.611*EXP(17.27*tmin/(tmin+237.3))
    ed=(eamax+eamin)/200*RHmean
    Rnl=2.45E-9*(0.9*nact/Np+0.1)*(0.34-0.14*ed**0.5)*((tmax+273)**4+(tmin+273)**4)
    rn=rns-rnl

    IF (tmean1.eq.0) THEN
        g=0.
    ELSE
        g=0.38*(tmean-tmean1)
    ENDIF
    tmean1=tmean
    amta=2.501-(2.361e-3)*tmean
    u2=4.87*uh/Log(67.8*hu-5.42)
    gama=0.00163*ato_p/lamta
    et0=(0.408*deta*(rn-G)+gama*900/(273+tmean)*u2*(ea-ed))/(deta+gama*(1+0.34*u2))
    et0=max(et0,0.0)    ! modified by MW 2016-09-26 17:36

    END SUBROUTINE Refet


text:
C    @(#)bus_hash.f	20.3 2/13/96
      integer function bus_hash (busname, basekv)

      character*(*)  busname
      real           basekv

c     compute hash value of busname

      include 'ipfinc/parametr.inc'

      include 'ipfinc/bushasht.inc'

      integer        k, i

      k = 0
      do i = 1, 8
         k = k + k + ichar (busname(i:i))
      end do
      k = k + int (basekv)
      bus_hash = mod (k, BUS_HASHSIZE) + 1

      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp