text:
!
! Copyright (c) Stanford University, The Regents of the University of
!               California, and others.
!
! All Rights Reserved.
!
! See Copyright-SimVascular.txt for additional details.
!
! Permission is hereby granted, free of charge, to any person obtaining
! a copy of this software and associated documentation files (the
! "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish,
! distribute, sublicense, and/or sell copies of the Software, and to
! permit persons to whom the Software is furnished to do so, subject
! to the following conditions:
!
! The above copyright notice and this permission notice shall be included
! in all copies or substantial portions of the Software.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!--------------------------------------------------------------------
!
!     Here data (mesh, solution) input/output is handled interfacing
!     with fortran-based VTK module.
!
!--------------------------------------------------------------------

      SUBROUTINE READVTU(lM, fName)
      USE COMMOD
      USE LISTMOD
      USE ALLFUN
      USE vtkXMLMod
      IMPLICIT NONE

      TYPE(mshType), INTENT(INOUT) :: lM
      CHARACTER(LEN=stdL) :: fName

      TYPE(vtkXMLType) :: vtu
      INTEGER :: iStat
      REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: tmpX

      iStat = 0
      std = " <VTK XML Parser> Loading file <"//TRIM(fName)//">"
      CALL loadVTK(vtu, fName, iStat)
      IF (iStat .LT. 0) err = "VTU file read error (init)"

      CALL getVTK_numPoints(vtu, lM%gnNo, iStat)
      IF (iStat .LT. 0) err = "VTU file read error (num points)"

      CALL getVTK_numElems(vtu, lM%gnEl, iStat)
      IF (iStat .LT. 0) err = "VTU file read error (num cells)"

      CALL getVTK_nodesPerElem(vtu, lM%eNoN, iStat)
      IF (iStat .LT. 0) err = "VTU file read error (nodes per cell)"

      ALLOCATE(lM%x(nsd,lM%gnNo),tmpX(maxNSD,lM%gnNo))
      CALL getVTK_pointCoords(vtu, tmpX, iStat)
      IF (iStat .LT. 0) err = "VTU file read error (coords)"
      lM%x(:,:) = tmpX(1:nsd,:)
      DEALLOCATE(tmpX)

      ALLOCATE(lM%gIEN(lM%eNoN,lM%gnEl))
      CALL getVTK_elemIEN(vtu, lM%gIEN, iStat)
      IF (iStat .LT. 0) err = "VTU file read error (ien)"
      lM%gIEN = lM%gIEN + 1

      CALL flushVTK(vtu)

      RETURN
      END SUBROUTINE READVTU
!--------------------------------------------------------------------
      SUBROUTINE READVTP(lFa, fName)
      USE COMMOD
      USE LISTMOD
      USE ALLFUN
      USE vtkXMLMod
      IMPLICIT NONE

      TYPE(faceType), INTENT(INOUT) :: lFa
      CHARACTER(LEN=stdL) :: fName

      TYPE(vtkXMLType) :: vtp
      INTEGER :: iStat, a, e, Ac

      REAL(KIND=8), ALLOCATABLE :: tmpX(:,:)

      iStat = 0
      std = " <VTK XML Parser> Loading file <"//TRIM(fName)//">"
      CALL loadVTK(vtp, fName, iStat)
      IF (iStat .LT. 0) err = "VTP file read error (init)"

      CALL getVTK_numPoints(vtp, lFa%nNo, iStat)
      IF (iStat .LT. 0) err = "VTP file read error (num points)"

      CALL getVTK_numElems(vtp, lFa%nEl, iStat)
      IF (iStat .LT. 0) err = "VTP file read error (num cells)"

      CALL getVTK_nodesPerElem(vtp, lFa%eNoN, iStat)
      IF (iStat .LT. 0) err = "VTP file read error (nodes per cell)"

      ALLOCATE(lFa%x(nsd,lFa%nNo), tmpX(maxNSD,lFa%nNo))
      CALL getVTK_pointCoords(vtp, tmpX, iStat)
      IF (iStat .LT. 0) err = "VTP file read error (coords)"
      lFa%x(:,:) = tmpX(1:nsd,:)
      DEALLOCATE(tmpX)

      ALLOCATE(lFa%IEN(lFa%eNoN,lFa%nEl))
      CALL getVTK_elemIEN(vtp, lFa%IEN, iStat)
      IF (iStat .LT. 0) err = "VTP file read error (ien)"

      ALLOCATE(lFa%gN(lFa%nNo))
      CALL getVTK_pointData(vtp, "GlobalNodeID", lFa%gN, iStat)
      IF (iStat .LT. 0) THEN
         DEALLOCATE(lFa%gN)
      ELSE
         DO e=1, lFa%nEl
            DO a=1, lFa%eNoN
               Ac = lFa%IEN(a,e)+1
               Ac = lFa%gN(Ac)
               lFa%IEN(a,e) = Ac
            END DO
         END DO
      END IF

      ALLOCATE(lFa%gE(lFa%nEl))
      CALL getVTK_elemData(vtp, "GlobalElementID", lFa%gE, iStat)
      IF (iStat .LT. 0) THEN
         DEALLOCATE(lFa%gE)
      ELSE
         lFa%gnEl = lFa%nEl
         ALLOCATE(lFa%gebc(1+lFa%eNoN,lFa%gnEl))
         lFa%gebc(1,:) = lFa%gE(:)
         lFa%gebc(2:1+lFa%eNoN,:) = lFa%IEN(:,:)
      END IF

      CALL flushVTK(vtp)

      RETURN
      END SUBROUTINE READVTP
!--------------------------------------------------------------------
      SUBROUTINE READVTUS(lA, lY, lD, fName)
      USE COMMOD
      USE LISTMOD
      USE ALLFUN
      USE vtkXMLMod
      IMPLICIT NONE

      REAL(KIND=8), INTENT(INOUT) :: lA(tDof,gtnNo), lY(tDof,gtnNo),
     2   lD(tDof,gtnNo)
      CHARACTER(LEN=STDL) :: fName

      TYPE(vtkXMLType) :: vtu

      INTEGER :: iStat, iEq, iOut, iM, l, s, e, a, b, Ac, nNo, oGrp
      CHARACTER(LEN=stdL) :: varName
      REAL(KIND=8), ALLOCATABLE :: tmpS(:,:), tmpGS(:,:)

      iStat = 0
      CALL loadVTK(vtu, fName, iStat)
      IF (iStat .LT. 0) err = "VTU file read error (init)"

      CALL getVTK_numPoints(vtu, nNo, iStat)
      IF (iStat .LT. 0) err = "VTU file read error (num points)"
      IF (nNo .NE. SUM(msh(:)%gnNo))
     2   err = "Incompatible mesh and "//TRIM(fName)

      DO iEq=1, nEq
         DO iOut=1, eq(iEq)%nOutput
            IF (.NOT.eq(iEq)%output(iOut)%wtn(1)) CYCLE
            l = eq(iEq)%output(iOut)%l
            s = eq(iEq)%s + eq(iEq)%output(iOut)%o
            e = s + l - 1
            varName = eq(iEq)%sym//"_"//TRIM(eq(iEq)%output(iOut)%name)

            oGrp = eq(iEq)%output(iOut)%grp
            SELECT CASE(oGrp)
            CASE (outGrp_A, outGrp_Y, outGrp_D)
               IF (l .GT. 1) THEN
                  ALLOCATE(tmpGS(maxNSD,nNo))
               ELSE
                  ALLOCATE(tmpGS(1,nNo))
               END IF
               CALL getVTK_pointData(vtu, TRIM(varName), tmpGS, iStat)
               IF (iStat .LT. 0) err ="VTU file read error (point data)"

               IF (nNo .NE. gtnNo) THEN
                  IF (l .GT. 1) THEN
                     ALLOCATE(tmpS(maxNSD,gtnNo))
                  ELSE
                     ALLOCATE(tmpS(1,gtnNo))
                  END IF
                  IF (.NOT.ALLOCATED(msh(1)%gpN)) err = "Unexpected "//
     2               "behavior. VTU file read error"
                  b = 0
                  DO iM=1, nMsh
                     DO a=1, msh(iM)%gnNo
                        b = b + 1
                        Ac = msh(iM)%gpN(a)
                        tmpS(:,Ac) = tmpGS(:,b)
                     END DO
                  END DO
                  DEALLOCATE(tmpGS)
                  ALLOCATE(tmpGS(SIZE(tmpS,1),gtnNo))
                  tmpGS = tmpS
                  DEALLOCATE(tmpS)
               END IF
            CASE DEFAULT
               CYCLE
            END SELECT

            SELECT CASE(oGrp)
            CASE (outGrp_A)
               lA(s:e,:) = tmpGS(1:l,:)
            CASE (outGrp_Y)
               lY(s:e,:) = tmpGS(1:l,:)
            CASE (outGrp_D)
               lD(s:e,:) = tmpGS(1:l,:)
            END SELECT
            DEALLOCATE(tmpGS)
         END DO
      END DO

      CALL flushVTK(vtu)

      RETURN
      END SUBROUTINE READVTUS
!--------------------------------------------------------------------
!     Read a particular point dataset from a vtu file. Point data will
!     be read and stored in lM%x array
      SUBROUTINE READVTUPDATA(lM, fName, kwrd, m, idx)
      USE COMMOD
      USE LISTMOD
      USE ALLFUN
      USE vtkXMLMod
      IMPLICIT NONE

      TYPE(mshType), INTENT(INOUT) :: lM
      CHARACTER(LEN=*) :: fName, kwrd
      INTEGER, INTENT(IN) :: m, idx

      TYPE(vtkXMLType) :: vtu
      INTEGER :: iStat, a
      REAL(KIND=8), ALLOCATABLE :: tmpR(:,:)

      iStat = 0
      std = " <VTK XML Parser> Loading file <"//TRIM(fName)//">"
      CALL loadVTK(vtu, fName, iStat)
      IF (iStat .LT. 0) err = "VTU file read error (init)"

      CALL getVTK_numPoints(vtu, a, iStat)
      IF (a .NE. lM%gnNo) err = "Mismatch in num points for "//
     2   TRIM(kwrd)

      IF (m .EQ. nsd) THEN
         ALLOCATE(tmpR(maxNSD,lM%gnNo))
         tmpR = 0D0
         CALL getVTK_pointData(vtu, TRIM(kwrd), tmpR, iStat)
         IF (iStat .LT. 0) err = "VTU file read error "//TRIM(kwrd)
         DO a=1, lM%gnNo
            lM%x((idx-1)*nsd+1:idx*nsd,a) = tmpR(1:nsd,a)
         END DO
         DEALLOCATE(tmpR)
      ELSE
         CALL getVTK_pointData(vtu, TRIM(kwrd), lM%x, iStat)
         IF (iStat .LT. 0) err = "VTU file read error "//TRIM(kwrd)
      END IF

      CALL flushVTK(vtu)

      RETURN
      END SUBROUTINE READVTUPDATA
!--------------------------------------------------------------------
!     Read a particular point dataset from a vtp file. Point data will
!     be read and stored in lFa%x array
      SUBROUTINE READVTPPDATA(lFa, fName, kwrd, m, idx)
      USE COMMOD
      USE LISTMOD
      USE ALLFUN
      USE vtkXMLMod
      IMPLICIT NONE

      TYPE(faceType), INTENT(INOUT) :: lFa
      CHARACTER(LEN=*) :: fName, kwrd
      INTEGER, INTENT(IN) :: m, idx

      TYPE(vtkXMLType) :: vtp
      INTEGER :: iStat, a

      REAL(KIND=8), ALLOCATABLE :: tmpR(:,:)

      iStat = 0
      std = " <VTK XML Parser> Loading file <"//TRIM(fName)//">"
      CALL loadVTK(vtp, fName, iStat)
      IF (iStat .LT. 0) err = "VTP file read error (init)"

      CALL getVTK_numPoints(vtp, a, iStat)
      IF (a .NE. lFa%nNo) err = "Mismatch in num points for "//
     2   TRIM(kwrd)

      IF (m .EQ. nsd) THEN
         ALLOCATE(tmpR(maxNSD,lFa%nNo))
         tmpR = 0D0
         CALL getVTK_pointData(vtp, TRIM(kwrd), tmpR, iStat)
         IF (iStat .LT. 0) err = "VTP file read error "//TRIM(kwrd)
         DO a=1, lFa%nNo
            lFa%x((idx-1)*nsd+1:idx*nsd,a) = tmpR(1:nsd,a)
         END DO
         DEALLOCATE(tmpR)
      ELSE
         CALL getVTK_pointData(vtp, TRIM(kwrd), lFa%x, iStat)
         IF (iStat .LT. 0) err = "VTU file read error "//TRIM(kwrd)
      END IF

      CALL flushVTK(vtp)

      RETURN
      END SUBROUTINE READVTPPDATA
!####################################################################
      SUBROUTINE WRITEVTU(lM, fName)
      USE COMMOD
      USE ALLFUN
      USE vtkXMLMod
      IMPLICIT NONE

      TYPE(mshType), INTENT(IN) :: lM
      CHARACTER(LEN=stdL), INTENT(IN) :: fName

      TYPE(vtkXMLType) :: vtu
      INTEGER :: iStat

      CALL vtkInitWriter(vtu, TRIM(fName), iStat)
      IF (iStat .LT. 0) err = "VTU file write error (init)"

      CALL putVTK_pointCoords(vtu, lM%x, iStat)
      IF (iStat .LT. 0) err = "VTU file write error (coords)"

      CALL putVTK_elemIEN(vtu, lM%gIEN, lM%vtkType, iStat)
      IF (iStat .LT. 0) err = "VTU file write error (ien)"

      CALL vtkWriteToFile(vtu, iStat)
      IF (iStat .LT. 0) err = "VTU file write error"
      CALL flushVTK(vtu)

      RETURN
      END SUBROUTINE WRITEVTU
!--------------------------------------------------------------------
      SUBROUTINE WRITEVTP(lFa, fName)
      USE COMMOD
      USE ALLFUN
      USE vtkXMLMod
      IMPLICIT NONE

      TYPE(faceType), INTENT(IN) :: lFa
      CHARACTER(LEN=stdL), INTENT(IN) :: fName

      TYPE(vtkXMLType) :: vtp
      INTEGER :: iStat, vtkType

      IF (nsd .EQ. 3) THEN
         SELECT CASE (lFa%eNoN)
         CASE (3)
            vtkTYpe = 5
         CASE (4)
            vtkType = 9
         END SELECT
      ELSE
         SELECT CASE (lFa%eNoN)
         CASE (2)
            vtkType = 3
         CASE (3)
            vtkType = 21
         END SELECT
      END IF

      CALL vtkInitWriter(vtp, TRIM(fName), iStat)
      IF (iStat .LT. 0) err = "VTP file write error (init)"

      CALL putVTK_pointCoords(vtp, lFa%x, iStat)
      IF (iStat .LT. 0) err = "VTP file write error (coords)"

      CALL putVTK_elemIEN(vtp, lFa%IEN, vtkType, iStat)
      IF (iStat .LT. 0) err = "VTP file write error (ien)"

      IF (ALLOCATED(lFa%gN)) THEN
         CALL putVTK_pointData(vtp, "GlobalNodeID", lFa%gN, iStat)
         IF (iStat .LT. 0) err = "VTP file write error (point data)"
      END IF

      IF (ALLOCATED(lFa%gE)) THEN
         CALL putVTK_elemData(vtp, "GlobalElementID", lFa%gE, iStat)
         IF (iStat .LT. 0) err = "VTP file write error (cell data)"
      END IF

      CALL vtkWriteToFile(vtp, iStat)
      IF (iStat .LT. 0) err = "VTP file write error"
      CALL flushVTK(vtp)

      RETURN
      END SUBROUTINE WRITEVTP
!--------------------------------------------------------------------
      SUBROUTINE WRITEVTUS(lA, lY, lD)
      USE COMMOD
      USE ALLFUN
      USE vtkXMLMod
      IMPLICIT NONE

      REAL(KIND=8), INTENT(IN) :: lA(tDof,tnNo), lY(tDof,tnNo),
     2   lD(tDof,tnNo)

      TYPE(dataType) :: d(nMsh)
      TYPE(vtkXMLType) :: vtu

      CHARACTER(LEN=stdL) :: fName
      INTEGER :: iStat, iEq, iOut, iM, a, e, Ac, Ec, nNo, nEl, s, l, ie,
     2   is, nSh, oGrp, outDof, nOut, cOut, itmp, ne, iFn, nFn
      LOGICAL :: l1, l2, l3

      CHARACTER(LEN=stdL), ALLOCATABLE :: outNames(:)
      INTEGER, ALLOCATABLE :: outS(:), tmpI(:,:)
      REAL(KIND=8), ALLOCATABLE :: tmpV(:,:)

      l1 = .FALSE.
      itmp = SUM(iblank(:))
      itmp = cm%reduce(itmp)
      IF (itmp .GT. 0) l1 = .TRUE.

      l2 = .FALSE.
      itmp = SUM(ighost(:))
      itmp = cm%reduce(itmp)
      IF (itmp .GT. 0) l2 = .TRUE.

      l3 = .FALSE.
      DO iEq=1, nEq
         IF (eq(iEq)%phys .EQ. phys_CMM .AND. ALLOCATED(Dinit)) THEN
            l3 = .TRUE.
            EXIT
         END IF
      END DO

      nOut   = 1
      outDof = nsd
      DO iEq=1, nEq
         DO iOut=1, eq(iEq)%nOutput
            IF (.NOT.eq(iEq)%output(iOut)%wtn(1)) CYCLE
            oGrp   = eq(iEq)%output(iOut)%grp
            IF (oGrp .EQ. outGrp_fN) THEN
               nFn = 1
               DO iM=1, nMsh
                  nFn = MAX(nFn, msh(iM)%nFn)
               END DO
               nOut   = nOut + nFn
               outDof = outDof + eq(iEq)%output(iOut)%l * nFn
            ELSE
               nOut   = nOut + 1
               outDof = outDof + eq(iEq)%output(iOut)%l
            END IF
         END DO
      END DO

!     iblank array for immersed solids
      IF (l1) THEN
         nOut   = nOut + 1
         outDof = outDof + 1
      END IF

!     ighost array for immersed solids/shells
      IF (l2) THEN
         nOut   = nOut + 1
         outDof = outDof + 1
      END IF

!     Initial displacements for CMM equation
      IF (l3) THEN
         nOut = nOut + 1
         outDof = outDof + nsd
      END IF

      ALLOCATE(outNames(nOut), outS(nOut+1))

!     Prepare all solultions in to dataType d
      DO iM=1, nMsh
         IF (ALLOCATED(tmpV)) DEALLOCATE(tmpV)
         ALLOCATE(d(iM)%x(outDof,msh(iM)%nNo), tmpV(maxnsd,msh(iM)%nNo))
         cOut           = 1
         outS(cOut)     = 1
         outS(cOut+1)   = nsd + 1
         outNames(cOut) = ""
         DO a=1, msh(iM)%nNo
            Ac = msh(iM)%gN(a)
            d(iM)%x(1:nsd,a) = x(:,Ac)
         END DO

         IF (l3) THEN
            cOut           = cOut + 1
            is             = outS(cOut)
            ie             = is + nsd - 1
            outS(cOut+1)   = ie + 1
            outNames(cOut) = "Initial_displacement"
            DO a=1, msh(iM)%nNo
               Ac = msh(iM)%gN(a)
               d(iM)%x(is:ie,a) = Dinit(:,Ac)
            END DO
         END IF

         DO iEq=1, nEq
            DO iOut=1, eq(iEq)%nOutput
               IF (.NOT.eq(iEq)%output(iOut)%wtn(1)) CYCLE
               l = eq(iEq)%output(iOut)%l
               s = eq(iEq)%s + eq(iEq)%output(iOut)%o
               e = s + l - 1

               cOut = cOut + 1
               is   = outS(cOut)
               ie   = is + l - 1
               outS(cOut+1) = ie + 1

               oGrp = eq(iEq)%output(iOut)%grp
               outNames(cOut) = TRIM(eq(iEq)%output(iOut)%name)

               SELECT CASE (oGrp)
                  CASE (outGrp_NA)
                  err = "Undefined output grp in VTK"
               CASE (outGrp_A)
                  DO a=1, msh(iM)%nNo
                     Ac = msh(iM)%gN(a)
                     d(iM)%x(is:ie,a) = lA(s:e,Ac)
                  END DO
               CASE (outGrp_Y)
                  DO a=1, msh(iM)%nNo
                     Ac = msh(iM)%gN(a)
                     d(iM)%x(is:ie,a) = lY(s:e,Ac)
                  END DO
               CASE (outGrp_D)
                  DO a=1, msh(iM)%nNo
                     Ac = msh(iM)%gN(a)
                     d(iM)%x(is:ie,a) = lD(s:e,Ac)
                  END DO
               CASE (outGrp_WSS)
                  CALL BPOST(msh(iM), tmpV, lY, lD, oGrp)
                  DO a=1, msh(iM)%nNo
                     d(iM)%x(is:ie,a) = tmpV(1:l,a)
                  END DO
               CASE (outGrp_trac)
                  CALL BPOST(msh(iM), tmpV, lY, lD, oGrp)
                  DO a=1, msh(iM)%nNo
                     d(iM)%x(is:ie,a) = tmpV(1:l,a)
                  END DO
               CASE (outGrp_stress)
                  IF (ALLOCATED(tmpV)) DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(nstd,msh(iM)%nNo))
                  tmpV = 0D0
                  IF (pstEq) THEN
                     DO a=1, msh(iM)%nNo
                        Ac = msh(iM)%gN(a)
                        tmpV(:,a) = pS0(:,Ac)
                     END DO
                  END IF
                  IF (.NOT.cmmInit) THEN
                     CALL TPOST(msh(iM), l, tmpV, lD, iEq, oGrp)
                  END IF
                  DO a=1, msh(iM)%nNo
                     d(iM)%x(is:ie,a) = tmpV(:,a)
                  END DO
                  DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(maxnsd,msh(iM)%nNo))
               CASE (outGrp_vort, outGrp_eFlx, outGrp_hFlx,
     2            outGrp_stInv, outGrp_vortex, outGrp_Visc)
                  CALL POST(msh(iM), tmpV, lY, lD, oGrp, iEq)
                  DO a=1, msh(iM)%nNo
                     d(iM)%x(is:ie,a) = tmpV(1:l,a)
                  END DO
               CASE (outGrp_absV)
                  DO a=1, msh(iM)%nNo
                     Ac = msh(iM)%gN(a)
                     d(iM)%x(is:ie,a) = lY(1:nsd,Ac)
     2                                - lY(nsd+2:2*nsd+1,Ac)
                  END DO
               CASE (outGrp_fN)
                  cOut = cOut - 1
                  IF (ALLOCATED(tmpV)) DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(nFn*nsd,msh(iM)%nNo))
                  tmpV = 0D0
                  IF (msh(iM)%nFn .NE. 0)
     2               CALL FIBDIRPOST(msh(iM), nFn, tmpV, lD, iEq)
                  DO iFn=1, nFn
                     cOut = cOut + 1
                     is   = outS(cOut)
                     ie   = is + l - 1
                     outS(cOut+1)   = ie + 1
                     outNames(cOut) = eq(iEq)%sym//"_"//
     2                  TRIM(eq(iEq)%output(iOut)%name)//STR(iFn)
                     DO a=1, msh(iM)%nNo
                        d(iM)%x(is:ie,a) =
     2                     tmpV((iFn-1)*nsd+1:iFn*nsd,a)
                     END DO
                  END DO
                  DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(maxnsd,msh(iM)%nNo))
               CASE (outGrp_fA)
                  IF (ALLOCATED(tmpV)) DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(1,msh(iM)%nNo))
                  tmpV = 0D0
                  IF (msh(iM)%nFn .EQ. 2)
     2               CALL FIBALGNPOST(msh(iM), tmpV, lD, iEq)
                  DO a=1, msh(iM)%nNo
                     d(iM)%x(is,a) = tmpV(1,a)
                  END DO
                  DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(maxnsd,msh(iM)%nNo))
               CASE (outGrp_J)
                  IF (ALLOCATED(tmpV)) DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(1,msh(iM)%nNo))
                  tmpV = 0D0
                  CALL TPOST(msh(iM), l, tmpV, lD, iEq, oGrp)
                  DO a=1, msh(iM)%nNo
                     d(iM)%x(is,a) = tmpV(1,a)
                  END DO
                  DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(maxnsd,msh(iM)%nNo))
               CASE (outGrp_F)
                  IF (ALLOCATED(tmpV)) DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(nsd*nsd,msh(iM)%nNo))
                  tmpV = 0D0
                  CALL TPOST(msh(iM), l, tmpV, lD, iEq, oGrp)
                  DO a=1, msh(iM)%nNo
                     d(iM)%x(is:ie,a) = tmpV(:,a)
                  END DO
                  DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(maxnsd,msh(iM)%nNo))
               CASE (outGrp_divV)
                  IF (ALLOCATED(tmpV)) DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(1,msh(iM)%nNo))
                  tmpV = 0D0
                  CALL DIVPOST(msh(iM), tmpV, lY, lD, iEq)
                  DO a=1, msh(iM)%nNo
                     d(iM)%x(is,a) = tmpV(1,a)
                  END DO
                  DEALLOCATE(tmpV)
                  ALLOCATE(tmpV(maxnsd,msh(iM)%nNo))
               CASE DEFAULT
                  err = "Undefined output"
               END SELECT
            END DO
         END DO

         IF (l1) THEN
            cOut = cOut + 1
            is   = outS(cOut)
            ie   = is
            outS(cOut+1)   = ie + 1
            outNames(cOut) = "IBLANK"
            DO a=1, msh(iM)%nNo
               Ac = msh(iM)%gN(a)
               d(iM)%x(is:ie,a) = REAL(iblank(Ac), KIND=8)
            END DO
         END IF

         IF (l2) THEN
            cOut = cOut + 1
            is   = outS(cOut)
            ie   = is
            outS(cOut+1)   = ie + 1
            outNames(cOut) = "IGHOST"
            DO a=1, msh(iM)%nNo
               Ac = msh(iM)%gN(a)
               d(iM)%x(is:ie,a) = REAL(ighost(Ac), KIND=8)
            END DO
         END IF
      END DO

!     Integrate data from all processors
      nNo = 0
      nEl = 0
      DO iM=1, nMsh
         IF (msh(iM)%eType .EQ. eType_NRB) THEN
            CALL INTNRBDATA(msh(iM), d(iM), outDof)
         ELSE
            CALL INTMSHDATA(msh(iM), d(iM), outDof)
         END IF
         nNo = nNo +  d(iM)%nNo
         nEl = nEl +  d(iM)%nEl
      END DO
      IF (cm%slv()) THEN
         savedOnce = .TRUE.
         RETURN
      END IF

      DEALLOCATE(tmpV)
      ALLOCATE(tmpV(maxnsd, nNo))

!     Writing to vtu file (master only)
      IF (cTS .GE. 1000) THEN
         fName = STR(cTS)
      ELSE
         WRITE(fName,'(I3.3)') cTS
      END IF

      fName = TRIM(saveName)//"_"//TRIM(ADJUSTL(fName))//".vtu"
      dbg = "Writing VTU"

      CALL vtkInitWriter(vtu, TRIM(fName), iStat)
      IF (iStat .LT. 0) err = "VTU file write error (init)"

!     Writing the position data
      iOut = 1
      s   = outS(iOut)
      e   = outS(iOut+1)-1
      nSh  = 0
      tmpV = 0D0
      DO iM=1, nMsh
         DO a=1, d(iM)%nNo
            tmpV(1:nsd,a+nSh) = d(iM)%gx(s:e,a)
         END DO
         nSh = nSh + d(iM)%nNo
      END DO
      CALL putVTK_pointCoords(vtu, tmpV(1:nsd,:), iStat)
      IF (iStat .LT. 0) err = "VTU file write error (coords)"

!     Writing the connectivity data
      nSh = -1
      DO iM=1, nMsh
         ALLOCATE(tmpI(d(iM)%eNoN, d(iM)%nEl))
         DO e=1, d(iM)%nEl
            tmpI(:,e) = d(iM)%IEN(:,e) + nSh
         END DO
         CALL putVTK_elemIEN(vtu, tmpI, d(iM)%vtkType, iStat)
         IF (iStat .LT. 0) err = "VTU file write error (ien)"
         DEALLOCATE(tmpI)
         nSh = nSh + d(iM)%nNo
      END DO

!     Writing all solutions
      DO iOut=2, nOut
         IF (ALLOCATED(tmpV)) DEALLOCATE(tmpV)
         s = outS(iOut)
         e = outS(iOut+1) - 1
         l = e - s + 1
         ALLOCATE(tmpV(l, nNo))
         nSh = 0
         DO iM=1, nMsh
            DO a=1, d(iM)%nNo
               tmpV(:,a+nSh) = d(iM)%gx(s:e,a)
            END DO
            nSh = nSh + d(iM)%nNo
         END DO
         CALL putVTK_pointData(vtu, outNames(iOut), tmpV, iStat)
         IF (iStat .LT. 0) err = "VTU file write error (point data)"
      END DO

!     Write element-based variables
      ne = 1
      IF (.NOT.savedOnce .OR. nMsh.GT.1) THEN
         ALLOCATE(tmpI(1,nEl))
         IF (.NOT.savedOnce) THEN
            savedOnce = .TRUE.
            ne = ne + 2

!        Write partition data
            IF (.NOT.cm%seq()) THEN
               Ec = 0
               DO iM=1, nMsh
                  DO e=1, d(iM)%nEl
                     Ec = Ec + 1
                     tmpI(1,Ec) = d(iM)%xe(e,2)
                  END DO
               END DO
               CALL putVTK_elemData(vtu, 'Proc_ID', tmpI, iStat)
               IF (iStat .LT. 0) err = "VTU file write error (proc id)"
            END IF
         ELSE
            ne = ne + 1
         END IF

!     Write the domain ID
         IF (ALLOCATED(dmnID)) THEN
            Ec = 0
            DO iM=1, nMsh
               DO e=1, d(iM)%nEl
                  Ec = Ec + 1
                  tmpI(1,Ec) = d(iM)%xe(e,1)
               END DO
            END DO
            CALL putVTK_elemData(vtu, 'Domain_ID', tmpI, iStat)
            IF (iStat .LT. 0) err = "VTU file write error (dom id)"
         END IF

!     Write the mesh ID
         IF (nMsh .GT. 1) THEN
            Ec = 0
            DO iM=1, nMsh
               DO e=1, d(iM)%nEl
                  Ec = Ec + 1
                  tmpI(1,Ec) = iM
               END DO
            END DO
            CALL putVTK_elemData(vtu, 'Mesh_ID', tmpI, iStat)
            IF (iStat .LT. 0) err = "VTU file write error (mesh id)"
         END IF
         DEALLOCATE(tmpI)
      END IF

!     Write element ghost cells if necessary
      IF (l1 .OR. l2) THEN
         ALLOCATE(tmpI(1,nEl))
         Ec = 0
         DO iM=1, nMsh
            DO e=1, d(iM)%nEl
               Ec = Ec + 1
               tmpI(1,Ec) = d(iM)%xe(e,ne)
            END DO
         END DO
         CALL putVTK_elemData(vtu, 'EGHOST', tmpI, iStat)
         IF (iStat .LT. 0) err = "VTU file write error (EGHOST)"
         DEALLOCATE(tmpI)
      END IF

      DO iM=1, nMsh
         CALL DESTROY(d(iM))
      END DO

      CALL vtkWriteToFile(vtu, iStat)
      IF (iStat .LT. 0) err = "VTU file write error"

      CALL flushVTK(vtu)

      RETURN
      END SUBROUTINE WRITEVTUS
!####################################################################
!     This routine prepares data array of a regular mesh
      SUBROUTINE INTMSHDATA(lM, d, outDof)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE

      TYPE(mshType), INTENT(IN) :: lM
      TYPE(dataType), INTENT(INOUT) :: d
      INTEGER, INTENT(IN) :: outDof

      INTEGER :: e, i, ierr, Ec, m
      INTEGER, ALLOCATABLE :: sCount(:)

      d%eNoN    = lM%eNoN
      d%vtkType = lM%vtkType
      IF (cm%mas()) THEN
         d%nNo  = lM%gnNo
         d%nEl  = lM%gnEl
         d%eNoN = lM%eNoN
      ELSE
         d%nNo  = 0
         d%nEl  = 0
      END IF
      ALLOCATE(d%IEN(d%eNoN,d%nEl), d%gx(outDof,d%nNo))

      d%gx = GLOBAL(lM, d%x)
      DEALLOCATE(d%x)

      IF (cm%mas()) THEN
         DO e=1, lM%gnEl
            Ec = lM%otnIEN(e)
            d%IEN(:,e) = lM%gIEN(:,Ec)
         END DO
      END IF

      m = 1
      IF (.NOT.savedOnce .OR. nMsh.GT.1) THEN
         IF (.NOT.savedOnce) THEN
            m = m + 2
         ELSE
            m = m + 1
         END IF
         ALLOCATE(d%xe(d%nEl,m))
         IF (ALLOCATED(dmnId)) THEN
            ALLOCATE(sCount(cm%np()))
            DO i=1, cm%np()
               sCount(i) = lM%eDist(i) - lM%eDist(i-1)
            END DO
            CALL MPI_GATHERV(lM%eId, lM%nEl, mpint, d%xe(:,1), sCount,
     2         lM%eDist, mpint, master, cm%com(), ierr)

            IF (cm%mas()) THEN
               DEALLOCATE(sCount)
               ALLOCATE(sCount(lM%gnEl))
               sCount = d%xe(:,1)
               DO e=1, lM%gnEl
                  Ec = lM%otnIEN(e)
                  d%xe(e,1) = sCount(Ec)
               END DO
            END IF
            DEALLOCATE(sCount)
         ELSE
            d%xe(:,1) = 1
         END IF

         IF (.NOT.savedOnce) THEN
            IF (cm%mas()) THEN
               IF (.NOT.cm%seq()) THEN
                  i = 0
                  DO e=1, lM%gnEl
                     IF (e .GT. lM%eDist(i)) THEN
                        DO
                           i = i + 1
                           IF (e .LE. lM%eDist(i)) EXIT
                        END DO
                     END IF
                     d%xe(e,2) = i
                  END DO
                  ALLOCATE(sCount(lM%gnEl))
                  sCount = d%xe(:,2)
                  DO e=1, lM%gnEl
                     Ec = lM%otnIEN(e)
                     d%xe(e,2) = sCount(Ec)
                  END DO
                  DEALLOCATE(sCount)
               END IF
            END IF
         END IF
      END IF

      IF (ALLOCATED(lM%iGC)) THEN
         IF (.NOT.ALLOCATED(d%xe)) ALLOCATE(d%xe(d%nEl,m))
         ALLOCATE(sCount(cm%np()))
         DO i=1, cm%np()
            sCount(i) = lM%eDist(i) - lM%eDist(i-1)
         END DO
         CALL MPI_GATHERV(lM%iGC, lM%nEl, mpint, d%xe(:,m), sCount,
     2      lM%eDist, mpint, master, cm%com(), ierr)

         IF (cm%mas()) THEN
            DEALLOCATE(sCount)
            ALLOCATE(sCount(lM%gnEl))
            sCount = d%xe(:,m)
            DO e=1, lM%gnEl
               Ec = lM%otnien(e)
               d%xe(e,m) = sCount(Ec)
            END DO
         END IF
         DEALLOCATE(sCount)
      END IF

      RETURN
      END SUBROUTINE INTMSHDATA
!--------------------------------------------------------------------
!     This routine will interpolates NURBS into a BIL/BRK mesh, so you
!     can write it into a VTK file.
      SUBROUTINE INTNRBDATA(lM, d, outDof)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE

      TYPE(mshType), INTENT(IN) :: lM
      TYPE(dataType), INTENT(INOUT) :: d
      INTEGER, INTENT(IN) :: outDof

      LOGICAL flag, clcDmn
      INTEGER ie, e, ex, ey, ez, s, sx, sy, sz, ia, Ac, a, ax, ay, az,
     2   jx, jy, jz, ierr, i, nShft, insd
      INTEGER, ALLOCATABLE :: IEN(:,:), tDmnId(:), sCe(:), dise(:),
     2   sCn(:), disn(:)
      REAL(KIND=8), ALLOCATABLE :: N(:,:), tmpX(:,:)

      insd = nsd
      IF (lM%lShl) insd = nsd - 1

      jy = lM%bs(2)%nEl*(lM%bs(2)%nSl - 1) + 1
      IF (insd .EQ. 2) THEN
         jz     = 0
         d%eNoN = 4
      ELSE IF (insd .EQ. 3) THEN
         jz     = lM%bs(3)%nEl*(lM%bs(3)%nSl - 1) + 1
         d%eNoN = 8
      END IF
      ALLOCATE(sCe(cm%np()), dise(0:cm%np()), sCn(cm%np()),
     2   disn(0:cm%np()))

!     Preparing sCn, sCe for the future communications
      dise = 0
      disn = 0
      DO i=1, cm%np()
         e = lM%bs(1)%nEl*(lM%eDist(i) - lM%eDist(i-1))/lM%gnEl
         jx = e*(lM%bs(1)%nSl - 1) + 1
         IF (insd .EQ. 2) jz = 2
         sCe(i) = (jx - 1)*(jy - 1)*(jz - 1)

         IF (insd .EQ. 2) jz = 1
         IF (lM%eDist(i).NE.lM%gnEl .OR. e.EQ.0) jx = jx - 1
         sCn(i) = jx*jy*jz

         dise(i) = dise(i-1) + sCe(i)
         disn(i) = disn(i-1) + sCn(i)
      END DO
      IF (cm%mas()) THEN
         d%nEl = dise(cm%np())
         d%nNo = disn(cm%np())
      ELSE
         d%nEl = 0
         d%nNo = 0
      END IF
      ALLOCATE(d%IEN(d%eNoN,d%nEl), d%gx(outDof,d%nNo),
     2   IEN(d%eNoN,sCe(cm%tF())), tmpX(outDof,sCn(cm%tF())),
     3   N(lM%eNoN,lM%nSl))

!     I am calculating the last column only if this is the last
!     partition that has a section of this mesh. This is to prevent
!     repetition associated with the node on the boundary of partitions
      flag = .FALSE.
      IF (lM%eDist(cm%tF()).NE.lM%gnEl .OR. sCe(cm%tF()).EQ.0) THEN
         flag = .TRUE.
      END IF
!     This would be the number of nodes in the lower procesors
      nShft = disn(cm%id())

      clcDmn = .FALSE.
      IF (.NOT.savedOnce) THEN
         ALLOCATE(d%xe(d%nEl,2))
         IF (cm%mas()) THEN
!     Based on the two posible fields
            IF (.NOT.cm%seq()) THEN
               i = 0
               DO e=1, d%nEl
                  IF (e .GT. dise(i)) THEN
                     DO
                        i = i + 1
                        IF (e .LE. dise(i)) EXIT
                     END DO
                  END IF
                  d%xe(e,1) = i
               END DO
            END IF
         END IF
         IF (ALLOCATED(dmnId)) THEN
            ALLOCATE(tDmnId(sCe(cm%tF())))
            IF (ALLOCATED(lM%eId)) THEN
               clcDmn = .TRUE.
            ELSE
               tDmnId = 1
            END IF
         END IF
      END IF

      ie = 0
      IF (insd .EQ. 2) THEN
!     I am using d%eType = eType_BIL
         d%vtkType = 9
         DO e=1, lM%nEl
            ex =    (e-1)/lM%bs(2)%nEl  + 1
            ey = MOD(e-1, lM%bs(2)%nEl) + 1
            ax = (ex - 1)*(lM%bs(1)%nSl - 1)
            DO sx=1, lM%bs(1)%nSl-1
               ax = ax + 1
               ay = (ey - 1)*(lM%bs(2)%nSl - 1)
               DO sy=1, lM%bs(2)%nSl-1
                  ie = ie + 1
                  ay = ay + 1
                  ia = (ax - 1)*jy + ay + nShft

                  IEN(1,ie) = ia + jy + 1
                  IEN(2,ie) = ia      + 1
                  IEN(3,ie) = ia
                  IEN(4,ie) = ia + jy
                  IF (clcDmn) tDmnId(ie) = lM%eId(e)
               END DO
            END DO

            s  = 0
            CALL NRBNNS(lM, N, e)
            ax = (ex - 1)*(lM%bs(1)%nSl - 1)
            DO sx=1, lM%bs(1)%nSl
               IF (flag .AND. sx.EQ.lM%bs(1)%nSl) CYCLE
               ax = ax + 1
               ay = (ey - 1)*(lM%bs(2)%nSl - 1)
               DO sy=1, lM%bs(2)%nSl
                  s  = s + 1
                  ay = ay + 1
                  ia = (ax - 1)*jy + ay
                  tmpX(:,ia) = 0D0
                  DO a=1, lM%eNoN
                     Ac = lM%IEN(a,e)
                     Ac = lM%lN(Ac)
                     tmpX(:,ia) = tmpX(:,ia) + d%x(:,Ac)*N(a,s)
                  END DO
               END DO
            END DO
         END DO
      ELSE
!     I am using d%eType = eType_BRK
         d%vtkType = 12
         DO e=1, lM%nEl
            ex =   (e-1)/(lM%bs(2)%nEl*lM%bs(3)%nEl) + 1
            ey = (MOD(e-1,lM%bs(2)%nEl*lM%bs(3)%nEl))/lM%bs(3)%nEl + 1
            ez = MOD(e-1, lM%bs(3)%nEl) + 1
            ax = (ex - 1)*(lM%bs(1)%nSl - 1)
            DO sx=1, lM%bs(1)%nSl-1
               ax = ax + 1
               ay = (ey - 1)*(lM%bs(2)%nSl - 1)
               DO sy=1, lM%bs(2)%nSl-1
                  ay = ay + 1
                  az = (ez - 1)*(lM%bs(3)%nSl - 1)
                  DO sz=1, lM%bs(3)%nSl-1
                     az = az + 1
                     ie = ie + 1
                     ia = ((ax - 1)*jy + ay - 1)*jz + az + nShft

                     IEN(1,ie) = ia + jy*jz + jz + 1
                     IEN(2,ie) = ia         + jz + 1
                     IEN(3,ie) = ia         + jz
                     IEN(4,ie) = ia + jy*jz + jz
                     IEN(5,ie) = ia + jy*jz      + 1
                     IEN(6,ie) = ia              + 1
                     IEN(7,ie) = ia
                     IEN(8,ie) = ia + jy*jz
                     IF (clcDmn) tDmnId(ie) = lM%eId(e)
                  END DO
               END DO
            END DO

            s  = 0
            CALL NRBNNS(lM, N, e)
            ax = (ex - 1)*(lM%bs(1)%nSl - 1)
            DO sx=1, lM%bs(1)%nSl
               IF (flag .AND. sx.EQ.lM%bs(1)%nSl) CYCLE
               ax = ax + 1
               ay = (ey - 1)*(lM%bs(2)%nSl - 1)
               DO sy=1, lM%bs(2)%nSl
                  ay = ay + 1
                  az = (ez - 1)*(lM%bs(3)%nSl - 1)
                  DO sz=1, lM%bs(3)%nSl
                     az = az + 1
                     s  = s + 1
                     ia = ((ax - 1)*jy + ay - 1)*jz + az
                     tmpX(:,ia) = 0D0
                     DO a=1, lM%eNoN
                        Ac = lM%IEN(a,e)
                        Ac = lM%lN(Ac)
                        tmpX(:,ia) = tmpX(:,ia) + d%x(:,Ac)*N(a,s)
                     END DO
                  END DO
               END DO
            END DO
         END DO
      END IF
      DEALLOCATE(d%x)

      IF (cm%seq()) THEN
         d%gx  = tmpX
         d%IEN = IEN
      ELSE
         IF (ALLOCATED(tDmnId)) CALL MPI_GATHERV(tDmnId, sCe(cm%tF()),
     2      mpint, d%xe(:,2), sCe, dise, mpint, master, cm%com(), ierr)

!     First collecting all the elements connectivity
         sCe  = sCe*d%eNoN
         sCn  = sCn*outDof
         dise = dise*d%eNoN
         disn = disn*outDof
         CALL MPI_GATHERV(IEN, sCe(cm%tF()), mpint, d%IEN, sCe, dise,
     2      mpint, master, cm%com(), ierr)
!     Now collecting the solutions
         CALL MPI_GATHERV(tmpX, sCn(cm%tF()), mpreal, d%gx, sCn, disn,
     2      mpreal, master, cm%com(), ierr)
      END IF

      RETURN
      END SUBROUTINE INTNRBDATA
!####################################################################
!     This routine calculates average of VTU files.
      SUBROUTINE CALCAVE
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE

      LOGICAL :: lExist
      INTEGER :: fid, n, lTS, fTS, iTS
      CHARACTER(LEN=stdL) fName
      REAL(KIND=8) :: cntr

      REAL(KIND=8), ALLOCATABLE :: tmpA(:,:), tmpY(:,:), tmpD(:,:),
     2   Ag(:,:), Yg(:,:), Dg(:,:)

      fid  = 1
      fTS  = nTS
      lTS  = 0
      cntr = 0D0
      iTS  = 0
      IF (zeroAve) iTS = rsTS
      IF (cm%mas()) THEN
         ALLOCATE(tmpA(tDof,gtnNo), tmpY(tDof,gtnNo),
     2      tmpD(tDof,gtnNo))
      ELSE
         ALLOCATE(tmpA(0,0), tmpY(0,0), tmpD(0,0))
      END IF
      ALLOCATE(Ag(tDof,tnNo), Yg(tDof,tnNo), Dg(tDof,tnNo))

      Ag   = 0D0
      Yg   = 0D0
      Dg   = 0D0
      std = " Computing average quantities"
      DO n=iTS, nTS, saveIncr
         IF (n .LT. saveATS) CYCLE
         IF (n .LT. fTS) fTS = n
         IF (n .GT. lTS) lTS = n
         cntr = cntr + 1D0
         IF (cm%mas()) THEN
            IF (n .GE. 1000) THEN
               fName = STR(n)
            ELSE
               WRITE(fName,'(I3.3)') n
            END IF
            fName = TRIM(saveName)//"_"//TRIM(ADJUSTL(fName))//".vtu"
            INQUIRE (FILE=TRIM(fName), EXIST=lExist)
            IF (lExist) THEN
               CALL READVTUS(tmpA, tmpY, tmpD, fName)
            ELSE
               err = "File "//TRIM(fName)//" is missing"
            END IF
         END IF
         Ag = Ag + LOCAL(tmpA)
         Yg = Yg + LOCAL(tmpY)
         Dg = Dg + LOCAL(tmpD)
      END DO
      n = NINT(cntr)
      IF (n .EQ. 0) RETURN
      Ag = Ag/cntr
      Yg = Yg/cntr
      Dg = Dg/cntr
      cTS = lTS
      fName = saveName
      saveName = TRIM(saveName)//"_ave_"//STR(fTS)//"_"//
     2   STR(saveIncr)
      CALL WRITEVTUS(Ag, Yg, Dg)
      saveName = fName

      RETURN
      END SUBROUTINE CALCAVE
!####################################################################


text:
!-----------------------------------------------------------------------
!     Program:       FIELDLINES
!     Authors:       S. Lazerson
!     Date:          02/21/2012
!     Description:   The VMEC2V690 code reads a VMEC2000 WOUT file from
!                    the current LIBSTELL library and writes out that
!                    file in v6.90 format.
!     References:
!-----------------------------------------------------------------------
      PROGRAM VMEC2V690
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE read_wout_mod
      USE safe_open_mod
!-----------------------------------------------------------------------
!     Local Variables
!          numargs      Number of input arguments
!          i            Index
!          arg_len      Length of input strings
!          arg1         Input file
!          args         Input arguments
!-----------------------------------------------------------------------
      IMPLICIT NONE
      integer                                      :: numargs,i,ier,&
                                                      istat, js, mn, &
                                                      fid
      integer, parameter                           :: arg_len =256
      character*(arg_len)                          :: wout_file
      character*(arg_len)                          :: input_ext
      character*(arg_len),allocatable,dimension(:) :: args
!-----------------------------------------------------------------------
!     Begin Program
!-----------------------------------------------------------------------
      ! First Handle the input arguments
      CALL GETCARG(1, wout_file, numargs)
      wout_file = TRIM(wout_file)
      CALL read_wout_file(wout_file,ier)
      IF (ier /=0) THEN
         WRITE(*,*) 'Error Reading File: ',TRIM(wout_file)
         STOP
      END IF
      nextcur = 0
      IF (ALLOCATED(extcur)) nextcur = SIZE(extcur)
      wout_file = 'wout.' // TRIM(input_extension) // '_v690'
      fid = 32
      call safe_open(fid, istat, wout_file, 'replace', 'formatted')
      WRITE(fid,'(a15,a)') 'VMEC VERSION = ','6.90'
      WRITE(fid,*) wb, wp, gamma, pfac, rmax_surf, rmin_surf, zmax_surf
      WRITE(fid,*) nfp, ns, mpol, ntor, mnmax, itfsq, niter, iasym, &
                   ireconstruct, ierr_vmec
      WRITE(fid,*) imse, itse, nbsets, nobd, nextcur, nstore_seq
      IF (nbsets .gt. 0) WRITE(fid,*) (nbfld(i),i=1,nbsets)
      WRITE(fid,'(a)') mgrid_file
      DO js = 1, ns
         DO mn = 1, mnmax
            IF (js == 1) WRITE(fid,*) INT(xm(mn)),INT(xn(mn))
            WRITE(fid,*) rmnc(mn,js),zmns(mn,js),lmns(mn,js),&
                         bmnc(mn,js),gmnc(mn,js),bsubumnc(mn,js),&
                         bsubvmnc(mn,js),bsubsmns(mn,js),&
                         bsupumnc(mn,js),bsupvmnc(mn,js),&
                         currvmnc(mn,js)
            IF (lasym) WRITE(fid,*) rmns(mn,js),zmnc(mn,js),lmnc(mn,js)
         END DO
      END DO
      WRITE(fid,*) (iotas(js), mass(js), pres(js), &
         beta_vol(js), phip(js), buco(js), bvco(js), phi(js), vp(js), &
         overr(js), jcuru(js), jcurv(js), specw(js),js=2,ns)
      WRITE(fid,*) aspect, betatot, betapol, betator, betaxis, b0
      WRITE(fid,*) isigng
      WRITE(fid,'(a)') input_extension
      WRITE(fid,*) IonLarmor, VolAvgB, rbtor0, rbtor, Itor, &
        Aminor, Rmajor, Volume
      WRITE(fid,*) (Dmerc(js), Dshear(js), Dwell(js), Dcurr(js), &
             Dgeod(js), equif(js), js=2,ns-1)
      IF (nextcur .gt. 0) then
         WRITE(fid,*) (extcur(i),i=1,nextcur)
         WRITE(fid,*) (curlabel(i),i=1,nextcur)
      ENDIF
      WRITE(fid,*) (fsqt(i),wdot(i),i=1,nstore_seq)
      WRITE(fid,*) (jdotb(js),bdotgradv(js),js=1,ns)
      CLOSE(fid)
!-----------------------------------------------------------------------
!     End Program
!-----------------------------------------------------------------------
      END PROGRAM VMEC2V690

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
#include "assert.h"

#ifdef DEBUG
#define pure
#define elemental
#endif

module equation_of_state

contains

!> \brief Calculate equation of state variables
!! - AF is the logarithm of the variable f that parametrizes the el. degeneracy
!! - AT is the logarithm of the temperature
!! - Composition variables are passed using a module
pure subroutine statef ( af, at, xa, abund, eos )
   use real_kind
   use constants
   use settings
   use opacity
   use neutrinos
   use atomic_data
   use eostate_types

   implicit none
   real(double), intent(in) :: af, at
   real(double), intent(in) :: xa(9)
   type(eostate), intent(out) :: eos
   type(abundance), intent(out) :: abund
   real(double), parameter :: epsilon = 1.0d-30

   integer :: i,j
   real(double) :: f,wf,psi,g,se,sef,set,ue,sp,up,spt,spf
   real(double) :: xip,paf,puf,pwf,ff,ti,de,dc,dvt,dvf,dpa,dpat,dpaf
   real(double) :: dsat,dsaf,dua,dv,nef,si,sif,sit,ui,vm,sha
   real(double) :: dsa,net,vx,sjha,scha,nx,nxf,nxt,dh2ti,expdh2ti
   real(double) :: d1ti,d2ti,d3ti,zet,dzh2t,dzh2tt,zh2,zh2t,zh2tt,zh2s
   real(double) :: h2a,h2bt,h2at,qa,qb,qc,hg,hi,np,h2,ni,qd,qaf,qat,qf,qbf,qbt
   real(double) :: hgf,hgt,hif,hit,h2f,h2t,db,dl,dpb,dpbt,dpbf,dsb,dsbt,dsbf,dub
   real(double) :: tcr,p0,pi,t2,t4,b,dsf,dst,q,scv,gamma,tst,frho,tf
   real(double) :: llambda,mu_plasma,mu_rad,nu
   real(double) :: dvt_over_dvf,scha_over_sha,sjha_over_sha
   real(double) :: inv_ne
   real(double) :: ha(26), va(26)

   real(double) :: ap, arho, u, u1, u2, u3, p, rho, fk, T, sf, st, zt, grada
   real(double) :: scp, rf, rt, xhi, s, pr, pg, pf, pt, uf, ut, en
   real(double) :: wmu, delta, phi, fkt, fkr, prandtl

   real(double) :: fdi(16)
   real(double) :: re, pe, qe, ret, pet, qet, ref, pef, qef, xtt, xft, xff, xttt, xftt, xfft, xfff
   real(double) :: rp, pp, qp, rpt, ppt, qpt, rpf, ppf, qpf

   real(double) :: na(9), neo, nio, nzz, avm, ne, xai(9, 26)
   real(double) :: lna(9)           ! log(na), used several times
   real(double) :: mu(9)            ! Chemical potentials (in units of kT)
   real(double) :: zeta(9)          ! Average charge for atomic species i
   real(double) :: dlrdni(9)        ! Derivative of density, -(dlnrho/dNi)_P,T
   real(double) :: dsdnip(9)        ! Derivative of entropy, (dS/dNi)_p
   real(double) :: dsdnif(9)        ! Derivative of entropy, (dS/dNi)_f
   real(double) :: nif(9)           ! dNi/dlnf, for atomic species i
   real(double) :: nit(9)           ! dNi/dlnT, for atomic species i

   i = 0  ! Always defined

   ! Calculate F and G parameters, needed to calculate Fermi-Dirac integrals
   f = exp(af)
   T = exp(at)
   !T = MAX(1000.0, T)         ! Don't trust T < 1000K
   uf = f/(1.0d0 + f)
   wf = sqrt(1.0d0 + f)
   psi = af + 2.0d0*(wf - log(1.0d0 + wf))
   g = cte*T*wf

   ! Evaluate Fermi-Dirac integrals according to Eggleton, Faulkner &
   ! Flannery (1973): re, pe, se and ue correspond to rho*, p*, s* and u*.
   ! psi is the usual electron-degeneracy parameter
   call fdirac ( f, g, fdi )
   re   = fdi(1); pe   = fdi(2); qe   = fdi(3)
   ret  = fdi(4); pet  = fdi(5); qet  = fdi(6)
   ref  = fdi(7); pef  = fdi(8); qef  = fdi(9)
   xtt  = fdi(10); xft  = fdi(11); xff  = fdi(12)
   xttt = fdi(13); xftt = fdi(14); xfft = fdi(15)
   xfff = fdi(16)

   pe = g*pe
   pet = pet + 1.0d0
   pef = pef + 0.5d0*uf
   qe = qe/(re*wf)
   se = qe + 2.0d0*wf - psi
   sef = qe*(qef - ref - 0.5d0*uf) - 1.0d0/wf
   set = qe*(qet - ret)
   ue = se + psi - pe/(re*cte*T)

   ! Contributions from positrons
   rp = 0.0d0
   pp = 0.0d0
   sp = 0.0d0
   up = 0.0d0
   rpt = 0.0d0
   ppt = 0.0d0
   spt = 0.0d0
   rpf = 0.0d0
   ppf = 0.0d0
   spf = 0.0d0
   xip = -psi - 2.0 / (cte*T)
   if (eos_include_pairproduction .and. (xip > -15.0)) then
      ! Calculate effective degeneracy parameter for positrons
      ! See also Cox & Giuli ?24.9
      ! NB: PF and PG are later returned as dP/dlogf and Pgas
      paf = find_positron_f(xip)
      pf = exp(paf)
      puf = pf/(1.0d0 + pf)
      pwf = sqrt(1.0d0 + pf)
      pg = cte*T*pwf

      call positron_fdirac(pf, pg, fdi)
      rp   = fdi(1); pp   = fdi(2); qp   = fdi(3)
      rpt  = fdi(4); ppt  = fdi(5); qpt  = fdi(6)
      rpf  = fdi(7); ppf  = fdi(8); qpf  = fdi(9)

      ! Derivatives are now known with respect to the *positron* degeneracy
      ! parameter pf, but we need to know them with respect to the *electron*
      ! degeneracy parameter f. The conversion factor is easy, however,
      ! FF = dlog pf / dlog f = - sqrt( (1+f)/(1+pf) )
      ! (because dxi/dxip = -1)
      if (rp > 0.0d0) then
         ff = -wf/pwf

         rpf = rpf * ff
         ppf = ppf * ff
         qpf = qpf * ff

         ! Derivatives with respect to temperature get an extra term because
         ! the effective positron degeneracy parameter depends on the temperature
         ! as well as on the electron degeneracy parameter.
         rpt = rpt - 2.0d0 * rpf / ( cte*T*g )
         ppt = ppt - 2.0d0 * ppf / ( cte*T*g )
         qpt = qpt - 2.0d0 * qpf / ( cte*T*g )

         pp = pg*pp
         ppt = ppt + 1.0d0  + puf / pg
         ppf = ppf + 0.5d0 * puf * ff
         qp = qp/(rp*pwf)
         sp = qp + 2.0d0*pwf - xip
         spf = qp*(qpf - rpf - 0.5d0*puf) + 1.0d0/wf
         spt = qp*(qpt - rpt)
         up = sp + xip - pp/(rp*cte*T)
      end if
   end if

   ! Evaluate some quantities that do not depend on the state of ionization:
   ! the na are the element number densities (per baryon); avm is the average
   ! mass per baryon (in amu); neo and nio are the numbers of electrons and
   ! ions assuming complete ionization.
   na(1:9) = xa(1:9)/cbn(1:9)
   nio = sum(na(1:9))
   neo = dot_product(dkzn(1:9), na(1:9))
   nzz = dot_product(dkzn(1:9)**2, na(1:9))
   avm = dot_product(can(1:9), na(1:9))

   ! Pressi gives a crude model for pressure ionization and
   ! a model for coulomb interactions, returning corrections to the electron
   ! chemical potential, pressure, entropy and internal energy.
   ! ti is 1ev/kt, de is 1 amu * number of electrons/cm3
   ti = cevb/T
   de = re*cd
   call pressi ( 1, ti, psi, neo, nio, nzz, de, ref, ret, xtt, xft, xff, xftt, xfft, xfff, f, &
                  dc, dvt, dvf, dpa, dpat, dpaf, dsa, dsat, dsaf, dua )
   de = max(re-rp, 1.0d-50)*cd
   dv = dc - psi
   dvf = dvf - wf

   lna = log(na + epsilon)
   nef = 0.0d0
   net = 0.0d0
   si = 0.0d0
   sif = 0.0d0
   sit = 0.0d0
   ui = 0.0d0
   vm = 0.0d0
   mu = 0.0d0
   nif = 0.0d0
   nit = 0.0d0
   dsdnip = 0.0d0
   dsdnif = 0.0d0
   ! Contributions of the completely ionized species
   ne = dot_product(dkzn(kion+1:9), na(kion+1:9))
   si = -dot_product(na(kion+1:9), lna(kion+1:9) - lcan(kion+1:9))
   dsdnif(kion+1:9) = -(lna(kion+1:9) - lcan(kion+1:9))
   mu(kion+1:9) = -dsdnif(kion+1:9) + chi(1, kion+1:9)*ti
   zeta(kion+1:9) = dkzn(kion+1:9)

   ! We don't really care about the ionisation state of all other elements,
   ! but we need it for some calculations (gravitational settling,
   ! radiative levitation), so we need to calculate what the ionisation
   ! fractions would be anyway. We use a flag to optionally skip this part
   ! of the calculation.
   if (.false.) then
      do i = kion+1, 9
         va(1) = -chi(1,i)*ti + dv
         ha(1) = fxp(va(1))*com(kzn(i))/com(kzn(i)+1)
         sha = 1.0d0 + ha(1)
         do j = 2, kzn(i)
            va(j) = -chi(j,i)*ti + j*dv
            ha(j) = ha(j-1)*fxp(va(j) - va(j-1)) * com(kzn(i)+1-j)/com(kzn(i)+2-j)
            sha = sha + ha(j)
         end do
         vx = na(i)/sha          ! Number density of ground state per baryon
         do j = 1, kzn(i)
            xai(i,j) = ha(j)*vx  ! Number density of state J per baryon
         end do
      end do
   else
      forall (i = kion+1:9)
         xai(i, kzn(i)-1) = 0.0d0
         xai(i, kzn(i)) = na(i)
      end forall
   end if
   ! Calculate ionization of the first kion elements, in reverse order so
   ! hydrogen is last (the H calculation is continued after the loop ends)
   dvt_over_dvf = dvt / dvf
   do i = kion, 1, -1
      ! Compute potentials va and number ratios ha of ionization state j
      ! relative to the ground state 1
      zeta(i) = 0.0d0
      va(1) = -chi(1,i)*ti + dv
      ha(1) = fxp(va(1))*com(kzn(i))/com(kzn(i)+1)    ! Boltzmann factor
      sha = 1.0d0 + ha(1)                             ! Partition function
      sjha = ha(1)                                    ! Derivative of partition function wrt electron chemical potential
      scha = chi(1,i)*ti*ha(1)                        ! Derivative of partition function wrt temperature
      do j = 2, kzn(i)
         va(j) = -chi(j,i)*ti + j*dv
         ha(j) = ha(j-1)*fxp(va(j) - va(j-1)) * com(kzn(i)+1-j) / com(kzn(i)+2-j)
         sha = sha + ha(j)
         sjha = sjha + j*ha(j)
         scha = scha + chi(j,i)*ti*ha(j)
         zeta(i) = zeta(i) + j * ha(j)
      end do
      scha_over_sha = scha / sha
      sjha_over_sha = sjha / sha
      si = si + na(i)*lcan(i)
      dsdnif(i) = dsdnif(i) + lcan(i)
      mu(i) = mu(i) - lcan(i)
      zeta(i) = zeta(i) / ( sha + epsilon )
      if ( i > 1 ) then
         ! Contributions to electron number density ne, entropy si and
         ! internal energy ui for helium and heavier
         vx = na(i)/sha                               ! Occupation of ground state (neutral) = na / partition function
         si = si - vx*log(vx/com(kzn(i)+1) + epsilon)

         do j = 1, kzn(i)
            nx = ha(j)*vx        ! Number in state J per baryon
            xai(i,j) = nx        ! Store
            nxf = nx*dvf*(j - sjha_over_sha)
            nxt = nxf*dvt_over_dvf + nx*(chi(j,i)*ti - scha_over_sha)
            ne = ne + j*nx
            nef = nef + j*nxf
            net = net + j*nxt
            sif = sif - va(j)*nxf
            sit = sit - va(j)*nxt
            si = si - nx*log(max(nx/com(kzn(i) + 1 - j), epsilon))
            dsdnif(i) = dsdnif(i) - ha(j)/sha * log(max(nx/com(kzn(i) + 1 - j), epsilon))
            mu(i) = mu(i) + ha(j)/sha * ( log(ha(j)/sha + epsilon) )
            ui = ui + chi(j,i)*nx
            nif(i) = nxf
            nit(i) = nxt
         end do
      end if
   end do

   ! Set derivatives to 0 for species that are not present.
   ! Not entirely correct, but the easiest way to avoid numerical problems
   where (xa <= 0.0d0)
      dsdnif = 0.0d0
      dsdnip = 0.0d0
   end where

   ! Ionization and molecular dissciation of hydrogen.
   ! Partition function for H2 from Vardya (1960), Webbink (1975)
   vm = can(1)*sqrt(can(1))
   dh2ti = ch2(1)*ti
   expdh2ti = exp(-min(dh2ti, 100.0d0))
   d1ti = ch2(2)*ti
   d2ti = (ch2(3)*ti)**2
   d3ti = (ch2(4)*ti)**3
   zet = 1.0d0 - (1.0d0 + dh2ti)*expdh2ti
   dzh2t = -dh2ti**2*expdh2ti/zet
   dzh2tt = (dh2ti - 2.0d0 - dzh2t)*dzh2t
   zh2 = 6608.8d0*zet*dh2ti**(-2.5d0)*exp(-d1ti - d2ti - d3ti)
   zh2 = min(zh2, 1.0d100)       ! Avoid overflows later
   zh2t = 2.5d0 + d1ti + 2.0d0*d2ti + 3.0d0*d3ti + dzh2t
   zh2tt = -d1ti - 4.0d0*d2ti - 9.0d0*d3ti + dzh2tt
   zh2s = zh2*sqrt(8.0d0)/vm
   h2a = cen*(zh2s/4.0d0)*de/(T*sqrt(T))/expdh2ti
   h2bt = dh2ti + 1.5d0 - zh2t
   h2at = ret - h2bt

   ! Solve for densities of H+, H, and H2
   qa = 2.0d0*h2a + ha(1)*(1.0d0 + ha(1))
   qb = ne + ha(1)*(ne - na(1))
   qc = na(1)*ne
   hg = 2.0d0*qc/(sqrt(qb*qb + 4.0d0*qa*qc) + qb)
   hi = ha(1)*hg + tiny(0.0_dbl)    ! Make sure it's never 0
   xai(1,1) = hi                    ! Number density of protons (ionised H)
   ne = ne + hi                     ! Ionisation electrons / baryon
   np = rp / abs(re - rp) * ne     ! Positrons / baryon
   inv_ne = 1.0d0/ne
   h2 = h2a*hg*hg*inv_ne
   ni = nio - h2
   mu(1) = mu(1) + hi/(hi+hg) * ( log(hi/com(1) + epsilon) - chi(1, 1)*ti )
   mu(1) = mu(1) + hg/(hi+hg) * ( log(hg/com(2) + epsilon) )
   dsdnif(1) = dsdnif(1) - hi/(hi+hg) * log(hi/com(1) + epsilon)
   dsdnif(1) = dsdnif(1) - hg/(hi+hg) * log(hg/com(2) + epsilon)
   zeta(1) = hi / (hi + hg + h2 + epsilon)

   ! Derivatives w.r.T. f and T
   qa = ne + 4.0d0*hg*h2a
   qb = ha(1)*(ne - 2.0d0*h2)
   qd = 1.0d0/(qa + qb)
   qc = 2.0d0*h2*qd
   qaf = (nef - ne*ref )*qc
   qat = (net - ne*h2at)*qc
   qf = hg*qd
   qbf = dvf*qf
   qbt = (chi(1,1)*ti + dvt)*qf
   hgf = qaf - qb*qbf
   hgt = qat - qb*qbt
   hif = ha(1)*(qaf + qa*qbf)
   hit = ha(1)*(qat + qa*qbt)
   nef = nef + hif
   net = net + hit
   h2f = h2*ref  + inv_ne*(2.0d0*h2a*hg*hgf - h2*nef)
   h2t = h2*h2at + inv_ne*(2.0d0*h2a*hg*hgt - h2*net)

   ! Hydrogen contribution to entropy, internal energy
   sif = sif - va(1)*hif - h2bt*h2f
   sit = sit - va(1)*hit - h2bt*h2t + h2*(zh2t + zh2tt)

   ! Avoid overflow problems when HI, HG, H2 -> 0
   si = si - hi*log(abs(hi/com(1)) + epsilon)              &
           - hg*log(abs(hg/com(2)) + epsilon)              &
           - h2*(log(abs(h2/zh2s) + epsilon) - zh2t)
   ui = ui + chi(1,1)*hi + 0.5*ch2(1)*(hi + hg)

   ! db is 1 amu * number of baryons/cm3; rho is the mass density in g/cm3
   db = inv_ne*de
   dl = log(db)
   rho = db*avm
   arho = log(rho)
   rt = (re*ret + rp*rpt)/(re + rp) - inv_ne*net
   rf = (re*ref + rp*rpf)/(re + rp) - inv_ne*nef

   ! Second call to pressi compensates for spurious pressure and entropy terms in the limit of complete ionisation
   de = db*neo
   !wmu = 1.0d0/(neo+ni+h2)
   wmu = 1.0d0/(ne+nio)
   call pressi ( 0, ti, psi, neo, nio, nzz, de, rf, rt, xtt, xft, xff, xftt, xfft, xfff, f, &
                  dc, dvt, dvf, dpb, dpbt, dpbf, dsb, dsbt, dsbf, dub )

   ! Pressure terms:
   pe = cb*pe              ! Partial pressure due to electrons
   pp = cb*pp              ! Partial pressure due to positrons
   tcr = T*cr
   p0 = tcr*db
   pi = ni*p0
   t2 = T*T
   t4 = t2*t2
   pr = ca*t4/3.0d0
   b = 4.0d0*pr/p0
   ! When positrons are present, ignore pressure ionisation - everything is fully ionised anyway.
   if (pp > 0.0d0) then
      dpa = 0.0
      dpb = 0.0
      dpaf = 0.0
      dpbf = 0.0
      dpat = 0.0
      dpbt = 0.0
      dsaf = 0.0
      dsbf = 0.0
      dsat = 0.0
      dsbt = 0.0
      dsa = 0.0
      dsb = 0.0
   end if

   ! Finalise calculation of chemical potentials
   mu = mu + log(cen) + dl - 1.5d0*at

   ! Derivative of the density with respect to the composition
   vm     =  pi + (pe*pef + tcr*(dpaf-dpbf)) / ref
   dlrdni = (p0 + (pe*pef + tcr*(dpaf-dpbf)) / ref * zeta / ne) / vm

   ! Finalise calculation of entropy derivatives
   ! Common terms are collected in dS/dNi_f, then we add
   ! the terms proportional to dNe/dNi to dS/dNi_f and terms
   ! proportional to dlnrho/dNi_T,p and dNe/dNi to dS/dNi_p
   dsdnif = CR*(dsdnif + (1.5d0*at - dl + 1.5d0 - log(cen))) / avm * wmu
   dsdnip = dsdnif
   dsdnif = dsdnif + CR*((b + na + ne*se + ne*dsa - neo*dsb)*zeta/ne) / (cbn*avm)
   dsdnip = dsdnip + (CR*(b + na)/(ne*cbn) - (ne*sef + ne*dsaf-neo*dsbf)/ref) * dlrdni/avm &
                   + (CR*ne*se/avm + (ne*sef + ne*dsaf-neo*dsbf)/ref * (pi/vm - 1.0d0) ) * zeta/ne/avm

   ! Convert from derivative with respect to number fraction per unit mass to derivative of mass fraction
   dsdnif = dsdnif*cbn*wmu

   ! Pressure, in dyne/cm^2
   pg = pe + pp + pi + tcr*(dpa - dpb)
   p = max(pg + pr, 1.0d-200)
   pf = (pe*pef + pp*ppf      + pi*rf - h2f*p0 + tcr*(dpaf - dpbf))/p
   pt = (pe*pet + pp*ppt + pi + pi*rt - h2t*p0 + tcr*(dpat - dpbt)+pr*4.d0)/p
   ap = log(p)

   ! Entropy, in erg/g/K:
   dsf = nef*dsa + ne*dsaf - neo*dsbf - rf*b
   dst = net*dsa + ne*dsat - neo*dsbt - (rt - 3.0d0)*b
   sf = cr*(-ni*rf           + nef*se + ne*sef + np*spf + sif + dsf)/avm
   st = cr*( ni*(1.5d0 - rt) + net*se + ne*set + np*spt + sit + dst)/avm
   s = cr*(se*ne + sp*np + dsa*ne - dsb*neo + b + (1.5d0*at - dl + 2.5d0  - log(cen))*ni + si)/avm

   ! Internal energy, in erg/g:
   !          electr? posit?  elect?   electron  ion        pressure
   u1 = tcr*(ue*ne + up*np + dua*ne - dub*neo + 1.5d0*ni  + 0.75d0*b)/avm      ! Heat, pressure
   u2 = tcr*(ti*ui)/avm                                                        ! Recombination
   u3 = tcr*(zh2t*h2)/avm                                                      ! H2 association
   u  = u1 + u2 + u3                                                           ! Total internal energy (erg/g)
   uf = T*sf + p/rho * rf + CR*T*( dot_product(mu, nif) - dv*nef )
   ut = T*st + p/rho * rt + CR*T*( dot_product(mu, nit) - dv*net )

   ! Other thermodynamic quantities:
   q = min(pt*sf - pf*st, 1.0d300)
   scp = -q/pf             ! Specific heat at constant pressure
   scv = st - sf*rt/rf     ! Specific heat at constant volume
   grada = sf/q
   gamma = q/(rt*sf-rf*st)
   zt = sqrt(abs((ne*ref/wf + nzz)/ni))

   ! Coefficients of thermal expansion (K&W (6.6))
   delta = rf*pt/pf - rt
   phi   = -dot_product(xa*dlrdni, -(ne+nio)/(1.0d0+zeta))

   ! tst ought to be zero, if all the above programming is correct
   tst = sf/cr - p*(rt*pf-rf*pt)/(tcr*rho)

   ! Copy abundance results to output struct
   abund%xa = xa
   abund%na = na
   abund%neo = neo
   abund%nio = nio
   abund%nzz = nzz
   abund%avm = avm
   abund%ne = ne
   abund%xai = xai

   !*** End of thermodynamic calculation. beginning of table-based calculations
   frho = arho/cln
   tf = at/cln

   ! Cap values to allowed range for opacity tables
   frho = min(10.0d0, max(frho,-12.0d0))
   tf = min(9.3d0, max(tf,  3.0d0))
   fk = get_opacity(frho, tf, abund)
   xhi = 4.0d0*cl*pr/(fk*rho*rho*scp*T)   ! cm**2/s
   ! Plasma and photon viscosity
   ! The plasma viscosity is from Spitzer (1962)
   ! Photon viscosity from Thomas (1930) and Weinberg (1971) with the
   ! correction factor (10/9) of van den Horn & van Weert (1981).
   llambda = log( 2.0d0/(3.0d0*echar**3) * sqrt(AMU*wmu*(boltzm*T)**3/(cpi*rho*zt**5)) )
   mu_plasma = 0.406*sqrt( AMU*wmu*(boltzm*T)**5 ) / ( (zt*echar)**4*llambda )
   mu_rad = 8.0d0*pr/(3.0d0*cl*fk*rho)    ! gram/cm s
   nu = ( mu_plasma + mu_rad ) / rho
   ! Prandtl number
   prandtl = nu/xhi

   ! Derivatives of opacity: (dlogk/dlogT)_rho and (dlogk/dlogrho)_T
   !> \todo FIXME: these are the results for Kramer's opacity, should take
   !! derivative of spline instead to get value from tables.
   !<
   fkr =  1.0d0
   fkt = -3.5d0

   ! Neutrino-loss rates from Itoh et al (1983-1992) or Itoh et al (1996)
   call get_neutrino_rate(tf, frho, abund, en)

   ! Copy all EoS results to output struct
   eos%at = at
   eos%af = af
   eos%ap    = ap;    eos%arho    = arho; eos%u     = u;     eos%p    = p
   eos%u1    = u1;    eos%u2      = u2;   eos%u3    = u3;
   eos%rho   = rho;   eos%fk      = fk;   eos%T     = T;     eos%sf   = sf
   eos%st    = st;    eos%zt      = zt;   eos%grada = grada; eos%scp  = scp
   eos%rf    = rf;    eos%rt      = rt;   eos%xhi   = xhi;   eos%s    = s
   eos%pr    = pr;    eos%pg      = pg;   eos%pf    = pf;    eos%pt   = pt
   eos%en    = en;    eos%wmu     = wmu;  eos%gamma1 = gamma
   eos%delta = delta; eos%phi     = phi;  eos%fkt   = fkt
   eos%fkr   = fkr;   eos%nu      = nu;   eos%prandtl = prandtl
   eos%dv    = dv
   eos%alpha = rf/pf
   eos%scv   = scv
   return

contains

  !> \brief Fudged exponential function, used to avoid too large or too small numbers in the ionisation state of the Saha equation.
  !! - Uses precomputed values of the limiting exponentials.
   elemental function fxp(x)
      use real_kind

      implicit none
      real(double) :: fxp
      real(double), intent(in) :: x
      real(double), parameter :: fxp_low = -50.0d0
      real(double), parameter :: fxp_high = 50.0d0
      real(double), parameter :: low = 1.928749847963917820563444d-22
      real(double), parameter :: high = 5.184705528587072045056000d+21
      if (x>fxp_high) then
         fxp = high
         return
      else if (x>fxp_low) then
         fxp = exp(x)
         return
      end if
      fxp = low
   end function fxp
end subroutine statef

! ------------------------------------------------------------------------------
! FIDRAC
! Calculate the Fermi-Dirac integrals, as parametrised by functions f and g.
! ------------------------------------------------------------------------------
!  Returns:
!   FDI(9) - array of Fermi-Dirac integrals.
! ------------------------------------------------------------------------------
pure subroutine fdirac ( f, g, fdi )
   use real_kind

   implicit none
   real(double), intent(in) :: f,g
   real(double), intent(out) :: fdi(16)

   integer :: i,ij,ik,il,im,ind
   real(double) :: ff(4), gg(4), vw(4,4), vx(4,4)
   real(double) :: vf,vg,uf,ug,fdf,wv,ww

   real(double) :: d(3,3)
   real(double) :: xtt, xft, xff, xttt, xftt, xfft, xfff

   real(double), parameter :: c(4,4,3) = reshape(                     &
        (/ 2.315472d0,  7.128660d0,  7.504998d0,  2.665350d0,   &
           7.837752d0, 23.507934d0, 23.311317d0,  7.987465d0,   &
           9.215560d0, 26.834068d0, 25.082745d0,  8.020509d0,   &
           3.693280d0, 10.333176d0,  9.168960d0,  2.668248d0,   &
           2.315472d0,  6.748104d0,  6.564912d0,  2.132280d0,   &
           7.837752d0, 21.439740d0, 19.080088d0,  5.478100d0,   &
           9.215560d0, 23.551504d0, 19.015888d0,  4.679944d0,   &
           3.693280d0,  8.859868d0,  6.500712d0,  1.334124d0,   &
           1.157736d0,  3.770676d0,  4.015224d0,  1.402284d0,   &
           8.283420d0, 26.184486d0, 28.211372d0, 10.310306d0,   &
           14.755480d0, 45.031658d0, 46.909420d0, 16.633242d0,  &
           7.386560d0, 22.159680d0, 22.438048d0,  7.664928d0/), &
           (/ 4, 4, 3 /) )

   ! Evaluate Fermi-Dirac integrals (Eggleton, Faulkner & Flannery 1973).
   ! Matrix D contains rho*, P* and Q* and 1st logarithmic derivatives w.r.t.
   ! T, f. XTT etc are 2nd and 3rd log derivs of rho*
   vf = 1.0d0/(1.0d0 + f)
   vg = 1.0d0/(1.0d0 + g)
   uf = f*vf
   ug = g*vg
   fdf = g*g + g
   fdf = uf*fdf*sqrt(fdf)
   ff(1) = 1.0d0
   gg(1) = 1.0d0
   xtt = 0.0d0
   xft = 0.0d0
   xff = 0.0d0
   xttt = 0.0d0
   xftt = 0.0d0
   xfft = 0.0d0
   xfff = 0.0d0
   do i = 2, 4
      ff(i) = f*ff(i - 1)
      gg(i) = g*gg(i - 1)
      fdf = fdf*vf*vg
   end do
   ind = 4
   do i = 1, 3
      vx(1:ind, 1:ind) = 0.0d0
      do ij = 1, 4
         do ik = 1, 4
            vw(1, 1) = c(ik, ij, i)*gg(ij)*ff(ik)
            do il = 1, ind - 1
               vw(il + 1, 1) = (ij - 1)*vw(il, 1)
               do im = 1, ind - il
                  vw(il, im + 1) = (ik - 1)*vw(il, im)
               end do
            end do
            vx(1:ind, 1:ind) = vx(1:ind, 1:ind) + vw(1:ind, 1:ind)
         end do
      end do

      wv = 1.0d0/vx(1, 1)
      vx(1, 2:ind) = vx(1, 2:ind)*wv
      vx(2:ind, 1:ind-1) = vx(2:ind, 1:ind-1)*wv
      d(i, 1) = fdf*vx(1, 1)
      d(i, 2) = vx(2, 1) + 1.5d0 - 1.5d0*ug
      ww = 0.5d0*d(i, 2) - 4.0d0
      d(i, 3) = vx(1, 2) + 1.0d0 + ww*uf

      if ( i == 1 ) then
         ! Second- and third-order derivatives of density, needed in pressi:
         xtt = vx(3, 1) - vx(2, 1)**2 - 1.5d0*ug*vg
         xft = vx(2, 2) - vx(2, 1)*vx(1, 2) + 0.5d0*uf*xtt
         xff = vx(1, 3) - vx(1, 2)**2 + uf*(xft + ww*vf - 0.25d0*xtt*uf)
         xttt = vx(4, 1) + vx(2, 1)*(2.0d0*vx(2, 1)**2 - 3.0d0*vx(3, 1))   &
              - 1.5d0*(1.0-g)*ug*vg**2
         xftt = vx(3, 2) - 2.0d0*vx(2, 1)*vx(2, 2) - vx(1, 2)*(vx(3, 1)    &
              - 2.0d0*vx(2, 1)**2) + 0.5d0*uf*xttt
         xfft = vx(2, 3) - 2.0d0*vx(2, 2)*vx(1, 2) - vx(2, 1)*(vx(1, 3)    &
              - 2.0d0*vx(1,2)**2) + uf*(xftt + 0.5d0*vf*xtt - 0.25d0*uf*xttt)
         xfff = vx(1, 4) + vx(1, 2)*(2.0d0*vx(1, 2)**2 - 3.0d0*vx(1, 3))   &
              + uf*(1.5d0*(xfft + vf*xft) - uf*(0.75d0*(xftt + vf*xtt) -   &
              0.125d0*uf*xttt) + ww*(1.0d0 - f)*vf**2)
      end if
      ind = 2
   end do

   ! Copy output
   fdi = (/ d(1, 1), d(2, 1), d(3, 1), d(1, 2), d(2, 2), d(3, 2),          &
            d(1, 3), d(2, 3), d(3, 3), xtt, xft, xff, xttt, xftt, xfft, xfff /)
end subroutine fdirac

! --------------------------------------------------------------------------
! PRESSI
! --------------------------------------------------------------------------
! Non-ideal corrections to the equation-of-state: pressure ionisation and
! Coulomb interactions.
! As explained in Pols&al. (1995), the effect of pressure ionisation is
! calculated twice: once for the actual number of electrons Ne and once for
! the total number of electrons Ne0, to make sure the corrections to
! pressure and entropy vanish in the limit of complete ionisation.
!
! Input parameters:
!  IPR      - 1 or 0, if 1 also calculates Coloumb interactions
!  TI       - 1eV/kT, reciprocal temperature
!  PSI      - electron degeneracy parameter
!  NEO      - electron number fraction
!  NIO      - ion number fraction
!  NZZ      - average charge squared, ~electron density * average ion charge
!  XI       - electron number density in funny units, 1 amu * electrons/cm3
!  XF       - logarithmic derivative of electron density to f, dln rhoe/dln f
!  XT       - logarithmic derivative of electron density to T, dln rhoe/dln T
!  xtt      - higher order logarithmic derivative of electron density
!  xft      - higher order logarithmic derivative of electron density
!  xff      - higher order logarithmic derivative of electron density
! (xttt     - higher order logarithmic derivative of electron density)
!  xftt     - higher order logarithmic derivative of electron density
!  xfft     - higher order logarithmic derivative of electron density
!  xfff     - higher order logarithmic derivative of electron density
!  F        - Eggleton F parameter (alternative degeneracy parameter)
! Output parameters:
!  DC       - Correction for electron chemical potential (sign flipped)
!  DCT      - Derivative of above with respect to log T
!  DCF      - Derivative of above with respect to log f
!  DP       - Correction for the pressure
!  DPT      - Derivative of above with respect to log T
!  DPF      - Derivative of above with respect to log f
!  DS       - Correction for the entropy
!  DST      - Derivative of above with respect to log T
!  DSF      - Derivative of above with respect to log f
!  DU       - Correction for the internal energy
! The output parameters are calculated by taking derivatives of the free
! energy, which requires high order derivatives of the electron density.
! --------------------------------------------------------------------------
elemental subroutine pressi(ipr, ti, psi, neo, nio, nzz, xi, xf, xt, xtt, xft, xff, xftt, xfft, xfff, f, &
     dc, dct, dcf, dp, dpt, dpf, ds, dst, dsf, du)
   use real_kind
   use cuberoot
   use constants
   ! Computes effect (`pressure ionisation') on EoS of a free energy contribution
   ! delta(f) = -R.T.ne.g(x,y), where x = xi = Ne/v = ne, and y = yi = xih/(R.T).
   ! also does Coulomb interaction.
   ! delta(el. chem. pot.) = -dc = 1/(R.T) df/dNe
   ! delta(p) = R.T.dp = -dF/dv
   ! delta(s) = R.ne.ds = -dF/dT.
   ! Works firstly with independent variables x, y (called xi, yi), which are
   ! effectively Ne, v and T, but then has to get derivatives of dc, dp, ds
   ! with respect to independent variables f and T.

   implicit none
   integer, intent(in) :: ipr
   real(double), intent(in) :: ti, psi, neo, nio, nzz
   real(double), intent(in) :: xi, xf, xt, f, xtt, xft, xff, xftt, xfft, xfff
   real(double), intent(out) :: dc, dct, dcf, dp, dpt, dpf, ds, dst, dsf, du

   real(double) :: yi,ee,wx,cc,bb,gi,ff1,aa,th,ff2,af,thf,tht,rxf,thx,thy,psx,psy
   real(double) :: w2,psxx,psxy,psyy,eex,bx,by,bxx,bxy,byy,cx,cxx
   real(double) :: dgdx,dgdy,dgdxx,dgdxy,dgdyy,aff,thff,thft,thtt,tof
   real(double) :: wxx,wxy,thxx,thxy,thyy,thc,ww,ge,eeg,beg,begg,geg,dgdg
   real(double) :: gam,rbe,dgdgg,rthc,wy,gamx,gamy,gamxx,gamxy,gamyy,wt

   real(double), parameter :: ca1 = 0.89752d0
   real(double), parameter :: ca2 = 0.768d0
   real(double), parameter :: ca3 = 0.208d0
   real(double), parameter :: cp1 = 3.0d0
   real(double), parameter :: cp2 = 0.35d0
   real(double), parameter :: cp3 = 2.0d0
   real(double), parameter :: cp4 = 3.0d-2

   ! Pressure ionization:
   yi = 13.595d0*ti
   ee = 1.0d0/(1.0d0 + xi/cp4)
   wx = min((cp1/xi)**cp2, 300.0d0)
   cc = exp(-wx)
   bb = yi + psi - cp3*log(ee)

   ! extra free energy is -R.T.Ne.GI(X,Y)
   ! GI = exp{-WX(X)}{Y + psi(X,Y) + 2ln(1+X/const)}
   ! It OUGHT to have 3psi/5 for psi, but works worse!
   gi = cc*bb
   ff1 = 1.0d0/(f + 1.0d0)
   ! AA = dlnf/dpsi; THeta = dlnne/dpsi; also needed in Coulomb corrections
   ! f, T derivs, then X, Y.
   aa = sqrt(ff1)
   th = aa*xf
   ff2 = -0.5d0*f*ff1
   af = aa*ff2
   thf = af*xf + aa*xff
   tht = aa*xft
   rxf = 1.0d0/xf
   thx = thf*rxf
   thy = thx*xt - tht
   ! first and second derivatives dPSI/dlnX ... d2PSI/dlnY2
   psx = 1.0d0/th
   psy = xt*psx
   w2 = psx*psx
   psxx = -thx*w2
   psxy = -thy*w2
   psyy = psxy*xt + (xft*xt*rxf - xtt)*psx
   ! derivative -dlnEE/dlnX; -d2lnEE/dlnX2 = -EE*dlnEE/dlnX
   eex = cp3*(xi/cp4)*ee
   ! derivatives of BB
   bx = psx + eex
   by = yi + psy
   bxx = psxx + ee*eex
   bxy = psxy
   byy = yi + psyy
   ! derivatives of CC
   cx = cp2*wx*cc
   cxx = cp2*cx*(wx - 1.0d0)
   ! derivatives of GI
   dgdx = cc*bx + cx*bb
   dgdy = cc*by
   dgdxx = cc*bxx + 2.0d0*cx*bx + cxx*bb
   dgdxy = cc*bxy + cx*by
   dgdyy = cc*byy
   if ( ipr == 1 ) then
      ! Coulomb interaction.
      ! further derivatives of AA, THeta
      aff = 3.0d0*af*ff2 + af
      thff = aff*xf + 2.0d0*af*xff + aa*xfff
      thft = af*xft + aa*xfft
      thtt = aa*xftt
      ! d2TH/dlnX2, etc
      tof = xt*rxf
      wxx = thff - thx*xff
      wxy = thx*xft - thft
      thxx = wxx*rxf*rxf
      thxy = wxy*rxf + thxx*xt
      thyy = tof*(tof*wxx + 2.0d0*wxy) + thtt - thx*xtt
      ! gam is the plasma interaction parameter. note that thc = zt**2*nio/neo
      thc = th + abs(nzz/neo)
      ww = 1.0d0/ca2
      gam = cbrt(xi*(cpl*neo/nio)**2*c3rd)*ti/cevb*thc
      ! new BB and EE, and their derivates
      bb = (ca1*sqrt(3.0d0/gam))**ww
      ee = (gam/(gam + ca3))**ww
      rbe = 1.0d0/(ee + bb)
      ! further addition GE to free en; adds to previous GI.
      ! GE = GE(GAMma), GAM = const. * X**0.33 * Y * (THeta(X,Y) + const)
      ge = (nio/neo)*ca1*gam*rbe**ca2
      eeg = ca3/(gam + ca3)
      beg = (eeg*ee - 0.5d0*bb)*rbe
      begg = (eeg*eeg*(1.0d0 - gam*ca2/ca3)*ee + 0.25d0*bb)*rbe
      geg = 1.0d0 - beg
      dgdg = ge*geg
      dgdgg = ge*(geg*geg + (beg*beg - begg)*ww)
      rthc = 1.0d0/thc
      wx = thx*rthc
      wy = thy*rthc
      ! dlnGAM/dlnX, etc
      gamx = c3rd + wx
      gamy = 1.0d0 + wy
      gamxx = thxx*rthc - wx*wx
      gamxy = thxy*rthc - wx*wy
      gamyy = thyy*rthc - wy*wy
      gi = gi + ge
      ! derivatives w.r.t. X, Y; in effect w.r.t. Ne, V, and T, since X = Ne/V
      dgdx = dgdx + dgdg*gamx
      dgdy = dgdy + dgdg*gamy
      dgdxx = dgdxx + dgdgg*gamx**2 + dgdg*gamxx
      dgdxy = dgdxy + dgdgg*gamx*gamy + dgdg*gamxy
      dgdyy = dgdyy + dgdgg*gamy**2 + dgdg*gamyy
   end if
   ! Evaluate changes to el. chem. potential (-dc), pressure (R.T.dp), and
   ! entropy (R.ne.ds), and their derivatives w.r.t. log(f) and log(T)
   dc = dgdx + gi
   ww = dgdxx + dgdx
   wt = ww*xt
   dcf = ww*xf
   dct = wt - dgdxy - dgdy
   dp = -xi*dgdx
   dpf = -xi*dcf
   dpt = xi*(dgdxy - wt) + dp
   ds = gi - dgdy
   dst = dgdx - dgdxy
   dsf = dst*xf
   dst = dst*xt - dgdy + dgdyy
   du = -dgdy
end subroutine pressi

! Compute rates of (at present) 20 nuclear rections, and the corresponding
! energy and neutrino release
elemental subroutine nucrat ( tl, abund, eos )
   use real_kind
   use cuberoot
   use settings, only: cnuc
   use constants
   use reaction_rate_data
   use eostate_types

   implicit none
   real(double), intent(in) :: tl
   type(abundance), intent(in) :: abund
   type(eostate), intent(inout) :: eos
   integer :: it
   real(double) :: rhb,wc,wb,wa,xb,vl,za,zb,zc,zd,tf,tt,tu
   real(double) :: scrn(21),strn(21),dstr(21)
   real(double) :: fpng,rpn,fccg,rcc,f34,pp2,pp3,qpp,qnpp

   real(double) :: rtt(22)
   real(double), parameter :: csa = 0.624
   real(double), parameter :: csb = 0.316
   real(double), parameter :: csc = 0.460
   real(double), parameter :: csd = 0.38
   real(double), parameter :: cxd = 0.86

   real(double) :: rrt(22)
   real(double) :: rpp, r33, r34,  rbe, rbp, rpc, rpna,  rpo,  r3a,  rac
   real(double) :: ran, rao, rane, rcca, rco, roo, rgne, rgmg, rccg, rpng
   real(double) :: ramg

   real(double) :: n1, n4, n12, n14, n16, n20, n24, n28, n56
   real(double) :: ne, ni, avm

   ! Not always defined:
   dstr = 0.0_dbl
   strn = 0.0_dbl

   ni = abund%nio
   ne = abund%neo
   n1 = abund%na(1)
   n4 = abund%na(2)
   n12 = abund%na(3)
   n14 = abund%na(4)
   n16 = abund%na(5)
   n20 = abund%na(6)
   n24 = abund%na(7)
   n28 = abund%na(8)
   n56 = abund%na(9)
   avm = abund%avm

   ! rhb is 'baryon density': 1 amu * number of baryons per cm3
   rhb = eos%rho/avm
   ! Electron-screening theory from Graboske, Dewitt, Grossman & Cooper (1973),
   ! for strong (za, zb, zc) are intermediate screening (zd). The reaction
   ! dependent charge parameters are stored in cza ... czd.
   wc = dot_product((/n1, n4, n12, n14, n16, n20, n24, n28, n56/), vz(1:9))
   wc = wc/ni
   wb = ne/ni
   wa = eos%zt*eos%zt/(wb*wb)
   xb = cbrt(abs(wb))
   vl = cpl*sqrt(abs(ni)*rhb*exp(-3.0d0*tl))
   za = csa*xb*vl**(2.0/3.0)
   zb = csb*xb*za
   zc = csc/(xb*xb)
   zd = csd*wc*wa*(vl/(wa*eos%zt))**cxd

   ! Reaction rates interpolated in T, mostly from Caughlan & Fowler (1988):
   tf = tl/cln
   rrt(:) = 0.0d0
   rtt(:) = 0.0d0
   tt = 50.0d0*(tf - 6.0d0) + 1.0d0
   if ( tt >= 1.0d0 ) then
      it = max0(1, min0(199, int(tt)))
      tt = tt - it
      tu = 1.0d0 - tt
      rrt(2:20) = tu*crt(it, 1:19) + tt*crt(it + 1, 1:19)
      rrt(21)   = tu*nucsyn_crt(it, 55) + tt*nucsyn_crt(it + 1, 55)
      scrn(1:19) = zd*czd(1:19)
      strn(1:19) = za*cza(1:19) + zb*czb(1:19)
      dstr(1:19) = zc*czc(1:19)
      scrn(20) = zd*czd(55)
      strn(20) = za*cza(55) + zb*czb(55)
      dstr(20) = zc*czc(55)
      where (dstr < 0.29d0*strn) scrn = min(scrn, strn - dstr)
      rrt(2:21) = exp(min(cln*rrt(2:21) + scrn(1:20), 300.d0))

      ! Derivative (dR/dlogT)_rho, by differentiation of the interpolation function:
      rtt(2:21) = 50.0/cln * (crt(it + 1, 1:20) - crt(it, 1:20))
   end if

   ! modulate with cnuc
   rrt = rrt * cnuc

   ! Copy rates to sensible variable names
   !> \todo FIXME: we really shouldn't need this, do we?
   rpp = rrt(2)
   r33 = rrt(3)
   r34 = rrt(4)
   rBe = rrt(5)
   rBp = rrt(6)
   rpc = rrt(7)
   rpN = rrt(8)
   rpO = rrt(9)
   r3a = rrt(10)
   raC = rrt(11)
   raN = rrt(12)
   raO = rrt(13)
   raNe = rrt(14)
   rCC = rrt(15)
   rCO = rrt(16)
   rOO = rrt(17)
   rgNe = rrt(18)
   rgMg = rrt(19)
   rCCg = rrt(20)
   ramg = rrt(21)

   ! Multiply with density and abundances to get rates per baryon per second,
   ! note that abundances of He3 and Be7 are not needed in equilibrium
   rpp = rhb*n1*n1*rpp/2.0d0
   r33 = rhb*r33/2.0d0
   r34 = rhb*n4*r34
   rbe = rhb*ne*rbe
   rbp = rhb*n1*rbp
   rpc = rhb*n1*n12*rpc
   rpn = rhb*n1*n14*rpn
   rpo = rhb*n1*n16*rpo
   r3a = rhb*rhb*n4*n4*n4*r3a/6.0d0
   rac = rhb*n4*n12*rac
   ran = rhb*n4*n14*ran
   rao = rhb*n4*n16*rao
   rane = rhb*n4*n20*rane
   rcc = rhb*n12*n12*rcc/2.0d0
   rco = rhb*n12*n16*rco
   roo = rhb*n16*n16*roo/2.0d0
   rgne = n20*rgne
   rgmg = n24*rgmg
   ramg = rhb*n4*n24*ramg
   ! Branching of pn and cc reactions
   fpng = 8.0d-4
   rpna = (1.0d0 - fpng)*rpn
   rpng = fpng*rpn
   rpn = rpna
   fccg = rccg
   rcca = (1.0d0 - fccg)*rcc
   rccg = fccg*rcc
   rcc = rcca
   ! pp chain in equilibrium, rpp becomes effective rate of 2 h1 -> 0.5 he4
   f34 = 0.0d0
   rpp = min(rpp, 1.0d303)
   r33 = min(r33, 1.0d303)
   r34 = min(r34, 1.0d303)
   if (r34 > 1.0d-20)  f34 = 2.0d0/(1.0d0 + sqrt(1.0d0 + 8.0d0*rpp*r33/(r34*r34)))
   rpp = rpp*(1.0d0 + f34)
   pp2 = 1.0d0
   if (rbe + rbp > 1.0d-20) pp2 = rbe/(rbe + rbp)
   pp3 = 1.0d0 - pp2
   qpp = qrt(1) + 0.5d0*qrt(2)
   qnpp = (qnt(1) + f34*(qnt(4)*pp2 + qnt(5)*pp3))/(1.0d0 + f34)

   ! Copy rates back to rrt array for loop below
   rrt(2:21) = (/ rpp, r33, r34,  rBe, rBp, rpc, rpN,  rpO,  r3a,  raC,&
                  raN, raO, raNe, rCC, rCO, rOO, rgNe, rgMg, rCCg, rpng /)

   ! Calculate energy release and neutrino loss, in erg/gram/sec
   eos%ex = qpp*rpp
   eos%enx = qnpp*rpp
   eos%ext = qpp*rpp*rtt(2)
   eos%ex = eos%ex + dot_product(qrt(6:20), rrt(7:21)) + 9.98414d0*ramg
   eos%enx = eos%enx + dot_product(qnt(6:20), rrt(7:21))
   eos%ext = eos%ext + dot_product(qrt(6:20), rrt(7:21)**2) + 9.98414d0*ramg**2
   eos%ext = eos%ext / max(eos%ex, 1.0d-50)
   ! Guard against overflows
   eos%ex = min(cme*eos%ex/avm, 1.0d303)
   eos%enx = -min(cme*eos%enx/avm, 1.0d303)

   ! Copy reaction rates to output struct
   eos%rpp  = rpp;  eos%r33  = r33;
   eos%r34  = r34;  eos%rbe  = rbe;
   eos%rbp  = rbp;  eos%rpc  = rpc;
   eos%rpna = rpna; eos%rpo  = rpo;
   eos%r3a  = r3a;  eos%rac  = rac;
   eos%ran  = ran;  eos%rao  = rao;
   eos%rane = rane; eos%rcca = rcca;
   eos%rco  = rco;  eos%roo  = roo;
   eos%rgne = rgne; eos%rgmg = rgmg;
   eos%rccg = rccg; eos%rpng = rpng;
   eos%ramg = ramg
end subroutine nucrat

!> \brief  Helper function: calculate degeneracy parameter xi given lnf
elemental function calculate_degeneracy_parameter(lnf)
   use real_kind

   implicit none
   real(double) :: calculate_degeneracy_parameter
   real(double), intent(in) :: lnf
   real(double) :: uf, wf, xi

   !   if (auf < -100.0) then
   !      xi = auf
   !   else if (auf > 100.0) then
   !      uf = exp(auf)
   !      xi = 2.0d0*sqrt(uf)
   !   else
   uf = exp(lnf)
   wf = sqrt(1.0d0 + uf)
   xi = log(uf) + 2.0d0*(wf - log(1.0d0 + wf))
   !   end if

   calculate_degeneracy_parameter = xi
end function calculate_degeneracy_parameter

!  Find the parameter "log f" from the positron degeneracy parameter XI
elemental function find_positron_f(xi)
   use real_kind

   implicit none
   real(double), intent(in) :: xi
   real(double) :: find_positron_f

   !  Asymptotic expansions
   if (xi > 15.0) then
      find_positron_f = log(0.25*xi**2 + 1)
      return
   else if (xi < -5) then
      find_positron_f = xi - 2.0 + log(4.0)
      return
   end if
   !  If logf is below -100.0 we really don't care enough to determine it
   !  more accurately; all Fermi functions go to 0 as exp(xi) in this region
   !  anyway.
   find_positron_f = solve_f_interval(xi,   -5.d0, 15.d0)
end function find_positron_f

!  Calculate Fermi-Dirac integrals for positrons, passing in the appropriate
!  values of F and G for Positrons
pure subroutine positron_fdirac ( f, g, dp )
   use real_kind

   implicit none
   real(double), intent (in) :: f, g
   real(double), intent (out) :: dp(16)

   !  If the degeneracy parameter for positrons is very small, we don't care
   !  because all contributions are going to be essentially 0 anyway.
   if (f < 1.0d-300) then
      dp(:) = 0.0d0
      return
   end if

   !  Compute positron FD functions
   call fdirac(f, g, dp)
end subroutine positron_fdirac

!  Solve function: find the point x in the interval (x1, x2) where the
!  function value equals y, using bisection.
!  Returns x, or x1 - 1 if the function did not converge.
!  Code based on the Numerical Recipes routine RTSEC
pure function solve_f_interval(y,x1,x2)
   use real_kind

   implicit none
   real(double), parameter :: epsilon = 1.0d-8
   integer, parameter :: maxiter = 10
   real(double), intent(in) :: y, x1, x2
   real(double) :: solve_f_interval
   real(double) :: dx, fh, fl, xl, xh
   integer :: i

   fl = calculate_degeneracy_parameter(x1) - y
   fh = calculate_degeneracy_parameter(x2) - y
   if (abs(fl) < abs(fh)) then
      !     Swap highest/lowest variables
      xh = x1
      xl = x2
      !     Swap function values
      dx = fl
      fl = fh
      fh = dx
   else
      xl = x1
      xh = x2
   end if
   !  Loop until converged to the required accuracy
   do i = 1, maxiter
      !     Calculate correction
      dx = (xl - xh) * fh / (fh - fl)
      !     On each iteration, flip boundaries
      xl = xh
      fl = fh
      xh = xh + dx
      fh = calculate_degeneracy_parameter(xh) - y
      !     Converged?
      if (abs(dx) < epsilon .or. fh == 0.0d0) exit
   end do
   solve_f_interval = xh
end function solve_f_interval

end module equation_of_state

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
    module ArrayUtils
    use MiscUtils
    use MpiUtils
    implicit none

    INTERFACE Reallocate
    module procedure Realloc_R, Realloc_D, Realloc_I
    END INTERFACE Reallocate

    contains

    subroutine realloc_R(arr, new_size, keep)
    real, allocatable, intent(inout) :: arr(:)
    integer, intent(in) :: new_size
    logical, intent(in), optional :: keep
    integer sz
    real, allocatable :: tmp(:)

    if (.not. allocated(arr)) then
        allocate(arr(new_size))
        return
    end if
    if (LBOUND(arr,1)/=1) call MpiStop('Realloc only works on arrays starting at 1')
    if (new_size/=size(arr)) then
        allocate(tmp(new_size))
        if (DefaultTrue(keep)) then
            sz = min(new_size, size(arr))
            tmp(:sz) = arr(:sz)
        end if
        call move_alloc(tmp, arr)
    end if

    end subroutine realloc_R

    subroutine realloc_D(arr, new_size, keep)
    double precision, allocatable, intent(inout) :: arr(:)
    integer, intent(in) :: new_size
    logical, intent(in), optional :: keep
    integer sz
    double precision, allocatable :: tmp(:)

    if (.not. allocated(arr)) then
        allocate(arr(new_size))
        return
    end if
    if (LBOUND(arr,1)/=1) call MpiStop('Realloc only works on arrays starting at 1')
    if (new_size/=size(arr)) then
        allocate(tmp(new_size))
        if (DefaultTrue(keep)) then
            sz = min(new_size, size(arr))
            tmp(:sz) = arr(:sz)
        end if
        call move_alloc(tmp, arr)
    end if

    end subroutine realloc_D

    subroutine realloc_I(arr, new_size, keep)
    integer, allocatable, intent(inout) :: arr(:)
    integer, intent(in) :: new_size
    logical, intent(in), optional :: keep
    integer sz
    integer, allocatable :: tmp(:)

    if (.not. allocated(arr)) then
        allocate(arr(new_size))
        return
    end if
    if (LBOUND(arr,1)/=1) call MpiStop('Realloc only works on arrays starting at 1')
    if (new_size/=size(arr)) then
        allocate(tmp(new_size))
        if (DefaultTrue(keep)) then
            sz = min(new_size, size(arr))
            tmp(:sz) = arr(:sz)
        end if
        call move_alloc(tmp, arr)
    end if

    end subroutine realloc_I

    function IndexOf(aval,arr, n)
    integer, intent(in) :: n, arr(n), aval
    integer IndexOf, i

    do i=1,n
        if (arr(i)==aval) then
            IndexOf= i
            return
        end if
    end do
    IndexOf = 0

    end function IndexOf

    function MaxIndex(arr, n)
    integer, intent(in) :: n
    real, intent(in) :: arr(n)
    integer locs(1:1), MaxIndex

    locs = maxloc(arr(1:n))
    MaxIndex = locs(1)

    end function MaxIndex

    function MinIndex(arr, n)
    integer, intent(in) :: n
    real, intent(in) :: arr(n)
    integer locs(1:1), MinIndex

    locs = minloc(arr(1:n))
    MinIndex = locs(1)

    end function MinIndex

    end module ArrayUtils

max_stars_repo_name: LBJ-Wade/forutils

max_issues_repo_name: LBJ-Wade/forutils

max_forks_repo_name: LBJ-Wade/forutils


text:
   10 program main
         integer   :: a,b
         logical   :: c,d
         real      :: e,f

*        This is only valid code for f77, not f90 (or later)
         real*8    :: g,h

         character ::i,j

*        This is only valid code for f77, not f90 (or later)
         complex*8 :: k,l

         complex   :: m,n
         double precision :: o
   20 end program main

ext: f

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
	SUBROUTINE AF_PLYR  ( cdtp, islyr, ielyr, nlyr, iret )
C************************************************************************
C* AF_PLYR								*
C*									*
C* This subroutine separates the turbulence, sky cover, or icing data	*
C* from a PIREP report into layers that can be individually decoded.	*
C* A layer of data is separated from the other layers by a "like-type"	*
C* group containing a "/" character.  The data must have already been	*
C* broken up into "like-type" groups using subroutine AF_BKGP.		*
C* Output consists of arrays containing the indicees of the "like-type"	*
C* groups which contain the start and the end of the data for each of	*
C* the layers.								*
C*									*
C* AF_PLYR  ( CDTP, ISLYR, IELYR, NLYR, IRET )				*
C*									*
C* Input parameters:							*
C*	CDTP		CHAR*		Tag describing type of data	*
C*									*
C* Output parameters:							*
C*	ISLYR (NLYR)	INTEGER		Indicees of "like-type" groups	*
C*					marking start of each layer	*
C*	IELYR (NLYR)	INTEGER		Indicees of "like-type" groups	*
C*					marking end of each layer	*
C*	NLYR		INTEGER		Number of layers 		*
C*	IRET		INTEGER		Return code 			*
C*					  0 = normal return 		*
C*									*
C**									*
C* Log:									*
C* J. Ator/NP12		09/96						*
C* J. Ator/NP12		08/97	New interface format, style changes	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'afcmn.cmn'
C*
	INTEGER		islyr (*), ielyr (*)
	CHARACTER*(*)	cdtp
C*
	LOGICAL		slash
C-----------------------------------------------------------------------
	iret = 0
	nlyr = 1
	islyr ( nlyr ) = 1
C
	IF  ( nflds .gt. 2 )  THEN
	    DO ii = 2, ( nflds - 1 )
		IF  ( itypsf (ii) .eq. NALNMR )  THEN
C
C*		    Does this non-alphanumeric "like-type"  group
C*		    contain a "/" character?
C
		    slash = .false.
		    DO jj = 1, lensf (ii)
			IF  ( fields (ii) (jj:jj)  .eq. '/' )  THEN
			  slash = .true.
			END IF
		    END DO
		    IF  ( slash )  THEN
C
C*			This "like-type" group contains a "/" character,
C*			so it separates two layers of the data.
C
			ielyr ( nlyr ) = ii - 1
			IF  ( nlyr .lt. MXNLYR )  THEN
			  nlyr = nlyr + 1
			  islyr ( nlyr ) = ii + 1
			ELSE
			  WRITE ( UNIT = logmsg, FMT = '( I2, 3A )' )
     +			    MXNLYR, ' ', cdtp, ' layers'
			  CALL DC_WLOG  ( 2, 'AF', 3, logmsg, ierwlg )
			  RETURN
			END IF
		    END IF
		END IF
	    END DO
	END IF
C
	ielyr ( nlyr ) = nflds
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


!***********************************/
!	Name:Sutherland's lawに基づき粘度を求めるプログラム
!	Alias:UGetLaminarViscosity
!	Description:
!	Type:
!	Input:
!	Output:
!	Note:
!	Author:Akitaka Toyota
!	Date:2018.11.13
!	Update:-
!	Other:
!***********************************/
subroutine UGetLaminarViscosity(UConf, UG,UCC)
    use StructVar_Mod
    use LoopVar_Mod
    use FrequentOperation
    use ConstantVar_Mod, gamma => SpecificOfHeatRatio
    implicit none
    type(Configulation), intent(in) :: UConf
    type(UnstructuredGrid), intent(in) :: UG
    type(CellCenter), intent(inout) :: UCC
    double precision :: Mu0, STS, SC1, Mach2!, Tref

    if(UConf%UseSutherlandLaw == 0) then
        UCC%LaminarViscosity = 1.0d0

    else
        Mu0 = ReferenceViscosity_Mu0
        !Tref = ReferenceTemperature_Tref
        STS = SutherlandTemperature_S
        SC1 = SutherlandCoefficient1
        Mach2 = MachNumber ** 2

        do iCell = 1, UG%GI%RealCells
            UCC%Temparature(iCell, 1, 1) = gamma * UCC%PrimitiveVariable(5, iCell, 1, 1) / UCC%PrimitiveVariable(1, iCell, 1, 1) * Mach2
            UCC%LaminarViscosity(iCell, 1, 1) = SC1 * UCC%Temparature(iCell, 1, 1) ** 1.5d0 / (UCC%Temparature(iCell, 1, 1) + STS)
        end do

    end if

    return
end subroutine UGetLaminarViscosity


text:
!/*-----------------------------------------------------------------------------
! *  get_acc_jerk_pot_coll  --  calculates accelerations and jerks, and as side
! *                             effects also calculates potential energy and
! *                             the time scale coll_time for significant changes
! *                             in local configurations to occur.
! *                                                  __                     __
! *                                                 |          -->  -->       |
! *               M                           M     |           r  . v        |
! *   -->          j    -->       -->          j    | -->        ji   ji -->  |
! *    a   ==  --------  r    ;    j   ==  -------- |  v   - 3 ---------  r   |
! *     ji     |-->  |3   ji        ji     |-->  |3 |   ji      |-->  |2   ji |
! *            | r   |                     | r   |  |           | r   |       |
! *            |  ji |                     |  ji |  |__         |  ji |     __|
! *
! *  note: it would be cleaner to calculate potential energy and collision time
! *        in a separate function.  However, the current function is by far the
! *        most time consuming part of the whole program, with a double loop
! *        over all particles that is executed every time step.  Splitting off
! *        some of the work to another function would significantly increase
! *        the total computer time (by an amount close to a factor two).
! *
! *  We determine the values of all four quantities of interest by walking
! *  through the system in a double {i,j} loop.  The first three, acceleration,
! *  jerk, and potential energy, are calculated by adding successive terms;
! *  the last, the estimate for the collision time, is found by determining the
! *  minimum value over all particle pairs and over the two choices of collision
! *  time, position/velocity and sqrt(position/acceleration), where position and
! *  velocity indicate their relative values between the two particles, while
! *  acceleration indicates their pairwise acceleration.  At the start, the
! *  first three quantities are set to zero, to prepare for accumulation, while
! *  the last one is set to a very large number, to prepare for minimization.
! *       The integration loops only over half of the pairs, with j > i, since
! *  the contributions to the acceleration and jerk of particle j on particle i
! *  is the same as those of particle i on particle j, apart from a minus sign
! *  and a different mass factor.
! *-----------------------------------------------------------------------------
! */

subroutine get_acc_jerk_pot_coll(mass, pos, &
                                 vel,  acc, &
                                 jerk, n, epot, &
                                 coll_time)
    use numz
    use dim
    implicit none
    real(b8) mass(0:), pos(0:,0:)
    real(b8) vel(0:,0:),  acc(0:,0:)
    real(b8) jerk(0:,0:)
    integer  n
    real(b8) epot
    real(b8) coll_time
    real(b8) rji(0:ndim-1)
    real(b8) vji(0:ndim-1)
    real r2,v2,rv_r2,r,r3,da(0:2),dj(0:2),da2
    real(b8) coll_time_q ! collision time to 4th power
    real(b8) coll_est_q ! collision time scale estimate to 4th power (quartic)
    real(b8), parameter :: vln =1.0e300_b8
    real(b16) mij
    integer i,k,j

    acc=0.0_b8
    jerk=0.0_b8
    epot=0.0_b8
    coll_time_q = vln      ! collision time to 4th power
    do i=0,n-1
      do j=i+1, n-1                              !! rji() is the vector from
!            real rji(NDIM)                        !! particle i to particle j
!            real vji(NDIM)                        !! vji() = d rji() / d t
            do k=0, ndim-1
                rji(k) = pos(j,k) - pos(i,k)
                vji(k) = vel(j,k) - vel(i,k)
            enddo
            r2 = 0.0_b8                           !! | rji |^2
            v2 = 0.0_b8                           !! | vji |^2
            rv_r2 = 0.0_b8                        !! ( rij . vij ) / | rji |^2
            do k=0, ndim-1
                r2=r2+ rji(k) * rji(k)
                v2=v2+ vji(k) * vji(k)
                rv_r2= rv_r2+ rji(k) * vji(k)
            enddo
            rv_r2 = rv_r2 / r2
            r = sqrt(r2)                     !! | rji |
            r3 = r * r2                      !! | rji |^3

!! add the {i,j} contribution to the total potential energy for the system:

            epot = epot - mass(i) * mass(j) / r

!! add the {j (i)} contribution to the {i (j)} values of acceleration and jerk:

!           real da(3)                            !! main terms in pairwise
!           real dj(3)                            !! acceleration and jerk
            do k=0, ndim-1
                da(k) = rji(k) / r3                           !! see equations
                dj(k) = (vji(k) - 3 * rv_r2 * rji(k)) / r3    !! in the header
            enddo
            do k=0, ndim-1
                acc(i,k) = acc(i,k) + mass(j) * da(k)                 !! using symmetry
                acc(j,k) = acc(j,k) - mass(i) * da(k)                 !! find pairwise
                jerk(i,k) = jerk(i,k) + mass(j) * dj(k)               !! acceleration
                jerk(j,k) = jerk(j,k) - mass(i) * dj(k)               !! and jerk
            enddo

!! first collision time estimate, based on unaccelerated linear motion:

            coll_est_q = (r2*r2) / (v2*v2)
            if (coll_time_q .gt. coll_est_q) coll_time_q = coll_est_q

!! second collision time estimate, based on free fall:

            da2 = 0.0_b8                                  !! da2 becomes the
            do k=0, ndim-1                              !! square of the
                da2= da2+ da(k) * da(k)
            enddo                                         !! pair-wise accel-
            mij = mass(i) + mass(j)                       !! eration between
            da2=da2* mij * mij                            !! particles i and j

            coll_est_q = r2/da2
            if (coll_time_q .gt. coll_est_q) coll_time_q = coll_est_q
        enddo
    enddo                                              !! from q for quartic back
    coll_time = sqrt(sqrt(coll_time_q))                !! to linear collision time
end subroutine get_acc_jerk_pot_coll

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
      BLOCK DATA
C
C     BLOCK DATA 7 - APRIL 1981
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON/ACFILT/ACPARM(20),IACFLT(20)
C
      COMMON/CONDMP/IDPHLD
C
      COMMON/CRASTR/VRASTR(36)
C
      COMMON/CSECBD/VSB(46)
C
      COMMON/DPCONS/DPPARM(200),IDPCS(50)
C
      COMMON/DPCSIC/DPIC(50),DPUP(50),DPDN(50)
C
      COMMON/GPCONS/GPPARM(200),IGPCS(50)
C
      COMMON/GPINCS/GPICS(50),GPUP(50),GPDN(50)
C
      COMMON/ISECBD/IVSB(10)
C
C
      COMMON/SECICS/VSBIC(18),SBUP(2),SBDN(2)
C
      COMMON/SBINTF/VSBK(12)
C
      DATA ACPARM/20*0.0D0/
      DATA IACFLT/20*0/
C
      DATA IDPHLD/0/
C
      DATA VRASTR/36*0.0D0/
C
      DATA VSB/46*0.0D0/
C
      DATA DPPARM/200*0.0D0/
      DATA IDPCS/50*0/
C
      DATA DPIC/50*0.0D0/
      DATA DPUP/50*1.0D-3/
      DATA DPDN/50*1.0D-5/
C
      DATA GPPARM/200*0.0D0/
      DATA IGPCS/50*0/
C
      DATA GPICS/50*0.0D0/
      DATA GPUP/50*1.0D-3/
      DATA GPDN/50*1.0D-5/
C
      DATA IVSB/10*0/
C
C
      DATA VSBIC/18*0.0D0/
      DATA SBUP/2*1.0D-3/
      DATA SBDN/2*1.0D-5/
C
      DATA VSBK/12*0.0D0/
C
      END

ext: for