A library of subroutines is available for CALLing from GROUND coding when arithmetic operations are to be performed on EARTH arrays.
These subroutines operate on STOREd or SOLVEd phi-variables, COefficients, VALues, and auxiliary variables, or on parts of the F-array for which storage has been created by MAKE, NFUSER or GXMAKE.
In the following argument lists:
The variable Y is calculated over the area of the 'PATCH' which is currently active. The limits of this patch are set by COVAL in the SATELLITE, when the function is responding to the setting of the coefficient or value of COVAL to GRND.
Otherwise it is calculated over the whole of the constant-z slab, unless IXF, IXL, IYF and IYL are explicitly set before the function CALL.
The functions are (with recently-introduced ones presented first:-
FNZERO_IN_SLD_VAC Y = 0.0 in cells occupied by solid
or having zero volume porosity
FNQUOT(Y,X1,X2) Y = X1 / X2
FN0(Y,X) Y = X
FN1(Y,A) Y = A
FN2(Y,X,A,B) Y = A + B*X
FN3(Y,X,A,B,C) Y = A + B*X + C*X**2
FN4(Y,X,A,B,C,D) Y = A + B*X + C*X**2 + D*X**3
FN5(Y,X,A,B,C,D,E) Y = A + B*X + C*X**2 + D*X**3+E*X**4
FN6(Y,X,A,B) Y = 1./(A + B*X)
FN7(Y,X,A,B,C) Y = (A + B*X)/(1.0 + C*X)
FN8(Y,X,A,B,C,D) Y = A*(X + B)**C + D
FN9(Y,X,A,B,C,D) Y = A*((X + B)**C)/(1.0 + D*X)
FN10(Y,X1,X2,A,B1,B2) Y = A + B1*X1 + B2*X2
FN11(Y,X1,X2,A,B1,B2,C1,C2,C12)
Y = A + B1*X1 + B2*X2 + C1*X1**2 + C2*X2**2 + C12*X1*X2
FN12(Y,X1,X2,X3,A,B1,B2,B3) Y = A + B1*X1 + B2*X2 + B3*X3
FN13(Y,X1,X2,X3,X4,A,B1,B2,B3,B4)
Y = A + B1*X1 + B2*X2 + B3*X3 + B4*X4
FN14(Y,X1,X2,X3,X4,X5,A,B1,B2,B3,B4,B5)
Y = A + B1*X1 + B2*X2 + B3*X3 + B4*X4 + B5*X5
FN15(Y,X1,X2,A,B1) Y = (A + B1*X1)/X2
FN16(Y,X1,X2,A,B1,B2) Y = (A + B1*X1)/(1.0 + B2*X2)
FN17(Y,X1,X2,X3,A,B1,B2,B3) Y = (A + B1*X1)/(1.0 + B2*X2 + B3*X3)
FN18(Y,X1,X2,X3,X4,A,B1,B2,B3,B4)
Y = (A + B1*X1 + B2*X2)/(1.0 + B3*X3 + B4*X4)
FN19(Y,X1,X2,X3,X4,A,B,C) Y = A*X1*X2*X3*EXP(B/(C + X4))
FN20(A,X1,X2,X3) A=PATCH-WISE SUM AT SLAB of X1*X2*X3
FN21(Y,X1,X2,A,B) Y = A + B*X1*X2
FN22(Y,A) Y = AMAX1(Y,A)
FN23(Y,A) Y = AMIN1(Y,A)
FN24(Y,X1,X2,A,B,C) Y = A * X1 + B * AMAX1(0.0,X2 + C)
FN25(Y,A) Y = A * Y
FN26(Y,X) Y = Y * X
FN026(Y,X1,X2): Y = X1 * X2
FN27(Y,X) Y = Y / X
FN27D(Y,X1,X2): Y = Y / (X1 * X2)
FN28(Y,X,A) Y = A / X
FN28P2(Y,X1,X2,A,B): Y = A / (X1 + B * X2)
FN29(Y,IPHASE,IDIR,IU,IV,IW) Y = resultant velocity of the phase
along wall, the normal to which is specified by IDIR as follows:
IDIR=1 for x direction;
IDIR=2 for y direction; and
IDIR=3 for z direction.
The last 3 arguments must be set as follows:
IU=1 if STORE(U1)=.true., else IU=0;
IV=1 if STORE(V1)=.true., else IV=0; and,
IW=1 if STORE(W1)=.true., else IW=0.
The second parameter IPHASE = 1 for first-phase velocities,
U1, V1, W1; and
IPHASE = 2 for second-phase velocities,
U2, V2, W2.
FN30(Y) Y = SQRT(Y)
FN30a(Y,X) Y = SQRT(X)
FN31(Y,X1,X2,A,B1,B2) Y = A * (X1**B1) * (X2**B2)
FN32(Y,X,A,B) Y = Y * exp(A/(B + X))
FN33(Y,A) Y = Y + A
FN34(Y,X,A) Y = Y + A * X
FN35(Y,X,A,B) Y = Y + A * X**B
FN36(Y,X,A,B) Y = Y + A * exp(B * X)
FN37(Y,X,A) Y = Y * X**A
FN38(Y,X,A) Y = Y * A * exp(B * X)
FN39(Y,X1,X2,X3,X4) Y = X1 * X2 + X3 * X4
FN40(Y) Y = ABS(Y)
FN41(Y,X,A,B) Y = SIN(A * X + B)
FN42(Y,X,A,B) Y = COS(A * X + B)
FN43(Y,X,A,B) Y = ALOG(A * X + B)
FN44(Y,X,A,B,C) Y = A * X / (1.0 + B * X ** C)
FN45(Y,X,A) Y = Y * SIN(X +A)
FN46(Y,X,A,B) Y = Y * (A + B * X)
FN47(Y,A,B) Y = AMAX1(A,AMIN1(B,Y))
FN48(Y,X1,X2,X3,X4,A,B) Y = A*X1*X2 + B*X3*X4
FN49(Y,X) Y = Y + X**2
FN50(Y,I) Y = Y**I
FN51(Y,A) Y = Y**A
FN52(Y,X,I) Y = Y * X**I
FN53(Y,X1,X2,A) Y = Y + A*X1*X2
FN54(Y,X1,X2,A) Y = Y + A*X1/X2
FN55(Y,X1,X2,X3,A) Y = A*X1*X2*X3
FN56(Y,X1,X2,X3,A) Y = A*X1*X2/X3
FN57(Y,X1,X2,X3,A) Y = Y + A*X1*X2*X3
FN58(Y,X) Y = X where X > 0.0
FN59(Y,X) Y = 0.0 where X = 0.0
FN60(Y,X) Y = Y + X
FN61(Y,X,A) Y = A * X**3
FN62(Y,X,A) Y = A * Y/X**3
FN63(Y,X,A) Y = A * X**4
FN64(Y,X,A) Y = Y + A * X**4
FN65(Y1,Y2,A) Y1 = A, Y2 = A
FN66(Y,X,A,B) Y = AMAX1(0.0, A + B*X)
FN67(Y,X1,X2,X3,B1,B2) Y = (B1*X1 + B2*X2)/X3
FN68(Y,X1,X2,X3,A,B) Y = (A + B*X1)/(X2*X3)
FN69(Y,X,A,B) Y = Y * (A + B/X)
FN70(Y,X,A,B) Y = B/X when X >> A
FN71(Y,X1,A,B) Y = A * X1**B
FN72(Y,X1,A,I) Y = A * X1**I
FN73(Y,X1,A,I) Y = A * Y * X1**I
FN74(Y,X1,A,B,C)
Y = A * ((X1+C)**3 + (X1+C)**2*(B+C) + (X1+C)*(B+C)**2 + (B+C)**3
FN75(X,X1,A,B,C)
Y = A * abs(X1 - B)**C
FN76(Y,X1,X2) Y = Y * X1 / X2
FN77(Y,X1,X2) Y = Y * X1 * X2
FN78(Y,X1,X2,X3,A) Y = Y + A*X1 * X2 / X3
FN79(A,X1,X2,X3,X4) A = Patchwise sum of X1*X2*X3*X4
FN80(A,X1,X2) A = Patchwise sum of X1*X2
FN82(Y,X1,A,B) Y = AMAX1(Y , A + B*X1)
FN83(Y,X1,A,B) Y = AMIN1(Y , A + B*X1)
FN84(Y,X1,X2,A,B) Y = A / (X1 + B*X2)
FN85(Y,X1,X2,A) Y = Y * X1 * (A - X2)
FN95(Y,X1,A,B,C) Y = A * (X1 - B)**C
FN96(Y,X1,X2,A) Y = Y / (X1 * X2)**A
FN97(Y,X1,X2,A,B,C) Y = A * X1 * (X2 - B)**C
FN98(Y,CM,A,B,C) superseded by GXMDOT
FN99(Y,X,CF,A,B,C,D) superseded by GXCFIP
FN101(Y,X,A) Y is a solution of x = (y - 0.5*sin(2.0*y)) / pi.
This is used for determining the half-angle y subtended at
the centre of a circular pipe by the surface of a liquid layer,
given the volume fraction, X, of the liquid.
The solution is obtained by iteration, which ceases after 20
cycles or when Y changes by less than A, whichever happens first.
FN102(Y,X,A,B,C) y=a*b*(c**3)*(.6666*(sinx)**3-cosx*(x-0.5*sin2x))
This calculates the hydrostatic force on the segment of the
cross-section of a circular-sectioned pipe occupied by the
heavier fluid.
The half-angle at each cross-section is X, supplied by FN101 .
A is the acceleration due to gravity.
B is the density difference between the heavier and lighter fluids.
C is the radius of the pipe.
FN103(Y,X,IDIR) Y = X(next)-X
where IDIR =1,2,3 or for north, south, east or west.
FN104(Y,X,IPLUS) Y = -X and X = Y(next)
This subroutine puts -X into Y, and then adds X;
here 'next' refers to the north when IPLUS=1 and
east when IPLUS=NY.
FN105(Y,X1,X2,IPLUS) Y = X1
unless X2 is negative, in which case Y = X1 shifted IPLUS cells
from the current one.
For example, if X2 refers to v velocities and IPLUS is 1,
Y = X1 in the North cell when v is negative.
FN107(Y,X1,X2,X3) Y = a weighted average of X1 and X2,
the weighting factor depending on the ratio of the values of
convection flux X3 to the time flux, ie. on the Courant number.
FN108(Y,SUM) which returns the patchwise sum of all Y's at the
current Z slab in SUM.
FN108D(X1,X2,SUM1,SUM2) which returns the patchwise sum of X1 in SUM1,
and that X2 in SUM2
FN109(Y,X1,X2,X3) Y = X1*(X2-X3)
FN110(Y,X1,X2) for X2 >= 0 , Y = Y + X2
for X2 >0 , X1="X1" X2
FN111(Y,X) Y = Y/SQRT(ABS(Y)) X = X/ABS(Y)
FN112(Y,X1,X2,X3) Y=Y+X1*(X2-X3)
FN117(Y,X1,X2,IPLUS) Source term for fully-implicit
central- differencing scheme with non-uniform grid
FNDUDX(Y,U1) for y=d(u1)/dx ; FNDUDX(Y,U2) for y=d(u2)/dx.
FNDUDY(Y,U1) for y=d(u1)/dy ; FNDUDY(Y,U2) for y=d(u2)/dy.
FNDUDZ(Y,U1) for y=d(u1)/dz ; FNDUDZ(Y,U2) for y=d(u2)/dz.
FNDVDX(Y,V1) for y=d(v1)/dx ; FNDVDX(Y,V2) for y=d(v2)/dx.
FNDVDY(Y,V1) for y=d(v1)/dy ; FNDVDY(Y,V2) for y=d(v2)/dy.
FNDVDZ(Y,V1) for y=d(v1)/dz ; FNDVDZ(Y,V2) for y=d(v2)/dz.
FNDWDX(Y,W1) for y=d(w1)/dx ; FNDWDX(Y,W2) for y=d(w2)/dx.
FNDWDY(Y,W1) for y=d(w1)/dy ; FNDWDY(Y,W2) for y=d(w2)/dy.
FNDWDZ(Y,W1) for y=d(w1)/dz ; FNDWDZ(Y,W2) for y=d(w2)/dz.
FNDSDX(Y,SCAL) for y=d(scal)/dx
FNDSDY(Y,SCAL) for y=d(scal)/dy
FNDSDZ(Y,SCAL) for y=d(scal)/dz
(scal is any EARTH index for a scalar dependent variable)
FNGENK(Y,1) for phase-1 , &
FNGENK(Y,2) for phase-2.
Y = 2*( (du/dx)**2 + (dv/dy)**2 + (dw/dz)**2 ) +
(du/dy + dv/dx)**2 + (du/dz + dw/dx)**2 +
(dw/dy + dv/dz)**2
Y = the turbulence-energy-generation rate.
x, y & z = distances. u,v & w = velocities.
For BFC calculations, the subroutine GXGENB uses the
cartesian velocity components UCRT, VCRT AND WCRT.
It is called from GXGENK.
FNSYM(Y,IFACE) Set symmetry values at boundaries
FNVLSQ(Y,1) for phase-1, and
FNVLSQ(Y,2) for phase-2, where Y = U**2 + V**2 + W**2 ,
for use in calculations involving the kinetic energy
of the mean motion.
Only velocities flowing into the cells are taken into account.
FNAV(Y,X,'next') Y = the average of X and X(next),
where 'next' is 'N', 'S', 'E', 'W', 'H' or 'L' .
FNIFDV(Y,X,A,LOGIC):
for LOGIC = .TRUE. , Y = Y / X
for LOGIC = .FALSE. , Y = Y / A
FNIFML(Y,X,A,LOGIC):
for LOGIC = .TRUE. , Y = Y * X
FNMAX(Y,X) Y = AMAX1(Y,X)
FNGBFC(Y,iref,IZ)
See SUBROU
section on BFC subroutines.
FNVSLP(Y,A) Y = the maximum of A and
the absolute interphase slip velocity into Y.
If STORE(VREL) is present in Q1, Y will be copied into
LBNAME('VREL').
See the Encyclopaedia entry SUBROU for further service subroutines which can be used in GROUND subroutines.
The following function subroutines are available for use in GROUND coding. Uses of many of them are exemplified in GREX3 and GX... They interact with the EARTH module of PHOENICS by accessing the F-array (See the corresponding PHENC entry for an explanation of how they are to be used).
WARNING: integer and logical functions must be declared as such in subroutines in which they are used, eg:
INTEGER HIGH LOGICAL WHF
The functions are:-
ANYZ(block-location index,any IZ)
AUX(SWV block-location index)
EAST(INDVAR)
HIGH(INDVAR)
INDONE(INDVAR,IYY,IXX)
INDYX(INDVAR)
INDZ(INDVRZ,IZZ)
L0F(INDVAR)
LOW(INDVAR)
LBNAME('variable name')
LBZ(LB,IZZ)
OLD(INDVAR)
NORTH(INDVAR)
SOUTH(INDVAR)
WEST(INDVAR)
GETONE(INDVAR,VALUE,IYY,IXX)
GETYX(LB,2D array name,NYDIM,NXDIM)
GETZ(INDVAR,1D array name,GZ,NZDIM)
PRNYX('name',array name,NYDIM,NXDIM)
SETYX(LB,2D array name,NYDIM,NXDIM)
VARONE(INDVAR,IYY,IXX)
VARYX(INDVAR)
VARZ(INDVRZ,IZZ)
tests whether the variable has its built-in convection terms activated.
tests whether the variable has its built-in diffusion terms activated.
tests whether the variable is equal to 0.0 .
tests whether the variable belongs to the first phase fluid or not.
tests whether the variable is greater than or equal to 0.0 .
tests whether the variable equals GRND.
tests whether the variable has been set equal to GRNDnumber in the SATELLITE.
tests whether the variable is greater than 0.0 .
tests whether the variable is to have its diffusion terms represented by harmonic averaging.
tests whether the variable has its inter-phase transport terms activated.
tests whether the variable is less than or equal to 0.0
tests whether the variable is less than 0.0 .
tests whether the variable is not equal to 0.0 .
tests whether the variable is to have its linear equations solved Jacobi point-by-point.
tests whether the variable is to have its field tabulated.
tests whether the variable is to have its linear-equation corrections monitored.
tests whether the variable is to have its whole-field residuals tabulated/plotted.
tests whether the variable is to have its slab-wise residuals printed (at IZMON).
tests whether the variable is to have its spot-values tabulated plotted.
tests whether the variable is to have its whole-field residuals printed.
tests whether real variable 1 is equal to variable 2 .
tests whether real variable 1 is greater than or equal to variable 2.
tests whether real variable 1 is greater than variable 2 .
tests whether real variable 1 is less than or equal to variable 2 .
tests whether real variable 1 is less than variable 2 .
tests whether real variable 1 is not equal to variable 2 .
tests whether the variable is the dependent variable of a differential equation which is being solved.
tests whether the variable has its built-in source terms activated.
tests whether storage has been provided for the variable in the SATELLITE.
tests whether the variable has its built-in transient term activated
tests whether the variable is to have its linear equations solved in a whole-field manner.
tests whether the variable is to be represented by explicit finite-domain equations.
wbs