XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX''"> ULIBFOR Contents - Example User-Supplied Subroutines
C############################################################################
C #
C This file contains dummy subroutines that can be altered by the user #
C and included with the QTRAN thermal analysis module. These #
C routines are called at key entry points in the calculation #
C loop as well as during initialization phases, reading of the #
C input data file, etc. The routines included are: #
C #
C ROUTINES: #
C #
C #
C UAFLOW --> Called from AFLOW when a user material property has #
C been specified for an advection conductor #
C #
C UCCPAC --> Called from CCAPAC when a user material property has #
C beeb defined for either the capacitor density or #
C specific heat. Even if the user properties are dummies #
C to flag a user defined capacitance determination, the #
C UPROP value must pass a value that will allow a #
C capacitor evaluation prior to the user defining it. #
C #
C UCNDUC --> Called from CONDUC when a user material property has #
C been defined for the thermal conductivity. Even if the #
C user properties are dummies to flag a user defined #
C conductance determination, the UPROP value must pass a #
C value that will allow conductor evaluation prior to the #
C user defining it. #
C #
C UEXITQ --> Called whenever EXITQ is called. This is just #
C before the standard output routine qtran.dat is closed. #
C Is called during error conditions as well as normal #
C exits. The UFNSHD routine will be call just before #
C this but only for normal exits. #
C #
C UFNSHD --> Called at the completion of a solution cycle when #
C crunch is going to be exited. #
C #
C UHVAL --> Called from QTRAN subroutine CONV0, this #
C subroutine is called whenever the user specifies #
C a convection configuration type of 1000 or #
C greater. This allows the user to build and #
C maintain custom convection configurations. #
C #
C UINIT1 --> Called from QTRAN subroutine INIT1, this #
C subroutine may be used to initialize any system #
C arrays provided by the user. UINIT1 is called #
C prior to reading in any input data. #
C #
C UINIT2 --> Called from QTRAN subroutine INIT2, this #
C subroutine may be used to set up pointer tables #
C or anything else that may need to be done after #
C reading in the input data file(s). #
C #
C UINPUT --> Called from QTRAN subroutine INPUT, this #
C subroutine may be used to read in customized #
C data not normally found in QTRAN's input data #
C file. It is called after QTRAN's normal input #
C data has been read in. #
C #
C ULOOP1 --> Called from QTRAN subroutine RESET1, this #
C subroutine may be used to perform any necessary #
C housekeeping or calculations prior to beginning #
C a steady state run or a transient time step. #
C #
C ULOOP2 --> Called from QTRAN subroutine RESET2, this #
C subroutine may be used to perform any necessary #
C housekeeping or calculations prior to beginning #
C an iteration for either steady state or #
C transient runs. #
C #
C ULOOP3 --> Called from QTRAN subroutine CRUNCH, this #
C subroutine is called at the beginning of #
C calculations for each node for each iteration. #
C Custom resistor types are one example of things #
C that might be appropriate for ULOOP3. #
C #
C ULOOP4 --> Called from QTRAN subroutine CRUNCH, this #
C subroutine is called at the end of calculations #
C for each node for each iteration immediately #
C after the temperature for that node has been #
C updated. #
C #
C ULOOP5 --> Called from QTRAN subroutine CRUNCH, this #
C subroutine is called at the end of each #
C steady state or transient iteration. This #
C routine can be used to perform auxiliary #
C calculations that must be performed in parallel #
C with the thermal calculations, e.g., fluid flow #
C calculations, mass transport calculations, or #
C anything else that is appropriate. #
C #
C ULOOP6 --> Called from QTRAN subroutine CRUNCH, this #
C subroutine is called only after a steady state #
C calculation has converged. #
C #
C ULOOP7 --> Called from QTRAN subroutine CRUNCH, this #
C subroutine is called after each transient time #
C step. #
C #
C UMCRPR --> Same as UMICRO call except the material property #
C parameters are passed as well. This enables the user #
C directly use material properties without having to #
C provide wrappers to properly define the dimensions #
C #
C UMICRO --> Called from QTRAN subroutine FLIB, this #
C subroutine is called whenever the user specifies #
C a microfunction option of 1000 or greater. This #
C allows the user to build and maintain custom #
C heat source or temperature boundary condition #
C functions that are too exotic to be covered by #
C any of the existing QTRAN microfunctions. #
C #
C UOUTPT --> Called from QTRAN subroutine QFLOW, this #
C subroutine may be used to print out customized #
C data immediately after the temperature data is #
C printed out for either transient or steady state #
C runs. UOUTPT data is printed out prior to the #
C resistor data. #
C #
C UPCPAC --> Called from CCAPAC when a user material property has #
C beeb defined for either the capacitor density or #
C specific heat. Even if the user properties are dummies #
C to flag a user defined capacitance determination, the #
C UPROP value must pass a value that will allow a #
C capacitor evaluation prior to the user defining it. #
C This routine is for users to be able to define and #
C control the phase change energy. #
C #
C UPLOT --> Called from CRUNCH after every converged #
C transient or steady state solution. #
C Specific plot information could be created and #
C written to a file here. #
C #
C UPRNTC --> Called from QTRAN subroutine TPSET, this #
C subroutine will inform the users if a print had #
C occured on the previous time step. It can be #
C used by the user to force a print at the current #
C time step. #
C #
C UPROP --> Called from QTRAN subroutine PROPS, this #
C subroutine is called whenever the user specifies #
C a material property evaluation option (IEVAL) of #
C U (User-Coded). This allows the user to build #
C and maintain a library of custom material #
C property subroutines for very exotic materials #
C or applications not covered by any of QTRAN's #
C pre-coded material property evaluation options. #
C #
C URADAT --> Called from RADIAT. Allows the user to modify #
C radiation conductors; however, they must have some #
C way to identify the conductor to change. Either by #
C type, node IDs, material properties or some #
C combination of them. #
C #
C URSTRT --> Called from QTRAN subroutine GETRST if QTRAN #
C is resuming execution from a restart file. #
C #
C USOL --> Called from QTRAN subroutine CRUNCH if the #
C QTRAN solution option SOL = 1000. You may #
C then invoke any solution system you wish. #
C QTRAN's solvers will be bypassed. Execution #
C passes through CRUNCH to RESET1, then USOL, #
C and then SSTATE or TRANS, depending upon whether #
C the problem being run is steady state or #
C transient. #
C #
C UWAVER --> Called from WAVER. Allows the user to modify #
C radiation conductors; however, they must have some #
C way to identify the conductor to change. Either by #
C type, node IDs, material properties or some #
C combination of them. Same as URADAT except this is for #
C wavelength dependent radiation. #
C #
C #
C NOTE: PATRAN CUSTOMERS MAY MODIFY THE ROUTINES IN THIS FILE #
C IN ANY MANNER THEY SEE FIT. HOWEVER, IT IS UP TO THE #
C CUSTOMERS TO MAINTAIN ANY SUCH MODIFICATIONS. PDA #
C ENGINEERING WILL NOT BE RESPONSIBLE FOR THE #
C CONSEQUENCES INCURRED DUE TO THE MODIFICATION OF THESE #
C ROUTINES. #
C #
C############################################################################
C
C############################################################################
C
C############################################################################
C #
C S U B R O U T I N E A F L O W #
C #
C############################################################################
C
SUBROUTINE UAFLOW( RMDOT, ICPFLO, COEFF, EXPO )
C
C############################################################################
C
C This user subroutine can be used by the user to compute the
C mass flow resistor heat flows.
C
C nrf is the heat flow resistor i.d. number.
C rmdot is the vector that stores the mass flow rate data.
C icpflo is the vector that stores the specific heat material property
C i.d. number for the flow resistor.
C coeff and expo are the material property arrays for subroutine props.
C t1 is the temperature of node #1 of the flow resistor.
C t2 is the temperature of node #2 of the flow resistor.
C q12 is the heat flow from node #1 to node #2 of the resistor.
C q21 is the heat flow from node #2 to node #1 of the resistor.
C since "upwind" differencing is used, one of the two heat
C flows will always be zero (the upstream node will never
C have heat flowing to it from the downstream node.)
C rmflow is the computed mass flow rate.
C
C
C##############################################################################
C
C > > > > > MODIFICATION HISTORY < < < < <
C
C Modification: Calculate the energy transferred as the difference in
C the energy state at the two end points rather than the average
C movement between the two points.
C By: Haddock Date: 5 February 1991
C
C Modification: Energy state at the end points must be integrated
C from the base reference state. If we are only interested in
C energy movement, then we can integrate between the two points.
C By: Haddock Date: 13 February 1992
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
double_precision T1
double_precision T2
double_precision CP
double_precision Q12
double_precision Q21
double_precision GVAL
double_precision RMFLOW
C
INTEGER NRF
C
COMMON / IB39 / NRF
C
COMMON / RB36 / T1
COMMON / RB37 / T2
COMMON / RB44 / CP
COMMON / RB46 / Q12
COMMON / RB47 / Q21
COMMON / RB48 / GVAL
COMMON / RB62 / RMFLOW
C
C############################################################################
C define call list variables
C
INTEGER ICPFLO( J30, 2 )
C
double_precision RMDOT(*), COEFF( J3, J4 ), EXPO( J3, J4 )
C
C############################################################################
C Declare the local variables.
C
double_precision RM1, RM2, DIFFT
double_precision CPEVAL
double_precision PROPS
C
C############################################################################
C Declare external functions
C
EXTERNAL CPEVAL
EXTERNAL PROPS
C
C############################################################################
C
C Check the flow direction, and then calculate the appropriate heat
C flow from the upstream node to the downstream node as
C follows:
C (1) compute the specific heat value CP
C (2) compute the conductance GVAL = mass flow rate * CP
C (3) compute the heat flow from the upstream node to the
C downstream node, and then set the heat flow from the
C downstream node to the upstream node to zero
C
C############################################################################
C
!C CP = CPEVAL( T1, T2, COEFF, EXPO, ICPFLO(NRF,1) )
C
!C RMFLOW = RMDOT(NRF)
C
!C IF(ICPFLO(NRF,2).NE.0) THEN
C
!C RM1 = RMFLOW * PROPS( COEFF, EXPO, T1, ICPFLO(NRF,2) )
!C RM2 = RMFLOW * PROPS( COEFF, EXPO, T2, ICPFLO(NRF,2) )
!C RMFLOW = ( RM1 + RM2 ) * 0.5D+00
!C ENDIF
C
!C GVAL = ABS( RMFLOW * CP )
C
!C DIFFT = ( T1 - T2 )
C
!C IF(RMFLOW.GE.0.D+00) THEN
!C Q12 = GVAL * DIFFT
!C Q21 = 0.0D+00
!C ELSE
!C Q21 = GVAL * (-DIFFT)
!C Q12 = 0.0D+00
!C ENDIF
C
C############################################################################
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U C C P A C #
C #
C############################################################################
C
SUBROUTINE UCCPAC( CT, CSUM, C, CP, CM,
1 ICPNT, CNPNT, CPROP, CRHO, CVOL,
2 ICP, IRHO, COEFF, EXPO, IPERTR,
3 OT, TNODE, TP, TM,
4 T, TPLUSP, TMNUSP,
5 TDIFF, TDIFFP, TDIFFM )
C
C############################################################################
C #
C This subroutine performs all capacitance when a user defined #
C property has be defined. #
C #
C############################################################################
C
C
C declare the subroutine arguments, where:
c ct --> thermal capacity of this capacitor.
c csum --> sum of all capacitances to this node.
c c --> capacitance energy at new temperature.
c cp --> capacitance energy at plus perturbed temperature.
c cm --> capacitance energy at minu sperturbed temperature.
c icpnt --> index to capacitor in cnpnt array.
c cnpnt --> array containing the capacitor i.d. numbers for each
c node. all zeros are packed out.
c cprop --> array of specific heat mid numbers for each capacitor.
c crho --> array of density mid numbers for each capacitor.
c cvol --> array containing capacitor volumes.
c icp --> specific heat material property type.
c irho --> density material property type.
c coeff --> material property data array.
c expo --> material property data array.
c ipertr --> perturbation flag.
c ot --> old temperature value.
c tnode --> weighted new temperature.
c tp --> weighted new temperature and plus perturbation.
c tm --> weighted new temperature and minus perturbation.
c t --> new temperature value.
c tplusp --> new temperature value plus perturb.
c tmnusp --> new temperature value minus perturb.
c tdiff --> temperature change.
c tdiffp --> temperature change plus perturb.
c tdiffm --> temperature change minus perturb.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Declare the arrays.
C
DOUBLE PRECISION CP
DOUBLE PRECISION TNODE
DOUBLE PRECISION COEFF( J3, J4 ), EXPO( J3, J4 )
DOUBLE PRECISION CVOL(*)
C
INTEGER IPERTR
INTEGER ICPNT
INTEGER CNPNT(*), CPROP(*), CRHO(*)
C
C############################################################################
C
C Some common blocks
C
C*C*RLH INTEGER ICMPID
C*C*RLH INTEGER IRMPID
C*C*RLH INTEGER MPID
C*C*RLH INTEGER MID
C
C*C*RLH COMMON / IA45 / MID(1)
C
C
C############################################################################
C
C Declare the local variables.
C
DOUBLE PRECISION CSUM, C, TP, TM, TDIFF
DOUBLE PRECISION TDIFFP, TDIFFM, CT, OT, T
DOUBLE PRECISION CM
DOUBLE PRECISION TPLUSP, TMNUSP
DOUBLE PRECISION CPEVAL, PROPS
DOUBLE PRECISION LCP
DOUBLE PRECISION LRHO
C
INTEGER ICP, IRHO
C
C*C*RLH DOUBLE PRECISION LOCT
C
C############################################################################
C Declare external functions
C
C*C*RLH EXTERNAL CPEVAL, PROPS
C
C#######################################################################
C
C
C-----------------------------------------------------------------------
C
! Note: If a user defined heat capacity or density has been defined
! as a user defined material property, then the uprops routine has
! been called and the heat capacity was determined with the
! following procedure. If the uprop for the given uprop is a
! dummy only to flag a user defined heat capacity, then the user
! can recalculate the heat capacity or mark appropriate modifications.
C
C#######################################################################
C
C this section is for non-phase change behavior. compute
C the nodal capacitance energy, where:
C c = cp * rho * vol * dT. then do the same thing
C for the perturbed temperature values. note that
C the C, CP, and CM variables will store the sum of
C all capacitance heat flows for the node. the
C variable csum will store the sum of all nodal
C capacitances, and will be used for computing
C explicit stable time steps.
C
C*C*RLH ICMPID = MID( ICP )
C*C*RLH IRMPID = MID( IRHO )
C
C*C*RLH IF( ICMPID .EQ. 303105 .OR. ICMPID .EQ. 370705 .OR.
C*C*RLH 1 IRMPID .EQ. 303104 .OR. IRMPID .EQ. 370704 ) THEN
C
C*C*RLH IF( IPERTR .EQ. 0 ) THEN
C
C............................................................................
! Evaluate at temperature T
!
C*C*RLH LOCT = T
C
C*C*RLH ELSE IF( IPERTR .EQ. 1 ) THEN
C
C............................................................................
! Evaluate at temperature plus perturbed temperature.
!
C*C*RLH LOCT = TPLUSP
C
C*C*RLH ELSE IF( IPERTR .LE. -1 ) THEN
C
C............................................................................
! Evaluate at temperature minus perturbed temperature.
C
C*C*RLH LOCT = TMNUSP
C
C*C*RLH END IF
C
C*C*RLH IF( ICMPID .EQ. 303105 .OR. ICMPID .EQ. 370705 ) THEN
C
C*C*RLH LCP = EXPO( ICP, 2 )
C
C*C*RLH ELSE
C
C*C*RLH LCP = CPEVAL( OT, LOCT, COEFF, EXPO, ICP )
C
C*C*RLH END IF
C
C*C*RLH IF( IRMPID .EQ. 303104 .OR. IRMPID .EQ. 370704 ) THEN
C
C*C*RLH LRHO = EXPO( IRHO, 2 )
C
C*C*RLH ELSE
C
C*C*RLH LRHO = PROPS( COEFF, EXPO, TNODE, IRHO )
C
C*C*RLH END IF
C
C*C*RLH IF( IPERTR .EQ. 0 ) THEN
C
C............................................................................
! Evaluate at temperature T
!
C*C*RLH CT = LCP * CVOL(ICPNT) * LRHO
C
C*C*RLH ELSE IF( IPERTR .EQ. 1 ) THEN
C
C............................................................................
! Evaluate at temperature plus perturbed temperature.
!
C*C*RLH CT = LCP * CVOL(ICPNT) * LRHO
C
C*C*RLH ELSE IF( IPERTR .LE. -1 ) THEN
C
C............................................................................
! Evaluate at temperature minus perturbed temperature.
C
C*C*RLH CT = LCP * CVOL(ICPNT) * LRHO
C
C*C*RLH END IF
C
C*C*RLH END IF
C
C############################################################################
C
RETURN
END
C
C , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
C
C compute q conducted from node 1 to node 2 (qc).
C
C*C*RLH QC = GVAL * ( T1 - T2 )
C
C*C*RLH END IF
C
C Other options can be defined for different sections of the code.
C
C****************************************************************************
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U E X I T Q #
C #
C############################################################################
C
SUBROUTINE UEXITQ
C
C............................................................................
C
C UEXITQ --> Called whenever EXITQ is called. This is just
C before the standard output routine qtran.dat is closed.
C Is called during error conditions as well as normal
C exits. The UFNSHD routine will be call just before
C this but only for normal exits.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C............................................................................
C
RETURN
END
C
C
C############################################################################
C #
C S U B R O U T I N E U F N S H D #
C #
C############################################################################
C
SUBROUTINE UFNSHD
C
C............................................................................
C
C UFNSHD --> Called at the completion of a solution cycle when
C crunch is going to be exited.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C............................................................................
C
RETURN
END
C
C
C
C############################################################################
C #
C PATRAN Customers May Modify This Routine In Any Manner They See Fit. #
C However, It Is The Customer's Responsibility To Maintain Any Such #
C Modifications. PDA Engineering Will Not Be Responsible For The #
C Consequences Incurred Due To The Modification Of This Routine. #
C #
C############################################################################
C #
C S U B R O U T I N E U H V A L #
C #
C############################################################################
C
SUBROUTINE UHVAL( ICFIG, IRESIS, COEFF, EXPO, JPROP, GP, T1, T2,
$ GVALH, Q, LOGP, J1, J2, J3, J4, J6 )
C
C############################################################################
C
C This subroutine is meant to be used as an example routine for those
C hardy souls who feel the urge to write their own User-Coded
C convection configuration library.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C ARGUMENTS:
C
C icfig --> configuration type of resistor
C iresis --> convective resistor i.d. number
C coeff --> material property data array
C expo --> material property data array
C jprop --> list of MPID numbers assigned to convective resistors
C gp --> list of geometric property data assigned to
C convective resistors
C t1 --> temperature of node #1 of the convective resistor
C (in degrees ICCALC)
C t2 --> temperature of node #2 of the convective resistor
C (also in degrees ICCALC)
C gvalh --> "conductance" of the resistor ( product of h * Area )
C q --> heat flow from node 1 to node 2
C logp --> resistor data print flag. if .true., resistor data
C is requested to be printed. if .false., no request
C has been made.
C j1-j6 --> array dimensions
C
C############################################################################
C
C Declare the arrays.
C
INTEGER J1, J2, J3, J4, J6
C
INTEGER JPROP(J1,J2)
C
DOUBLE PRECISION COEFF( J3, J4 ), EXPO( J3, J4 ), GP( J1, J6 )
C
C############################################################################
C
C Declare the other arguments.
C
LOGICAL LOGP
C
INTEGER ICFIG, IRESIS
C
DOUBLE PRECISION T1, T2, GVALH, Q
C
C############################################################################
C
C Declare the logical unit number and its common block for the output
C file.
C
INTEGER IO
C
COMMON / IB33 / IO
C
C############################################################################
C
C Declare the common block for the film coefficient H. You must be sure
C to calculate the value of the convective film coefficient and store it
C in the common block variable H. This value will then be used when
C computing the average nodal H value which is placed in the nodal
C results files and which can be used for post-processing.
C
DOUBLE PRECISION H
C
COMMON / RB69 / H
C
C############################################################################
C #
C The preceeding is standard boiler plate for subroutine UHVAL and #
C should ALWAYS be included exactly as is. #
C #
C############################################################################
C #
C The following shows how to set up 3 seperate configurations, and is #
C entirely optional. The following guidelines do apply, #
C however: #
C (1) ALWAYS return GVALH as the product of h * Area. #
C Failure to do this will result in driving the #
C SNPSOR equation solver off the deep end. #
C (2) ALWAYS return Q as the heat transfer from node 1 #
C to node 2, i.e., Q = h * Area * ( T1 - T2 ) #
C (3) ALWAYS calculate and store the film coefficient #
C in the common block variable H. #
C #
C############################################################################
C
C Declare the local variables.
C
C*C INTEGER CONFIG
C
C*C DOUBLE PRECISION AREA, RHO, MU, K, CP, EXPAN, TBAR, TFILM,
C*C $ DTEMP, G, L, PR, RA, HLOW, RAL, HHIGH
C*C DOUBLE PRECISION PROPS
C*C
C*C EXTERNAL PROPS
C
C############################################################################
C
C Set up the logic to choose between user-supplied configurations
C 1000, 1001, and 1002. The computed go-to's are frequently
C significantly faster than nested if-then-else structures
C for large lists. The computed go-to structure below may
C easily be expanded to handle many more configurations. The
C logic assumes that you have decided to use user-supplied
C configurations 1000, 1001, and 1002. The first thing done
C is to subtract 999 from the ICFIG configuration value. Thus
C if ICFIG = 1000, CONFIG will be set to 1. The go-to will then
C branch to stmt number 1000. The other configurations are
C handled identically, as is obvious from the following coding.
C
C Note that while all coding is currently in this routine, it is easy
C (and better practice) to modularize each configuration into
C its own subroutine, e.g., UHVAL1, UHVAL2, etc., which would
C then by called by this master routine UHVAL. UHVAL would then
C contain nothing but the selection logical (the computed go-to)
C and calls to your list of routines.
C
C*C CONFIG = ICFIG - 999
C
C*C GO TO (1000,1001,1002),CONFIG
C
C############################################################################
C
C Set up the target statement label for configuration 1000.
C
C*C 1000 CONTINUE
C
C This configuration is going to be identical to QTRAN library
C configuration number 31, "Generic Constant H Value."
C
C Get the resistor surface area. This we have decided to store as
C Geometric Property number 1. Geometric properties for this
C resistor (number IRESIS) will be stored in the IRESIS'th row
C of the GP array. GP #1 will be in GP(IRESIS,1), GP #2 will
C be in GP(IRESIS,2), etc.
C
C*C AREA = GP(IRESIS,1)
C
C............................................................................
C
C Get the resistor h value. This we have decided to store as
C Geometric Property number 2.
C
C*C H = GP(IRESIS,2)
C
C............................................................................
C
C Compute the conductance. This is a requirement and its value will
C be returned to the calling routine.
C
C*C GVALH = H * AREA
C
C............................................................................
C
C Compute the heat transfer from node 1 to node 2. This is a requirement
C and its value will be returned to the calling routine.
C
C*C Q = GVALH * ( T1 - T2 )
C
C............................................................................
C
C Check to see if resistor data is being printed out. If not, return.
C If so, write out the h value and area.
C
C*C IF(.NOT.LOGP) RETURN
C
C............................................................................
C
C The following two QTRAN subroutines print out the h value and
C resistor area very easily.
C
C*C CALL HOUT( H )
C*C CALL AREAO( AREA )
C
C............................................................................
C
C Return to the calling routine.
C
C*C RETURN
C
C############################################################################
C
C Set up the target statement label for configuration 1001.
C
C*C 1001 CONTINUE
C
C This configuration will be identical to configuration number 29,
C "Generic Variable H value, H = H(TBAR)." For this property,
C it is assumed that the h value is to be taken directly from
C the material property specified by JPROP #1.
C
C First, get the TBAR value.
C
C*C TBAR = ( T1 + T2 ) / 2.0D+00
C
C............................................................................
C
C Second, get the JPROP #1 number. This will be stored in the JPROP
C array as JPROP(IRESIS,1). JPROP #2 would be stored in
C JPROP(IRESIS,2), and so on. Normally, MPID numbers are stored
C in the JPROP array. We will therefore call JPROP #1 by the
C name MPID.
C
C*C MPID = JPROP(IRESIS,1)
C
C............................................................................
C
C Third, get the h value. This is done by a call to QTRAN subroutine
C PROPS. PROPS is used to evaluate all of QTRAN's material
C property values. NOTE: If MPID is negative, subroutine PROPS
C will ignore it and will use the current TIME value as the
C evaluation argument. Since you are doing the coding,
C you may of course use any argument that you wish.
C
C*C H = PROPS( COEFF, EXPO, TBAR, MPID )
C
C............................................................................
C
C Fourth, get the surface area for the resistor. This is assumed to
C have been entered as GP #1, although you of course may enter
C it as any GP value that you wish (or get it from anywhere
C else that you are able to).
C
C*C AREA = GP(IRESIS,1)
C
C............................................................................
C
C Compute the conductance ( h * Area ).
C
C*C GVALH = H * AREA
C
C............................................................................
C
C Compute the heat transferred from node 1 to node 2.
C
C*C Q = GVALH * ( T1 - T2 )
C
C............................................................................
C
C Check to see if resistor data is being printed out. If not, return.
C If so, write out the h value and area.
C
C*C IF(.NOT.LOGP) RETURN
C
C............................................................................
C
C The following two QTRAN subroutines print out the h value and
C resistor area very easily.
C
C*C CALL HOUT( H )
C*C CALL AREAO( AREA )
C
C............................................................................
C
C Return to the calling routine.
C
C*C RETURN
C
C############################################################################
C
C Set up the target statement label for configuration 1002.
C
C*C 1002 CONTINUE
C
C Now that we're done with the warm-up exercises, let's do a semi-serious
C convection configuration. Let's assume that:
C
C h = GP#5 { Rayleigh < 1.5E+05 }C
C h = (GP#2) * Rayleigh**(GP#3) * Prandtl ** (GP#4) { otherwise }C
C First, compute the necessary material properties.
C
C Assume: density = jprop #1
C viscosity = jprop #2
C coeffecient of expansion = jprop #3
C specific heat = jprop #4
C conductivity = jprop #5
C
C Compute the film temperature.
C
C*C TFILM = ( T1 + T2 ) * 0.5D+00
C
C Calculate the fluid properties. Since negative values for
C material properties would likely make our calculations
C blow up (as would a zero value), we will protect ourselves
C just a little bit and take the absolute value of the material
C properties and then fudge them by 1.0D-10 (the D exponent
C implies double precision).
C
C*C AREA = GP( IRESIS, 1 )
C
C*C RHO = DABS( PROPS( COEFF, EXPO, TFILM, JPROP(IRESIS,1)) ) +
C*C $ 1.D-10
C*C MU = DABS( PROPS( COEFF, EXPO, TFILM, JPROP(IRESIS,2)) ) +
C*C $ 1.D-10
C*C EXPAN = DABS( PROPS( COEFF, EXPO, TFILM, JPROP(IRESIS,3)) ) +
C*C $ 1.D-10
C*C CP = DABS( PROPS( COEFF, EXPO, TFILM, JPROP(IRESIS,4)) ) +
C*C $ 1.D-10
C*C K = DABS( PROPS( COEFF, EXPO, TFILM, JPROP(IRESIS,5)) ) +
C*C $ 1.D-10
C
C............................................................................
C
C Calculate the Prandtl number.
C
C*C PR = MU * CP / K
C
C............................................................................
C
C Calculate the Rayleigh number.
C
C First, get the gravitational constant entered as GP #6. Then get
C the characteristic length L used by the Rayleigh number.
C GP #7 will be the shortest distance between the resistor area
C and the surface edge whose boundary layer thickness is zero,
C and GP #8 will be the longest distance between the resistor
C area and the same edge. Note that adding GP #7 and GP #8 is
C one way of placing the resistor surface at the center of an
C artificially constructed effective characteristic length.
C This assumes that h(L) varies linearly from h(0) to h(L),
C and is a first order approximation only. Although not a
C particularly tremendous approximation, it does limit correctly
C when GP #7 = 0 and GP #8 = L and it is a whole lot better than
C simply assuming that h is constant across the whole convective
C surface. There are better ways, but this way is especially
C easy.
C
C*C G = GP(IRESIS,6)
C*C L = GP(IRESIS,7) + GP(IRESIS,8)
C
C Compute the temperature difference between the surface and the free
C stream fluid.
C
C*C DTEMP = DABS( T1 - T2 ) + 1.0D-10
C
C*C RA = DABS( G * EXPAN * ( RHO * RHO ) * ( L * L * L ) * DTEMP /
C*C $ ( MU * MU ) * PR ) + 1.0D-10
C
C............................................................................
C
C We now have a Rayleigh number. Check to see if we are in range of
C the first correlation or the second.
C NOTE: Most iterative codes (QTRAN included) do not appreciate step
C function discontinuities such as would be generated by a
C transition NUMBER. We will therefore use an arbitrary
C transition RANGE. What this means is that instead of using
C the transition number of 1.5E+05, we will smear the transition
C arbitrarily over the range of 1.3E+05 to 1.7E+05. Since step
C function discontinuities do not really occur anyway (laminar
C to turbulent transitions really do occur over a range) there is
C even some physical justification for this.
C
C First, calculate both the low range and the high range h values.
C
C*C HLOW = GP(IRESIS,5)
C
C , , , , , , , , , , , , , , , , , , , , , , , , ,
C
C Limit the RA value to the range of the high correlation. This is done
C for two reasons. One, it limits the correlation to a valid
C Rayleigh number range. Two, it provides us with an
C interpolation point for the transition range. If the lower
C correlation were not a constant, we would have performed the
C same type of limiting operation with it.
C
C*C RAL = DMAX1( RA, 1.7D+05)
C
C*C HHIGH = GP(IRESIS,2) * RAL ** GP(IRESIS,3) * PR ** GP(IRESIS,4)
C
C , , , , , , , , , , , , , , , , , , , , , , , , ,
C
C Next, check to see if we are in the low range, high range, or
C transition regime. If we are in the transition regime,
C we will interpolate. This can all be done transparently by
C simply calling QTRAN subroutine INTERP. The calling sequence
C for INTERP is:
C
C CALL INTERP( low_transition_range, low_h_value,
C $ high_transition_range, high_h_value,
C $ Rayleigh (or whatever), returned_h_value )
C
C INTERP does more than just interpolate. If the Rayleigh
C number is below the low_transition_range, INTERP will return
C the low_h_value. If above the high_transition_range, INTERP
C will return the high_h_value. INTERP can obviously be used
C in the same manner for Reynolds number, Prandtl number, and
C other transition range calculations.
C
C*C CALL INTERP( 1.3D+05, HLOW, 1.7D+05, HHIGH, RA, H )
C
C............................................................................
C
C Compute the conductance.
C
C*C GVALH = H * GP(IRESIS,1)
C
C............................................................................
C
C Compute the heat transferred from node #1 to node #2.
C
C*C Q = GVALH * ( T1 - T2 )
C
C............................................................................
C
C Check to see if resistor data is being printed out. If not, return.
C If so, write out the h value and area.
C
C*C IF( .NOT.LOGP ) RETURN
C
C............................................................................
C
C The two QTRAN subroutines HOUT and AREAO print out the h value and
C resistor area very easily. You may also write out any other
C data that you wish, such as Rayleigh number, Prandtl number,
C etc.
C
C*C CALL HOUT( H )
C
C*C WRITE(IO,1) MU, CP, K, PR
C*C 1 FORMAT(' Viscosity:',T20,1PD20.10/' Specific Heat:',T20,D20.10/C*C $' Conductivity:',T20,D20.10 / 'Prandtl No.:', T20, D20.10 )
C
C*C WRITE(IO,2) RHO, TFILM, G, EXPAN, RA
C*C 2 FORMAT(' Density:',T20,1PD20.10/' Film Temp:',T20,D20.10/C*C $' Grav Constant:',T20,D20.10/' Coeff Expan:',T20,D20.10 /
C*C $' Rayleigh:', T20, D20.10 )
C
C*C CALL AREAO( AREA )
C
C............................................................................
C
C Return to the calling routine.
C
RETURN
C
C############################################################################
C
END
C
C############################################################################
C #
C S U B R O U T I N E U I N I T 1 #
C #
C############################################################################
C
SUBROUTINE UINIT1
C
C............................................................................
C
C UINIT1 --> Called from QTRAN subroutine INIT1, this
C subroutine may be used to initialize any system
C arrays provided by the user. UINIT1 is called
C prior to reading in any input data.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C
C............................................................................
C
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U I N I T 2 #
C #
C############################################################################
C
SUBROUTINE UINIT2
C
C............................................................................
C
C UINIT2 --> Called from QTRAN subroutine INIT2, this
C subroutine may be used to set up pointer tables
C or anything else that may need to be done after
C reading in the input data file(s).
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Load the dimension definitions - ( from common.dims )
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C............................................................................
C
C
C............................................................................
C
C
RETURN
END
C############################################################################
C #
C This file contains dummy subroutines that can be altered by the user #
C and included with the QTRAN thermal analysis module. These #
C routines are called at key entry points in the calculation #
C loop as well as during initialization phases, reading of the #
C input data file, etc. The routines included are: #
C #
C ROUTINES: #
C #
C UINPUT --> Called from QTRAN subroutine INPUT, this #
C subroutine may be used to read in customized #
C data not normally found in QTRAN's input data #
C file. It is called after QTRAN's normal input #
C data has been read in. #
C #
C UOUTPT --> Called from QTRAN subroutine QFLOW, this #
C subroutine may be used to print out customized #
C data immediately after the temperature data is #
C printed out for either transient or steady state #
C runs. UOUTPT data is printed out prior to the #
C resistor data. #
C #
C UINIT1 --> Called from QTRAN subroutine INIT1, this #
C subroutine may be used to initialize any system #
C arrays provided by the user. UINIT1 is called #
C prior to reading in any input data. #
C #
C UINIT2 --> Called from QTRAN subroutine INIT2, this #
C subroutine may be used to set up pointer tables #
C or anything else that may need to be done after #
C reading in the input data file(s). #
C #
C ULOOP1 --> Called from QTRAN subroutine RESET1, this #
C subroutine may be used to perform any necessary #
C housekeeping or calculations prior to beginning #
C a steady state run or a transient time step. #
C #
C ULOOP2 --> Called from QTRAN subroutine RESET2, this #
C subroutine may be used to perform any necessary #
C housekeeping or calculations prior to beginning #
C an iteration for either steady state or #
C transient runs. #
C #
C ULOOP3 --> Called from QTRAN subroutine CRUNCH, this #
C subroutine is called at the beginning of #
C calculations for each node for each iteration. #
C Custom resistor types are one example of things #
C that might be appropriate for ULOOP3. #
C #
C ULOOP4 --> Called from QTRAN subroutine CRUNCH, this #
C subroutine is called at the end of calculations #
C for each node for each iteration immediately #
C after the temperature for that node has been #
C updated. #
C #
C ULOOP5 --> Called from QTRAN subroutine CRUNCH, this #
C subroutine is called at the end of each #
C steady state or transient iteration. This #
C routine can be used to perform auxiliary #
C calculations that must be performed in parallel #
C with the thermal calculations, e.g., fluid flow #
C calculations, mass transport calculations, or #
C anything else that is appropriate. #
C #
C ULOOP6 --> Called from QTRAN subroutine CRUNCH, this #
C subroutine is called only after a steady state #
C calculation has converged. #
C #
C ULOOP7 --> Called from QTRAN subroutine CRUNCH, this #
C subroutine is called after each transient time #
C step. #
C #
C UHVAL --> Called from QTRAN subroutine CONV0, this #
C subroutine is called whenever the user specifies #
C a convection configuration type of 1000 or #
C greater. This allows the user to build and #
C maintain custom convection configurations. #
C #
C UMICRO --> Called from QTRAN subroutine FLIB, this #
C subroutine is called whenever the user specifies #
C a microfunction option of 1000 or greater. This #
C allows the user to build and maintain custom #
C heat source or temperature boundary condition #
C functions that are too exotic to be covered by #
C any of the existing QTRAN microfunctions. #
C #
C UPROP --> Called from QTRAN subroutine PROPS, this #
C subroutine is called whenever the user specifies #
C a material property evaluation option (IEVAL) of #
C U (User-Coded). This allows the user to build #
C and maintain a library of custom material #
C property subroutines for very exotic materials #
C or applications not covered by any of QTRAN's #
C pre-coded material property evaluation options. #
C #
C URSTRT --> Called from QTRAN subroutine GETRST if QTRAN #
C is resuming execution from a restart file. #
C #
C USOL --> Called from QTRAN subroutine CRUNCH if the #
C QTRAN solution option SOL = 1000. You may #
C then invoke any solution system you wish. #
C QTRAN's solvers will be bypassed. Execution #
C passes through CRUNCH to RESET1, then USOL, #
C and then SSTATE or TRANS, depending upon whether #
C the problem being run is steady state or #
C transient. #
C *
C UPLOT --> Called from CRUNCH after every converged *
C transient or steady state solution. *
C Specific plot information could be created and *
C written to a file here. *
C #
C NOTE: PATRAN CUSTOMERS MAY MODIFY THE ROUTINES IN THIS FILE #
C IN ANY MANNER THEY SEE FIT. HOWEVER, IT IS UP TO THE #
C CUSTOMERS TO MAINTAIN ANY SUCH MODIFICATIONS. PDA #
C ENGINEERING WILL NOT BE RESPONSIBLE FOR THE #
C CONSEQUENCES INCURRED DUE TO THE MODIFICATION OF THESE #
C ROUTINES. #
C #
C############################################################################
C############################################################################
C############################################################################
C #
C S U B R O U T I N E U I N P U T #
C #
C############################################################################
C
SUBROUTINE UINPUT( IINU, IOU, LNUM )
C
C............................................................................
C
C UINPUT --> Called from QTRAN subroutine INPUT, this
C subroutine may be used to read in customized
C data not normally found in QTRAN's input data
C file. It is called after QTRAN's normal input
C data has been read in.
C
C ARGUMENTS:
C IINU --> logical unit number for QTRAN's input data file.
C IOU --> logical unit number for QTRAN's output data file.
C LNUM --> line number of QTRAN's input data file when
C this routine was called. It is suggested that
C the user keep track of line numbers in the input
C data file so that when an error is encountered
C you can write a diagnostic message to the output
C and/or status files informing the user which line
C number the error occurred on.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
INTEGER IINU, IOU, LNUM
C
C............................................................................
C
C Any read statements that you wish to use may be included in this
C section of the routine in ANSI standard Fortran.
C
C You will have to set up your own common blocks to store this data in,
C since no data can be passed back to QTRAN or any of your
C user-supplied routines through the argument list of this
C subroutine. For example, suppose that you wanted to read data
C into a group of arrays called UA1, UA2, and UA3. The method
C of doing this is as follows:
C
C {1} Set up common blocks for each array.C {2} Read in the data into each of these arrays.C {3} Return to the calling routine.C
C Example code for this is as follows. Note that the code is
C commented out. Also, assume for the moment that each array
C is 100 elements long.
C
C............................................................................
C
C Declare the common blocks and arrays.
C
C INTEGER MAXBLK
C
C PARAMETER ( MAXBLK = 100 )
C
C DOUBLE PRECISION UA1, UA2, UA3
C
C COMMON / UB1 / UA1( MAXBLK )
C COMMON / UB2 / UA2( MAXBLK )
C COMMON / UB3 / UA3( MAXBLK )
C
C INTEGER I
C
C , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
C
C Read in MAXBLK lines of UA1, UA2, and UA3 data.
C
C DO 1 I=1,MAXBLK
C READ( IINU, 2, ERR=900, END=999 ) UA1(I), UA2(I), UA3(I)
C 2 FORMAT( 3G20.10 )
C LNUM = LNUM + 1
C 1 CONTINUE
C
C Return to the calling routine.
C
C RETURN
C
C , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
C
C 900 CONTINUE
C
C This section would be for when a format error was encountered.
C Write out an error message and quit.
C
C WRITE( IOU, 901 ) LNUM
C
C 901 FORMAT(//' ***>>> ERROR <<<***'//
C $' An error has occurred while reading in data for User-Supplied'/
C $' subroutine UINPUT. The error occurred on line number',I10,'.'/
C $/' ***>>> EXECUTION TERMINATING <<<***'//)
C
C CALL EXITQ
C
C , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
C
C 999 CONTINUE
C
C This section would be for when an unexpected end-of-file occurred
C while trying to read in the data. Write out an error message
C and quit.
C
C WRITE( IOU, 998 ) LNUM
C WRITE( IOU, 998 ) LNUM
C
C 998 FORMAT(//' ***>>> ERROR <<<***'//
C $' An unexpected END-OF-FILE occurred while reading in data in'/
C $' User-Supplied subroutine UINPUT. The error occurred on line'/
C $' number',I10,'.'//' ***>>> EXECUTION TERMINATING <<<***'//)
C
C CALL EXITQ
C
C............................................................................
C
END
C
C############################################################################
C #
C S U B R O U T I N E U L O O P 1 #
C #
C############################################################################
C
SUBROUTINE ULOOP1( COEFF, EXPO )
C
C............................................................................
C
C ULOOP1 --> Called from QTRAN subroutine RESET1, this
C subroutine may be used to perform any necessary
C housekeeping or calculations prior to beginning
C a steady state run or a transient time step.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Declare the arguments
C
DOUBLE PRECISION COEFF( J3, J4 ), EXPO( J3, J4 )
C
C############################################################################
C
C
C............................................................................
C
C
C............................................................................
C
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U L O O P 2 #
C #
C############################################################################
C
SUBROUTINE ULOOP2( COEFF, EXPO )
C
C............................................................................
C
C ULOOP2 --> Called from QTRAN subroutine RESET2, this
C subroutine may be used to perform any necessary
C housekeeping or calculations prior to beginning
C an iteration for either steady state or
C transient runs.
C
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Declare the arguments
C
DOUBLE PRECISION COEFF( J3, J4 ), EXPO( J3, J4 )
C
C############################################################################
C
C............................................................................
C
C
C............................................................................
C
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U L O O P 3 #
C #
C############################################################################
C
SUBROUTINE ULOOP3( COEFF, EXPO )
C
C............................................................................
C
C ULOOP3 --> Called from QTRAN subroutine CRUNCH, this
C subroutine is called at the beginning of
C calculations for each node for each iteration.
C Custom resistor types are one example of things
C that might be appropriate for ULOOP3.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Declare the arguments
C
DOUBLE PRECISION COEFF( J3, J4 ), EXPO( J3, J4 )
C
C############################################################################
C
C
C............................................................................
C
C
C............................................................................
C
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U L O O P 4 #
C #
C############################################################################
C
SUBROUTINE ULOOP4( COEFF, EXPO )
C
C............................................................................
C
C ULOOP4 --> Called from QTRAN subroutine CRUNCH, this
C subroutine is called at the end of calculations
C for each node for each iteration immediately
C after the temperature for that node has been
C updated.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Declare the arguments
C
DOUBLE PRECISION COEFF( J3, J4 ), EXPO( J3, J4 )
C
C############################################################################
C
C
C............................................................................
C
C
C............................................................................
C
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U L O O P 5 #
C #
C############################################################################
C
SUBROUTINE ULOOP5( COEFF, EXPO )
C
C............................................................................
C
C ULOOP5 --> Called from QTRAN subroutine CRUNCH, this
C subroutine is called at the end of each
C steady state or transient iteration. This
C routine can be used to perform auxiliary
C calculations that must be performed in parallel
C with the thermal calculations, e.g., fluid flow
C calculations, mass transport calculations, or
C anything else that is appropriate.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Declare the arguments
C
DOUBLE PRECISION COEFF( J3, J4 ), EXPO( J3, J4 )
C
C############################################################################
C
C
C............................................................................
C
C
C............................................................................
C
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U L O O P 6 #
C #
C############################################################################
C
SUBROUTINE ULOOP6( COEFF, EXPO )
C
C............................................................................
C
C ULOOP6 --> Called from QTRAN subroutine CRUNCH, this
C subroutine is called only after a steady state
C calculation has converged.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Declare the arguments
C
DOUBLE PRECISION COEFF( J3, J4 ), EXPO( J3, J4 )
C
C############################################################################
C
C
C............................................................................
C
C
C............................................................................
C
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U L O O P 7 #
C #
C############################################################################
C
SUBROUTINE ULOOP7( COEFF, EXPO )
C
C............................................................................
C
C ULOOP7 --> Called from QTRAN subroutine CRUNCH, this
C subroutine is called after each transient time
C step.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Declare the arguments
C
DOUBLE PRECISION COEFF( J3, J4 ), EXPO( J3, J4 )
C
C############################################################################
C
C
C............................................................................
C
C............................................................................
C
RETURN
END
C If greater, adjust TMAX30.
C
C*C TMAX30 = MAX( TEMPS(NODE), TMAX30 )
C
C If the maximum temperature is greater than or equal to 1000,
C set the microfunction value to 723.15.
C
C*C IF( TMAX30 .GE. 1000.0 ) VAL = 723.15
C
C*C ENDIF
C
C############################################################################
C
RETURN
END
C
C
C############################################################################
C #
C S U B R O U T I N E U M C R P R #
C #
C############################################################################
C
SUBROUTINE UMCRPR( X, MICRO, IFUNC, VAL, P, TX, TY, COEFF, EXPO )
C
C############################################################################
C
C this subroutine calculates the microfunction value and returns it
C in argument "val" for microfunction "micro". The evaluation
C option is passed to UMICRO as IFUNC. X is the independent
C variable specified by the calling macro function.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare the arguments.
C
INCLUDE 'common.dims'
C
INTEGER MICRO, IFUNC
C
DOUBLE PRECISION VAL, X
DOUBLE PRECISION P( J7, J8 ), TX( J9, J10 ), TY( J9, J10 )
DOUBLE PRECISION COEFF( J3, J4 ), EXPO( J3, J4 )
C*C DOUBLE PRECISION PROPS
C*C
C*C EXTERNAL PROPS
C
C############################################################################
C
C This routine is the same as UMICRO except the arguments have been
C added that allow the user to access material properties in
C addition to the micro functions. The following is an example of
C a call to properties. The user will have to figure out how to
C acquire the proper MPID's.
C
C PROPS is used to evaluate all of QTRAN's material
C property values. NOTE: If MPID is negative, subroutine PROPS
C will ignore it and will use the current TIME value as the
C evaluation argument. Since you are doing the coding,
C you may of course use any argument that you wish.
C
C*C SOMTHG = PROPS( COEFF, EXPO, TBAR, MPID )
C
C You may fill this area in with any computational algorithm that you
C desire. If you wish access to QTRAN's arrays for some reason,
C I suggest that you look at QTRAN's main module and reference
C these arrays through the common blocks set up in the main
C module. This should prove rather trivial. Note that the node
C number that this microfunction is going to be applied to is
C available as:
C
C Node_number = IALIAS(NODE),
C where Node_number will be the number that you used in the input
C data file. The value of NODE is actually QTRAN's internal
C reference to the packed node numbers stored in IALIAS.
C For example, suppose that the 23rd node number defined was
C node number 10070. Then IALIAS(23) will contain the value
C of 10070. For internal computational purposes for QTRAN,
C this is NODE 23. The value of the NODE being processed when
C this microfunction is called is available through the loop
C parameters common block above.
C
C In a similar manner, the microfunction I.D. number is stored in
C MFID(MICRO).
C
C If you need to use local variables whose values must be saved between
C calls to this routine, I suggest that you place these variables
C in a local common block, e.g., suppose that you need an array
C for storing maximum and minimum temperature values. The
C following common block could be used (without the "C" in
C column 1, of course):
C
C COMMON / LOCAL1 / TEMPMX(maxt), TEMPMN(maxt)
C
C where "maxt" is available from the "DIMS" common block above
C and happens to be the maximum number of nodes for which the
C problem is currently dimensioned. Exact values of "maxt"
C or other dimensions are available from the main program
C module.
C
C############################################################################
C
C Consider the following hypothetical microfunction application. You
C have a need to apply a heat source to two nodes if their maximum
C temperature exceeds 1000 degrees in whatever units you specified
C for temperature with the ICCALC variable. The nodes that are
C to have heat sources applied are nodes 15 and 30. It is further
C specified that the heat sources for the two nodes have different
C values. The following programming shows how to do this.
C
C CAUTION: If your Fortran compiler does not automatically
C initialize variables to zero, you may need to
C add initialization code for your "common"
C variables in either a BLOCK DATA subroutine or
C in the MAIN program module.
C
C............................................................................
C
C Declare the common blocks for the IALIAS and TEMPS arrays.
C
C Dimensions for following array are defined in QTRAN.FOR
C
C COMMON / IA17 / IFLIST( J17, J16 )
C COMMON / IA19 / ITLIST( J18, J19 )
C
C*C DOUBLE PRECISION TEMPS
C*C
C*C INTEGER IFLIST, ITLIST, IALIAS, NODE, KQMAC, KTMAC
C*C
C*C COMMON / IA17 / IFLIST( 10, 5 )
C*C COMMON / IA19 / ITLIST( 7, 5 )
C*C COMMON / IA25 / IALIAS(1)
C*C COMMON / IB21 / NODE
C*C COMMON / IB73 / KQMAC
C*C COMMON / IB74 / KTMAC
C*C COMMON / RA16 / TEMPS(1)
C
C CAUTION: You should really dimension IALIAS and TEMPS correctly
C instead of using the dummy dimensions of "1". However,
C since this is purely for demonstration purposes....
C Correct dimensions are available in the MAIN module.
C
C............................................................................
C
C Set up a common block to save the maximum temperatures for nodes 15 &
C 30.
C
C*C DOUBLE PRECISION TMAX15, TMAX30
C*C
C*C INTEGER NODEA, NODE1, NODE2
C*C
C*C COMMON / LOCAL1 / TMAX15, TMAX30
C
C............................................................................
C
C Get the node number that the microfunction is applied to.
C
C*C TMAX15 = 15.0
C*C TMAX30 = 30.0
C*C NODNUM = ABS( IALIAS( NODE ) )
C*C IF( KTMAC .GT. 0 ) THEN
C
C TEMPERATURE MICRO FUNCTION IS DEFINED, GET NODES FROM THE IFLIST ARRAY
C
C*C NODEA = IFLIST( KTMAC, 1)
C*C NODE1 = IFLIST( KTMAC, 3)
C*C NODE2 = IFLIST( KTMAC, 4)
C*C ELSE IF( KQMAC .GT. 0 ) THEN
C*C NODEA = ITLIST( KTMAC, 1)
C*C NODE1 = ITLIST( KTMAC, 3)
C*C NODE2 = ITLIST( KTMAC, 4)
C*C END IF
C
C............................................................................
C
C Check for node 15.
C
C*C IF (NODEA .LE. 15 ) THEN
C
C Check the nodes temperature against it's maximum temperature.
C If greater, adjust TMAX15.
C
C*C TMAX15 = MAX( TEMPS (NODE ), TMAX15 )
C
C If the maximum temperature is greater than or equal to 1000,
C set the microfunction value to 100.
C
C*C IF( TMAX15 .GE. 1000.0 ) VAL = 100.0
C
C
C , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
C
C Check for node 30.
C
C*C ELSE IF( NODEA .LE. 30 ) THEN
C
C Check the node's temperature against it's maximum temperature.
C If greater, adjust TMAX30.
C
C*C TMAX30 = MAX( TEMPS(NODE), TMAX30 )
C
C If the maximum temperature is greater than or equal to 1000,
C set the microfunction value to 723.15.
C
C*C IF( TMAX30 .GE. 1000.0 ) VAL = 723.15
C
C*C ENDIF
C
C############################################################################
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U M I C R O #
C #
C############################################################################
C
SUBROUTINE UMICRO( X, MICRO, IFUNC, VAL )
C
C############################################################################
C
C this subroutine calculates the microfunction value and returns it
C in argument "val" for microfunction "micro". The evaluation
C option is passed to UMICRO as IFUNC. X is the independent
C variable specified by the calling macro function.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare the arguments.
C
INTEGER MICRO, IFUNC
C
DOUBLE PRECISION VAL, X
C
C############################################################################
C
C You may fill this area in with any computational algorithm that you
C desire. If you wish access to QTRAN's arrays for some reason,
C I suggest that you look at QTRAN's main module and reference
C these arrays through the common blocks set up in the main
C module. This should prove rather trivial. Note that the node
C number that this microfunction is going to be applied to is
C available as:
C
C Node_number = IALIAS(NODE),
C where Node_number will be the number that you used in the input
C data file. The value of NODE is actually QTRAN's internal
C reference to the packed node numbers stored in IALIAS.
C For example, suppose that the 23rd node number defined was
C node number 10070. Then IALIAS(23) will contain the value
C of 10070. For internal computational purposes for QTRAN,
C this is NODE 23. The value of the NODE being processed when
C this microfunction is called is available through the loop
C parameters common block above.
C
C In a similar manner, the microfunction I.D. number is stored in
C MFID(MICRO).
C
C If you need to use local variables whose values must be saved between
C calls to this routine, I suggest that you place these variables
C in a local common block, e.g., suppose that you need an array
C for storing maximum and minimum temperature values. The
C following common block could be used (without the "C" in
C column 1, of course):
C
C COMMON / LOCAL1 / TEMPMX(maxt), TEMPMN(maxt)
C
C where "maxt" is available from the "DIMS" common block above
C and happens to be the maximum number of nodes for which the
C problem is currently dimensioned. Exact values of "maxt"
C or other dimensions are available from the main program
C module.
C
C############################################################################
C
C Consider the following hypothetical microfunction application. You
C have a need to apply a heat source to two nodes if their maximum
C temperature exceeds 1000 degrees in whatever units you specified
C for temperature with the ICCALC variable. The nodes that are
C to have heat sources applied are nodes 15 and 30. It is further
C specified that the heat sources for the two nodes have different
C values. The following programming shows how to do this.
C
C CAUTION: If your Fortran compiler does not automatically
C initialize variables to zero, you may need to
C add initialization code for your "common"
C variables in either a BLOCK DATA subroutine or
C in the MAIN program module.
C
C............................................................................
C
C Declare the common blocks for the IALIAS and TEMPS arrays.
C
C Dimensions for following array are defined in QTRAN.FOR
C
C COMMON / IA17 / IFLIST( J17, J16 )
C COMMON / IA19 / ITLIST( J18, J19 )
C
C*C DOUBLE PRECISION TEMPS
C*C
C*C INTEGER IFLIST, ITLIST, IALIAS, NODE, KQMAC, KTMAC
C*C
C*C COMMON / IA17 / IFLIST( 10, 5 )
C*C COMMON / IA19 / ITLIST( 7, 5 )
C*C COMMON / IA25 / IALIAS(1)
C*C COMMON / IB21 / NODE
C*C COMMON / IB73 / KQMAC
C*C COMMON / IB74 / KTMAC
C*C COMMON / RA16 / TEMPS(1)
C
C CAUTION: You should really dimension IALIAS and TEMPS correctly
C instead of using the dummy dimensions of "1". However,
C since this is purely for demonstration purposes....
C Correct dimensions are available in the MAIN module.
C
C............................................................................
C
C Set up a common block to save the maximum temperatures for nodes 15 &
C 30.
C
C*C DOUBLE PRECISION TMAX15, TMAX30
C*C
C*C INTEGER NODEA, NODE1, NODE2
C*C
C*C COMMON / LOCAL1 / TMAX15, TMAX30
C
C............................................................................
C
C Get the node number that the microfunction is applied to.
C
C*C TMAX15 = 15.0
C*C TMAX30 = 30.0
C*C NODNUM = ABS( IALIAS( NODE ) )
C*C IF( KTMAC .GT. 0 ) THEN
C
C TEMPERATURE MICRO FUNCTION IS DEFINED, GET NODES FROM THE IFLIST ARRAY
C
C*C NODEA = IFLIST( KTMAC, 1)
C*C NODE1 = IFLIST( KTMAC, 3)
C*C NODE2 = IFLIST( KTMAC, 4)
C*C ELSE IF( KQMAC .GT. 0 ) THEN
C*C NODEA = ITLIST( KTMAC, 1)
C*C NODE1 = ITLIST( KTMAC, 3)
C*C NODE2 = ITLIST( KTMAC, 4)
C*C END IF
C
C............................................................................
C
C Check for node 15.
C
C*C IF (NODEA .LE. 15 ) THEN
C
C Check the nodes temperature against it's maximum temperature.
C If greater, adjust TMAX15.
C
C*C TMAX15 = MAX( TEMPS (NODE ), TMAX15 )
C
C If the maximum temperature is greater than or equal to 1000,
C set the microfunction value to 100.
C
C*C IF( TMAX15 .GE. 1000.0 ) VAL = 100.0
C
C
C , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
C
C Check for node 30.
C
C*C ELSE IF( NODEA .LE. 30 ) THEN
C
C Check the node's temperature against it's maximum temperature.
C If greater, adjust TMAX30.
C
C*C TMAX30 = MAX( TEMPS(NODE), TMAX30 )
C
C If the maximum temperature is greater than or equal to 1000,
C set the microfunction value to 723.15.
C
C*C IF( TMAX30 .GE. 1000.0 ) VAL = 723.15
C
C*C ENDIF
C
C############################################################################
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U O U T P T #
C #
C############################################################################
C
SUBROUTINE UOUTPT( IO )
C
C............................................................................
C
C UOUTPT --> Called from QTRAN subroutine QFLOW, this
C subroutine may be used to print out customized
C data immediately after the temperature data is
C printed out for either transient or steady state
C runs. UOUTPT data is printed out prior to the
C resistor data.
C
C ARGUMENT:
C IO --> logical unit number for QTRAN's output data file.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
INTEGER IO
C
C............................................................................
C
C Any write statements that you wish to use may be included in this
C section of the routine in ANSI standard Fortran.
C
C As with many of the user-supplied routines, the main method of
C communicating between your user-supplied routines is through
C user-specified common blocks. For example, suppose the
C array UA1 contains results data from other user-coded
C subroutines. You will have to have built a common block
C containing UA1 already in your other routines. Assuming
C that this common block is called UB1 (similar to the example
C common block for subroutine UINPUT), the remainder of this
C example subroutine (which is commented out) shows an example
C of how your might go about printing the array data out.
C
C............................................................................
C
C Declare the example common block and array.
C
C INTEGER MAXBLK
C
C PARAMETER ( MAXBLK = 100 )
C
C DOUBLE PRECISION UA1
C
C COMMON / UB1 / UA1( MAXBLK )
C
C INTEGER I
C
C , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
C
C Write out MAXBLK lines of UA1 data.
C
C WRITE( IO, 1 )
C 1 FORMAT( '1USER-SUPPLIED ARRAY UA1 DATA' /
C $ ' ----------------------------' //
C $11X, 'I', 14X, 'UA1(I)' /
C $11X, '-', 14X, '------' )
C
C DO 2 I = 1, MAXBLK
C WRITE( IO, 3 ) I, UA1( I )
C 3 FORMAT( 1X, I10, 1PD20.10 )
C 2 CONTINUE
C
C , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
C
C Return to the calling routine.
C
C RETURN
C
C............................................................................
C
END
C
C############################################################################
C #
C S U B R O U T I N E U P C P A C #
C #
C############################################################################
C
SUBROUTINE UPCPAC( ICPNT, OT, T, TPHASE, QPHASE, COEFF,
$ EXPO, CVOL, ICP, IRHO, C, IPHASE, QPHRHV )
C
C............................................................................
C
C this subroutine performs any necessary phase change energy balance
C calculations.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C > > > > > MODIFICATION HISTORY < < < < <
C
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Declare variable in call list
C
INTEGER IPHASE
INTEGER ICPNT
INTEGER ICP
INTEGER IRHO
C
double_precision TPHASE
double_precision OT
double_precision T
double_precision QPHASE
double_precision C
C
C real arrays.
C
double_precision CVOL(*)
double_precision COEFF(J3,J4)
double_precision EXPO(J3,J4)
double_precision QPHRHV
C
C****************************************************************************
C
C Declare external functions
C
double_precision PROPS, CPEVAL
EXTERNAL PROPS, CPEVAL
C
C############################################################################
C
C local variables.
C
C*C*RLH double_precision TDIFF
C*C*RLH double_precision T12
C*C*RLH double_precision T23
C*C*RLH double_precision Q1
C*C*RLH double_precision Q2
C*C*RLH double_precision Q3
C*C*RLH double_precision TP2
C*C*RLH double_precision TP3
C*C*RLH double_precision FRAC
C*C*RLH double_precision TBAR
C*C*RLH double_precision TPCBM
C*C*RLH double_precision TPCBP
C
C############################################################################
C
C *---------------------------*
C | +|
C | + |
C | + |
C t | + |
C e | +++++++++++++ |
C m | +############# 3 |
C p | + ############# |
C | + #### #### |
C | + #### 2 #### |
C | + #### #### |
C | + 1 ############# |
C | + ############# |
C *---------------------------*
C
C energy
C
C----------------------------------------------------------------------------
CC
C Define the temperature at the beginning and end of the phase change
C
C*C*RLH TPCBM = TPHASE - PCBAND
C*C*RLH TPCBP = TPHASE + PCBAND
C
C check to see if going up or going down the t/e curve.
C if going down (iphase < 0), if going up (iphase >= 0).
C
C*C*RLH IF( IPHASE .LT. 0 ) THEN
C
C............................................................................
C
C section 100 is for phase change going down the t/e curve. it assumes
C that ot is in region 2 or 3, and that t is in region 1 or 2
C of the t/e curve.
C
C begin the calculations for this section by calculating Q3,
C where Q3 is the capacitor stored energy between
C ot and tphase.
C
C*C*RLH TDIFF = TPCBP - OT
C
C*C*RLH IF( TDIFF .GT. 0.D+00 ) THEN
C
C*C*RLH TDIFF = 0.D+00
C*C*RLH Q3 = 0.0D+00
!
C*C*RLH ELSE
C
C*C*RLH T23 = OT + TDIFF * 0.5D+00
C*C*RLH Q3 = TDIFF * CVOL(ICPNT) *
C*C*RLH $ CPEVAL( OT, TPCBP, COEFF, EXPO, ICP ) *
C*C*RLH $ PROPS( COEFF, EXPO, T23, IRHO )
!
C*C*RLH END IF
C
C............................................................................
C
C calculate energy q2, which is the amount of latent heat energy
C being released through "freezing".
C
C*C*RLH TP2 = MIN( TPCBP, OT )
C*C*RLH TP3 = MAX( TPCBM, T )
C*C*RLH FRAC = ( TP3 - TP2 ) / ( PCBAND * 2.D+00 )
C*C*RLH QPHRHV = QPHASE * CVOL(ICPNT) *
C*C*RLH 1 PROPS( COEFF, EXPO, TPHASE, IRHO )
C*C*RLH Q2 = FRAC * QPHRHV
C
C Add the normal capacitance effects to Q2.
C
C*C*RLH TBAR = ( TP2 + TP3 ) * 0.5D+00
C*C*RLH Q2 = Q2 + ( TP3 - TP2 ) * CVOL(ICPNT) *
C*C*RLH $ CPEVAL( TP2, TP3, COEFF, EXPO, ICP ) *
C*C*RLH $ PROPS( COEFF, EXPO, TBAR, IRHO )
C
C............................................................................
C
C calculate the energy q1, which is the amount of energy stored
C between t and tphase.
C
C*C*RLH TDIFF = T - TPCBM
C
C*C*RLH IF( TDIFF .GT. 0.D+00 ) THEN
C
C*C*RLH TDIFF = 0.D+00
C*C*RLH Q1 = 0.0D+00
!
C*C*RLH ELSE
C
C*C*RLH T12 = T - TDIFF * 0.5D+00
C*C*RLH Q1 = TDIFF * CVOL(ICPNT) *
C*C*RLH $ CPEVAL( TPCBM, T, COEFF, EXPO, ICP ) *
C*C*RLH $ PROPS( COEFF, EXPO, T12, IRHO )
C
C*C*RLH END IF
C
C............................................................................
C
C all the potential capacitance effects have now been computed.
C add them to the total nodal capacitance computed so
C far, and then return to the calling routine.
C
C*C*RLH C = C + Q1 + Q2 + Q3
C
C*C*RLH RETURN
C
C============================================================================
C
C*C*RLH ELSE
C
C This section is for phase change going up the t/e curve.
C
C begin the calculations by calculating Q1, where Q1 is the
C capacitor energy stored between ot and tphase.
C
C*C*RLH TDIFF = TPCBM - OT
C
C*C*RLH IF( TDIFF .LT. 0.D+00 ) THEN
C
C*C*RLH TDIFF = 0.D+00
C*C*RLH Q1 = 0.0D+00
!
C*C*RLH ELSE
C
C*C*RLH T12 = OT + TDIFF * 0.5D+00
C*C*RLH Q1 = TDIFF * CVOL(ICPNT) *
C*C*RLH $ CPEVAL( OT, TPCBM, COEFF, EXPO, ICP ) *
C*C*RLH $ PROPS( COEFF, EXPO, T12, IRHO )
!
C*C*RLH END IF
C
C............................................................................
C
C calculate energy q2, which is the amount of latent heat being
C absorbed by "melting".
C
C*C*RLH TP2 = MAX( TPCBM, OT )
C*C*RLH TP3 = MIN( TPCBP, T )
C*C*RLH FRAC = ( TP3 - TP2 ) / ( PCBAND * 2.D+00 )
C*C*RLH QPHRHV = QPHASE * CVOL(ICPNT) *
C*C*RLH 1 PROPS( COEFF, EXPO, TPHASE, IRHO )
C*C*RLH Q2 = FRAC * QPHRHV
C
C add the normal capacitance effects to Q2.
C
C*C*RLH TBAR = ( TP2 + TP3 ) * 0.5D+00
C
C*C*RLH Q2 = Q2 + ( TP3 - TP2 ) * CVOL(ICPNT) *
C*C*RLH $ CPEVAL( TP2, TP3, COEFF, EXPO, ICP ) *
C*C*RLH $ PROPS( COEFF, EXPO, TBAR, IRHO )
C
C............................................................................
C
C calculate energy q23, which is the amount of energy stored
C in the capacitor between t and tphase.
C
C*C*RLH TDIFF = T - TPCBP
C
C*C*RLH IF( TDIFF .LT. 0.D+00 ) THEN
C
C*C*RLH TDIFF = 0.D+00
C*C*RLH Q3 = 0.0D+00
!
C*C*RLH ELSE
C
C*C*RLH T23 = T - TDIFF * 0.5D+00
C*C*RLH Q3 = TDIFF * CVOL(ICPNT) *
C*C*RLH $ CPEVAL( TPCBP, T, COEFF, EXPO, ICP ) *
C*C*RLH $ PROPS( COEFF, EXPO, T23, IRHO )
!
C*C*RLH END IF
C
C............................................................................
C
C all the potential capacitance effects have now been computed.
C add them to the total nodal capacitance computed so far,
C and then return to the calling routine.
C
C*C*RLH C = C + Q1 + Q2 + Q3
C
C*C*RLH RETURN
C
C*C*RLH ENDIF
C
C############################################################################
C
C end of subroutine phasor
C
RETURN
C
END
C
C############################################################################
C
C############################################################################
C #
C S U B R O U T I N E U P L O T #
C #
C############################################################################
C
SUBROUTINE UPLOT( TEMPS )
C
C............................................................................
C
C UPLOT --> Called from OUTPLT subroutine, this
C subroutine is called after each converged solution.
C
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Define dimension statement, included in common block
C
C MAXT -> Maximum number of nodes
C
C............................................................................
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.blk file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Define variables in call list
C
C TEMPS -> Current temperature of nodes
C
DOUBLE PRECISION TEMPS( MAXT )
C
C............................................................................
C
C User operations
C
C
C............................................................................
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U P R N T C #
C #
C############################################################################
C
SUBROUTINE UPRNTC( TEMPS, IALIAS )
C
C############################################################################
C #
C This routine give the user the ability to set a special print flag
C if special conditions are meet.
C
C############################################################################
C
C > > > > > MODIFICATION HISTORY < < < < <
C
C Modification: Add special print conditions
C By: Haddock Date: 22 October 1996
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C declare the arrays.
C
INTEGER FPRTFL, JPRTFL
C
COMMON / IB107 / FPRTFL
COMMON / IB108 / JPRTFL
C
INTEGER IALIAS(*)
C
C real arrays
C
DOUBLE PRECISION TEMPS(*)
C
C############################################################################
C declare the local variables.
C
C If JPRTFL is 1 then a printout occurred on the previous time step.
C To force printout, set FPRTFL to 1.
C
C############################################################################
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U P R O P #
C #
C############################################################################
C
SUBROUTINE UPROP( PROP, X, COEFF, EXPO, IPMPID )
C
C############################################################################
C
C This subroutine can be modified by the user for the express purpose
C of supplying special-coded user-supplied material properties
C that may be too exotic to be handled by QTRAN's normal
C material property library functions. This routine is called
C by QTRAN whenever a material property has been defined with
C an evaluation option of U (IEVAL parameter in Section 5.3.1 of
C the input data file).
C
C Arguments:
C
C PROP --> value of the material property to be returned
C to QTRAN for use in calculations.
C
C X --> Temperature normally used to evaluate
C the material properties, e.g., temperature of
C a node or average temperature of a conductive
C resistor.
C
C COEFF()--> Array of MDATA1 values. For this material
C property, the MATA1(i) values are stored as
C COEFF(IPMPID,i+1). COEFF(IPMPID,1) contains
C the number of MDATA1/MDATA2 data pairs.
C
C EXPO() --> Array of MDATA2 values. For this material
C property, the MATA2(i) values are stored as
C EXPO(IPMPID,i+1). EXPO(IPMPID,1) stores an
C evaluation code, which is used to determine
C whether the material property is a constant,
C a linear data table, etc. (from subroutine
C UPROP, this code always specifies "user-supplied").
C
C IPMPID --> packed material property i.d. This is not the same
C as the MPID number normally entered in Section 5.3.1
C of the input data file. IPMPID here corresponds to
C the IPMPID'th material property that was defined in
C Section 5.3.1. To retrieve the MPID number that you
C gave in Section 5.3.1 for this material property,
C you must reference the MID array in the following
C manner:
C
C COMMON / IA45 / MID(1)
C MPID = MID( IPMPID )
C
C You will then have the correct MPID number for
C the IPMPID'th material property. It should be noted
C that the MPID number is normally not needed for
C anything here, although you might wish to have
C access to it in the event that you set up some
C error message routines.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.blk file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Declare the arguments and local variables.
C
DOUBLE PRECISION PROP, X, COEFF( J3, J4 ), EXPO( J3, J4 )
C
INTEGER IPMPID
C
C*C*RLH INTEGER MPID
C*C*RLH INTEGER MID
C
C*C*RLH DOUBLE PRECISION TIME
C
C*C*RLH COMMON / IA45 / MID(1)
C
C*C*RLH COMMON / RB1 / TIME
C
C############################################################################
C
C The remainder of this subroutine shows an EXAMPLE ONLY of what you
C might do if you wished one or more custom material properties.
C
C Let us suppose that the material property that you wish to define is a
C mix of time and temperature.
C
C############################################################################
C
C This property section shows how to build a material property
C evaluation that is a mix of both temperature and time.
C Suppose that the form of the material property to be
C defined is as follows:
C
C Property(Temp,Time) = { A + B * Temp } * { C + D * Time }C
C The following coding accomplishes this. The assumption is that:
C
C A = 1.0
C B = 2.0
C C = 3.0
C D = -2.0
C
C*C PROP = ( 1.0D+00 + 2.0D+00 * X ) *
C*C $ ( 3.0D+00 - 2.0D+00 * TIME )
C
C NOTE: The value of TIME is passed to this routine via
C a named COMMON block which is inserted at the
C beginning of this routine. The block's name is:
C
C############################################################################
C
! This is an example of how to use user routines to modify resistor values
! User material properties have to be used to flag which resistor or
! capacitors are to be modified by the user.
!
C*C*RLH MPID = MID( IPMPID )
!
C*C*RLH IF( MPID .EQ. 203101 ) THEN
! Silver conductivity MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 203104 ) THEN
! Silver density MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 203105 ) THEN
! Silver specific heat MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 203117 ) THEN
! Silver Emissivity
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 270701 ) THEN
! Flint Glass Conductivity MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 270704 ) THEN
! Flint Glass Density MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 270705 ) THEN
! Flint Glass Specific Heat MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 270717 ) THEN
! Flint Glass Emissivity
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 303101 ) THEN
! Silver conductivity MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 303104 ) THEN
! Silver density MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 303105 ) THEN
! Silver specific heat MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 303117 ) THEN
! Silver Emissivity
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 370701 ) THEN
! Flint Glass Conductivity MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 370704 ) THEN
! Flint Glass Density MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 370705 ) THEN
! Flint Glass Specific Heat MKS units
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH ELSE IF( MPID .EQ. 370717 ) THEN
! Flint Glass Emissivity
C*C*RLH PROP = COEFF( IPMPID, 2 )
C*C*RLH END IF
!
RETURN
C
C############################################################################
C
C As you can see, probably the hardest part of this whole exercise was
C wading through my verbage. Good Luck!
C
C############################################################################
C
END
C
C############################################################################
C #
C S U B R O U T I N E U R A D A T #
C #
C############################################################################
C
SUBROUTINE URADAT( NRAD, TT1, TT2, TT3, T1SAVE, T2SAVE,
1 COEFF, EXPO, GVAL, QGRAY )
C
C############################################################################
C
C This subroutine modifies or calculates the radiant interchange between
C two gray radiosity nodes as the user desires. The user can define
C determine the resistor type and other information necessary by
C defining different common blocks that contain the radiation resistor
C information. This routine is called three times for each resistor
C evaluation - once with current temeratures and then with the plus
C and minus perturbed values. One should keep this in mind if they are
C doing integration or other such evaluations during the calculation
C cycle.
C
C Some variables that might be used defined below:
C
C NRAD = resistor for which the heat flow is to be calculated.
C COEFF & EXPO = the system material property arrays which are
C used to evaluate emissivities and transmissivities.
C TT1 & TT2 = Node temperature in absolute units
C TT3, T1SAVE & T2SAVE are node temperatures in ICCALC units
C
C Other variables one might want
C
C temps = system temperature array.
C qgray = heat transmitted from node 1 to node 2 by the gray resistor.
C igtype = 1 => r = (1-e)/(e*a)
C 2 => r = 1/(f*a*tau)
C 3 => r = 1/( f*a*(1-tau) )
C 4 => r = 1/(f*a)
C 5 => r = 1/(f)
C 6 => r = (1-e)/(e*a), e=constant
C 7 => r = 1/(f*a*tau), tau from extinction coefficient
C 8 => r = 1/( f*a*(1-tau) ), tau from extinction coefficient
C where tau is the transmissivity of the participative media
C and e is the surface emissivity.
C (r is the value of the gray radiative resistor).
C 9 => r = 1/(f*a*tau)
C 10 => r = 1/( f*a*(1-tau) )
C 11 => r = 1/(f*a*tau), tau from extinction coefficient
C 12 => r = 1/( f*a*(1-tau) ), tau from extinction coefficient
C
C Note: Subtypes 9-12 have F and A pre-multiplied and
C stored in the GSFACT array.
C
C igprop = material property index for emissivity, transmissivity, or
C extinction coefficient.
C gsfact is the shape factor.
C garea is the surface area.
C sbc is the stephan-boltzman constant.
C logp = logical print / no-print variable for resistor data.
C gdist = view factor distance, used with extinction coefficients.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C define call list variables
C
INTEGER NRAD
C
DOUBLE PRECISION COEFF( J3, J4 ), EXPO( J3, J4 )
DOUBLE PRECISION GVAL
DOUBLE PRECISION QGRAY
DOUBLE PRECISION TT1
DOUBLE PRECISION TT2
DOUBLE PRECISION TT3
DOUBLE PRECISION T1SAVE
DOUBLE PRECISION T2SAVE
C
C############################################################################
C
C Declare the common block variables
C
C
C############################################################################
C
C*C*RLH INTEGER MID
C*C*RLH INTEGER IO
C*C*RLH INTEGER IGTYPE
C*C*RLH INTEGER IGPROP
!
C*C*RLH LOGICAL LOGP
!
C*C*RLH DOUBLE PRECISION SBC
C*C*RLH DOUBLE PRECISION EMISS
C*C*RLH DOUBLE PRECISION TAU
C*C*RLH DOUBLE PRECISION GAREA
!!
C*C*RLH COMMON / IA14 / IGTYPE( 1 )
C*C*RLH COMMON / IA15 / IGPROP( 1 )
C*C*RLH COMMON / IA45 / MID( 1 )
C
C*C*RLH COMMON / IB33 / IO
!
C*C*RLH COMMON / LB3 / LOGP
!
C*C*RLH COMMON / RA23 / GAREA( 1 )
!
C*C*RLH COMMON / RB22 / SBC
C*C*RLH COMMON / RB51 / EMISS
C*C*RLH COMMON / RB52 / TAU
C
C############################################################################
C Declare the local variables.
C
C*C*RLH DOUBLE PRECISION PROPS
C*C*RLH DOUBLE PRECISION F12
C
C*C*RLH INTEGER N3
C*C*RLH INTEGER LIPROP
C*C*RLH INTEGER MPID
C
C############################################################################
C Declare external functions
C
C*C*RLH EXTERNAL PROPS
C
C############################################################################
C
C*C*RLH LIPROP = ABS( IGPROP( NRAD ) )
C*C*RLH MPID = MID( LIPROP )
C
C*C*RLH IF( IGTYPE( NRAD ) .EQ. 13 .AND. ( MPID .EQ. 303117 .OR.
C*C*RLH 1 MPID .EQ. 370717 ) ) THEN
C
C this section is for radiative resistors with one surface node.
C
C calculate the emissivity.
C
C*C*RLH EMISS = EXPO( LIPROP, 2 )
C
C............................................................................
C
C check to see that emiss doesn't blow up the calculation.
C
C*C*RLH IF( EMISS .GT. 0.9999D+00 ) THEN
C*C*RLH EMISS = 0.9999D+00
C*C*RLH ELSE IF( EMISS .LE. 0.D+00 ) THEN
C*C*RLH QGRAY = 0.D+00
C*C*RLH GOTO 9000
C*C*RLH ENDIF
C
C............................................................................
C
C calculate the heat transmitted from node 1 to 2.
C
C*C*RLH GVAL = SBC * ( EMISS * GAREA(NRAD) ) *
C*C*RLH $ ( TT1 * TT1 + TT2 * TT2 ) * ABS( TT1 + TT2 )
C
C*C*RLH QGRAY = GVAL * ( TT1 - TT2 )
C
C*C*RLH END IF
C
C############################################################################
C
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C*C*RLH 9000 CONTINUE
C
RETURN
C
C############################################################################
C
END
C
C############################################################################
C #
C S U B R O U T I N E U R S T R T #
C #
C############################################################################
C
SUBROUTINE URSTRT
C
C............................................................................
C
C URSTRT --> Called from QTRAN subroutine GETRST, this
C subroutine is called if QTRAN is resuming
C execution from a restart file.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C
C............................................................................
C
C
C............................................................................
C
C
RETURN
END
C
C############################################################################
C #
C S U B R O U T I N E U S O L #
C #
C############################################################################
C
SUBROUTINE USOL( T1RH, T2RH, T3RH, JPROP, JTYPE, RHNPNT, T1RC,
$ T2RC, RCPROP, RCNPNT, T1RR, T2RR, T3RR, IGTYPE, IGPROP, GSFACT,
$ GAREA, RRNPNT, QCARD, IFLIST, QIPPNT, ITLIST, FIXNUM, FIXVAL,
$ TFIX, CPROP, CRHO, RCN, RHN, RRN, QIP, ALPHA, IALIAS, GP, COEFF,
$ EXPO, RCAL, P, TX, TY, SETTIM, PIDSET, PIDPAR, PIDID, CVOL,
$ TEMPS, QVECT, OTEMPS, T1RW, T2RW, T3RW, IWTYPE, IWPROP, RWNPNT,
$ RWN, WSPROP, WLPROP, SFACTR, WAREA, T1RF, T2RF, ICPFLO, RFN,
$ RFNPNT, RMDOT, DTMAXA, PRINTA, T1CAPS, CNPNT, CN, MID, PID,
$ MFID, QMFACT, TMFACT, QBASE, GDIST, WDIST,
$ TERROR, GVALCA, QINMAC, AVNODH, FVAR, GUMTRX, ICNTRL,
$ INDRLX, IPRLXC, RLXTBS, RLXTBT, RELAXV, EFACTB, RELAXM, RERROR,
$ IRRLXC,
$ T1RP, T2RP, PTYPE, PPROP, RPN, RPNPNT, IMLIST, MIP,
$ MIPPNT, IPLIST, PFIX, PALIAS, PIALAS, HIALAS,
$ PGP, MMFACT, PMFACT, MDBASE, PRHOE, VELOCP, PRHO,
$ PRESS, OPRESS, MDOTP, MDOTND, HYCCE, DIFHED, QMDOTP, QBASEF )
C
C############################################################################
C #
C This subroutine contains solution option 1000. #
C #
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.blk file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C
C Define variables
C
INTEGER ICNTRL
C
C integer arrays.
C
INTEGER JPROP( J1, J2 )
INTEGER T1RH(*)
INTEGER T2RH(*)
INTEGER T3RH(*)
INTEGER JTYPE(*)
INTEGER RCPROP(*)
INTEGER T1RC(*)
INTEGER T2RC(*)
INTEGER RCNPNT(*)
INTEGER RHNPNT(*)
INTEGER T1RR(*)
INTEGER T2RR(*)
INTEGER T3RR(*)
INTEGER QCARD( J14, J15 )
INTEGER IFLIST( J17, J16 )
INTEGER QIPPNT(*)
INTEGER ITLIST( J18, J19 )
INTEGER RRNPNT(*)
INTEGER FIXNUM(*)
INTEGER FIXVAL(*)
INTEGER T1RW(*)
INTEGER T2RW(*)
INTEGER T3RW(*)
INTEGER IWTYPE(*)
INTEGER IWPROP(*)
INTEGER RWNPNT(*)
INTEGER RWN(*)
INTEGER TFIX(*)
INTEGER CPROP(*)
INTEGER CRHO(*)
INTEGER RCN(*)
INTEGER RRN(*)
INTEGER RHN(*)
INTEGER QIP(*)
INTEGER T1RF(J30)
INTEGER T2RF(J30)
INTEGER ICPFLO( J30, 2 )
INTEGER RFN(*)
INTEGER RFNPNT(*)
INTEGER T1CAPS(*)
INTEGER CNPNT(*)
INTEGER CN(*)
INTEGER PIDPAR(*)
INTEGER PIDID(*)
INTEGER IALIAS(*)
INTEGER MID(*)
INTEGER PID(*)
INTEGER MFID(*)
INTEGER IGTYPE(*)
INTEGER IGPROP(*)
INTEGER INDRLX( MAXT )
INTEGER IPRLXC( MAXT )
INTEGER IRRLXC( MAXT )
C
C Hydraulic variables
C
INTEGER T1RP ( * )
INTEGER T2RP ( * )
INTEGER PTYPE ( * )
INTEGER PPROP ( J44, J46 )
INTEGER RPN ( * )
INTEGER RPNPNT( * )
INTEGER IMLIST( J48, J49 )
INTEGER MIP ( * )
INTEGER MIPPNT( * )
INTEGER IPLIST( J50, J51 )
INTEGER PFIX ( * )
INTEGER PALIAS( * )
INTEGER PIALAS( * )
INTEGER HIALAS( * )
C
C............................................................................
C
C real arrays.
C
DOUBLE PRECISION ALPHA(*)
DOUBLE PRECISION GP( J1, J6 )
DOUBLE PRECISION COEFF( J3, J4 )
DOUBLE PRECISION EXPO( J3, J4 )
DOUBLE PRECISION RCAL(*)
DOUBLE PRECISION P(J7,J8)
DOUBLE PRECISION TX( J9, J10 )
DOUBLE PRECISION TY( J9, J10 )
DOUBLE PRECISION GSFACT(*)
DOUBLE PRECISION GAREA(*)
DOUBLE PRECISION SETTIM(*)
DOUBLE PRECISION WSPROP(*)
DOUBLE PRECISION WLPROP(*)
DOUBLE PRECISION SFACTR(*)
DOUBLE PRECISION WAREA(*)
DOUBLE PRECISION CVOL(*)
DOUBLE PRECISION TEMPS(*)
DOUBLE PRECISION QVECT(*)
DOUBLE PRECISION OTEMPS(*)
DOUBLE PRECISION RMDOT(*)
DOUBLE PRECISION DTMAXA( J32, J33 )
DOUBLE PRECISION PRINTA( J34, J35 )
DOUBLE PRECISION PIDSET( J22, J23 )
DOUBLE PRECISION QMFACT( J17 )
DOUBLE PRECISION TMFACT( J18 )
DOUBLE PRECISION QBASE( MAXT )
DOUBLE PRECISION GDIST(*)
DOUBLE PRECISION WDIST(*)
DOUBLE PRECISION TERROR (MAXT, J27 )
DOUBLE PRECISION GVALCA(*)
DOUBLE PRECISION QINMAC( MAXT )
DOUBLE PRECISION AVNODH( MAXT )
DOUBLE PRECISION FVAR( MAXT, * )
DOUBLE PRECISION GUMTRX( M1, M1 )
DOUBLE PRECISION RLXTBS( J39, J40)
DOUBLE PRECISION RLXTBT( J39, J40)
DOUBLE PRECISION RELAXV( MAXT )
DOUBLE PRECISION EFACTB( MAXT )
DOUBLE PRECISION RELAXM( MAXT )
DOUBLE PRECISION RERROR( MAXT, J41 )
DOUBLE PRECISION QBASEF(MAXT)
C
C hydraulic variables
C
DOUBLE PRECISION PGP ( J44, J45 )
DOUBLE PRECISION MMFACT( * )
DOUBLE PRECISION PMFACT( * )
DOUBLE PRECISION MDBASE( * )
DOUBLE PRECISION PRHOE ( * )
DOUBLE PRECISION VELOCP( * )
DOUBLE PRECISION PRHO ( * )
DOUBLE PRECISION OPRESS( * )
DOUBLE PRECISION PRESS ( * )
DOUBLE PRECISION MDOTP ( * )
DOUBLE PRECISION MDOTND(*)
DOUBLE PRECISION HYCCE(*)
DOUBLE PRECISION DIFHED(*)
DOUBLE PRECISION QMDOTP(*)
C
C############################################################################
C
C At this point, you must build your solution matrix in whatever form
C necessary and solve the resultant system.
C
C############################################################################
C
RETURN
END
C
C############################################################################
C
SUBROUTINE UWAVER( NWRAD, TT1, TT2, TT3, T1SAVE, T2SAVE, F1, F2,
$ COEFF, EXPO, GVAL, QWAVE )
C
C############################################################################
C
C This subroutine calculates the wave length dependent radiation
C as defined by the user.
C
C NWRAD = radiant resistor number
C TT1 = temperature of node 1 in absolute units
C TT2 = temperature of node 2 in absolute units
C TT3 = temperature of node 3 in iccalc units
C T1SAVE = temperature of node 1 in iccalc units
C T2SAVE = temperature of node 2 in iccalc units
C F1 = fraction of energy leaving node 1 that falls in the
C wave band between ws and wl.
C F2 = fraction of energy leaving node 2 that falls in the
C wave band between ws and wl.
C COEFF & EXPO = the system material property arrays which are
C used to evaluate emissivities and transmissivities.
C GVAL = conductor value
C QWAVE = heat flux through this radiation resistor
C
C IWPROP = internal property ID
C SFACTS = view factor value
C WDIST = distance
C WAREA = area
C
C The property, view factor value, and distance are array to hold that the
C user wishes to pass through to their radiation calculations. IWPROP
C is defined in the field material property and could be information
C that is passed through to reference other properties or property type
C information. SFACTS and WDIST could be arrays that hold spatial
C information. WAREA is the area associated with node 1 and will be
C supplied by PATRAN.
C
C The Gval and QGRAY values must be calculated and passed back to the
C calling routine.
C
C############################################################################
C
IMPLICIT NONE
C
C############################################################################
C
C Declare problem dimensions and common block where they are defined.
C This is taken from the common.dims file.
C
INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
COMMON /DIMS/ J1, J2, J3, J4, J5, J6, J7, J8, J9,
$ J10, J11, J12, J13, J14, J15, J16, J17, J18,
$ J19, J20, J21, J22, J23, J24, J25, J26, J27,
$ J28, J29, J30, J31, J32, J33, J34, J35, J36,
$ J37, J38, J39, J40, J41, J42, J43,
$ J44, J45, J46, J47, J48, J49, J50, J51, J52,
$ MAXT, MAXTNC, MAXTNH, MAXTNR, MAXTNW, MAXTQN,
$ MAXTNF, M1, MAXP
C
C############################################################################
C define call list variables
C
DOUBLE PRECISION TT1
DOUBLE PRECISION TT2
DOUBLE PRECISION TT3
DOUBLE PRECISION T1SAVE
DOUBLE PRECISION T2SAVE
DOUBLE PRECISION F1
DOUBLE PRECISION F2
DOUBLE PRECISION GVAL
DOUBLE PRECISION QWAVE
DOUBLE PRECISION COEFF( J3, J4 )
DOUBLE PRECISION EXPO( J3, J4 )
C
INTEGER NWRAD
C
C############################################################################
C
C Declare the common variables
C
C*C*RLH DOUBLE PRECISION WAREA
C
C*C*RLH DOUBLE PRECISION SBC
C
C*C*RLH INTEGER IWTYPE
C*C*RLH INTEGER IWPROP
C
C*C*RLH COMMON / IA29 / IWTYPE( 1 )
C*C*RLH COMMON / IA30 / IWPROP( 1 )
C
C*C*RLH COMMON / RA14 / WAREA( 1 )
C
C*C*RLH COMMON / RB22 / SBC
C
C############################################################################
C Declare the local variables.
C
C*C*RLH DOUBLE PRECISION EMISS
C
C*C*RLH DOUBLE PRECISION PROPS
C
C
C############################################################################
C Declare external functions
C
C*C*RLH EXTERNAL PROPS
C
C############################################################################
C
C*C*RLH IF( IWTYPE( NWRAD ) .EQ. 1 .AND. IWPROP( NWRAD ) .EQ. 100103 ) THEN
C
C Now the user can make the necessary calculation.
C The example shown is what could be done to evaluate the function
C ( e * A ) / ( 1 - e ) between a surface and radiosity node.
C
C
C calculate the emissivity.
C
C*C*RLH EMISS = PROPS( COEFF, EXPO, T1SAVE, IWPROP(NWRAD) )
C
C............................................................................
C
C check to see that emiss doesn't blow up the calculation.
C
C*C*RLH IF( EMISS .GT. 0.9999D+00 ) THEN
C*C*RLH EMISS = 0.9999D+00
C*C*RLH ELSE IF( EMISS .LE. 0.00001D+00 ) THEN
C*C*RLH QWAVE = 0.D+00
C*C*RLH GVAL = 0.D+00
C*C*RLH RETURN
C*C*RLH ENDIF
C
C............................................................................
C
C calculate the heat transmitted from node 1 to 2.
C
C*C*RLH GVAL = SBC * ( EMISS * WAREA(NWRAD) ) / ( 1.D+00 - EMISS ) *
C*C*RLH $ ( F1 * F1 * TT1 * TT1 + F2 * F2 * TT2 * TT2 ) *
C*C*RLH $ ABS( F1 * TT1 + F2 * TT2 )
C
C*C*RLH QWAVE = GVAL * ( F1 * TT1 - F2 * TT2 )
C
C*C*RLH ENDIF
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
RETURN
C
C############################################################################
C
END