C
C                     FISOLV User Subroutine Templates
C
#ifdef SOLVER_MODULE
#undef SOLVER_MODULE
#endif
#ifndef PREPOST_MODULE
#define PREPOST_MODULE
#endif
#include "FIDAP_DATA_TYPES.h"
C
C     When using a user subroutine (USRxxx) the user inserts
C     their commands after the dimension statements and either
C     deletes or comments out the call to ERMSGS.
C
      SUBROUTINE USRACP (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR,UTNG)
C
C     USER DEFINED A COEFFICIENT FOR PERMEABILITY
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = A COEFFICIENT FOR PERMEABILITY EQUATION
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C     UTNG  = NORM OF VELOCITY (SPEED)
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS),UTNG(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (8301, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRBCF (VAL,NODE,IDF,TIME,SOL,NUMEQA,NDOF,NUMNP,LDOFU,
     1                   CONSTR,nodepr,xyz,iflag)
C
C     USER DEFINED FLUX COEFFICIENTS FOR APPLIED FLUX BOUNDARY
C     CONDITIONS
C
C     VAL    = COMPUTED (SPECIFIED) FLUX COEFFICIENT
C     SOL    = GLOBAL SOLUTION VECTOR
C     NUMEQA = GLOBAL EQUATION NUMBER ARRAY
C     NDOF   = ACTIVE NUMBER OF DEGREES OF FREEDOM
C     NUMNP  = NUMBER OF NODAL POINTS
C     NODE   = NODE NUMBER OF B.C.
C     IDF    = DEGREE OF FREEDOM FOR NODE
C     TIME   = TIME
C     LDOFU  = ACTIVE DEGREE OF FREEDOM ARRAY
C     CONSTR = ARRAY OF SPECIFIED NONZERO BOUNDARY CONDITIONS
C     XYZ    = nodal coordinates
C     iflag  = flag for user to set (not equal to 0) if coordinates
C              are updated
C     NODEPR = reverse permutation array
C              node(external) = NODEPR(NODE) where NODE = internal node no.
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION SOL(*),NUMEQA(NUMNP,NDOF),CONSTR(*),LDOFU(*)
      DIMENSION NODEPR(*),XYZ(NUMNP,*)
      ZRO = 0.D0
C
      CALL ERMSGS (8902, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRBCP (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR,UTNG)
C
C     USER DEFINED B COEFFICIENT FOR PERMEABILITY
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = B COEFFICIENT FOR PERMEABILITY EQUATION
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C     UTNG  = NORM OF VELOCITY (SPEED)
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS),UTNG(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (8401, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRBCC (VAL,NODE,TIME,NDP,VARN,X,Y,Z,XNORM,YNORM,
     1                   ZNORM,SRFTNS)
C
C     USER DEFINED BOUNDARY CONDITION FOR CONTACT ANGLE
C
C     VAL    = COMPUTED (SPECIFIED) BOUNDARY CONDITION (IN DEGREES)
C     NODE   = NODE NUMBER OF B.C.
C     TIME   = TIME
C     NDP    = NUMBER OF NODES IN THE ELEMENT
C     VARN   = ARRAY CONTAINING SOLUTION VARIABLES IN THIS ELEMENT
C     X      = X-COORDINATE OF THE NODE
C     Y      = Y-COORDINATE OF THE NODE
C     Z      = Z-COORDINATE OF THE NODE
C     XNORM  = X-COMPONENT OF THE LOCAL NORMAL VECTOR
C     YNORM  = Y-COMPONENT OF THE LOCAL NORMAL VECTOR
C     ZNORM  = Z-COMPONENT OF THE LOCAL NORMAL VECTOR
C     SRFTNS = VALUE OF SURFACE TENSION AT THIS NODE
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION VARN(NDP,*)
      ZRO = 0.D0
C
      CALL ERMSGS (8903, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRBCN (VAL,NODE,IDF,TIME,SOL,NUMEQA,NDOF,NUMNP,LDOFU,
     1                   CONSTR,NODEPR,XYZ,iflag)
C
C     USER DEFINED BOUNDARY CONDITIONS
C
C     VAL    = COMPUTED (SPECIFIED) BOUNDARY CONDITION
C     SOL    = GLOBAL SOLUTION VECTOR
C     NUMEQA = GLOBAL EQUATION NUMBER ARRAY
C     NDOF   = ACTIVE NUMBER OF DEGREES OF FREEDOM
C     NUMNP  = NUMBER OF NODAL POINTS
C     NODE   = NODE NUMBER OF B.C.
C     IDF    = DEGREE OF FREEDOM FOR NODE
C     TIME   = TIME
C     LDOFU  = ACTIVE DEGREE OF FREEDOM ARRAY
C     CONSTR = ARRAY OF SPECIFIED NONZERO BOUNDARY CONDITIONS
C     XYZ    = nodal coordinates
C     iflag  = flag for user to set (not equal to 0) if coordinates
C              are updated
C     NODEPR = reverse permutation array
C              node(external) = NODEPR(NODE) where NODE = internal node no.
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION SOL(*),NUMEQA(NUMNP,NDOF),CONSTR(*),LDOFU(*)
      DIMENSION NODEPR(*),XYZ(NUMNP,*)
      ZRO = 0.D0
C
      CALL ERMSGS (8901, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRBDY (NELT,NE,NG,BDYF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR,IOPT)
C
C     USER DEFINED BODY FORCES
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     BDYF  = BODY FORCES
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C     IOPT  = 0 - BODY FORCE
C     IOPT  = 1 - LORENTZ FORCE
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION BDYF(3,NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      IP1 = IOPT + 1
      GO TO (10,20) IP1
   10 CALL ERMSGS (7101, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
   20 CALL ERMSGS (7102, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRBTA (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR,DCOEF,
     2                   IFLAG)
C
C     USER DEFINED VOLUME EXPANSION COEFFICIENTS
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = VOLUME EXPANSION COEFFICIENT
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
C     IFLAG = 1 AND DCOEF ARE ONLY USED FOR DENSITY VARIATIONS IN
C             THE CONTINUITY EQUATION
C
C     IFLAG = 1 , COMPUTE THE GRADIENTS OF BETA W.R.T. TEMPERATURE AND
C                 SPECIES
C           = 0 , DON'T COMPUTE GRADIENTS
C     IFLAG IS SET IN FISOLV
C
C     DCOEF = ARRAY OF THE GRADIENTS OF BETA WHERE DCOEF(J,I) IS THE
C             (I-1)TH SPECIES (I-1=0 IS TEMPERATURE) FOR THE JTH POINT
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS),DCOEF(NPTS,*)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (7701, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRCAP (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED SPECIES CAPACITY
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = SPECIES CAPACITY
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (8601, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRCND (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED CONDUCTIVITY
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = CONDUCTIVITY
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (7501, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRCNV (NELT,NE,NG,HC,VARI,NDFCD,LDOFU,SHP,XYZL,
     1                   PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED CONVECTION HEAT TRANSFER COEFFICIENT
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     HC    = CONVECTION COEFFICIENT
C     VARI  = SOLUTION VARIABLE ARRAY AT THE INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION HC(NPTS),SHP(NPTS,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (7901, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRDIF (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED SPECIES DIFFUSIVITY
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = SPECIES DIFFUSIVITY
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (8701, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USREMS (NELT,NE,NG,HC,VARI,NDFCD,LDOFU,SHP,XYZL,
     1                   PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED EMISSIVITY COEFFICIENT FOR RADIATION
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     HC    = EMISSIVITY COEFFICIENT
C     VARI  = SOLUTION VARIABLE ARRAY AT THE INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION HC(NPTS),SHP(NPTS,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (7801, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRENT (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED ENTHALPY
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = ENTHALPY
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     PROP  = USER DEFINED PARAMETERS
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (7402, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      subroutine usreos (nelt,ne,nmgp,dens,uv,temp,pres,spec,nspecs,
     1                   xyzl,prop,npts,wtmol,ugasr,presrf,sos,ierr)
C
C     user defined equation of state
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     dens  = density
c     uv    = velocity components
c             u-velocity is in uv(npts,1)
c             v-velocity is in uv(npts,2)
c             w-velocity is in uv(npts,3)
C     TEMP  = TEMPERATURE
C     pres  = pressure
C     SPEC  = ARRAY OF SPECIES VALUES WHERE SPEC(J,I) IS THE ITH active
C                SPECIES AT THE JTH POINT
C     XYZL  = X,Y,Z COORDINATES
C     PROP  = USER DEFINED PARAMETERS
c     wtmol = array of molecular weights where wtmol(i) is the
c             molecular weight of the ith species and the last
c             entry in the array (i.e., wtmol(16)) contains
c             the molecular weight of the carrier fluid
c     ugasr = universal gas constant
c     presrf= reference pressure
c     sos   = local speed of sound
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      dimension dens(npts),uv(npts,*),temp(npts),pres(npts),spec(npts,*)
      dimension wtmol(16),sos(npts),prop(*),xyzl(npts,*)
      ZRO = 0.D0
C
      CALL ERMSGS (7007, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USREOSLP (VAL,NODE,IDF,TIME,SOL,NUMEQA,NDOF,NUMNP,
     1        LDOFU,CONSTR,DERIVATIVE,NDFCD,NODEPR,XYZ,iflag)
C
C     THE  ELECTRO-OSMOTIC SLIP BOUNDARY CONDITION IS IMPLEMENTED
C     HERE THROUGH USER SUBROUTINE. THE IMPLEMENTATION OF THE ABOVE
C     BOUNDARY CONDITION REQUIRES THE DERIVATIVES OF THE POTENTIAL
C     FIELD ON THE SOLID SURFACE. 
C
C     VAL    = COMPUTED (SPECIFIED) BOUNDARY CONDITION
C     SOL    = GLOBAL SOLUTION VECTOR
C     NUMEQA = GLOBAL EQUATION NUMBER ARRAY
C     NDOF   = ACTIVE NUMBER OF DEGREES OF FREEDOM
C     NUMNP  = NUMBER OF NODAL POINTS
C     NODE   = NODE NUMBER OF B.C.
C     NDFCD   = NODE NUMBER OF B.C.
C     IDF    = DEGREE OF FREEDOM FOR NODE
C     TIME   = TIME
C     LDOFU  = ACTIVE DEGREE OF FREEDOM ARRAY
C     CONSTR = ARRAY OF SPECIFIED NONZERO BOUNDARY CONDITIONS
C     DERIVATIVE = DERIVATIVES OF SPECIES (WHICH ACTS AS A POTENTIAL
C                  FIELD FOR THE EHD PROBLEMS) AT THE BOUNDARY NODES. 
C     XYZ    = nodal coordinates
C     iflag  = flag for user to set (not equal to 0) if coordinates
C              are updated
C     NODEPR = reverse permutation array
C              node(external) = NODEPR(NODE) where NODE = internal node no.
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION SOL(*),NUMEQA(NUMNP,NDOF),CONSTR(*),LDOFU(*)
      DIMENSION NODEPR(*),XYZ(NUMNP,*)
      DIMENSION DERIVATIVE(NDFCD)
      ZRO = 0.D0
C
      CALL ERMSGS (7008, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
C
CC      VAL = ZRO
CC      GASCONST =  8.3144D7      /* Gas Constant g-cm2/s2-K-mole*/
CC      TEMP     = 298.D0         /* Temperature K */
CC      VISCOSITY = 1.D-2         /* Viscosity gm/cm-s */
CC      Q = 1.E13                 /* Surface charge concentration /cm2 */
CC      FARADAY= 9.65D11
CC      DIELECTRIC=7.08D-5
CC      ELEMENTARYCHARGE=1.60199D-12
CC      CHARGE=1.0
CC      IEQ = NUMEQA(NODE,LDOFU(KDS+1))
CC      IF (IEQ.GT.0) THEN
CC        C1 = SOL(IEQ)
CC      ELSE IF (IEQ.EQ.0) THEN
CC        C1 = ZRO
CC      ELSE
CC        C1 = CONSTR(-IEQ)
CC      ENDIF
CC      V1 = 2.D0*DIELECTRIC*GASCONST*TEMP*C1
CC      IF (V1.LE.ZRO) RETURN
CCC    
CC      TEMPVAL = r_ASINH(Q*ELEMENTARYCHARGE/(2.D0*SQRT(V1)))
CC      IF ( (IDF.GE.KDU) .AND. (IDF.LE.KDW) ) THEN
CC              VAL = -2*DIELECTRIC*GASCONST*TEMP*
CC     1              DERIVATIVE(IDF)*TEMPVAL
CC     2              /(CHARGE*FARADAY*VISCOSITY)
CC      ENDIF
      RETURN
      END
C
      SUBROUTINE USRFNT (VARI ,DVARI,PCUX ,PCUY ,PCUZ ,PCLC ,PCLT ,PCLD,
     1                   FLDN ,FLVS ,FLSP ,FLKN ,FLDF ,PCDN ,PDNL ,PDNS,
     2                   DIAM ,DMST ,PHIC ,PHIS ,PSPV ,PSPL ,PSPS ,PSPH,
     3                   XYZL ,IDSPC,PCLTV,PCLTB,PMASS,PMAST,TIME ,ALPH,
     4                   LDOFU,NDFCD,IERR ,INDP ,IOPT ,USLH ,USRH)
C
C
C     COMPUTES ADDITIONAL USER DEFINED DRIVING TERMS FOR
C     PARTICULATE EQUATIONS - IN THE FOLLOWING SECTION THE
C     TERM "VOLATILE" REFERS TO THE EVAPORATING PORTION OF
C     A DRYING PARTICLE OR THE COMBUSTIBLE COMPONENT OF A
C     PARTICLE UNDERGOING COMBUSTION
C
C     MOMENTUM -  USER DEFINED FORCES          (INDP = 1,2,3)
C     MASS     -  USER DEFINED MASS FLOW RATE  (INDP = 4)
C                    USLH RETURNS DRDT
C                    URHS RETURNS DMDT
C     ENERGY   -  USER DEFINED HEAT TRANSFER   (INDP = 5)
C
C     OUTPUT
C            USRH      EXPLICIT CONTRIBUTION
C            USLH      IMPLICIT CONTRIBUTION
C
C     INPUT
C
C            IDSPC     SPECIE ASSOCIATED WITH THIS PARTICLE
C
C            PCUX      UX COMPONENT OF PARTICLE VELOCITY
C            PCUY      UY COMPONENT OF PARTICLE VELOCITY
C            PCUZ      UZ COMPONENT OF PARTICLE VELOCITY
C            PCLC      PARTICLE SURFACE CONCENTRATION
C            PCLT      PARTICLE TEMPERATURE
C            PCLD      PARTICLE DIAMETER
C
C            PDNL      PARTICLE DENSITY "VOLATILE"
C            PDNS      PARTICLE DENSITY "NON VOLATILE CORE"
C            PCDN      EQUIVALENT (VOLATILE + NON VOLATILE) DENSITY
C            NOTE(1)   THE DISTINCTIONS ARE ONLY RELEVANT FOR MIXED
C                      PARTICLES
C
C            PSPS      SPECIFIC HEAT OF "NON VOLATILE"
C            PSPL      SPECIFIC HEAT OF "VOLATILE"
C            PSPV      SPECIFIC HEAT OF VAPOR OF PARTICLE
C            PSPH      SPECIFIC HEAT OF PARTICLE IS THE MASS WEIGHTED
C                      AVERAGE  OF THE  VOLATILE AND NON  VOLATILE
C                      SUBSTANCES IN THE PARTICLE
C            NOTE(2)   THE DISTINCTIONS ARE ONLY RELEVANT FOR MIXED
C                      PARTICLES
C
C            DIAM      (NON VOLATILE) DIAMETER
C            DMST      INITIAL PARTICLE DIAMETER
C                      (VOLATILE + NON VOLATILE)
C            PHIC      CURRENT VOLATILE MASS FRACTION (TIME TSTRT)
C            PHIS      INITIAL VOLATILE MASS FRACTION
C            PMASS     CURRENT PARTICLE MASS
C            PMAST     INITIAL PARTICLE MASS
C
C            TIME      CURRENT TIME
C            ALPH      CURRENT TIME INCREMENT
C
C            FLDN      FLUID DENSITY
C            FLVS      FLUID VISCOSITY
C            FLSP      FLUID SPECIFIC HEAT
C            FLKN      FLUID CONDUCTIVITY
C            FLDF      DIFFUSIVITY COEFFICIENT OF SPECIE IDSPC
C
C            PCLTV     VAPORIZATION TEMPERATURE
C            PCLTB     BOILING TEMPERATURE
C
C            VARI      ARRAY OF SOLUTION VARIABLES AT PARTICLE LOCATION
C            DVARI     GRADIENTS OF SOLUTION VARIABLES AT PARTICLE
C                      LOCATION
C            LDOFU     POINTER ARRAY FOR ACCESSING VARI AND DVARI
C                      INFORMATION
C            XYZL      X,Y,Z COORDINATES
C
C            IOPT      USER DEFINED OPTION FLAG
C                      IF INDP = 1,2,3  IOPT = VALUE OF DRUSER KEYWORD
C                      IF INDP = 4  IOPT = VALUE OF TEMPUSER KEYWORD
C                      IF INDP = 5  IOPT = VALUE OF PHSUSER KEYWORD
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
#include "TAPES.COM"
      DIMENSION XYZL(NDFCD),VARI(*),DVARI(NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      IERR = 1
      IF (INDP.EQ.1.OR.INDP.EQ.2.OR.INDP.EQ.3) THEN
          CALL ERMSGP (7500,2,2,1,IOUT,BATCH,0,0,0,0,0,ZRO,ZRO,
     1                 'PARTICLE MOMENTUM (USRFNT)',' ',' ')
      ELSEIF (INDP.EQ.4) THEN
          CALL ERMSGP (7500,2,2,1,IOUT,BATCH,0,0,0,0,0,ZRO,ZRO,
     1                 'PARTICLE MASS TRANSFER (USRFNT)',' ',' ')
      ELSEIF (INDP.EQ.5) THEN
          CALL ERMSGP (7500,2,2,1,IOUT,BATCH,0,0,0,0,0,ZRO,ZRO,
     1                 'PARTICLE HEAT TRANSFER (USRFNT)',' ',' ')
      ENDIF
C
      RETURN
      END
C
      SUBROUTINE USRGAP (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR,ZIGN,
     2                   ZIGT,GAP,TEMP2)
C
C     USER DEFINED GAP ELEMENT HEAT TRANSFER COEFFICIENT
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = SPECIFIC HEAT
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C     ZIGN  = NORMAL STRESS AT INTEGRATION POINTS
C     ZIGT  = TANGENTIAL STRESS AT INTEGRATION POINTS
C     GAP   = GAP WIDTH ALONG ELEMENT SIDE
C     TEMP2 = WALL TEMPERATURES
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS),GAP(NPTS),TEMP2(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (7901, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRINI (INIVAR)
C
C     USER SET FLAGS FOR ENSURING OF THE COMPUTATION
C     OF VARI AND DVARI ARRAYS FOR VARIABLES CHOSEN BY THE
C     USER
C     INIVAR COMES IN SET TO ZERO, IF THE USER DESIRES A VARIABLE
C     TO BE COMPUTED AT THE INTEGRATION POINTS FOR A USER SUBROUTINE
C     THEN SETTING INIVAR TO A NONZERO VALUE WILL FORCE COMPUTATION
C     THE INDICES ARE  KDU = U-VELOCITY      KDV = V-VELOCITY
C                      KDW = W-VELOCITY      KDK = KINETIC ENERGY
C                      KDE = DISSIPATION     KDT = TEMPERATURE
C                      KDS+I = FOR SPECIES I (I = 1 to 15)
C                      KDP = PRESSURE
C                      KDX = X-DISPLACEMENT  KDY = Y-DISPLACEMENT
C                      KDZ = Z-DISPLACEMENT  
C     
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION INIVAR(*)
      RETURN
      END
C
      SUBROUTINE USRINT (IOPT,X,Y,Z,VEC,NUMNP)
C
C     USER DEFINED INITIAL CONDITIONS
C
C     IOPT  = DEGREE OF FREEDOM FLAG
C     X     = X COORDINATES
C     Y     = Y COORDINATES
C     Z     = Z COORDINATES
C     VEC   = INITIAL CONDITIONS
C             IF INITIAL CONDITIONS ARE FOR VELOCITY THEN
C                VEC IS 2(OR3) TIMES NUMNP
C             ELSE VEC IS NUMNP LONG
C     NUMNP = GLOGAL NUMBER OF NODAL POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION X(NUMNP),Y(NUMNP),Z(NUMNP),VEC(*)
      ZRO = 0.D0
C
      CALL ERMSGS (8201, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRJAC (NELT,NE,NG,STIF,MNDF,VARI,DVARI,NDFCD,LDOFU,
     1                   SHP,DSDX,XYZL,VARN,DETER,TIME,NPTS,NDP,MNDP)
C
C     USER DEFINED JACOBIAN MATRIX
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     STIF  = ELEMENT JACOBIAN MATRIX
C     MNDF  = DIMENSION OF STIFFNESS MATRIX
C     VARI  = SOLUTION VARIABLE ARRAY AT THE INTEGRATION POINTS
C     DVARI = GRADIENT OF SOLUTION VARIABLE ARRAY AT THE
C             INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES
C     VARN  = VALUES OF NODAL DEGREES OF FREEDOM AT NODAL
C             POINTS OF ELEMENT NELT
C     DETER = (JACOBIAN OF ISOPARAMETRIC TRANSFORMATION)*(INTGRATION
C             WEIGHTS)
C     TIME  = TIME
C     MNDF  = DIMENSION OF ELEMENT JACOBIAN MATRIX
C     NPTS  = NUMBER OF POINTS
C     NDP   = NUMBER OF NODAL POINTS
C
#include "IMPLCT.COM"
C     INCLUDE 'INTPRP.COM'
      DIMENSION STIF(MNDF,MNDF),SHP(NPTS,NDP),DSDX(NPTS,NDP,NDFCD)
      DIMENSION DETER(NPTS),varn(ndp,*),XYZL(NPTS,ndfcd),LDOFU(*)
      DIMENSION VARI(NPTS,*),DVARI(NPTS,NDFCD,*)
      ZRO = 0.D0
C
      RETURN
      END
C
      SUBROUTINE USRLAT (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED LATENT HEAT
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = LATENT HEAT
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "TAPES.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      IERR = 1
      CALL ERMSGS (7502, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRLHS (NELT,NE,NG,STIF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,DETER,VARN,MNDF,NDP,mndp,NPTS,TIME)
C
C     CALCULATE THE ELEMENT MATRIX CONTRIBUTION FOR A USER
C     SUPPLID IMPLICIT TERM IN THE MOMENTUM , ENERGY EQUATIONS
C     AND SPECIES EQUATIONS
C        COEF*PHI*PHI IS COMPUTED FOR EACH EQUATION AND ADDED
C        TO THE LOCAL ELEMENT STIFFNESS MATRIX
C
C     INPUT       - SHP,DSDX,XYZL,DETER,VARN,NDP,NPTS,MNDF
C     NELT        = GLOBAL ELEMENT NUMBER
C     NE          = LOCAL ELEMENT NUMBER
C     NG          = GROUP NUMBER
C     MNDP        = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME        = TIME
C     NPTS        = NUMBER OF POINTS
C
C     OUTPUT      - STIF
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION STIF(MNDF,MNDF),SHP(NPTS,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION VARI(NPTS,*),DVARI(NPTS,NDFCD,*),varn(ndp,*)
      DIMENSION DSDX(NPTS,NDFCD,MNDP),LDOFU(*)
      DIMENSION DETER(NPTS)
C
C     CALCULATE USER CONTRIBUTION TO THE STIFFNESS MATRIX
C
      RETURN
C
      END
C
      SUBROUTINE USRMBC (VAL,NODE,TIME,SOL,NUMEQA,NDOF,NUMNP,LDOFU,
     1                   CONSTR,NODEPR,XYZ,iflag,ierr)
c
c       Purpose:
c         This routine decides whether a given point (x,y,z)
c         is affected by the moving body constraint.
c
c       Input:
c         node            - internal node number
c         time            - value of time counter
c         XYZ             - nodal position
c         NODEPR          - reverse permutation array
c               nodext = NODEPR(NODE)
c
c       Output:
c         VAL             - array imposed values for dofs
c         iflag()         - array of flags for dofs,
c                             0 if not affected, 1 if affected
c
#include "IMPLCT.COM" 
#include "PARUSR.COM" 
        DIMENSION SOL(*),ID(NUMNP,NDOF),LDOFU(*),XYZ(NUMNP,*) 
        DIMENSION IFLAG(*),VAL(*),NODEPR(*),CONSTR(*) 
C 
C
      ZRO = 0.D0 
C 
C     RETRIEVE NODAL COORDINATES: 
C 
      X = XYZ(NODE,1) 
      Y = XYZ(NODE,2) 
C 
C     INITIAL CALCULATIONS: 
C 
      VV = SQRT(X**2+Y**2) 
      IF (VV .LT. 1.0D-8) RETURN 
      XV = X / VV 
      YV = Y / VV 
C 
      THICK  = 0.02D0 
      OMEGA  = 1.047D0 
      ANGLE = OMEGA * TIME 
      PI     = 3.14159265359D0 
C 
C     LOOP OVER THE FOUR BLADES: 
C 
      DO 10 I=1,4 
C 
        XANG = COS(ANGLE+FLOAT(I-1)*PI/2.0D0) 
        YANG = SIN(ANGLE+FLOAT(I-1)*PI/2.0D0) 
        DIST = SQRT( (XV-XANG)**2 + (YV-YANG)**2 ) * VV 
C 
C       IF THE NODE COINCIDES WITH THE BLADE: 
C 
        IF (DIST .LT. THICK) THEN 
C 
C         SET FLAGS AND VALUES FOR THE TWO 
C         VELOCITY COMPONENTS: 
C 
          IFLAG(KDU) = 1 
          IFLAG(KDV) = 1 
C 
          VAL(KDU) = - OMEGA * Y 
          VAL(KDV) =   OMEGA * X 
C 
        ENDIF 
C 
   10 CONTINUE 
C 
      RETURN 
      END
      RETURN
      END
C
      SUBROUTINE USRMLT (NELT,NE,NG,COEF,VARI,NDFCD,LDOFU,SHP,XYZL,
     1                   PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED LIQUIDUS
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = MELT LIQUIDUS
C     VARI  = SOLUTION VARIABLE ARRAY AT THE INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS),SHP(NPTS,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (8801, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRMU (NELT,NE,NG,VISC,VARI,DVARI,NDFCD,LDOFU,SHP,DSDX,
     1                  XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED NONNEWTONIAN VISCOSITY MODEL
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     VISC  = VISCOSITY
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION VISC(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (7401, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRMXL (NELT,NE,NG,VISC,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR,
     2                   DENS)
C
C     USER DEFINED MIXING LENGTH
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     VISC  = VISCOSITY
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION LDOFU(*)
      DIMENSION VISC(NPTS),DENS(NPTS),SHP(NPTS,MNDP),
     1          DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD),
     2          PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*)
      ZRO = 0.D0
C
      CALL ERMSGS (7301, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRPRT (VARI ,DVARI,PCUX ,PCUY ,PCUZ ,PCLC ,PCLT ,PCLD,
     1                   FLDN ,FLVS ,FLSP ,FLKN ,FLDF ,PCDN ,PSPH ,PSPV,
     2                   XYZL ,IDSPC,PCLTV,PCLTB,TIME ,ALPH ,IOPT ,IERR,
     3                   LDOFU,NDFCD,INDP ,COEF )
C
C     USER DEFINED DRAG COEFFICIENT, NUSSELT OR SHERWOOD NUMBER
C     FOR PARTICLE TRACKING
C
C     OUTPUT
C            DRAG COEFFICIENT (INDP = 1)
C            SHERWOOD NUMBER (INDP = 2)
C            NUSSELT NUMBER (INDP = 3)
C
C     INPUT
C            IDSPC     SPECIE ASSOCIATED WITH THIS PARTICLE
C            PCUX      UX COMPONENT OF PARTICLE VELOCITY
C            PCUY      UY COMPONENT OF PARTICLE VELOCITY
C            PCUZ      UZ COMPONENT OF PARTICLE VELOCITY
C            PCLC      PARTICLE SURFACE CONCENTRATION
C            PCLT      PARTICLE TEMPERATURE
C            PCLD      PARTICLE DIAMETER
C            PCDN      PARTICLE DENSITY
C            PSPH      PARTICLE SPECIFIC HEAT
C            PSPV      SPECIFIC HEAT OF DIFFUSING VAPOR
C            TIME      CURRENT TIME
C            ALPH      CURRENT TIME INCREMENT
C            FLDN      FLUID DENSITY
C            FLVS      FLUID VISCOSITY
C            FLSP      FLUID SPECIFIC HEAT
C            FLKN      FLUID CONDUCTIVITY
C            FLDF      DIFFUSIVITY COEFFICIENT OF SPECIE IDSPC
C            VARI      ARRAY OF SOLUTION VARIABLES AT PARTICLE LOCATION
C            DVARI     GRADIENTS OF SOLUTION VARIABLES AT PARTICLE LOCATION
C            LDOFU     POINTER ARRAY FOR ACCESSING VARI AND DVARI INFORMATION
C            XYZL      X,Y,Z COORDINATES
C            IOPT      USER OPTION FLAG
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
#include "TAPES.COM"
      DIMENSION XYZL(NDFCD),VARI(*),DVARI(NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      IERR = 1
C
      IF (INDP.EQ.1) THEN
          CALL ERMSGP (7500,2,2,1,IOUT,BATCH,0,0,0,0,0,ZRO,ZRO,
     1                 'DRAG COEFFICIENT (USRPRT)',' ',' ')
      ELSEIF (INDP.EQ.2) THEN
          CALL ERMSGP (7500,2,2,1,IOUT,BATCH,0,0,0,0,0,ZRO,ZRO,
     1                 'SHERWOOD NUMBER (USRPRT)',' ',' ')
      ELSEIF (INDP.EQ.3) THEN
          CALL ERMSGP (7500,2,2,1,IOUT,BATCH,0,0,0,0,0,ZRO,ZRO,
     1                 'NUSSELT NUMBER (USRPRT)',' ',' ')
      ENDIF
      RETURN
      END
C
      SUBROUTINE USRPSI (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR,IOPT)
C
C     USER DEFINED SECOND ORDER FLUID COEFFICIENTS
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = PSI(IOPT) FOR THE SECOND ORDER FLUID
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
C     IOPT  = 1 - FIRST COEFFICIENT FOR SECOND ORDER FLUID
C     IOPT  = 2 - SECOND COEFFICIENT FOR SECOND ORDER FLUID
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      if (iopt.eq.1) then
         CALL ERMSGS (7202, 0, 3,
     1                0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      elseif (iopt.eq.2) then
         CALL ERMSGS (7203, 0, 3,
     1                0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      endif
C
      RETURN
      END
C
      SUBROUTINE USRREF (NELT,NE,NG,REFENV,VARI,NDFCD,LDOFU,
     1                   XYZL,PROP,TIME,NPTS,IPRP,IERR)
C
C     USER DEFINED REFERENCE ENVIRONMENTAL CONDITION
C     FOR BOUNDARY CONDITIONS WITH VARYING ENVIRONMENTS
C
C     NELT    = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG      = GROUP NUMBER
C     REFENV  = USER DEFINED ENVIRONMENTAL REFERENCE
C     VARI    = SOLUTION VARIABLE ARRAY AT THE INTEGRATION POINTS
C     XYZL    = X,Y,Z COORDINATES
C     PROP    = VALUE OF ENVIRONMENTAL REFERENCE
C     TIME    = TIME
C     NPTS    = NUMBER OF POINTS
C     IPRP    = 0 - TEMPERATURE FOR CONVECTION
C              -1 - TEMPERATURE FOR EMISSIVITY (RADIATION)
C               I - SPECIES I
C              -2 - PRESSURE
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION REFENV(NPTS),VARI(NPTS,*),XYZL(NPTS,NDFCD)
      DIMENSION LDOFU(*)
      ZRO = 0.D0
C
      if ( iprp.eq.-2 ) then
C
C        user defined pressure for free surfaces
C
         CALL ERMSGS (7905, 0, 3,
     1                0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      elseif ( iprp.eq.-1 ) then
C
C        user defined reference temperature for
C        radiation heat transfer
C
         CALL ERMSGS (7906, 0, 3,
     1                0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      elseif ( iprp.eq.0 ) then
C
C        user defined reference temperature for
C        convection heat transfer
C
         CALL ERMSGS (7906, 0, 3,
     1                0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      else
C
C        user defined reference species concentration
C
         CALL ERMSGS (7907, 0, 3,
     1                IPRP,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      endif
      RETURN
      END
C
      SUBROUTINE USRrgh (NELT,NE,NG,rough,xl,yl,zl,PROP,NPTS,IERR)
C
C     USER DEFINED surface roughness
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     rough  = surface roughness
c     xl     = x coordinates
c     yl     = y coordinates
c     zl     = z coordinates
C     PROP  = USER DEFINED PARAMETERS
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION rough(NPTS),xl(npts),yl(npts),zl(npts),PROP(*)
      ZRO = 0.D0
C
      CALL ERMSGS (8602, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
c
      SUBROUTINE USRRHS (NELT,NE,NG,BE,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,DETER,VARN,MNDF,NDP,mndp,NPTS,TIME)
C
C     CALCULATE RIGHT HAND SIDE VECTOR DUE TO USER CONTRIBUTION
C
C     INPUT       - SHP,DSDX,XYZL,DETER,VARN,NDP,NPTS,MNDF
C     NELT        = GLOBAL ELEMENT NUMBER
C     NE          = LOCAL ELEMENT NUMBER
C     NG          = GROUP NUMBER
C     MNDP        = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME        = TIME
C     NPTS        = NUMBER OF POINTS
C
C     OUTPUT      - BE
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION BE(MNDF),SHP(NPTS,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION VARI(NPTS,*),DVARI(NPTS,NDFCD,*),varn(ndp,*)
      DIMENSION DSDX(NPTS,NDFCD,MNDP),LDOFU(*)
      DIMENSION DETER(NPTS)
C
      RETURN
      END
C
      SUBROUTINE USRRXN (NELT,NE,NG,REAC,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR,IOPT)
C
C     USER DEFINED SOURCE FOR ENERGY OR SPECIES EQUATIONS
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     REAC  = REACTION TERM
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
C     IOPT = 0   TEMPERATURE REACTION TERM
C     IOPT = N   REACTION TERM FOR SPECIES N (0<N<16)
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION REAC(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
      IF (IOPT.EQ.0) THEN
         CALL ERMSGS (7003, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      ELSE
         CALL ERMSGS (7004, 0, 3,
     1             IOPT,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      ENDIF
      RETURN
      END
C
      SUBROUTINE USRXNJ (NELT,NE,NG,WV,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,ISRF,ISP,
     2                   STIF,DETER,MNDF,NSPEC)
C
C     JACOBIAN FOR USER-DEFINED SPECIES SOURCE TERMS
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     WV    = WORK VECTOR
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF INTEGRATION POINTS
c     STIF  = LOCAL STIFFNESS MATRIX
C     DETER = PRODUCT OF INTEGRATION WEIGHTS AND ISOPARAMETERIC JACOBIAN
C     NSPEC = NUMBER OF SPECIES
C
C     ISP   = 0   TEMPERATURE REACTION TERM
C     ISP   = N   REACTION TERM FOR SPECIES N (0<N<16)
C
C     ISRF  = 1  USE BULK PHASE REACTION MODEL
C     ISRF  = 2  USE SURFACE REACTION MODEL
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
C
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      DIMENSION STIF(MNDF,MNDF),DETER(*),WV(NPTS)
C
      ZRO = 0.D0
C
      RETURN
      END
C
      SUBROUTINE USRSKE (NELT,NE,NG,SORCE,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR,IOPT,
     2                   visct,viscl,den)
C
C     USER DEFINED SOURCES FOR KINETIC ENERGY OR DISSIPATION
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     SORCE = TURBULENT KINETIC ENERGY OR DISSIPATION (RETURNED VALUES)
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C     IOPT = 0   ENERGY EQUATION
C     IOPT = N   TRANSPORT EQUATION FOR SPECIES N (0<N<16)
C     VISCT = TURBULENT VISCOSITY
C     VISCL = LAMINAR VISCOSITY
C     IOPT = 0   KINETIC ENERGY SOURCES
C     IOPT = 1   DISSIPATION SOURCES
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION SORCE(NPTS),visct(npts),viscl(npts),den(npts)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
C
      ZRO = 0.D0
      IF (IOPT.EQ.0) THEN
         CALL ERMSGS (7005, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      ELSE
         CALL ERMSGS (7006, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      ENDIF
      RETURN
      END
C
      SUBROUTINE USRSLP (NELT,NE,NG,COEF,VARI,NDFCD,LDOFU,SHP,XYZL,PROP,
     1                   TIME,NPTS,ndp,MNDP,IERR,utng,sigman,sigmat)
C
C     USER DEFINED SLIP COEFFICIENT
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = SLIP COEFFICIENT
C     VARI  = SOLUTION VARIABLE ARRAY AT THE INTEGRATION POINTS
C     UTNG  = TANGENTIAL VELOCITY COMPONENT
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS),SHP(NPTS,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),LDOFU(*),utng(npts)
      DIMENSION sigman(npts),sigmat(npts)
      ZRO = 0.D0
C
      CALL ERMSGS (8001, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRSPH (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED SPECIFIC HEAT
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = SPECIFIC HEAT
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (7601, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRSPT (NELT,NE,NG,COEF,VARI,NDFCD,LDOFU,SHP,XYZL,
     1                   PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED SPECIES TRANSFER COEFFICIENT
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = SPECIES TRANSFER COEFFICIENT
C     VARI  = SOLUTION VARIABLE ARRAY AT THE INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS),SHP(NPTS,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (7902, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRSRC (NELT,NE,NG,SORCE,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR,IOPT)
C
C     USER DEFINED SOURCE FOR ENERGY OR SPECIES EQUATIONS
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     SORCE = HEAT OR SPECIES SOURCE (RETURNED VALUES)
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C     IOPT = 0   ENERGY EQUATION
C     IOPT = N   TRANSPORT EQUATION FOR SPECIES N (0<N<16)
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION SORCE(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
      IF (IOPT.EQ.0) THEN
         CALL ERMSGS (7001, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      ELSE
         CALL ERMSGS (7002, 0, 3,
     1             IOPT,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      ENDIF
      RETURN
      END
C
      SUBROUTINE USRSRF (NELT,NE,NG,COEF,VARI,NDFCD,LDOFU,SHP,XYZL,
     1                   PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED SURFACE TENSION
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = SURFACE TENSION
C     VARI  = SOLUTION VARIABLE ARRAY AT THE INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS),SHP(NPTS,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (8101, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRVAP (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED VAPOR PRESSURE
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = VAPOR PRESSURE
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "TAPES.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      IERR = 1
      CALL ERMSGS (7503, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRVSC (NELT,NE,NG,VISC,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED VISCOSITY
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     VISC  = VISCOSITY
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION VISC(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      CALL ERMSGS (7201, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRYNG (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)

C
C     USER DEFINED YOUNG MODULUS
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = YOUNG MODULUS
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      IERR = 1
      CALL ERMSGS (7503, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END

C
      SUBROUTINE USRPSN (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C
C     USER DEFINED POISSON COEFFICIENT
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = POISSON COEFFICIENT
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      IERR = 1
      CALL ERMSGS (7503, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRSPR (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C
C
C     USER DEFINED SPRING COEFFICIENT
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = SPRING COEFFICIENT
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      IERR = 1
      CALL ERMSGS (7503, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
      SUBROUTINE USRSEG (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C
C
C     USER DEFINED SOLUTE SEGREGATION COEFFICIENT
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = SOLUTE SEGREGATION COEFFICIENT
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS)
      DIMENSION SHP(NPTS,MNDP),DSDX(NPTS,NDFCD,MNDP),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),DVARI(NPTS,NDFCD,*),LDOFU(*)
      ZRO = 0.D0
C
      IERR = 1
      CALL ERMSGS (7503, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')

      RETURN
      END
C
      SUBROUTINE USRMOB (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED  MOBILITY COEFFICIENT
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  =  MOBILITY COEFFICIENT
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
C
      DIMENSION COEF(NPTS),SHP(MNDP,NPTS),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),LDOFU(*)
      DIMENSION DSDX(MNDP,NPTS,NDFCD)
      DIMENSION DVARI(NPTS,NDFCD,*)
      ZRO = 0.D0
C
      CALL ERMSGS (8702, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')
      RETURN
      END
C
C
      SUBROUTINE USRELD (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP,
     1                   DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP,IERR)
C
C     USER DEFINED  ELECTRICAL-CONDUCTIVITY COEFFICIENT
C
C     NELT  = GLOBAL ELEMENT NUMBER
C     NE    = LOCAL ELEMENT NUMBER
C     NG    = GROUP NUMBER
C     COEF  = ELECTRICAL-CONDUCTIVITY
C     VARI  = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS
C     DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS
C     LDOFU = pointer array for accessing vari and dvari information
C     XYZL  = X,Y,Z COORDINATES
C     SHP   = ELEMENT SHAPE FUNCTIONS
C     DSDX  = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION
C     PROP  = USER DEFINED PARAMETERS
C     MNDP  = FIRST DIMENSION OF SHAPE FUNCTION MATRICES
C     TIME  = TIME
C     NPTS  = NUMBER OF POINTS
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
      DIMENSION COEF(NPTS),SHP(MNDP,NPTS),XYZL(NPTS,NDFCD)
      DIMENSION PROP(*),VARI(NPTS,*),LDOFU(*)
      DIMENSION DSDX(MNDP,NPTS,NDFCD)
      DIMENSION DVARI(NPTS,NDFCD,*)
      ZRO = 0.D0
C
      CALL ERMSGS (8703, 0, 3,
     1             0,0,0,0,0,ZRO,ZRO,ZRO,' ',' ',' ')

      RETURN
      END
C
      SUBROUTINE USROPT (IOPT,DELTAS,IERR)
C
C     USER DEFINED FEED-BACK LOOP FOR MELT-INTERFACE PROBLEMS
C     RESTRICTIONS: 1. SEGREGRATED ALGORITHM
C                   2. STEADY-STATE
C
C     IOPT   = FEED-BACK LOOP FLAG
C              0 = NO FEED-BACK LOOP (OR) FEED-BACK LOOP CONVERGED
C                  THE CODE REVERTS BACK TO THE USUAL SEGREGRATED
C                  FREE-SURFACE ALGORITHM AND UPDATES THE FREE-SURFACE
C                  FOLLOWING A FREE-SURFACE ITERATION
C              1 = FEED-BACK LOOP IN PROGRESS
C                  DOES NOT UPDATE THE FREE-SURFACE BUT MODIFIES THE
C                  BCNO, BCFL, HTTRANS, OR ANY OTHER USER-DEFINED
C                  QUANTITY (ONLY ONE PARAMETER CAN BE OPTIMIZED THROUGH
C                  THE FEEDBACK LOOP)
C     DELTAS = TOTAL (AREA INTEGRATED) HEAT FLUX ACROSS THE MELT-INTERFACE
C     IERR   = ERROR FLAG
C
C     LOCAL VARIABLES
C     ===============
C     HCOEF  =  THE PARAMETER THAT IS BEING OPTIMIZED THROUGH THE
C               FEEDBACK LOOP. THIS VALUE IS PASSED ONTO THE USER
C               DEFINED USER-SUBROUTINE USRBCN, USRBCF, OR ANY OTHER
C               ROUTINE THROUGH THE COMMON BLOCK: COMMON/FDBACK/HCOEF.
C               THUS, THE USER NEEDS TO DO THE FOLLOWING:
C               INCLUDE THE COMMON BLOCK "COMMON/FDBACK/HCOEF"
C               IN THE USER SUBROUTINE AND ALSO NEED TO SET "VAL = HCOEF".
C
C     ALL THE REMAINING VARIABLES ARE LOCAL VARIABLES AND ARE SPECIFIC
C     TO A GIVEN FEED-BACK LOOP. SETTING UP THE APPROPRIATE FEEDBACK
C     LOOP STRATEGY IS THE RESPONSIBILITY OF THE USER.
C     ONE SAMPLE IMPLEMENTATION OF USER FEEDBACK IS
C     LEFT (COMMENTED OUT) FOR THE PURPOSE OF ILLUSTRATION.
C
#include "IMPLCT.COM"
#include "PARUSR.COM"
#include "TAPES.COM"
C
C THE DEFAULT USROPT RETURNS THE IOPT FLAG AS 0
C
      IOPT = 0
      RETURN
      END
C
C THE FOLLOWING IS A SAMPLE ILLUSTRATION OF USER FEEDBACK LOOP
C COMMENT THE ABOVE LINES AND PUT IN THE APPROPORAITE STRATEGY
C HERE
C     LOCAL VARIABLES (CONTD.)
C     ========================
C     DELTAH:   CONSTANT INCREMENT TO BE USED FOR INITIAL GUESS OF HCOEF
C     DELTAS:   change in surface position, computed in sgfupd
C     seps:     Error tolerance for convergence of surface position
C     ssum:     Cumulative sum of scoef for averaging
C     scoef:    Sensitivity coefficient (Jacobian deltah/deltas)
C     hrelax:   Relaxation to be applied to h changes
C     istep:    Current iteration of feedback loop
C     niter:    Maximum iteration number for feedback loop
C     inistep:  Initial number of steps over which DELTAH is applied for IC
C     iconv:    Convergence flag to be set when feedback is to be stopped
C               = 0 not converged, = 1 converged
C     iopt:     Flag to indicate whether optimization is to be activated or not
C               = 0, Not active, = 1, active. This may be needed for sgfupd.
C
C
C
C
C     COMMON/FDBACK/HCOEF
C     common/usrinp/istep,niter,inistep, iconv, iniflag,
C    1      ssum, seps, scoef, hrelax, deltah
C     F_INTEGER istep,niter,inistep, iconv, iniflag, iopt
C     F_DOUBLEPRECISION HCOEF, ssum, seps, scoef, hrelax, deltah
C
C the following are for initialization.  No need to change.
C     data istep/0/, iconv/0/, ssum/0.0/, iniflag/0/
C
C the following are user definable variables.
C     data inistep/5/, niter/100/
C     data hcoef/7000/, hrelax/1.0/, seps/1.0e-4/, deltah/100.0/
C
C     if (iniflag.eq.0) then
C         open(unit=99,file="xy")
C         write(99,*) 'iteration           FLX  aaaa          bbbb',
C    1'          ccc          '
C         write(99,*) 'crv1'
C         iniflag = 1
C         iopt = 1
C         return
C     endif
C  if converged or no optimization, return
C     if (iconv .eq. 1 ) then
C        iopt = 0
C        return
C     else
C        iopt = 1
C     endif
C
C  if initial call, turn on optimization, initialize variables
C     if (istep .eq. 0) then
C       iopt = 1
C       hcoef = hcoef + deltah
C       write(99,*) istep, hcoef
C       istep = istep + 1
C       return
C     endif
C  find out if jacobian claculations need to be done
C     if (istep .le. inistep) then
C  if initial guess stage, ramp up the h by deltah and compute scoef
C       hcoef = hcoef + deltah
C       ssum = ssum + deltah/deltas
C       scoef = ssum/istep
C       istep = istep + 1
C       write(99,*) istep, hcoef
Creturn
C     else
C
C  initial value of scoef is obtained, start feedback.
C
C      if (abs(deltas) .ge. seps) then
C  if feedback loop not converged, update hcoef and return
C           deltah = -scoef*deltas
C           hcoef = hcoef + deltah*hrelax
C           write(istat,100) istep, hcoef
C           write(99,*) istep, hcoef
C100   format('>> ITERATION: ',i4,', HEAT TRANSFER COEFFICIENT: ',e12.4)
C           istep = istep + 1
C if reached maximum iteration, set convergence flag and return
C           if (istep .gt. niter) then
C              write(istat,200)
C200   format('>> MAXIMUM NUMBER OF ITERATION REACHED FOR FEEDBACK LOOP',
C    1       /,'>> FEEDBACK LOOP TERMINATED')
C            iconv = 1
C            iopt = 0
C           endif
C           return
C      else
C  feedback converged, set convergence flag, turn off optimization and return
C          write(istat,300) hcoef
C300   format('>> FEEDBACK LOOP CONVERGED.  FINAL HT COEF.: ',E12.4)
C          iconv = 1
C          iopt = 0
C          return
C      endif
C     endif
C     RETURN
C     END
C
C
#undef PREPOST_MODULE
#define SOLVER_MODULE
C Patches would go here