text:
subroutine shape1(s,nnode,shape,deriv)

!-----------------------------------------------------------------------
!
! This routine evaluates shape functions and their first derivates
! for 1-d continuos with 2 3 4 & 5 nodes
!
!-----------------------------------------------------------------------

  use      typre
  implicit none
  integer(ip), intent(in)  :: nnode
  real(rp),    intent(in)  :: s
  real(rp),    intent(out) :: deriv(nnode), shape(nnode)
  real                     :: s1, s2, s3, s4, a, c

  if(nnode==2) then
     shape(1)= 0.5_rp*(1.0_rp-s)              !  1    2
     shape(2)= 0.5_rp*(1.0_rp+s)
     deriv(1)=-0.5_rp
     deriv(2)= 0.5_rp
  else if(nnode==3) then
     shape(1)= 0.5_rp*s*(s-1.0_rp)            !  1    3    2
     shape(2)= 0.5_rp*s*(s+1.0_rp)
     shape(3)=-(s+1.0_rp)*(s-1.0_rp)
     deriv(1)= s-0.5_rp
     deriv(2)= s+0.5_rp
     deriv(3)=-2.0_rp*s
  else if(nnode==4) then
     a = 9.0_rp/16.0_rp
     s1= s + 1.0_rp
     s2= s - 1.0_rp
     s3= s + 1.0_rp/3.0_rp
     s4= s - 1.0_rp/3.0_rp
     c = 3.0_rp*a
     shape(1)= -a*s2*s3*s4                    !  1    3    4    2
     shape(2)=  a*s1*s3*s4
     shape(3)=  c*s1*s2*s4
     shape(4)= -c*s1*s2*s3
     deriv(1)= - a*(s2*s3+s2*s4+s3*s4)
     deriv(2)=   a*(s1*s3+s1*s4+s3*s4)
     deriv(3)=   c*(s1*s2+s1*s4+s2*s4)
     deriv(4)= - c*(s1*s2+s1*s3+s2*s3)
! * * * * * * * * * * * * * * * * * *
! *Valores calculados para nnode=5  *
! * * * * * * * * * * * * * * * * * *
  else if(nnode==5) then
     s1 = s + 1.0_rp
     s2 = s - 1.0_rp
     s3 = s + 1.0_rp/2.0_rp
     s4 = s - 1.0_rp/2.0_rp
     shape(1) = 2.0_rp/3.0_rp*s*s3*s4*s2     !  1    3    4    5    2
     shape(2) =-8.0_rp/3.0_rp*s*s1*s4*s2
     shape(3) = 4.0_rp*s1*s3*s4*s2
     shape(4) =-8.0_rp/3.0_rp*s*s1*s3*s2
     shape(5) = 2.0_rp/3.0_rp*s*s1*s3*s4
     deriv(1) = 2.0_rp/3.0_rp*(s*s3*s4+s*s3*s2+s*s4*s2+s3*s4*s2)
     deriv(2) =-8.0_rp/3.0_rp*(s*s1*s4+s*s1*s2+s*s4*s2+s1*s4*s2)
     deriv(3) =        4.0_rp*(s1*s3*s4+s1*s3*s2+s1*s4*s2+s3*s4*s2)
     deriv(4) =-8.0_rp/3.0_rp*(s*s1*s3+s*s1*s2+s*s3*s2+s1*s3*s2)
     deriv(5) = 2.0_rp/3.0_rp*(s*s1*s3+s*s1*s4+s*s3*s4+s1*s3*s4)

  end if
!  a= 0.0_rp
!  c= 0.0_rp
!  a = a + shape(1)+ shape(2)+ shape(3)+ shape(4)
!  c = c + deriv(1)+ deriv(2)+ deriv(3)+ deriv(4)

end subroutine shape1

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
!##############################################################################
!# ****************************************************************************
!# <name> postprocessing </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains postprocessing routines for the cc2dminim2_method2
!# CC2D solver.
!#
!# The following routines can be found here:
!#
!# 1.) cc_postprocessing
!#     -> Evaluate the solution of the stationary solver, write GMV-file.
!# </purpose>
!##############################################################################

module postprocessing

  use fsystem
  use storage
  use linearsolver
  use boundary
  use bilinearformevaluation
  use linearformevaluation
  use cubature
  use matrixfilters
  use vectorfilters
  use bcassembly
  use triangulation
  use spatialdiscretisation
  use coarsegridcorrection
  use spdiscprojection
  use nonlinearsolver
  use paramlist
  use pprocerror
  use pprocgradients

  use collection
  use convection

  use ucd

  use pprocnavierstokes

  use basicstructures
  use optcanalysis
  use user_callback
  use spatialbcdef
  use spacematvecassembly

  use spacetimediscretisation

  implicit none

!<types>

!<typeblock>

  ! Postprocessing structure. Whenever the postprocessing calculation routines are
  ! called, they fills this structure using the current solution vector (and other
  ! information if necessary). The information in this structure can then be used
  ! for GMV output e.g.
  type t_c2d2postprocessing

    ! A discretisation structure that describes a piecewise constant discretisation
    ! (usually P0 or Q0).
    type(t_spatialDiscretisation) :: rdiscrConstant

    ! A discretisation structure that describes a piecewise linear discretisation
    ! (usually P1 or Q1).
    type(t_spatialDiscretisation) :: rdiscrLinear

    ! A discretisation structure that describes a piecewise quadratic discretisation
    ! (usually P2 or Q2).
    type(t_spatialDiscretisation) :: rdiscrQuadratic

    ! A vector that describes the X-velocity field in the vertices
    type(t_vectorScalar) :: rvectorVelX

    ! A vector that describes the Y-velocity field in the vertices
    type(t_vectorScalar) :: rvectorVelY

    ! A vector that describes the pressure field in the vertices
    type(t_vectorScalar) :: rvectorPressure

    ! A vector that describes the pressure field in the cells
    type(t_vectorScalar) :: rvectorPressureCells

    ! A vector that describes the streamfunction
    type(t_vectorScalar) :: rvectorStreamfunction

    ! A vector that describes the H1-error of the velocity field in the vertices
    type(t_vectorScalar) :: rvectorH1err

    ! A vector that describes the H1-error of the pressure in the cells
    type(t_vectorScalar) :: rvectorH1errCells

  end type

!</typeblock>

!</types>

contains

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

!<subroutine>

  subroutine cc_postprocessingStationary (rproblem,rvector)

!<description>
  ! Postprocessing of solutions of stationary simulations.
  ! Writes the solution into a GMV file, calculates forces,...
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT), target :: rproblem
!</inputoutput>

!<input>
  ! The solution vector which is to be evaluated by the postprocessing routines.
  type(t_vectorBlock), intent(IN) :: rvector
!</input>

!</subroutine>

  ! local variables
    integer :: ieltype

    ! We need some more variables for postprocessing - i.e. writing
    ! a GMV file.
    real(DP), dimension(:), pointer :: p_Ddata,p_Ddata2

    ! A pointer to the triangulation.
    type(t_triangulation), pointer :: p_rtriangulation

    ! A vector accepting Q1 data
    type(t_vectorBlock) :: rprjVector

    ! A discretisation structure for Q1
    type(t_blockDiscretisation) :: rprjDiscretisation

    ! Discrete boundary conditions for the output vector
    type(t_discreteBC), target :: rdiscreteBC
    type(t_discreteFBC), target :: rdiscreteFBC

    ! Output block for UCD output to GMV file
    type(t_ucdExport) :: rexport

    ! Forces on the object
    real(DP), dimension(NDIM2D) :: Dforces
    real(DP) :: df1,df2
    type(t_boundaryRegion) :: rregion

    ! Divergence
    type(t_vectorScalar), target :: rtempVector

    ! If we have a uniform discreisation, calculate the body forces on the
    ! 2nd boundary component - if it exists.
    if ((rvector%p_rblockDiscr%RspatialDiscr(1)% &
         ccomplexity .eq. SPDISC_UNIFORM) .and. &
        (boundary_igetNBoundComp(rproblem%rboundary) .ge. 2)) then

      ! Calculate drag-/lift coefficients on the 2nd boundary component.
      ! This is for the benchmark channel!
      call boundary_createRegion (rproblem%rboundary, &
          2, 0, rregion)
      rregion%iproperties = BDR_PROP_WITHSTART+BDR_PROP_WITHEND
      df1 = 1.0_DP/1000.0_DP
      df2 = 0.1_DP * 0.2_DP**2
      call ppns2D_bdforces_uniform (rvector,rregion,Dforces,CUB_G1_1D,df1,df2)

      call output_lbrk()
      call output_line ('Body forces real bd., bdc/horiz/vert')
      call output_line (' 2 / ' &
          //trim(sys_sdEP(Dforces(1),15,6)) // ' / '&
          //trim(sys_sdEP(Dforces(2),15,6)) )

    endif

    ! Print out the value of the optimal control functional.
    call output_lbrk ()
    call cc_printControlFunctionalStat (rproblem,0.0_DP,rvector)

    ! If we have a simple Q1~ discretisation, calculate the streamfunction.
    if (rvector%p_rblockDiscr%RspatialDiscr(1)% &
        ccomplexity .eq. SPDISC_UNIFORM) then

      ieltype = rvector%p_rblockDiscr%RspatialDiscr(1)% &
                RelementDistr(1)%celement

      if (elem_getPrimaryElement(ieltype) .eq. EL_Q1T) then

        ! Create a temporary vector
        call lsyssc_createVecByDiscr (rvector%RvectorBlock(3)%p_rspatialDiscr,&
            rtempVector,.true.)

        ! Calculate divergence = B1^T u1 + B2^T u2
        call lsyssc_scalarMatVec (&
            rproblem%RlevelInfo(rproblem%nlmax)%rstaticInfo%rmatrixD1, rvector%RvectorBlock(1), &
            rtempVector, 1.0_DP, 0.0_DP)
        call lsyssc_scalarMatVec (&
            rproblem%RlevelInfo(rproblem%nlmax)%rstaticInfo%rmatrixD2, rvector%RvectorBlock(2), &
            rtempVector, 1.0_DP, 1.0_DP)

        call output_lbrk()
        call output_line ('Divergence = ' &
            //trim(sys_sdEP(lsyssc_vectorNorm(rtempVector,LINALG_NORML2),15,6)) )

        call lsyssc_releaseVector (rtempVector)

      end if

    end if

    ! 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 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(rvector%p_rblockDiscr,rprjDiscretisation)

    call spdiscr_deriveSimpleDiscrSc (&
                 rvector%p_rblockDiscr%RspatialDiscr(1), &
                 EL_Q1, CUB_G2X2, &
                 rprjDiscretisation%RspatialDiscr(1))

    call spdiscr_deriveSimpleDiscrSc (&
                 rvector%p_rblockDiscr%RspatialDiscr(2), &
                 EL_Q1, CUB_G2X2, &
                 rprjDiscretisation%RspatialDiscr(2))

    ! Also use Q1 for the dual velocity field.
    call spdiscr_deriveSimpleDiscrSc (&
                 rvector%p_rblockDiscr%RspatialDiscr(4), &
                 EL_Q1, CUB_G2X2, &
                 rprjDiscretisation%RspatialDiscr(4))

    call spdiscr_deriveSimpleDiscrSc (&
                 rvector%p_rblockDiscr%RspatialDiscr(5), &
                 EL_Q1, CUB_G2X2, &
                 rprjDiscretisation%RspatialDiscr(5))

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

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

    ! Discretise the boundary conditions for this discretisation
    call bcasm_initDiscreteBC(rdiscreteBC)
    call cc_assembleBDconditions (rproblem,0.0_DP,rprjDiscretisation,CCDISCBC_PRIMALDUAL,&
      rdiscreteBC,rproblem%rcollection)

    call bcasm_initDiscreteFBC(rdiscreteFBC)
    call cc_assembleFBDconditions (rproblem,0.0_DP,rprjDiscretisation,CCDISCBC_PRIMALDUAL,&
      rdiscreteFBC,rproblem%rcollection)

    ! Connect the vector with the BC's
    call lsysbl_assignDiscreteBC (rprjVector,rdiscreteBC)
    call lsysbl_assignDiscreteFBC (rprjVector,rdiscreteFBC)

    ! Filter the solution vector to implement discrete BC's.
    call vecfil_discreteBCsol (rprjVector)

    ! Filter the solution vector to implement discrete BC's for fictitious
    ! boundary components.
    call vecfil_discreteFBCsol (rprjVector)

    ! Clean up the collection (as we are done with the assembly.
    call cc_doneCollectForAssembly (rproblem,rproblem%rcollection)

    ! Now we have a Q1/Q1/Q0 solution in rprjVector.
    !
    ! From the attached discretisation, get the underlying triangulation
    p_rtriangulation => &
      rvector%RvectorBlock(1)%p_rspatialDiscr%p_rtriangulation

    ! Start UCD export to GMV file:
    call ucd_startGMV (rexport,UCD_FLAG_STANDARD,p_rtriangulation,'gmv/u.gmv')

    ! Write the configuration of the application as comment block
    ! to the output file.
    call ucd_addCommentLine (rexport,'Configuration:')
    call ucd_addCommentLine (rexport,'---------------')
    call ucd_addParameterList (rexport,rproblem%rparamList)
    call ucd_addCommentLine (rexport,'---------------')

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

    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)

    ! Dual velocity field
    call lsyssc_getbase_double (rprjVector%RvectorBlock(4),p_Ddata)
    call lsyssc_getbase_double (rprjVector%RvectorBlock(5),p_Ddata2)
    call ucd_addVariableVertexBased (rexport,'X-vel-dual',UCD_VAR_STANDARD, p_Ddata)
    call ucd_addVariableVertexBased (rexport,'Y-vel-dual',UCD_VAR_STANDARD, p_Ddata2)

    ! Dual pressure
    call lsyssc_getbase_double (rprjVector%RvectorBlock(6),p_Ddata)
    call ucd_addVariableElementBased (rexport,'pressure-dual',UCD_VAR_STANDARD, p_Ddata)

    ! If we have a simple Q1~ discretisation, calculate the streamfunction.
    if (rvector%p_rblockDiscr%RspatialDiscr(1)% &
        ccomplexity .eq. SPDISC_UNIFORM) then

      ieltype = rvector%p_rblockDiscr%RspatialDiscr(1)% &
                RelementDistr(1)%celement

      if (elem_getPrimaryElement(ieltype) .eq. EL_Q1T) then

        call ppns2D_streamfct_uniform (rvector,rprjVector%RvectorBlock(1))

        call lsyssc_getbase_double (rprjVector%RvectorBlock(1),p_Ddata)
        call ucd_addVariableVertexBased (rexport,'streamfunction',&
            UCD_VAR_STANDARD, p_Ddata)

      end if

    end if

    ! Calculate the derivative of the pressure.
    ! First, project the pressure from Q0 to Q1 (if it's Q0).
    if (rvector%p_rblockDiscr%RspatialDiscr(3)% &
        ccomplexity .eq. SPDISC_UNIFORM) then

      ieltype = rvector%p_rblockDiscr%RspatialDiscr(3)% &
                RelementDistr(1)%celement

      if (elem_getPrimaryElement(ieltype) .eq. EL_Q0) then
        ! rprjVector%RvectorBlock(4) is already prepared to accept Q1 solutions
        call spdp_projectSolutionScalar (rvector%RvectorBlock(3),rprjVector%RvectorBlock(4))

        ! Calculate the derivative of the pressure into the first two
        ! subvectors -- overwriting the projected velocity field, which we don't need
        ! anymore.
        call ppgrd_calcGradient (rprjVector%RvectorBlock(4),rprjVector)

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

        call ucd_addVariableVertexBased (rexport,'X-deriv-pres',UCD_VAR_STANDARD, p_Ddata)
        call ucd_addVariableVertexBased (rexport,'Y-deriv-pres',UCD_VAR_STANDARD, p_Ddata2)

      end if

      ieltype = rvector%p_rblockDiscr%RspatialDiscr(6)% &
                RelementDistr(1)%celement

      if (elem_getPrimaryElement(ieltype) .eq. EL_Q0) then
        ! rprjVector%RvectorBlock(4) is already prepared to accept Q1 solutions
        call spdp_projectSolutionScalar (rvector%RvectorBlock(6),rprjVector%RvectorBlock(4))

        ! Calculate the derivative of the dual pressure into the first two
        ! subvectors -- overwriting the projected velocity field, which we don't need
        ! anymore.
        call ppgrd_calcGradient (rprjVector%RvectorBlock(4),rprjVector)

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

        call ucd_addVariableVertexBased (rexport,'X-deriv-dpres',UCD_VAR_STANDARD, p_Ddata)
        call ucd_addVariableVertexBased (rexport,'Y-deriv-dpres',UCD_VAR_STANDARD, p_Ddata2)

      end if

    end if

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

    ! Release the auxiliary vector
    call lsysbl_releaseVector (rprjVector)

    ! Release the discretisation strucutre
    call spdiscr_releaseBlockDiscr(rprjDiscretisation)

    ! Throw away the discrete BC's - not used anymore.
    call bcasm_releaseDiscreteBC (rdiscreteBC)
    call bcasm_releaseDiscreteFBC (rdiscreteFBC)

  end subroutine

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

!<subroutine>

  subroutine cc_postprocessingNonstat (rproblem,itimestep,dtime,rvector)

!<description>
  ! Postprocessing of solutions of stationary simulations.
  ! Writes the solution into a GMV file, calculates forces,...
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT), target :: rproblem
!</inputoutput>

!<input>
  ! The solution vector which is to be evaluated by the postprocessing routines.
  type(t_vectorBlock), intent(IN) :: rvector

  ! Number of the current timestep.
  integer, intent(in) :: itimestep

  ! Current simulation time.
  real(dp), intent(in) :: dtime
!</input>

!</subroutine>

  ! local variables
    integer :: ieltype

    ! We need some more variables for postprocessing - i.e. writing
    ! a GMV file.
    real(DP), dimension(:), pointer :: p_Ddata,p_Ddata2

    ! A pointer to the triangulation.
    type(t_triangulation), pointer :: p_rtriangulation

    ! A vector accepting Q1 data
    type(t_vectorBlock) :: rprjVector

    ! A discretisation structure for Q1
    type(t_blockDiscretisation) :: rprjDiscretisation

    ! Discrete boundary conditions for the output vector
    type(t_discreteBC), target :: rdiscreteBC
    type(t_discreteFBC), target :: rdiscreteFBC

    ! Output block for UCD output to GMV file
    type(t_ucdExport) :: rexport

    ! Forces on the object
    real(DP), dimension(NDIM2D) :: Dforces
    real(DP) :: df1,df2
    type(t_boundaryRegion) :: rregion

    ! Divergence
    type(t_vectorScalar), target :: rtempVector

    character(SYS_STRLEN) :: sgmvName,stemp

    ! If we have a uniform discretisation, calculate the body forces on the
    ! 2nd boundary component - if it exists.
    if ((rvector%p_rblockDiscr%RspatialDiscr(1)% &
         ccomplexity .eq. SPDISC_UNIFORM) .and. &
        (boundary_igetNBoundComp(rproblem%rboundary) .ge. 2)) then

      ! Calculate drag-/lift coefficients on the 2nd boundary component.
      ! This is for the benchmark channel!
      call boundary_createRegion (rproblem%rboundary, &
          2, 0, rregion)
      rregion%iproperties = BDR_PROP_WITHSTART+BDR_PROP_WITHEND
      df1 = 1.0_DP/1000.0_DP
      df2 = 0.1_DP * 0.2_DP**2
      call ppns2D_bdforces_uniform (rvector,rregion,Dforces,CUB_G1_1D,df1,df2)

      call output_lbrk()
      call output_line ('Body forces real bd., bdc/horiz/vert')
      call output_line (' 2 / ' &
          //trim(sys_sdEP(Dforces(1),15,6)) // ' / '&
          //trim(sys_sdEP(Dforces(2),15,6)) )

    endif

    ! If we have a simple Q1~ discretisation, calculate the streamfunction.
!    IF (rvector%p_rblockDiscretisation%RspatialDiscr(1)% &
!        ccomplexity .EQ. SPDISC_UNIFORM) THEN
!
!      ieltype = rvector%p_rblockDiscretisation%RspatialDiscr(1)% &
!                RelementDistr(1)%itrialElement
!
!      IF (elem_getPrimaryElement(ieltype) .EQ. EL_Q1T) THEN
!
!        ! Create a temporary vector
!        CALL lsyssc_createVecByDiscr (rvector%RvectorBlock(3)%p_rspatialDiscretisation,&
!            rtempVector,.TRUE.)
!
!        ! Calculate divergence = B1^T u1 + B2^T u2
!        CALL lsyssc_transposeMatrix (rproblem%RlevelInfo(rproblem%nlmax)%rmatrixB1,&
!            rBmatrix,LSYSSC_TR_VIRTUAL)
!        CALL lsyssc_scalarMatVec (&
!            rBmatrix, rvector%RvectorBlock(1), &
!            rtempVector, 1.0_DP, 0.0_DP)
!        CALL lsyssc_transposeMatrix (rproblem%RlevelInfo(rproblem%nlmax)%rmatrixB2,&
!            rBmatrix,LSYSSC_TR_VIRTUAL)
!        CALL lsyssc_scalarMatVec (&
!            rBmatrix, rvector%RvectorBlock(2), &
!            rtempVector, 1.0_DP, 1.0_DP)
!
!        CALL output_lbrk()
!        CALL output_line ('Divergence = ' &
!            //TRIM(sys_sdEP(lsyssc_vectorNorm(rtempVector,LINALG_NORML2),15,6)) )
!
!        CALL lsyssc_releaseVector (rtempVector)
!
!      END IF
!
!    END IF

    ! 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 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(rvector%p_rblockDiscr,&
        rprjDiscretisation)

    call spdiscr_deriveSimpleDiscrSc (&
                 rvector%p_rblockDiscr%RspatialDiscr(1), &
                 EL_Q1, CUB_G2X2, &
                 rprjDiscretisation%RspatialDiscr(1))

    call spdiscr_deriveSimpleDiscrSc (&
                 rvector%p_rblockDiscr%RspatialDiscr(2), &
                 EL_Q1, CUB_G2X2, &
                 rprjDiscretisation%RspatialDiscr(2))

    ! Also use Q1 for the dual velocity field.
    call spdiscr_deriveSimpleDiscrSc (&
                 rvector%p_rblockDiscr%RspatialDiscr(4), &
                 EL_Q1, CUB_G2X2, &
                 rprjDiscretisation%RspatialDiscr(4))

    call spdiscr_deriveSimpleDiscrSc (&
                 rvector%p_rblockDiscr%RspatialDiscr(5), &
                 EL_Q1, CUB_G2X2, &
                 rprjDiscretisation%RspatialDiscr(5))

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

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

    call cc_initCollectForAssembly(rproblem,dtime,rproblem%rcollection)

    ! Discretise the boundary conditions according to the Q1/Q1/Q0
    ! discretisation for implementing them into a solution vector.
    ! Discretise the boundary conditions for this discretisation
    call bcasm_initDiscreteBC(rdiscreteBC)
    call cc_assembleBDconditions (rproblem,dtime,rprjDiscretisation,CCDISCBC_PRIMALDUAL,&
      rdiscreteBC,rproblem%rcollection)

    call bcasm_initDiscreteFBC(rdiscreteFBC)
    call cc_assembleFBDconditions (rproblem,dtime,rprjDiscretisation,CCDISCBC_PRIMALDUAL,&
      rdiscreteFBC,rproblem%rcollection)

    call cc_doneCollectForAssembly(rproblem,rproblem%rcollection)

    ! Connect the vector with the BC's
    call lsysbl_assignDiscreteBC (rprjVector,rdiscreteBC)
    call lsysbl_assignDiscreteFBC (rprjVector,rdiscreteFBC)

    ! Filter the solution vector to implement discrete BC's.
    call vecfil_discreteBCsol (rprjVector)
    !DEBUG: CALL vecfil_discreteBCdef (rprjVector)

    ! Filter the solution vector to implement discrete BC's for fictitious
    ! boundary components.
    call vecfil_discreteFBCsol (rprjVector)

    ! Now we have a Q1/Q1/Q0 solution in rprjVector.
    !
    ! From the attached discretisation, get the underlying triangulation
    p_rtriangulation => &
      rvector%RvectorBlock(1)%p_rspatialDiscr%p_rtriangulation

    ! Get GMV filename
    call parlst_getvalue_string (rproblem%rparamList, 'TIME-POSTPROCESSING', &
                                 'sgmvFileName', stemp, 'gmv/u.gmv')
    read(stemp,*) sgmvName

    if (sgmvName .ne. '') then

      ! Start UCD export to GMV file:
      call output_lbrk ()
      call output_line ('Writing GMV file: ' &
          //trim(sgmvName)//'.'//sys_si0(itimeStep,5))

      call ucd_startGMV (rexport,UCD_FLAG_STANDARD,p_rtriangulation,&
          trim(sgmvName)//'.'//sys_si0(itimeStep,5))

      ! Write the configuration of the application as comment block
      ! to the output file.
      call ucd_addCommentLine (rexport,'Configuration:')
      call ucd_addCommentLine (rexport,'---------------')
      call ucd_addParameterList (rexport,rproblem%rparamList)
      call ucd_addCommentLine (rexport,'---------------')

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

      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)

      ! Dual velocity field
      call lsyssc_getbase_double (rprjVector%RvectorBlock(4),p_Ddata)
      call lsyssc_getbase_double (rprjVector%RvectorBlock(5),p_Ddata2)
      call ucd_addVariableVertexBased (rexport,'X-vel-dual',UCD_VAR_STANDARD, p_Ddata)
      call ucd_addVariableVertexBased (rexport,'Y-vel-dual',UCD_VAR_STANDARD, p_Ddata2)

      ! Dual pressure
      call lsyssc_getbase_double (rprjVector%RvectorBlock(6),p_Ddata)
      call ucd_addVariableElementBased (rexport,'pressure-dual',UCD_VAR_STANDARD, p_Ddata)

      ! Control u = P[min/max](-1/alpha lambda)
      call lsyssc_getbase_double (rprjVector%RvectorBlock(4),p_Ddata)
      call lsyssc_scaleVector (rprjVector%RvectorBlock(4),-1.0_DP/rproblem%roptControl%dalphaC)
      if (rproblem%roptControl%ccontrolConstraints .ne. 0) then
        call cc_projectControlTimestep (rprjVector%RvectorBlock(4),&
          rproblem%roptControl%dumin1,rproblem%roptControl%dumax1)
      end if
      call ucd_addVariableVertexBased (rexport,'X-control',UCD_VAR_STANDARD, p_Ddata)

      call lsyssc_getbase_double (rprjVector%RvectorBlock(5),p_Ddata2)
      call lsyssc_scaleVector (rprjVector%RvectorBlock(5),-1.0_DP/rproblem%roptControl%dalphaC)
      if (rproblem%roptControl%ccontrolConstraints .ne. 0) then
        call cc_projectControlTimestep (rprjVector%RvectorBlock(5),&
          rproblem%roptControl%dumin2,rproblem%roptControl%dumax2)
      end if
      call ucd_addVariableVertexBased (rexport,'Y-control',UCD_VAR_STANDARD, p_Ddata2)

      ! If we have a simple Q1~ discretisation, calculate the streamfunction.
      if (rvector%p_rblockDiscr%RspatialDiscr(1)% &
          ccomplexity .eq. SPDISC_UNIFORM) then

        ieltype = rvector%p_rblockDiscr%RspatialDiscr(1)% &
                  RelementDistr(1)%celement

        if (elem_getPrimaryElement(ieltype) .eq. EL_Q1T) then

          call ppns2D_streamfct_uniform (rvector,rprjVector%RvectorBlock(1))

          call lsyssc_getbase_double (rprjVector%RvectorBlock(1),p_Ddata)
          call ucd_addVariableVertexBased (rexport,'streamfunction',&
              UCD_VAR_STANDARD, p_Ddata)

        end if

      end if

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

    end if

    ! Release the auxiliary vector
    call lsysbl_releaseVector (rprjVector)

    ! Release the discretisations structure
    call spdiscr_releaseBlockDiscr (rprjDiscretisation)

    ! Throw away the discrete BC's - not used anymore.
    call bcasm_releaseDiscreteBC (rdiscreteBC)
    call bcasm_releaseDiscreteFBC (rdiscreteFBC)

  end subroutine

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

!<subroutine>

  subroutine cc_postprocSpaceTimeGMV (rproblem,rdiscr,rvector,sfilename)

!<description>
  ! For every sub-solution in the global space-time vector rvector,
  ! a GMV file is written to disc.
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT), target :: rproblem
!</inputoutput>

!<input>
  ! A space-time discretisation structure defining the discretisation of
  ! rvector.
  type(t_ccoptSpaceTimeDiscretisation), intent(IN) :: rdiscr

  ! A space-time vector. For every timestep, a GMV is written.
  type(t_spaceTimeVector), intent(IN) :: rvector

  ! A path + basic filename for the GMV-files. A number '.00000','.00001',...
  ! is appended for every timestep.
  character(LEN=*), intent(IN) :: sfilename
!</input>

!</subroutine>

    ! local variables
    type(t_vectorBlock) :: rvectorTmp
    integer :: i,ieltype
    real(dp) :: dtime

    ! We need some more variables for postprocessing - i.e. writing
    ! a GMV file.
    real(DP), dimension(:), pointer :: p_Ddata,p_Ddata2

    ! A pointer to the triangulation.
    type(t_triangulation), pointer :: p_rtriangulation

    ! A vector accepting Q1 data
    type(t_vectorBlock) :: rprjVector

    ! A discretisation structure for Q1
    type(t_blockDiscretisation) :: rprjDiscretisation

    ! Discrete boundary conditions for the output vector
    type(t_discreteBC), target :: rdiscreteBC
    type(t_discreteFBC), target :: rdiscreteFBC

    ! Output block for UCD output to GMV file
    type(t_ucdExport) :: rexport

    ! Create a temp vector
    call lsysbl_createVecBlockByDiscr (&
        rdiscr%p_rlevelInfo%rdiscretisation,&
        rvectorTmp,.true.)

    ! Postprocessing of all solution vectors.
    do i = 0,rdiscr%rtimeDiscr%nintervals

      dtime = rproblem%rtimedependence%dtimeInit + i*rdiscr%rtimeDiscr%dtstep

      ! Evaluate the space time function in rvector in the point
      ! in time dtime. Independent of the discretisation in time,
      ! this will give us a vector in space.
      !CALL sptivec_getTimestepData (rvector, i, rvectorTmp)
      call tmevl_evaluate(rvector,dtime,rvectorTmp)

      ! 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 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(rvectorTmp%p_rblockDiscr,&
          rprjDiscretisation)

      call spdiscr_deriveSimpleDiscrSc (&
                  rvectorTmp%p_rblockDiscr%RspatialDiscr(1), &
                  EL_Q1, CUB_G2X2, &
                  rprjDiscretisation%RspatialDiscr(1))

      call spdiscr_deriveSimpleDiscrSc (&
                  rvectorTmp%p_rblockDiscr%RspatialDiscr(2), &
                  EL_Q1, CUB_G2X2, &
                  rprjDiscretisation%RspatialDiscr(2))

      ! Also use Q1 for the dual velocity field.
      call spdiscr_deriveSimpleDiscrSc (&
                  rvectorTmp%p_rblockDiscr%RspatialDiscr(4), &
                  EL_Q1, CUB_G2X2, &
                  rprjDiscretisation%RspatialDiscr(4))

      call spdiscr_deriveSimpleDiscrSc (&
                  rvectorTmp%p_rblockDiscr%RspatialDiscr(5), &
                  EL_Q1, CUB_G2X2, &
                  rprjDiscretisation%RspatialDiscr(5))

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

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

      call cc_initCollectForAssembly(rproblem,dtime,rproblem%rcollection)

      ! Discretise the boundary conditions for this discretisation
      call bcasm_initDiscreteBC(rdiscreteBC)
      call cc_assembleBDconditions (rproblem,dtime,rprjDiscretisation,CCDISCBC_PRIMALDUAL,&
        rdiscreteBC,rproblem%rcollection)

      call bcasm_initDiscreteFBC(rdiscreteFBC)
      call cc_assembleFBDconditions (rproblem,dtime,rprjDiscretisation,CCDISCBC_PRIMALDUAL,&
        rdiscreteFBC,rproblem%rcollection)

      ! Connect the vector with the BC's
      call lsysbl_assignDiscreteBC (rprjVector,rdiscreteBC)
      call lsysbl_assignDiscreteFBC (rprjVector,rdiscreteFBC)

      call cc_doneCollectForAssembly(rproblem,rproblem%rcollection)

      ! Filter the solution vector to implement discrete BC's.
      call vecfil_discreteBCsol (rprjVector)
      !DEBUG: CALL vecfil_discreteBCdef (rprjVector)

      ! Filter the solution vector to implement discrete BC's for fictitious
      ! boundary components.
      call vecfil_discreteFBCsol (rprjVector)

      ! Now we have a Q1/Q1/Q0 solution in rprjVector.
      !
      ! From the attached discretisation, get the underlying triangulation
      p_rtriangulation => &
        rvectorTmp%rvectorBlock(1)%p_rspatialDiscr%p_rtriangulation

      ! Start UCD export to GMV file:
      call output_line ('Writing GMV file: ' &
          //trim(sfilename)//'.'//sys_si0(i,5))

      call ucd_startGMV (rexport,UCD_FLAG_STANDARD,p_rtriangulation,&
          trim(sfilename)//'.'//sys_si0(i,5))

      ! Write the configuration of the application as comment block
      ! to the output file.
      call ucd_addCommentLine (rexport,'Configuration:')
      call ucd_addCommentLine (rexport,'---------------')
      call ucd_addParameterList (rexport,rproblem%rparamList)
      call ucd_addCommentLine (rexport,'---------------')

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

      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)

      ! Dual velocity field
      call lsyssc_getbase_double (rprjVector%rvectorBlock(4),p_Ddata)
      call lsyssc_getbase_double (rprjVector%rvectorBlock(5),p_Ddata2)
      call ucd_addVariableVertexBased (rexport,'X-vel-dual',UCD_VAR_STANDARD, p_Ddata)
      call ucd_addVariableVertexBased (rexport,'Y-vel-dual',UCD_VAR_STANDARD, p_Ddata2)

      ! Dual pressure
      call lsyssc_getbase_double (rprjVector%rvectorBlock(6),p_Ddata)
      call ucd_addVariableElementBased (rexport,'pressure-dual',UCD_VAR_STANDARD, p_Ddata)

      ! If we have a simple Q1~ discretisation, calculate the streamfunction.
      if (rvectorTmp%p_rblockDiscr%RspatialDiscr(1)% &
          ccomplexity .eq. SPDISC_UNIFORM) then

        ieltype = rvectorTmp%p_rblockDiscr%RspatialDiscr(1)% &
                  RelementDistr(1)%celement

        if (elem_getPrimaryElement(ieltype) .eq. EL_Q1T) then

          call ppns2D_streamfct_uniform (rvectorTmp,rprjVector%rvectorBlock(1))

          call lsyssc_getbase_double (rprjVector%rvectorBlock(1),p_Ddata)
          call ucd_addVariableVertexBased (rexport,'streamfunction',&
              UCD_VAR_STANDARD, p_Ddata)

        end if

      end if

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

      ! Release the auxiliary vector
      call lsysbl_releaseVector (rprjVector)

      ! Release the discretisation structure
      call spdiscr_releaseBlockDiscr(rprjDiscretisation)

      ! Throw away the discrete BC's - not used anymore.
      call bcasm_releaseDiscreteBC (rdiscreteBC)
      call bcasm_releaseDiscreteFBC (rdiscreteFBC)

    end do

    ! Release memory.
    call lsysbl_releaseVector (rvectorTmp)

  end subroutine

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

!<subroutine>

  subroutine cc_initpostprocessing (rproblem,rpostprocessing)

!<description>
  ! Initialises the given postprocessing structure rpostprocessing
  ! according to the main problem rproblem. The structure can then be used
  ! to generate postprocessing data.
!</description>

!<input>
  ! A problem structure that describes the main problem to solve.
  type(t_problem), intent(IN),target :: rproblem
!</input>

!<output>
  type(t_c2d2postprocessing), intent(OUT) :: rpostprocessing
!</output>

!</subroutine>

  ! local variables
  type(t_blockDiscretisation), pointer :: p_rdiscr

    ! For postprocessing, we need discretisation structures in the Q0 and Q1 space,
    ! later perhaps in the Q2 space. For this purpose, derive the corresponding
    ! discretisation structure using the 'main' discretisation structure on the
    ! maximum level.
    !
    ! For simplicity, we use only the discretisation structure of the X-velocity
    ! to derive everything.

    p_rdiscr => rproblem%RlevelInfo(rproblem%NLMAX)%rdiscretisation

    ! Piecewise constant space:
    call spdiscr_deriveSimpleDiscrSc (&
                 p_rdiscr%RspatialDiscr(1), &
                 EL_Q0, CUB_G1X1, &
                 rpostprocessing%rdiscrConstant)

    ! Piecewise linear space:
    call spdiscr_deriveSimpleDiscrSc (&
                 p_rdiscr%RspatialDiscr(1), &
                 EL_Q1, CUB_G2X2, &
                 rpostprocessing%rdiscrLinear)

    ! Piecewise quadratic space:
    call spdiscr_deriveSimpleDiscrSc (&
                 p_rdiscr%RspatialDiscr(1), &
                 EL_Q2, CUB_G3X3, &
                 rpostprocessing%rdiscrQuadratic)

  end subroutine

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

!<subroutine>

  subroutine cc_clearpostprocessing (rpostprocessing)

!<description>
  ! Releases all calculated data in the given postprocessing structure
  ! so that it can be allocated again in the calculation routines.
  ! This routine must be called at the end of the postprocessing routines
  ! to release the temporary memory that was allocated for the vectors
  ! in the postprocessing structure.
!</description>

!<inputoutput>
  type(t_c2d2postprocessing), intent(INOUT) :: rpostprocessing
!</inputoutput>

!</subroutine>

    ! Release all vectors which might be allocated.
    call lsyssc_releaseVector (rpostprocessing%rvectorVelX)
    call lsyssc_releaseVector (rpostprocessing%rvectorVelY)
    call lsyssc_releaseVector (rpostprocessing%rvectorPressure)
    call lsyssc_releaseVector (rpostprocessing%rvectorPressureCells)
    call lsyssc_releaseVector (rpostprocessing%rvectorStreamfunction)
    call lsyssc_releaseVector (rpostprocessing%rvectorH1err)
    call lsyssc_releaseVector (rpostprocessing%rvectorH1errCells)

  end subroutine

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

!<subroutine>

  subroutine cc_donepostprocessing (rpostprocessing)

!<description>
  ! Releases a given problem structure. All allocated memory of this structure
  ! is released.
!</description>

!<inputoutput>
  type(t_c2d2postprocessing), intent(INOUT) :: rpostprocessing
!</inputoutput>

!</subroutine>

    ! Release all vectors -- if there are still some allocated
    call cc_clearpostprocessing (rpostprocessing)

    ! Release the discretisation structures allocated above.
    call spdiscr_releaseDiscr(rpostprocessing%rdiscrQuadratic)
    call spdiscr_releaseDiscr(rpostprocessing%rdiscrLinear)
    call spdiscr_releaseDiscr(rpostprocessing%rdiscrConstant)

  end subroutine

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

!<subroutine>

  subroutine cc_printControlFunctionalStat (rproblem,dtime,rvector)

!<description>
  ! Calculates and prints the value of the optimal control functional J(y,u)
  ! in the stationary case.
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT), target :: rproblem
!</inputoutput>

!<input>
  ! The solution vector which is to be evaluated by the postprocessing routines.
  type(t_vectorBlock), intent(IN) :: rvector

  ! Current simulation time.
  real(dp), intent(in) :: dtime
!</input>

!</subroutine>

    ! local variables
    real(DP), dimension(3) :: Derror
    real(DP) :: dalphaC

    ! Initialise the collection for the assembly process with callback routines.
    ! Basically, this stores the simulation time in the collection if the
    ! simulation is nonstationary.
    call cc_initCollectForAssembly (rproblem,dtime,rproblem%rcollection)

    ! Analyse the deviation from the target velocity field.
    call parlst_getvalue_double (rproblem%rparamList,'OPTIMALCONTROL',&
                                'dalphaC',dalphaC,0.1_DP)
    call cc_optc_stationaryFunctional (rvector,dalphaC,Derror,&
        rproblem%rcollection)

    call output_line ('||y-z||_L2: '//trim(sys_sdEL(Derror(1),2)))
    call output_line ('||u||_L2  : '//trim(sys_sdEL(Derror(2),2)))
    call output_line ('J(y,u)    : '//trim(sys_sdEL(Derror(3),2)))

    ! Release assembly-stuff from the collection
    call cc_doneCollectForAssembly (rproblem,rproblem%rcollection)

  end subroutine

!  !****************************************************************************
!
!!<subroutine>
!
!  SUBROUTINE ppns2D_fbdforces (rvector,rregion,Dforces,cformulation,CcubU,CcubP)
!
!!<description>
!  ! Calculates the drag-/lift-forces acting on a part of the real
!  ! boundary for a vector rvector with the solution of the 2D
!  ! (Navier-)Stokes equation. It's assumed that
!  !  rvector%rvectorBlock(1) = X-velocity,
!  !  rvector%rvectorBlock(2) = Y-velocity,
!  !  rvector%rvectorBlock(3) = pressure,
!  ! and that X- and Y-velocity is discretised with the same
!  ! finite element.
!  ! rregion specifies a boundary region where to calculate
!  ! the force integral. Dforces(1:NDIM2D) receives the forces in the
!  ! X- and Y-direction.
!  !
!  ! Double precision version
!!</description>
!
!!<input>
!  ! The FE solution vector.
!  TYPE(t_vectorBlock), INTENT(INOUT) :: rvector
!
!  ! Boundary region where to calculate the boundary forces.
!  ! Can be created e.g. by boundary_createRegion.
!  TYPE(t_boundaryRegion), INTENT(OUT) :: rregion
!
!  ! Type identifier that specifies the formulation of the Navier Stokes
!  ! equation that should be used for performing the integration.
!  ! ... = gradient formulation,
!  ! ... = deformation formulation
!  INTEGER, INTENT(IN) :: cformulation
!
!  ! OPTIONAL: Array with cubature formula identifiers when integrating
!  ! the velocity. If specified, there must be a cubature formula identifier
!  ! for each of the element distributions specified in the discretisation
!  ! structure that is associated to rvector.
!  ! If not specified, the standard cubature formula in the discretisation
!  ! structure for integrating linear forms is used.
!  INTEGER, DIMENSION(:), INTENT(IN) :: CcubU
!
!  ! OPTIONAL: Array with cubature formula identifiers when integrating
!  ! the pressure. If specified, there must be a cubature formula identifier
!  ! for each of the element distributions specified in the discretisation
!  ! structure that is associated to rvector.
!  ! If not specified, the standard cubature formula in the discretisation
!  ! structure for integrating linear forms is used.
!  INTEGER, DIMENSION(:), INTENT(IN) :: CcubP
!!</input>
!
!!<output>
!  ! Array receiving the forces acting on the boundary specified by rregion.
!  ! Note: These are the drag-/lift-FORCES, not the coefficients!!!
!  REAL(DP), DIMENSION(:), INTENT(OUT) :: Dforces
!!</output>
!
!!</subroutine>
!
!  ! Array to tell the element which derivatives to calculate.
!  LOGICAL, DIMENSION(EL_MAXNDER) :: BderU, BderP
!
!  ! Cubature point coordinates on the reference element.
!  REAL(DP), DIMENSION(CUB_MAXCUBP, NDIM3D) :: DxiU, DxiP
!
!  ! For every cubature point on the reference element,
!  ! the corresponding cubature weight
!  REAL(DP), DIMENSION(CUB_MAXCUBP) :: DomegaU, DomegaP
!
!  ! number of cubature points on the reference element
!  INTEGER :: ncubpU, ncubpP
!
!  ! Pointer to the vector entries
!  REAL(DP), DIMENSION(:), POINTER :: p_DdataUX, p_DdataUY, p_DdataP
!
!  ! An allocateable array accepting the DOF's of a set of elements.
!  integer, DIMENSION(:,:), ALLOCATABLE, TARGET :: IdofsTrialU, IdofsTrialP
!
!  ! Allocateable arrays for the values of the basis functions.
!  REAL(DP), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: DbasTrialU,DbasTrialP
!
!  ! Number of entries in the vector - for quicker access
!  INTEGER(I32) :: neqU, neqP
!
!  ! Arrays for saving Jacobian determinants and matrices
!  REAL(DP), DIMENSION(:,:), POINTER :: p_Ddetj
!  REAL(DP), DIMENSION(:,:,:), POINTER :: p_Djac
!
!  ! Pointer to KVERT of the triangulation
!  INTEGER(I32), DIMENSION(:,:), POINTER :: p_IverticesAtElement
!
!  ! Pointer to DCORVG of the triangulation
!  REAL(DP), DIMENSION(:,:), POINTER :: p_DvertexCoords
!
!  ! Current element distribution
!  TYPE(t_elementDistribution), POINTER :: p_elementDistribution
!
!  ! Number of elements in a block. Normally =BILF_NELEMSIM,
!  ! except if there are less elements in the discretisation.
!  INTEGER :: nelementsPerBlock
!
!
!  BderU = .FALSE.
!  BderP = .FALSE.
!
!  ! Get the vector data
!  neqU = rvector%RvectorBlock(1)%NEQ
!  neqP = rvector%RvectorBlock(3)%NEQ
!
!  IF (rvector%cdataType .NE. ST_DOUBLE) THEN
!    PRINT *,'ppns2D_fbdforces: Unsupported vector precision.'
!    STOP
!  END IF
!
!  ! We support only uniform and conformal discretisation structures.
!  IF (.NOT. ASSOCIATED(rvector%p_rblockDiscretisation)) THEN
!    PRINT *,'ppns2D_fbdforces: No discretisation structure!'
!    STOP
!  END IF
!
!  IF ((rvector%p_rblockDiscretisation%ccomplexity .NE. SPDISC_UNIFORM) .AND. &
!      (rvector%p_rblockDiscretisation%ccomplexity .NE. SPDISC_CONFORMAL)) THEN
!    PRINT *,'ppns2D_fbdforces: Discretisation too complex!'
!    STOP
!  END IF
!
!  ! Get pointers to the subvectors from the block vector
!  CALL lsyssc_getbase_double (rvector%RvectorBlock(1),p_DdataUX)
!  CALL lsyssc_getbase_double (rvector%RvectorBlock(2),p_DdataUY)
!  CALL lsyssc_getbase_double (rvector%RvectorBlock(3),p_DdataP)
!
!  ! Now loop over the different element distributions in the discretisation.
!  ! Each element distribution has a different element combination.
!  ! E.g. distribution 1 may describe P1~/P1~/P0 triangular elements, while
!  ! distrbution 2 may describe Q1~/Q1~/Q0.
!
!  !DO icurrentElementDistr = 1,rdiscretisation%inumFESpaces
!
!
!  !END DO ! icurrentElementDistr
!
!  Dforces = 0.0_DP
!
!  END SUBROUTINE
!
!
!  !****************************************************************************
!
!!<subroutine>
!
!  SUBROUTINE ppns2D_bdforces (rvector,rregion,Dforces,cformulation,CcubU,CcubP)
!
!!<description>
!  ! Calculates the drag-/lift-forces acting on a part of the real
!  ! boundary for a vector rvector with the solution of the 2D
!  ! (Navier-)Stokes equation. It's assumed that
!  !  rvector%rvectorBlock(1) = X-velocity,
!  !  rvector%rvectorBlock(2) = Y-velocity,
!  !  rvector%rvectorBlock(3) = pressure,
!  ! and that X- and Y-velocity is discretised with the same
!  ! finite element.
!  ! rregion specifies a boundary region where to calculate
!  ! the force integral. Dforces(1:NDIM2D) receives the forces in the
!  ! X- and Y-direction.
!  !
!  ! Double precision version
!!</description>
!
!!<input>
!  ! The FE solution vector.
!  TYPE(t_vectorBlock), INTENT(INOUT) :: rvector
!
!  ! Boundary region where to calculate the boundary forces.
!  ! Can be created e.g. by boundary_createRegion.
!  TYPE(t_boundaryRegion), INTENT(OUT) :: rregion
!
!  ! Type identifier that specifies the formulation of the Navier Stokes
!  ! equation that should be used for performing the integration.
!  ! ... = gradient formulation,
!  ! ... = deformation formulation
!  INTEGER, INTENT(IN) :: cformulation
!
!  ! 1D Cubature formula identifier to use for the line integration.
!  ! One of the CUB_xxxx_1D constants in the cubature.f90.
!  INTEGER, DIMENSION(:), INTENT(IN) :: CcubU
!!</input>
!
!!<output>
!  ! Array receiving the forces acting on the boundary specified by rregion.
!  ! Note: These are the drag-/lift-FORCES, not the coefficients!!!
!  REAL(DP), DIMENSION(:), INTENT(OUT) :: Dforces
!!</output>
!
!!</subroutine>
!
!  ! Spatial discretisation structure of velocity and pressure
!  TYPE(t_spatialDiscretisation), POINTER :: p_rdiscrU, p_rdiscrP
!
!  ! Element distribution of velocity and pressure
!  TYPE(t_elementDistribution), POINTER :: p_elemDistrU,p_elemDistrP
!
!
!  ! Array to tell the element which derivatives to calculate.
!  LOGICAL, DIMENSION(EL_MAXNDER) :: BderU, BderP
!
!  ! Cubature point coordinates on the reference element.
!  REAL(DP), DIMENSION(CUB_MAXCUBP, NDIM3D) :: Dxi1D, DXi2D
!
!  ! For every cubature point on the reference element,
!  ! the corresponding cubature weight
!  REAL(DP), DIMENSION(CUB_MAXCUBP) :: Domega
!
!  ! number of cubature points on the reference element
!  INTEGER :: ncubp
!
!  ! Pointer to the vector entries
!  REAL(DP), DIMENSION(:), POINTER :: p_DdataUX, p_DdataUY, p_DdataP
!
!  ! An allocateable array accepting the DOF's of a set of elements.
!  integer, DIMENSION(:,:), ALLOCATABLE, TARGET :: IdofsTrialU, IdofsTrialP
!
!  ! Allocateable arrays for the values of the basis functions.
!  REAL(DP), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: DbasTrialU,DbasTrialP
!
!  ! Number of entries in the vector - for quicker access
!  INTEGER(I32) :: neqU, neqP
!
!  ! Arrays for saving Jacobian determinants and matrices
!  REAL(DP), DIMENSION(:,:), POINTER :: p_Ddetj
!  REAL(DP), DIMENSION(:,:,:), POINTER :: p_Djac
!
!  ! Pointer to KVERT of the triangulation
!  INTEGER(I32), DIMENSION(:,:), POINTER :: p_IverticesAtElement
!
!  ! Pointer to DCORVG of the triangulation
!  REAL(DP), DIMENSION(:,:), POINTER :: p_DvertexCoords
!
!
!  ! Number of elements in a block. Normally =BILF_NELEMSIM,
!  ! except if there are less elements in the discretisation.
!  INTEGER :: nelementsPerBlock
!
!  BderU = .FALSE.
!  BderP = .FALSE.
!
!  ! Get the vector data
!  neqU = rvector%RvectorBlock(1)%NEQ
!  neqP = rvector%RvectorBlock(3)%NEQ
!
!  IF (rvector%cdataType .NE. ST_DOUBLE) THEN
!    PRINT *,'ppns2D_bdforces: Unsupported vector precision.'
!    STOP
!  END IF
!
!  ! We support only uniform and conformal discretisation structures.
!  IF (.NOT. ASSOCIATED(rvector%p_rblockDiscretisation)) THEN
!    PRINT *,'ppns2D_bdforces: No discretisation structure!'
!    STOP
!  END IF
!
!  IF ((rvector%p_rblockDiscretisation%ccomplexity .NE. SPDISC_UNIFORM) .AND. &
!      (rvector%p_rblockDiscretisation%ccomplexity .NE. SPDISC_CONFORMAL)) THEN
!    PRINT *,'ppns2D_bdforces: Discretisation too complex!'
!    STOP
!  END IF
!
!  ! Get pointers to the subvectors from the block vector
!  CALL lsyssc_getbase_double (rvector%RvectorBlock(1),p_DdataUX)
!  CALL lsyssc_getbase_double (rvector%RvectorBlock(2),p_DdataUY)
!  CALL lsyssc_getbase_double (rvector%RvectorBlock(3),p_DdataP)
!
!  ! Get pointers to the spatial discretisation structures of the
!  ! velocity and pressure
!  p_rdiscrU => rvector%RvectorBlock(1)%p_spatialDiscretisation
!  p_rdiscrP => rvector%RvectorBlock(3)%p_spatialDiscretisation
!
!  ! Now loop over the different element distributions in the discretisation.
!  ! Each element distribution has a different element combination.
!  ! E.g. distribution 1 may describe P1~/P1~/P0 triangular elements, while
!  ! distrbution 2 may describe Q1~/Q1~/Q0.
!  !
!  ! For every element distribution, we have to loop through all elements
!  ! along the boundary component, calculate their contribution to the
!  ! integral and add everything together.
!
!  DO icurrentElementDistr = 1,rdiscretisation%inumFESpaces
!
!    ! Get the element distribution structure of that FE-space comination
!    p_relemDistrU => p_rdiscrU%RelementDistr(RelementDistr)
!    p_relemDistrP => p_rdiscrP%RelementDistr(RelementDistr)
!
!    ! By that we know
!    ! - what's the element
!    ! - which transformation from the reference element to use
!
!
!
!
!
!  END DO ! icurrentElementDistr
!
!  END SUBROUTINE

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
      function dens(r,z)

      include 'commonblocks'

      if( idiskflag .eq. 1 ) then
         addens = appdiskdens(r,z)
         dens=totdens(r,z)
         dens = dens - addens
      else
         dens=totdens(r,z)
      endif
      return
      end

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
 !
 ! This file is part of the ExaHyPE project.
 ! Copyright (c) 2016  http://exahype.eu
 ! All rights reserved.
 !
 ! The project has received funding from the European Union's Horizon
 ! 2020 research and innovation programme under grant agreement
 ! No 671698. For copyrights and licensing, please consult the webpage.
 !
 ! Released under the BSD 3 Open Source License.
 ! For the full license text, see LICENSE.txt
 !

SUBROUTINE ADERRiemannSolverNonlinear(lQbndL,lFbndL,lQbndR,lFbndR,nv)
    USE typesDef

    USE, INTRINSIC :: ISO_C_BINDING
    IMPLICIT NONE
    ! Argument list
    REAL, INTENT(IN)     :: lQbndL(nVar,nDOF(2),nDOF(3))              ! left space-time degrees of freedom
    REAL, INTENT(IN)     :: lQbndR(nVar,nDOF(2),nDOF(3))              ! right space-time degrees of freedom
    REAL, INTENT(IN)     :: nv(d)                                     ! normal vector
    REAL, INTENT(INOUT)  :: lFbndL(nVar,nDOF(2),nDOF(3))              ! left flux
    REAL, INTENT(INOUT)  :: lFbndR(nVar,nDOF(2),nDOF(3))              ! right flux
    ! Local variables
    INTEGER           :: i,j,k,l
    REAL              :: aux(d), QavL(nVar), QavR(nVar), smax
    REAL              :: LL(nVar), LR(nVar), Qav(nVar), ncp(nVar), Bn(nVar,nVar)

    !OPEN(UNIT=12, FILE="aoutput_lQbndL.txt", ACTION="write", STATUS="replace")
    !WRITE(12, '(ES24.16,1x)') , lQbndL

    !OPEN(UNIT=12, FILE="aoutput_lQbndR.txt", ACTION="write", STATUS="replace")
    !WRITE(12, '(ES24.16,1x)') , lQbndR
    !CALL EXIT

 !   PRINT *, nv
 !   CALL EXIT

    !
    ! Compute the average states from the left and the right, which we need to compute the numerical dissipation
    QavL = 0.
    QavR = 0.
    DO k = 1, nDOF(3)
      DO j = 1, nDOF(2)
            aux = (/ 1.d0, wGPN(j), wGPN(k) /)
            QavL = QavL + PRODUCT(aux(1:nDim))*lQbndL(:,j,k)
            QavR = QavR + PRODUCT(aux(1:nDim))*lQbndR(:,j,k)
        ENDDO
    ENDDO
    !
    Qav = 0.5*(QavL+QavR)
    !
    ! Here, we implement a very simple Rusanov scheme with scalar dissipation (smax*Id).
    ! We can change this into a more sophisticated Osher or HLLEM Riemann solver whenever needed!
    !
    CALL PDEEigenvalues(LL,QavL,nv)
    CALL PDEEigenvalues(LR,QavR,nv)
    smax = MAX( MAXVAL(ABS(LL)), MAXVAL(ABS(LR)) )
    !
    ! Non-conservative product with path conservative schemes
    CALL PDEMatrixB(Bn,Qav,nv)
    !
    ! We now compute the numerical flux. Note that the scheme is at the moment written in
    ! CONSERVATION FORM => no fluctuations, but real fluxes.
    ! Later, this will be converted into the left and right fluctuations.
    !
    DO k = 1, nDOF(3)
      DO j = 1, nDOF(2)
            lFbndL(:,j,k) = 0.5*( lFbndR(:,j,k) + lFbndL(:,j,k) ) - 0.5*smax*( lQbndR(:,j,k) - lQbndL(:,j,k) )
            ncp = MATMUL( Bn, lQbndR(:,j,k) - lQbndL(:,j,k) )
            lFbndR(:,j,k) = lFbndL(:,j,k) - 0.5*ncp(:)
            lFbndL(:,j,k) = lFbndL(:,j,k) + 0.5*ncp(:)
        ENDDO
    ENDDO
    !

    !PRINT *, ' --------lFbndL-ADERRiemannSolver-------------------------------- '
    !PRINT *, lFbndL
    !PRINT *, ' --------lFbndL-ADERRiemannSolver-------------------------------- '

    !OPEN(UNIT=12, FILE="aoutput_lFbndL.txt", ACTION="write", STATUS="replace")
    !WRITE(12, '(ES24.16,1x)') , lFbndL
    !CALL EXIT
END SUBROUTINE ADERRiemannSolverNonlinear

max_stars_repo_name: linusseelinger/ExaHyPE-Tsunami

max_issues_repo_name: linusseelinger/ExaHyPE-Tsunami

max_forks_repo_name: linusseelinger/ExaHyPE-Tsunami


text:
      SUBROUTINE zbrak(fx,x1,x2,n,xb1,xb2,nb)
      INTEGER n,nb
      REAL x1,x2,xb1(nb),xb2(nb),fx
      EXTERNAL fx
      INTEGER i,nbb
      REAL dx,fc,fp,x
      nbb=0
      x=x1
      dx=(x2-x1)/n
      fp=fx(x)
      do 11 i=1,n
        x=x+dx
        fc=fx(x)
        if(fc*fp.le.0.) then
          nbb=nbb+1
          xb1(nbb)=x-dx
          xb2(nbb)=x
          if(nbb.eq.nb)goto 1
        endif
        fp=fc
11    continue
1     continue
      nb=nbb
      return
      END

ext: for

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


text:
c-----------------------------------------------------------------------
      subroutine read_re2_data(ifbswap, ifxyz, ifcur, ifbc)  ! .re2 reader

      include 'SIZE'
      include 'TOTAL'
      include 'RESTART'
      include 'CTIMER'

      logical ifbswap
      logical ifxyz, ifcur, ifbc
      integer idummy(100)

      common /nekmpi/ nidd,npp,nekcomm,nekgroup,nekreal

      etime0 = dnekclock_sync()

                  ibc = 2
      if (ifflow) ibc = 1

                  nfldt = 1
      if (ifheat) nfldt = 2+npscal
      if (ifmhd ) nfldt = 2+npscal+1

      ! first field to read
      if (param(33).gt.0) ibc = int(param(33))

      ! number of fields to read
      if (param(32).gt.0) then
        nfldt = ibc + int(param(32)) - 1
        nfldt = max(nfldt,1)
        if (nelgt.gt.nelgv) nfldt = max(nfldt,2)
      endif

      call blank(cbc,3*size(cbc))
      call rzero(bc ,size(bc))

#ifndef NOMPIIO
      call fgslib_crystal_setup(cr_re2,nekcomm,np)

      call byte_open_mpi(re2fle,fh_re2,.true.,ierr)
      call err_chk(ierr,' Cannot open .re2 file!$')

      call readp_re2_mesh (ifbswap,ifxyz)
      call readp_re2_curve(ifbswap,ifcur)
      do ifield = ibc,nfldt
        call readp_re2_bc(cbc(1,1,ifield),bc(1,1,1,ifield),
     &    ifbswap,ifbc)
      enddo

      call fgslib_crystal_free(cr_re2)
      call byte_close_mpi(fh_re2,ierr)
#else
      call byte_open(re2fle,ierr)
      call byte_read(idummy,21,ierr) ! skip hdr+endian code

      call bin_rd1_mesh (ifbswap)
      call bin_rd1_curve(ifbswap)
      do ifield = ibc,nfldt
         call bin_rd1_bc (cbc(1,1,ifield),bc(1,1,1,ifield),ifbswap)
      enddo

      call byte_close(ierr)
#endif

      etime_t = dnekclock_sync() - etime0
      if(nio.eq.0) write(6,'(A,1(1g9.2),A,/)')
     &                   ' done :: read .re2 file   ',
     &                   etime_t, ' sec'

      return
      end
c-----------------------------------------------------------------------
      subroutine readp_re2_mesh(ifbswap,ifread) ! version 2 of .re2 reader

      include 'SIZE'
      include 'TOTAL'

      logical ifbswap
      logical ifread

      parameter(nrmax = lelt)             ! maximum number of records
      parameter(lrs   = 1+ldim*(2**ldim)) ! record size: group x(:,c) ...
      parameter(li    = 2*lrs+2)

      integer         bufr(li-2,nrmax)
      common /scrns/  bufr

      integer         vi  (li  ,nrmax)
      common /ctmp1/  vi

      integer*8       lre2off_b,dtmp8
      integer*8       nrg

      nrg       = nelgt
      nr        = nelt
      irankoff  = igl_running_sum(nr) - nr
      dtmp8     = irankoff
      re2off_b  = 84 ! set initial offset (hdr + endian)
      lre2off_b = re2off_b + dtmp8*lrs*wdsizi
      lrs4      = lrs*wdsizi/4

      ! read coordinates from file
      nwds4r = nr*lrs4
      call byte_set_view(lre2off_b,fh_re2)
      call byte_read_mpi(bufr,nwds4r,-1,fh_re2,ierr)
      re2off_b = re2off_b + nrg*4*lrs4
      if (ierr.gt.0) goto 100

      if (.not.ifread) return

      if (nio.eq.0) write(6,*) 'reading mesh '

      ! pack buffer
      do i = 1,nr
         jj      = (i-1)*lrs4 + 1
         ielg    = irankoff + i ! elements are stored in global order
         vi(1,i) = gllnid(ielg)
         vi(2,i) = ielg
         call icopy(vi(3,i),bufr(jj,1),lrs4)
      enddo

      ! crystal route nr real items of size lrs to rank vi(key,1:nr)
      n   = nr
      key = 1
      call fgslib_crystal_tuple_transfer(cr_re2,n,nrmax,vi,li,
     &   vl,0,vr,0,key)

      ! unpack buffer
      ierr = 0
      if (n.gt.nrmax) then
         ierr = 1
         goto 100
      endif

      do i = 1,n
         iel = gllel(vi(2,i))
         call icopy     (bufr,vi(3,i),lrs4)
         call buf_to_xyz(bufr,iel,ifbswap,ierr)
      enddo

 100  call err_chk(ierr,'Error reading .re2 mesh$')

      return
      end
c-----------------------------------------------------------------------
      subroutine readp_re2_curve(ifbswap,ifread)

      include 'SIZE'
      include 'TOTAL'

      logical ifbswap
      logical ifread

      common /nekmpi/ nidd,npp,nekcomm,nekgroup,nekreal

      parameter(nrmax = 12*lelt) ! maximum number of records
      parameter(lrs   = 2+1+5)   ! record size: eg iside curve(5) ccurve
      parameter(li    = 2*lrs+1)

      integer         bufr(li-1,nrmax)
      common /scrns/  bufr

      integer         vi  (li  ,nrmax)
      common /ctmp1/  vi

      integer*8       lre2off_b,dtmp8
      integer*8       nrg
      integer*4       nrg4(2)

      integer*8       i8gl_running_sum

      ! read total number of records
      nwds4r    = 1*wdsizi/4
      lre2off_b = re2off_b
      call byte_set_view(lre2off_b,fh_re2)
      call byte_read_mpi(nrg4,nwds4r,-1,fh_re2,ierr)
      if(ierr.gt.0) goto 100

      if(wdsizi.eq.8) then
         if(ifbswap) call byte_reverse8(nrg4,nwds4r,ierr)
         call copy(dnrg,nrg4,1)
         nrg = dnrg
      else
         if(ifbswap) call byte_reverse (nrg4,nwds4r,ierr)
         nrg = nrg4(1)
      endif
      re2off_b = re2off_b + 4*nwds4r

      if(nrg.eq.0) return

      ! read data from file
      dtmp8 = np
      nr = nrg/dtmp8
      do i = 0,mod(nrg,dtmp8)-1
         if(i.eq.nid) nr = nr + 1
      enddo
      dtmp8     = i8gl_running_sum(int(nr,8)) - nr
      lre2off_b = re2off_b + dtmp8*lrs*wdsizi
      lrs4      = lrs*wdsizi/4

      re2off_b = re2off_b + nrg*4*lrs4

      if (.not.ifread) return
      if(nio.eq.0) write(6,*) 'reading curved sides '

      nwds4r = nr*lrs4
      call byte_set_view(lre2off_b,fh_re2)
      call byte_read_mpi(bufr,nwds4r,-1,fh_re2,ierr)
      if(ierr.gt.0) goto 100

      ! pack buffer
      do i = 1,nr
         jj = (i-1)*lrs4 + 1

         if(ifbswap) then
           lrs4s = lrs4 - wdsizi/4 ! words to swap (last is char)
           if(wdsizi.eq.8) call byte_reverse8(bufr(jj,1),lrs4s,ierr)
           if(wdsizi.eq.4) call byte_reverse (bufr(jj,1),lrs4s,ierr)
         endif

         ielg = bufr(jj,1)
         if(wdsizi.eq.8) call copyi4(ielg,bufr(jj,1),1)

         if(ielg.le.0 .or. ielg.gt.nelgt) goto 100
         vi(1,i) = gllnid(ielg)

         call icopy (vi(2,i),bufr(jj,1),lrs4)
      enddo

      ! crystal route nr real items of size lrs to rank vi(key,1:nr)
      n    = nr
      key  = 1
      call fgslib_crystal_tuple_transfer(cr_re2,n,nrmax,vi,li,vl,0,vr,0,
     &                                   key)

      ! unpack buffer
      if(n.gt.nrmax) goto 100
      do i = 1,n
         call icopy       (bufr,vi(2,i),lrs4)
         call buf_to_curve(bufr)
      enddo

      return

 100  ierr = 1
      call err_chk(ierr,'Error reading .re2 curved data$')

      end
c-----------------------------------------------------------------------
      subroutine readp_re2_bc(cbl,bl,ifbswap,ifread)

      include 'SIZE'
      include 'TOTAL'

      character*3  cbl(  6,lelt)
      real         bl (5,6,lelt)
      logical      ifbswap
      logical      ifread

      parameter(nrmax = 6*lelt) ! maximum number of records
      parameter(lrs   = 2+1+5)  ! record size: eg iside bl(5) cbl
      parameter(li    = 2*lrs+1)

      integer         bufr(li-1,nrmax)
      common /scrns/  bufr

      integer         vi  (li  ,nrmax)
      common /ctmp1/  vi

      integer*8       lre2off_b,dtmp8
      integer*8       nrg
      integer*4       nrg4(2)

      integer*8       i8gl_running_sum

      ! read total number of records
      nwds4r    = 1*wdsizi/4
      lre2off_b = re2off_b
      call byte_set_view(lre2off_b,fh_re2)
      call byte_read_mpi(nrg4,nwds4r,-1,fh_re2,ierr)
      if(ierr.gt.0) goto 100

      if(wdsizi.eq.8) then
         if(ifbswap) call byte_reverse8(nrg4,nwds4r,ierr)
         call copy(dnrg,nrg4,1)
         nrg = dnrg
      else
         if(ifbswap) call byte_reverse (nrg4,nwds4r,ierr)
         nrg = nrg4(1)
      endif
      re2off_b = re2off_b + 4*nwds4r

      if(nrg.eq.0) return

      ! read data from file
      dtmp8 = np
      nr = nrg/dtmp8
      do i = 0,mod(nrg,dtmp8)-1
         if(i.eq.nid) nr = nr + 1
      enddo
      dtmp8     = i8gl_running_sum(int(nr,8)) - nr
      lre2off_b = re2off_b + dtmp8*lrs*wdsizi
      lrs4      = lrs*wdsizi/4

      re2off_b = re2off_b + nrg*4*lrs4

      if (.not.ifread) return
      if(nio.eq.0) write(6,*) 'reading bc for ifld',ifield

      nwds4r = nr*lrs4
      call byte_set_view(lre2off_b,fh_re2)
      call byte_read_mpi(bufr,nwds4r,-1,fh_re2,ierr)
      if(ierr.gt.0) goto 100

      ! pack buffer
      do i = 1,nr
         jj = (i-1)*lrs4 + 1

         if(ifbswap) then
           lrs4s = lrs4 - wdsizi/4 ! words to swap (last is char)
           if(wdsizi.eq.8) call byte_reverse8(bufr(jj,1),lrs4s,ierr)
           if(wdsizi.eq.4) call byte_reverse (bufr(jj,1),lrs4s,ierr)
         endif

         ielg = bufr(jj,1)
         if(wdsizi.eq.8) call copyi4(ielg,bufr(jj,1),1)

         if(ielg.le.0 .or. ielg.gt.nelgt) goto 100
         vi(1,i) = gllnid(ielg)

         call icopy (vi(2,i),bufr(jj,1),lrs4)
      enddo

      ! crystal route nr real items of size lrs to rank vi(key,1:nr)
      n    = nr
      key  = 1

      call fgslib_crystal_tuple_transfer(cr_re2,n,nrmax,vi,li,vl,0,vr,0,
     &                                   key)

      ! fill up with default
      do iel=1,nelt
      do k=1,6
         cbl(k,iel) = 'E  '
      enddo
      enddo

      ! unpack buffer
      if(n.gt.nrmax) goto 100
      do i = 1,n
         call icopy    (bufr,vi(2,i),lrs4)
         call buf_to_bc(cbl,bl,bufr)
      enddo

      return

 100  ierr = 1
      call err_chk(ierr,'Error reading .re2 boundary data$')

      end
c-----------------------------------------------------------------------
      subroutine buf_to_xyz(buf,e,ifbswap,ierr)! version 1 of binary reader

      include 'SIZE'
      include 'TOTAL'
      logical ifbswap

c      integer e,eg,buf(0:49)
      integer e,eg,buf(0:49)

      nwds = (1 + ldim*(2**ldim))*(wdsizi/4) ! group + 2x4 for 2d, 3x8 for 3d

      if     (ifbswap.and.ierr.eq.0.and.wdsizi.eq.8) then
          call byte_reverse8(buf,nwds,ierr)
      elseif (ifbswap.and.ierr.eq.0.and.wdsizi.eq.4) then
          call byte_reverse (buf,nwds,ierr)
      endif
      if(ierr.ne.0) return

      if(wdsizi.eq.8) then
         call copyi4(igroup(e),buf(0),1) !0-1
         if (ldim.eq.3) then
            call copy  (xc(1,e),buf( 2),8) !2 --17
            call copy  (yc(1,e),buf(18),8) !18--33
            call copy  (zc(1,e),buf(34),8) !34--49
         else
            call copy  (xc(1,e),buf( 2),4) !2 --9
            call copy  (yc(1,e),buf(10),4) !10--17
          endif
      else
         igroup(e) = buf(0)
         if (if3d) then
            call copy4r(xc(1,e),buf( 1),8)
            call copy4r(yc(1,e),buf( 9),8)
            call copy4r(zc(1,e),buf(17),8)
         else
            call copy4r(xc(1,e),buf( 1),4)
            call copy4r(yc(1,e),buf( 5),4)
         endif
      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine buf_to_curve(buf)    ! version 1 of binary reader

      include 'SIZE'
      include 'TOTAL'

      integer e,eg,f,buf(30)

      if(wdsizi.eq.8) then
        call copyi4(eg,buf(1),1) !1-2
        e  = gllel(eg)

        call copyi4(f,buf(3),1) !3-4

        call copy  ( curve(1,f,e),buf(5) ,5) !5--14
        call chcopy(ccurve(  f,e),buf(15),1)!15
      else
        eg = buf(1)
        e  = gllel(eg)
        f  = buf(2)

        call copy4r( curve(1,f,e),buf(3),5)
        call chcopy(ccurve(f,e)  ,buf(8),1)
      endif

c     write(6,1) eg,e,f,(curve(k,f,e),k=1,5),ccurve(f,e)
c   1 format(2i7,i3,5f10.3,1x,a1,'ccurve')

      return
      end
c-----------------------------------------------------------------------
      subroutine buf_to_bc(cbl,bl,buf)    ! version 1 of binary reader

      include 'SIZE'
      include 'TOTAL'

      character*3 cbl(6,lelt)
      real         bl(5,6,lelt)

      integer e,eg,f,buf(30)

      if(wdsizi.eq.8) then
        call copyi4(eg,buf(1),1) !1-2
        e  = gllel(eg)

        call copyi4(f,buf(3),1) !3-4

        call copy  (bl(1,f,e),buf(5),5) !5--14
        call chcopy(cbl( f,e),buf(15),3)!15-16

        if(nelt.ge.1000000.and.cbl(f,e).eq.'P  ')
     $   call copyi4(bl(1,f,e),buf(5),1) !Integer assign connecting P element

      else
        eg = buf(1)
        e  = gllel(eg)
        f  = buf(2)

        call copy4r ( bl(1,f,e),buf(3),5)
        call chcopy (cbl(  f,e),buf(8),3)

        if (nelgt.ge.1000000.and.cbl(f,e).eq.'P  ')
     $     bl(1,f,e) = buf(3) ! Integer assign of connecting periodic element
      endif

c      write(6,1) eg,e,f,cbl(f,e),' CBC',nid
c  1   format(2i8,i4,2x,a3,a4,i8)

      return
      end
c-----------------------------------------------------------------------
      subroutine bin_rd1_mesh(ifbswap)    ! version 1 of binary reader

      include 'SIZE'
      include 'TOTAL'
      logical ifbswap

      integer e,eg,buf(55)

      if (nio.eq.0) write(6,*)    '  reading mesh '

      nwds = (1 + ldim*(2**ldim))*(wdsizi/4) ! group + 2x4 for 2d, 3x8 for 3d
      len  = 4*nwds                          ! 4 bytes / wd

      if (nwds.gt.55.or.isize.gt.4) then
         write(6,*) nid,' Error in bin_rd1_mesh: buf size',nwds,isize
         call exitt
      endif

      call nekgsync()

      niop = 10
      do k=1,8
         if (nelgt/niop .lt. 100) goto 10
         niop = niop*10
      enddo
   10 continue

      ierr  = 0
      ierr2 = 0
      len1  = 4
      do eg=1,nelgt             ! sync NOT needed here

         mid = gllnid(eg)
         e   = gllel (eg)
#ifdef DEBUG
         if (nio.eq.0.and.mod(eg,niop).eq.0) write(6,*) eg,' mesh read'
#endif
         if (mid.ne.nid.and.nid.eq.0) then              ! read & send

            if(ierr.eq.0) then
              call byte_read  (buf,nwds,ierr)
              call csend(e,ierr,len1,mid,0)
              if(ierr.eq.0) call csend(e,buf,len,mid,0)
            else
              call csend(e,ierr,len1,mid,0)
            endif

         elseif (mid.eq.nid.and.nid.ne.0) then          ! recv & process

            call crecv      (e,ierr,len1)
            if(ierr.eq.0) then
              call crecv      (e,buf,len)
              call buf_to_xyz (buf,e,ifbswap,ierr2)
            endif

         elseif (mid.eq.nid.and.nid.eq.0) then          ! read & process

            if(ierr.eq.0) then
              call byte_read  (buf,nwds,ierr)
              call buf_to_xyz (buf,e,ifbswap,ierr2)
            endif
         endif

      enddo
      ierr = ierr + ierr2
      call err_chk(ierr,'Error reading .re2 mesh. Abort. $')

      return
      end
c-----------------------------------------------------------------------
      subroutine bin_rd1_curve (ifbswap) ! v. 1 of curve side reader

      include 'SIZE'
      include 'TOTAL'
      logical ifbswap

      integer e,eg,buf(55)
      real rcurve

      nwds = (2 + 1 + 5)*(wdsizi/4) !eg+iside+ccurve+curve(6,:,:) !only 5 in rea
      len  = 4*nwds      ! 4 bytes / wd

      if (nwds.gt.55.or.isize.gt.4) then
         write(6,*)nid,' Error in bin_rd1_curve: buf size',nwds,isize
         call exitt
      endif

      call nekgsync()

      ierr = 0
      len1 = 4
      if (nid.eq.0) then  ! read & send/process

         if(wdsizi.eq.8) then
           call byte_read(rcurve,2,ierr)
           if (ifbswap) call byte_reverse8(rcurve,2,ierr)
           ncurve = rcurve
         else
           call byte_read(ncurve,1,ierr)
           if (ifbswap) call byte_reverse(ncurve,1,ierr)
         endif

         if(ncurve.ne.0) write(6,*) '  reading curved sides '
         do k=1,ncurve
           if(ierr.eq.0) then
              call byte_read(buf,nwds,ierr)
              if(wdsizi.eq.8) then
                if(ifbswap) call byte_reverse8(buf,nwds-2,ierr)
                call copyi4(eg,buf(1),1)  !1,2
              else
                if (ifbswap) call byte_reverse(buf,nwds-1,ierr) ! last is char
                eg  = buf(1)
              endif

              mid = gllnid(eg)
              if (mid.eq.0.and.ierr.eq.0) then
                 call buf_to_curve(buf)
              else
                 if(ierr.eq.0) then
                   call csend(mid,buf,len,mid,0)
                 else
                   goto 98
                 endif
              endif
           else
              goto 98
           endif
         enddo
  98     call buf_close_out  ! notify all procs: no more data

      else               ! wait for data from node 0

         ncurve_mx = 12*nelt
         do k=1,ncurve_mx+1   ! +1 to make certain we receive the close-out

            call crecv(nid,buf,len)
            if(wdsizi.eq.8) then
               call copyi4(ichk,buf(1),1)
               if(ichk.eq.0) goto 99
               call buf_to_curve(buf)
            elseif (buf(1).eq.0) then
               goto 99
            else
               call buf_to_curve(buf)
            endif

         enddo
   99    call buf_close_out

      endif
      call err_chk(ierr,'Error reading .re2 curved data. Abort.$')

      return
      end
c-----------------------------------------------------------------------
      subroutine bin_rd1_bc (cbl,bl,ifbswap) ! v. 1 of bc reader

      include 'SIZE'
      include 'TOTAL'
      logical ifbswap

      character*3 cbl(6,lelt)
      real         bl(5,6,lelt)

      integer e,eg,buf(55)
      real rbc_max

      nwds = (2 + 1 + 5)*(wdsizi/4)   ! eg + iside + cbc + bc(5,:,:)
      len  = 4*nwds      ! 4 bytes / wd

      if (nwds.gt.55.or.isize.gt.4) then
         write(6,*) nid,' Error in bin_rd1_bc: buf size',nwds,isize
         call exitt
      endif

      do e=1,nelt   ! fill up cbc w/ default
      do k=1,6
         cbl(k,e) = 'E  '
      enddo
      enddo

      call nekgsync()
      ierr=0
      len1=4
      if (nid.eq.0) then  ! read & send/process

         if(wdsizi.eq.8) then
           call byte_read(rbc_max,2,ierr)
           if (ifbswap) call byte_reverse8(rbc_max,2,ierr) ! last is char
           nbc_max = rbc_max
         else
           call byte_read(nbc_max,1,ierr)
           if (ifbswap) call byte_reverse(nbc_max,1,ierr) ! last is char
         endif

         if(nbc_max.ne.0) write(6,*) '  reading bc for ifld',ifield
         do k=1,nbc_max
c           write(6,*) k,' dobc1 ',nbc_max
            if(ierr.eq.0) then
               call byte_read(buf,nwds,ierr)
               if(wdsizi.eq.8) then
                 if (ifbswap) call byte_reverse8(buf,nwds-2,ierr)
                 call copyi4(eg,buf(1),1) !1&2 of buf
               else
                 if (ifbswap) call byte_reverse(buf,nwds-1,ierr) ! last is char
                 eg  = buf(1)
               endif
               mid = gllnid(eg)
c              write(6,*) k,' dobc3 ',eg,mid

               if (mid.eq.0.and.ierr.eq.0) then
                   call buf_to_bc(cbl,bl,buf)
               else
c                  write(6,*) mid,' sendbc1 ',eg
                   if(ierr.eq.0) then
                     call csend(mid,buf,len,mid,0)
                   else
                     goto 98
                   endif
c                  write(6,*) mid,' sendbc2 ',eg
               endif
c              write(6,*) k,' dobc2 ',nbc_max,eg
            else
               goto 98
            endif
         enddo
c        write(6,*) mid,' bclose ',eg,nbc_max
  98     call buf_close_outv ! notify all procs: no more data

      else               ! wait for data from node 0

         nbc_max = 2*ldim*nelt
         do k=1,nbc_max+1  ! Need one extra !

c           write(6,*) nid,' recvbc1',k
            call crecv(nid,buf,len)
c           write(6,*) nid,' recvbc2',k,buf(1)

            if(wdsizi.eq.8) then
               call copyi4(ichk,buf(1),1)
               if(ichk.eq.0) goto 99
               call buf_to_bc(cbl,bl,buf)
            elseif (buf(1).eq.0) then
                goto 99
            else
                call buf_to_bc(cbl,bl,buf)
            endif

         enddo
   99    call buf_close_outv

      endif

      call err_chk(ierr,'Error reading boundary data for re2. Abort.$')

      return
      end
c-----------------------------------------------------------------------
      subroutine buf_close_outv  ! this is the stupid O(P) formulation

      include 'SIZE'
      include 'PARALLEL'
      integer*4 zero
      real      rzero

      len   = wdsizi
      rzero = 0
      zero  = 0
c     write(6,*) nid,' bufclose'
      if (nid.eq.0) then
         do mid=1,np-1
            if(wdsizi.eq.8)call csend(mid,rzero,len,mid,0)
            if(wdsizi.eq.4)call csend(mid, zero,len,mid,0)
c           write(6,*) mid,' sendclose'
         enddo
      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine buf_close_out  ! this is the stupid O(P) formulation

      include 'SIZE'
      include 'PARALLEL'
      integer*4 zero
      real      rzero

c     len  = 4
      len   = wdsizi
      zero = 0
      rzero = 0
      if (nid.eq.0) then
         do mid=1,np-1
            if(wdsizi.eq.8)call csend(mid,rzero,len,mid,0)
            if(wdsizi.eq.4)call csend(mid, zero,len,mid,0)
         enddo
      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine read_re2_hdr(ifbswap, ifverbose) ! open file & chk for byteswap

      include 'SIZE'
      include 'TOTAL'

      logical ifbswap, ifverbose
      logical if_byte_swap_test

      integer fnami (33)
      character*132 fname
      equivalence (fname,fnami)

      character*132 hdr
      character*5 version
      real*4      test

      logical iffound

      ierr=0

      if (nid.eq.0) then
         if (ifverbose) write(6,'(A,A)') ' Reading ', re2fle
         call izero(fnami,33)
         m = indx2(re2fle,132,' ',1)-1
         call chcopy(fname,re2fle,m)

         inquire(file=fname, exist=iffound)
         if(.not.iffound) ierr = 1
      endif
      call err_chk(ierr,' Cannot find re2 file!$')

      if (nid.eq.0) then
         call byte_open(fname,ierr)
         if(ierr.ne.0) goto 100
         call byte_read(hdr,20,ierr)
         if(ierr.ne.0) goto 100

         read (hdr,1) version,nelgt,ldimr,nelgv
    1    format(a5,i9,i3,i9)

         wdsizi = 4
         if(version.eq.'#v002') wdsizi = 8
         if(version.eq.'#v003') then
           wdsizi = 8
           param(32)=1
         endif

         call byte_read(test,1,ierr)
         if(ierr.ne.0) goto 100
         ifbswap = if_byte_swap_test(test,ierr)
         if(ierr.ne.0) goto 100

         call byte_close(ierr)
      endif

 100  call err_chk(ierr,'Error reading re2 header$')

      call bcast(wdsizi, ISIZE)
      call bcast(ifbswap,LSIZE)
      call bcast(nelgv  ,ISIZE)
      call bcast(nelgt  ,ISIZE)
      call bcast(ldimr  ,ISIZE)
      call bcast(param(32),WDSIZE)

      if(wdsize.eq.4.and.wdsizi.eq.8)
     $   call exitti('wdsize=4 & wdsizi(re2)=8 not compatible$',wdsizi)

      return
      end

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


text:
module vegetables_example_m
    use vegetables_input_m, only: input_t

    implicit none
    private
    public :: example_t

    type :: example_t
        private
        class(input_t), allocatable :: input_
    contains
        private
        procedure, public :: input
    end type

    interface example_t
        module procedure constructor
    end interface
contains
    function constructor(input) result(example)
        class(input_t), intent(in) :: input
        type(example_t) :: example

        allocate(example%input_, source = input)
    end function

    function input(self)
        class(example_t), intent(in) :: self
        class(input_t), allocatable :: input

        allocate(input, source = self%input_)
    end function
end module

max_stars_repo_name: everythingfunctional/vegetables

max_issues_repo_name: everythingfunctional/vegetables

max_forks_repo_name: everythingfunctional/vegetables


text:
      SUBROUTINE magnetics_data
      USE vmec_main
      USE vacmod
      USE vsvd
      USE mgrid_mod
      IMPLICIT NONE
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: iobset, j, i, n, ibfld
      REAL(rprec) :: denwgt, dsitotal,
     1   chierr, wght0, abs_flux, chiwgt, raddeg, denbwgt, btotal,
     2   relerr_b, wghtb, chisqx, chibwgt
C-----------------------------------------------
!
!       ALGEBRAIC SUM OF FLUXES FOR GROUPS OF OBSERVER POSITIONS
!
      IF (nobd .ne. 0) THEN
         flmwgt = 0
         denwgt = 0
         WRITE (nthreed, 400)
  400    FORMAT(/' External Poloidal Flux Loop Measurements (Wb)',/,
     1      '   I     Connection          Coils       Plasma    ',
     2      '    Total     Measured   Chi-Sq Err        Sigma',
     3      '   Designation',/,1x,3('-'),1x,16('-'),7(1x,12('-')),'-')

         DO iobset = 1, nobd
            dsitotal = dsiext(iobset) + plflux(iobset)
            chierr = dsiobt(iobset) - dsitotal
            wght0 = cbig + one
            IF (ANY(iobset.eq.indxflx(:nflxs)))
     1          wght0 = sigma_flux(iobset)
            flmwgt = flmwgt + (chierr/wght0)**2
            denwgt = denwgt + (dsitotal/wght0)**2
            chierr = (chierr/wght0)**2
            IF (wght0 .lt. cbig) THEN
               WRITE (nthreed, 410) iobset, (iconnect(j,iobset),j=1,4),
     1            dsiext(iobset), plflux(iobset), dsitotal, dsiobt(
     2            iobset), chierr, wght0, dsilabel(iobset)
            ELSE
               WRITE (nthreed, 420) iobset, (iconnect(j,iobset),j=1,4),
     1            dsiext(iobset), plflux(iobset), dsitotal, dsiobt(
     2            iobset), chierr, dsilabel(iobset)
            ENDIF
         END DO

         WRITE (nthreed, 600)
         DO i = 1, nobser
            abs_flux = plflux(nobd+i)
            WRITE (nthreed, 610) i, xobser(i), zobser(i), psiext(i),
     1         abs_flux, psiext(i) + ABS_flux
         END DO
  600    FORMAT(//,' Poloidal Flux at Individual Coil Locations',/,
     1      '   I         R [M]        Z [M]        Coils       Plasma',
     2      '        Total',/,2x,3('-'),5(1x,12('-')))
  610    FORMAT(i4,1x,5f13.4)
      ENDIF

      nchisaddle = nflxs
      total_saddle_chi = 0.
      IF (nflxs .GT. 0) THEN
         chiwgt = flmwgt/(nflxs)
         flmwgt = SQRT(flmwgt/denwgt)
         total_saddle_chi = (nflxs)*chiwgt
         WRITE (nthreed, 430) chiwgt, flmwgt
      ENDIF
  410 FORMAT(i4,1x,4i4,6f13.4,7x,a8)
  420 FORMAT(i4,1x,4i4,5f13.4,5x,'No Match',7x,a8)
  430 FORMAT(/' MEAN CHI-SQ ERROR IN SADDLE LOOP MATCH: ',1p,e10.3,/,
     1   ' RMS ERROR IN SADDLE LOOP MATCH: ',e10.3)

!
!       BR, BZ FIELD COMPONENTS FOR MATCHING AT OBSERVER POSITIONS
!
      total_b_chi = 0.
      IF (nbfldn .ne. 0) THEN
  500    FORMAT(//' BFIELD GROUP: ',a,' // ',
     1      'External B Loop Magnetic Field Measurements (T)',/,
     2      '  I  Rcoil[M]  Zcoil[M]    Degree    B[coil]',
     3      '  B[plasma]   B[Total]  ','  B[Meas]    Sigma  Chi-Sq Err',
     4      /,1x,2('-'),3(1x,9('-')),4(1x,10('-')),1(1x,8('-')),1(1x,11(
     5      '-')))

         raddeg = 360.0_dp/twopi
         DO n = 1, nbsets                        ! over all allowed sets
            WRITE (nthreed, 500) bloopnames(n)
            bcwgt = 0.
            denbwgt = 0.
            b_chi(n) = 0.
            ibfld = 0
            DO j = 1, nbcoils(n)
               ibfld = ibfld + 1
               btotal = bcoil(j,n) + plbfld(j,n)
               relerr_b = btotal - bbc(j,n)
               wghtb = cbig + 1.0_dp
c                                                ! over NAMELIST sets
               IF (ANY(j.eq.indxbfld(:nbfld(n),n))) wghtb = sigma_b(j,n)
               bcwgt = bcwgt + (relerr_b/wghtb)**2
               denbwgt = denbwgt + (btotal/wghtb)**2
               chisqx = (relerr_b/wghtb)**2
               IF (wghtb .LT. cbig) THEN
                  WRITE (nthreed, 520) ibfld, rbcoil(j,n), zbcoil(j,n),
     1               abcoil(j,n)*raddeg, bcoil(j,n), plbfld(j,n), btotal
     2               , bbc(j,n), wghtb, chisqx
               ELSE
                  WRITE (nthreed, 530) ibfld, rbcoil(j,n), zbcoil(j,n),
     1               abcoil(j,n)*raddeg, bcoil(j,n), plbfld(j,n), btotal
     2               , bbc(j,n)
               ENDIF
            END DO
            IF (nbfld(n) .gt. 0) THEN
               chibwgt = bcwgt/(nbfld(n))
               bcwgt = SQRT(bcwgt/denbwgt)
               b_chi(n) = (nbfld(n))*chibwgt
               WRITE (nthreed, 540) bloopnames(n), chibwgt,
     1            bloopnames(n), bcwgt
               total_b_chi = total_b_chi + b_chi(n)
            ENDIF
         END DO

      ENDIF
  520 FORMAT(i3,3f10.3,4f11.3,f9.4,f12.4)
  530 FORMAT(i3,3f10.3,4f11.3,3x,'No Match')
  540 FORMAT(/' MEAN CHI-SQ ERROR IN B-COIL GROUP ',a,': ',1p,e10.3,/,
     1   ' RMS ERROR IN B-COIL GROUP ',a,': ',e10.3)

      END SUBROUTINE magnetics_data


text:
subroutine amrex_probinit (init,name,namlen,problo,probhi) bind(c)

  use amrex_error_module
  use eos_module
  use eos_type_module
  use network
  use probdata_module

  use amrex_fort_module, only : rt => amrex_real
  implicit none

  integer init, namlen
  integer name(namlen)
  real(rt)         problo(2), probhi(2)
  real(rt)         xn(nspec)

  integer untin,i

  type (eos_t) :: eos_state

  namelist /fortin/ p_l, u_l, rho_l, p_r, u_r, rho_r, T_l, T_r, frac, idir, &
       use_Tinit

  !
  !     Build "probin" filename -- the name of file containing fortin namelist.
  !
  integer maxlen
  parameter (maxlen=256)
  character probin*(maxlen)

  if (namlen .gt. maxlen) then
     call amrex_error("probin file name too long")
  end if

  do i = 1, namlen
     probin(i:i) = char(name(i))
  end do

  ! set namelist defaults

  p_l = 1.0               ! left pressure (erg/cc)
  u_l = 0.0               ! left velocity (cm/s)
  rho_l = 1.0             ! left density (g/cc)
  T_l = 1.0

  p_r = 0.1               ! right pressure (erg/cc)
  u_r = 0.0               ! right velocity (cm/s)
  rho_r = 0.125           ! right density (g/cc)
  T_r = 1.0

  idir = 1                ! direction across which to jump
  frac = 0.5              ! fraction of the domain for the interface

  use_Tinit = .false.     ! optionally use T_l/r instead of p_l/r for initialization

  !     Read namelists
  untin = 9
  open(untin,file=probin(1:namlen),form='formatted',status='old')
  read(untin,fortin)
  close(unit=untin)

  !     set local variable defaults -- the 'center' variables are the location of the
  !     interface
  split(1) = frac*(problo(1)+probhi(1))
  split(2) = frac*(problo(2)+probhi(2))

  !     compute the internal energy (erg/cc) for the left and right state
  xn(:) = 0.0e0_rt
  xn(1) = 1.0e0_rt

  if (use_Tinit) then

     eos_state%rho = rho_l
     eos_state%T = T_l
     eos_state%xn(:) = xn(:)

     call eos(eos_input_rt, eos_state)

     rhoe_l = rho_l*eos_state%e
     p_l = eos_state%p

     eos_state%rho = rho_r
     eos_state%T = T_r
     eos_state%xn(:) = xn(:)

     call eos(eos_input_rt, eos_state)

     rhoe_r = rho_r*eos_state%e
     p_r = eos_state%p

  else

     eos_state%rho = rho_l
     eos_state%p = p_l
     eos_state%T = 100000.e0_rt  ! initial guess
     eos_state%xn(:) = xn(:)

     call eos(eos_input_rp, eos_state)

     rhoe_l = rho_l*eos_state%e
     T_l = eos_state%T

     eos_state%rho = rho_r
     eos_state%p = p_r
     eos_state%T = 100000.e0_rt  ! initial guess
     eos_state%xn(:) = xn(:)

     call eos(eos_input_rp, eos_state)

     rhoe_r = rho_r*eos_state%e
     T_r = eos_state%T

  endif

end subroutine amrex_probinit

! ::: -----------------------------------------------------------
! ::: This routine is called at problem setup time and is used
! ::: to initialize data on each grid.
! :::
! ::: NOTE:  all arrays have one cell of ghost zones surrounding
! :::        the grid interior.  Values in these cells need not
! :::        be set here.
! :::
! ::: INPUTS/OUTPUTS:
! :::
! ::: level     => amr level of grid
! ::: time      => time at which to init data
! ::: lo,hi     => index limits of grid interior (cell centered)
! ::: nstate    => number of state components.  You should know
! :::		   this already!
! ::: state     <=  Scalar array
! ::: delta     => cell size
! ::: xlo,xhi   => physical locations of lower left and upper
! :::              right hand corner of grid.  (does not include
! :::		   ghost region).
! ::: -----------------------------------------------------------
subroutine ca_initdata(level,time,lo,hi,nscal, &
                       state,state_l1,state_l2,state_h1,state_h2, &
                       delta,xlo,xhi)

  use network, only: nspec
  use probdata_module
  use meth_params_module, only : NVAR, URHO, UMX, UMY, UEDEN, UEINT, UTEMP, UFS

  use amrex_error_module
  use amrex_fort_module, only : rt => amrex_real
  implicit none

  integer level, nscal
  integer lo(2), hi(2)
  integer state_l1,state_l2,state_h1,state_h2
  real(rt)         xlo(2), xhi(2), time, delta(2)
  real(rt)         state(state_l1:state_h1,state_l2:state_h2,NVAR)

  real(rt)         xcen,ycen
  integer i,j

  do j = lo(2), hi(2)
     ycen = xlo(2) + delta(2)*(float(j-lo(2)) + 0.5e0_rt)

     do i = lo(1), hi(1)
        xcen = xlo(1) + delta(1)*(float(i-lo(1)) + 0.5e0_rt)

        if (idir == 1) then
           if (xcen <= split(1)) then
              state(i,j,URHO) = rho_l
              state(i,j,UMX) = rho_l*u_l
              state(i,j,UMY) = 0.e0_rt

              state(i,j,UEDEN) = rhoe_l + 0.5*rho_l*u_l*u_l
              state(i,j,UEINT) = rhoe_l
              state(i,j,UTEMP) = T_l
           else
              state(i,j,URHO) = rho_r
              state(i,j,UMX) = rho_r*u_r
              state(i,j,UMY) = 0.e0_rt

              state(i,j,UEDEN) = rhoe_r + 0.5*rho_r*u_r*u_r
              state(i,j,UEINT) = rhoe_r
              state(i,j,UTEMP) = T_r
           endif

        else if (idir == 2) then
           if (ycen <= split(2)) then
              state(i,j,URHO) = rho_l
              state(i,j,UMX) = 0.e0_rt
              state(i,j,UMY) = rho_l*u_l

              state(i,j,UEDEN) = rhoe_l + 0.5*rho_l*u_l*u_l
              state(i,j,UEINT) = rhoe_l
              state(i,j,UTEMP) = T_l
           else
              state(i,j,URHO) = rho_r
              state(i,j,UMX) = 0.e0_rt
              state(i,j,UMY) = rho_r*u_r

              state(i,j,UEDEN) = rhoe_r + 0.5*rho_r*u_r*u_r
              state(i,j,UEINT) = rhoe_r
              state(i,j,UTEMP) = T_r
           endif

        else
           call amrex_abort('invalid idir')
        endif

        state(i,j,UFS:UFS-1+nspec) = 0.0e0_rt
        state(i,j,UFS  ) = state(i,j,URHO)

     enddo
  enddo

end subroutine ca_initdata

max_stars_repo_name: harpolea/Castro

max_issues_repo_name: harpolea/Castro

max_forks_repo_name: harpolea/Castro


text:
c to fix: interpolation in rcirc does not seem to work well. after
c that, hope that things come together!

! fip: changed central sig0 to be equal (except for a factor sigroz) to
! central z dispersion

      program diskdf
      parameter(nrmax=1000)
      parameter (pi=3.1415926535)
      dimension drat(0:nrmax),
     +          dz2rat(0:nrmax),
     +          fzrat(0:nrmax)
      dimension d0rat(0:nrmax),d1rat(0:nrmax),d2rat(0:nrmax),d3rat(0:nrmax),
     +          d4rat(0:nrmax)
      dimension fplot(0:10*nrmax),rplot(0:10*nrmax)
      dimension ibuf(100)
      character*60 filename

      common /fileroot/ filename
      common /potconstants/ apot(20,0:20000), frad(20,0:20000),
     +     dr, nr, lmax, potcor
      common /legendre/ plcon(0:40)
      common /gparameters/  a, b, c, v0, q, psi00,
     +                      psiout, rho1, sigbulge2,
     +                      rmdisk, rdisk, zdisk, outdisk, drtrunc,
     +                      potcor1
      common /moreconstants/ v02, v03, rdisk2, diskconst, bulgea
      common /diskpars/ sigr0, disksr, nrdisk
      common /splines/ rr(0:nrmax),fdrat(0:nrmax),drat2(0:nrmax),
     +                 fszrat(0:nrmax),szrat2(0:nrmax), nrspl

c      write(0,*) 'Root name for potential?'
c      read(*,'(a60)') filename
      write(0,*) 'Central velocity dispersion, scale length of sigR**2?'
      read(*,*) sigroz,disksr
      write(0,*) 'number of radial steps for correction fns (min. 6)?'
      read(*,*) nrspl
      write(0,*) 'number of iterations?'
      read(*,*) niter

      filename='dbh.dat'
      call readharmfile(filename,ibuf)

      psi00=pot(0.0,0.0)
      rmax = nrdisk*dr

      sigr0=sigroz*sqrt(3.1415*zdisk*diskdensf(0.,0.)*2.0*zdisk)

      drspl = (outdisk + 2*drtrunc)/nrspl
      do ir=0,nrspl
        rr(ir)=drspl*ir
        fdrat(ir)=1
        fszrat(ir)=1
      enddo
c      drspl=(outdisk-2*drtrunc)/(nrspl-4)
c      do ir=0,nrspl-4
c        rr(ir)=drspl*ir
c        fdrat(ir)=1
c        fszrat(ir)=1
c        enddo
c      do ir=3,0,-1
c        rr(nrspl-ir)=outdisk+(2-ir)*drtrunc
c        fdrat(nrspl-ir)=1
c        fszrat(nrspl-ir)=1
c     enddo
      call splined(rr(0),fdrat(0),nrspl+1,1.e32,1.e32,drat2(0))
      call splined(rr(0),fszrat(0),nrspl+1,1.e32,1.e32,szrat2(0))

      open(12,file='cordbh.dat',status='unknown')

c Make niter iterations on the correction functions for midplane density
c and vertical velocity dispersion.
      do icor=1,niter
       do ir=0,nrspl
        r=rr(ir)
        call omekap(r,fom,fka)
        sigr=sqrt(sigr2(r))
        sigz=sqrt(sigz2(r))
        vc=r*fom
        d0=0
        dz=0
        dz2=0
        dz3=0
        dsz=0
c integrate analytically over vr and vz, numerically over vphi
        dvr=0.1*sigr
        do ivt=1,101
           vt=vc+(ivt-51)*dvr
           df=diskdf5intez(vt,r,0.0)
           c=0.333333333333*(mod(ivt,2)*2+2)
           d0=d0+c*dvr*df
           df=diskdf5intez(vt,r,zdisk)
           dz2=dz2+c*dvr*df
           enddo
        drat(ir)=d0/diskdensf(r,0.)
        dz2rat(ir)=dz2/diskdensf(r,zdisk)
        fzrat(ir)=log(diskdensf(r,0.)/diskdensf(r,zdisk))/log(d0/dz2)
c        write(*,*) 'density=',d0
        enddo
       write(0,*)
       nrd=max(1,nrspl/10)
       write(0,'(''      r: '',11f6.3)') (rr(ir),ir=0,nrspl,nrd)
       write(0,'(f6.3,''  d'',11f6.3)') 0.,(drat(ir),ir=0,nrspl,nrd)
       write(0,'(f6.3,''  d'',11f6.3)') zdisk,(dz2rat(ir),ir=0,nrspl,nrd)
       do ir=0,nrspl
        fdrat(ir)=fdrat(ir)/drat(ir)
        fszrat(ir)=fszrat(ir)/fzrat(ir)
        enddo
       call splined(rr(0),fdrat(0),nrspl+1,1.e32,1.e32,drat2(0))
       call splined(rr(0),fszrat(0),nrspl+1,1.e32,1.e32,szrat2(0))
c       do ir=0,nrspl
c        r = rr(ir)
c        call splintd(rr(0),fdrat(0),drat2(0),nrspl+1,r,fd)
c        call splintd(rr(0),fszrat(0),szrat2(0),nrspl+1,r,fz)
c        write(11,'(x,5g12.5)') rr(ir),fdrat(ir), fd, fszrat(ir), fz
c        enddo
c       do ir=0,nrspl*10-1
c        r=rr(ir/10)+(rr(ir/10+1)-rr(ir/10))*0.1*mod(ir,10)
c        call splintd(rr(0),fdrat(0),drat2(0),nrspl+1,r,fd)
c        call splintd(rr(0),fszrat(0),szrat2(0),nrspl+1,r,fz)
c        write(11,'(x,3g12.5)') r,fd,fz
c        enddo
       print*,'sigr0,sigz0:', sigr0,sqrt(sigz2(0.))
       sigr0=sigroz*sqrt(sigz2(0.))
       enddo

      rfid = 2.5*rdisk
      call omekap(rfid,fom,fka)
      sigr = sqrt(sigr2(rfid))
      sigden = diskdensf(rfid,0.0)*2.0*zdisk
      print*,'central surfdens:',diskdensf(0.,0.0)*2.0*zdisk
c
c a close estimate of the disk surface density
c
      sigrcrit = 3.36*sigden/fka
      qtoomre = sigr/sigrcrit
      write(0,*) 'Toomre Q = ',qtoomre, ' at R = 2.5 R_d'
c first call initializes the spline fitting functions by the way
      call omekap(outdisk,fom,fka)

      write(0,'('' ==> Outer dispersion, truncation dispersion:'',2g10.3)')
     +  sqrt(sigr2(outdisk)),drtrunc*fka
      if (sqrt(sigr2(outdisk)).gt.drtrunc*fka) then
        write(0,*) '** WARNING **: you have asked for an outer velocity'
        write(0,*) '  dispersion which is larger than that used in the'
        write(0,*) '  disk truncation. Strange things may happen.'
        endif

c Finally, write out a more detailed table of residuals.
      do ir=0,nrspl
        r=rr(ir)
        call omekap(r,fom,fka)
        sigr=sqrt(sigr2(r))
        sigz=sqrt(sigz2(r))
        vc=r*fom
        d0=0
        d1=0
        d2=0
        d3=0
        d4=0
c integrate analytically over vr and vz, numerically over vphi
        dvr=0.1*sigr
        do ivt=1,101
           vt=vc+(ivt-51)*dvr
           df=diskdf5intez(vt,r,0.)
           c=0.333333333333*(mod(ivt,2)*2+2)
           d0=d0+c*dvr*df
           df=diskdf5intez(vt,r,0.5*zdisk)
           d1=d1+c*dvr*df
           df=diskdf5intez(vt,r,zdisk)
           d2=d2+c*dvr*df
           df=diskdf5intez(vt,r,1.5*zdisk)
           d3=d3+c*dvr*df
           df=diskdf5intez(vt,r,2*zdisk)
           d4=d4+c*dvr*df
           enddo
        d0rat(ir)=d0/diskdensf(r,0.)
        d1rat(ir)=d1/diskdensf(r,0.5*zdisk)
        d2rat(ir)=d2/diskdensf(r,zdisk)
        d3rat(ir)=d3/diskdensf(r,1.5*zdisk)
        d4rat(ir)=d4/diskdensf(r,2*zdisk)
        enddo
       write(0,*)
       nrd=max(1,(nrspl-1)/10+1)
       write(0,'(''      r: '',11f6.3)') (rr(ir),ir=0,nrspl,nrd)
       write(0,'(f6.3,''  d'',11f6.3)') 0.,(d0rat(ir),ir=0,nrspl,nrd)
       write(0,'(f6.3,''  d'',11f6.3)') 0.5*zdisk,(d1rat(ir),ir=0,nrspl,nrd)
       write(0,'(f6.3,''  d'',11f6.3)') zdisk,(d2rat(ir),ir=0,nrspl,nrd)
       write(0,'(f6.3,''  d'',11f6.3)') 1.5*zdisk,(d3rat(ir),ir=0,nrspl,nrd)
       write(0,'(f6.3,''  d'',11f6.3)') 2*zdisk,(d4rat(ir),ir=0,nrspl,nrd)
       write(12,'(''#'', 2g17.7,x,i4)') sigr0,disksr,nrspl
       do i=0,nrspl
          write(12,'(x,3g17.7)') rr(i),fdrat(i),fszrat(i)
       enddo
c       close(11)
       close(12)

c make plots of the correction functions.

       call pgbeg(0,'/ps',2,2)
       call pgsch(2.)
       write(*,*)

       do ir=0,nrspl*10-1
          rplot(ir)=rr(ir/10)+(rr(ir/10+1)-rr(ir/10))*0.1*mod(ir,10)
       enddo
       do ir=0,nrspl*10-1
          call splintd(rr(0),fdrat(0),drat2(0),nrspl+1,rplot(ir),fplot(ir))
       enddo

       call pgenv(0.,rplot(nrspl*10-1),0.,2.,0,0)
       call modstamp
       call pgline(nrspl*10,rplot(0),fplot(0))
       call pglabel('R',' ','Surfden Correction')

       do ir=0,nrspl*10-1
          call splintd(rr(0),fszrat(0),szrat2(0),nrspl+1,rplot(ir),fplot(ir))
       enddo

       call pgenv(0.,rplot(nrspl*10-1),0.,2.,0,0)
       call pgline(nrspl*10,rplot(0),fplot(0))
       call pglabel('R',' ','sigz**2 Correction')

       call modstamp

       call pgend

      end