SUBROUTINES in PIL
SUBROUTINE LIBRARY FOR GROUND
SUBROUTINES,
GROUND-ONLY
Subroutines may be appended to Q1 files. They are called by the statement
CALL <name>, where <name> is a character string of any length, of which only the first 20 characters
are significant.
Subroutine names must be upper-case, as must the SUBROUTINE statement itself.
A subroutine consists of a sequence of PIL commands beginning with a SUBROUTINE statement and ending with an ENDSUB statement. Execution of a Q1 file starts at the first line and is terminated by an ENDMAIN command. All subroutines must be below this MAIN program in the Q1 file. An example is:
Multiple ENDMAIN statements are permitted as alternative exit points.
Subroutines can be called to a maximum depth of 20. The maximum number of subroutines allowed in any Q1 file is 300. Subroutines can be called from within DO loops and other such constructs, which can also be used within subroutines.
Examples of Input-library cases which use subroutines are: 162, b532, b577, b579, f206, f209 and r220.
--------------
SUBROUtines which may be CALLed in Fortran coding sequences fall into three categories, namely:
For information about these see the PHENC entries:
The following subroutines can be called from GROUND, GREX3, and GX.. coding. They are grouped as:-
DATPRN(Y,N,Y,N,N,Y,Y,N,Y,Y,..........) A 24-argument subroutine called in GROUP 20 to activate the print-out of the input data set in the Satellite.
DUMP GROUND subroutine used to dump on to the disc all EARTH arrays for restarts. For information on use see entry for SAVE in Group 24.
GETONE(EARTH-variable index,name of variable,IY,IX) GROUND subroutine for getting the value of a variable in an EARTH array which pertains to a selected point (IX,IY) in the currently- visited slab.
GETCOV('patch name',variable index,coefficient,value) Returns coefficient & value for specified patch name & variable.
GETPTC('patch name',type,ixf,ixl,iyf,iyl,izf,izl,itf,itl) Returns patch parameters of specified patch name.
GETSPD(MODEL,VNAME,ITYPE,RVAL,IVAL,LVAL,CVAL,IERR) Reads the EARDAT file for lines starting MODEL VNAME , these lines being the result of including the PIL SPEDAT command in the Q1 file. It sets the value of:
GETSOR('patch name',variable index,net source) Returns net source of variable at specified patch, when called at the end of a sweep ( GROUND Group 19, section 7 ).
GETYX(EARTH-variable index,name of array,nydim,nxdim) GROUND subroutine for getting, for one constant-z slab, the values of an EARTH variable, which are then put into a 2D GROUND array, named and dimensioned by the user.
GETZ(EARTH-variable index,name of array,array dimension) GROUND subroutine for getting NZ values of those EARTH quantities that vary with IZ but not with IX or IY, such as ZGNZ,ZWNZ or WGRDNZ, whose indices are in COMMON/GRDLOC.
GRAPH(AA,NDA,NAV,AL,OA,NDO,NOP,OLA,ISA,APS,OPS,ISC,LU) GROUND subroutine used to produce line printer plots. The meanings of the arguments are:
| AA | array name of abscissa, |
| NDA | dimension of abscissa array (= first dimension of ordinate array) |
| NAV | number of abscissa values |
| AL | abscissa label (character constant or variable) |
| OA | array name of ordinate (2d array) |
| NDO | second dimension of ordinate array |
| NOP | number of ordinates to be plotted |
| OLA | ordinate label array (character) |
| ISA | scaling indicator array (elements set to 1) |
| APS | size of abscissa plot (0.5 for half page) |
| OPS | size of ordinate plot (0.5) |
| ISC | scaling control index (3) |
| LU | logical unit number (LUPR1 for RESULT file) |
MAKE(EARTH-array name) GROUND subroutine used in Group 1 Section 1 for activating storage in EARTH of variables for which conditional provision exists. The argument is one of the names EASP1......EASP10, which represent arrays of dimension NY,NX,which may be used for temporary variables, or the index names XG2D, XU2D, DXG2D, DXU2D, YG2D, YV2D, DYG2D, DYV2D, RG2D, RV2D, and LGEN1 (qv).
ONLYIF(lowest INDVAR,highest INDVAR,patch name) GROUND subroutine used to define conditions on INDVAR and patch name, under which the following calls to GROUND functions will be executed until another ONLYIF, a RETURN or an END.
PLOTL(line index, value of variable, minimum value of variable, maximum value of variable, number of columns, name of variable, logical-unit number) This subroutine prints a line as follows: index, value, : C :
<---------- no. of columns>
where C is the first character of the name, placed in a horizontal position indicating ( value - min. value)/(max. value - min. value)
PRN(character variable,EARTH-variable name) GROUND subroutine for printing an EARTH array. The character variable is a name of up to four characters to identify the printed array.
PRNYX(character variable,2D GROUND-array name,nydim,nxdim) GROUND subroutine for printing a 2D field array.
PRNYXZ(character variable,3D GROUND-array name,nydim,nxdim, NZ dimension,plane indicator,plane location) GROUND subroutine for printing a 3D field array. The sixth argument, ie. the plane indicator, is set as follows:
0 for a field tabulation at plane location IX;
1 " " " " " " " IY;
2 " " " " " " " IZ.
SETYX(EARTH-variable index,2D GROUND-array name,nydim,nxdim) GROUND subroutine for setting an EARTH variable to the values contained in a user-dimensioned array.
ZMOVE(location of j-th part of grid,velocity of grid,location of the grid at current time,location of the grid at old time, number of parts in the grid) GROUND subroutine for setting up a moving grid with N parts.
The following subroutines are for use in body-fitted grid cases only:
BFCROT(Y,XA,XB,angular speed of rotation,iflag) BFCROT is used to insert rotation source terms to simulate a rotating frame of reference, and is used in GXROTA. The second argument is an array of dimension 3 which specifies the cartesian coordinates of any point on the axis of rotation. The third argument so specifies any other point. The angular speed is taken as positive in the sense of a right-handed screw directed from point A to point B. The last argument is an integer which is normally set to zero, if set to 1 then a 'reduced' pressure system is used in which the centrifugal pressure variation is removed. The subroutine inserts the Coriolis and centrifugal body force terms for the velocity resolutes, over a PHASEM type patch representing the whole solution domain.
CORNER(index, K location, name of user array, nydim, nxdim) GROUND subroutine for use in BFC calculations only, to transfer the CORNER coordinates of all cell corners, for a specified slab-face location, K, into a user-array of dimension nydim, nxdim. In total there are (NY+1)*(NX+1) values of each cartesian coordinate to be returned at each K. The first argument, 'index' is set to 1 to return the XC coordinates, 2 for the YC coordinates or 3 for the ZC coordinates. The slab-face label, K, varies from a value of 1 at the LOW boundary of the domain to NZ+1 at the HIGH boundary.
FNGBFC(Y,iref,IZ) This subroutine puts into Y the BFC geometrical factor referenced by 'iref', a table for which is provided for the SUBROUtine GTIZYX. Y denotes an array index as in a FUNCtion subroutine.
GETCAR is a GROUND subroutine for use when BFC=T to get the cartesian components of the total velocity relative to the global grid coordinate system. In order to use GETCAR, storage for these quantities must be provided by means of the PIL command: STORE(UCRT,VCRT,WCRT). GETCAR is called automatically on the last sweep when this storage is reserved, so needs to be called by the user only when the cartesian components are needed during the computation.
GETPT(i,j,k,XC,YC,ZC) GROUND subroutine for use when BFC=T to get the cartesian coordinates XC, YC, ZC of the corner labelled i, j, k. If the coordinates of a large number of corners, at a given slab, k, are required, then the subroutine CORNER provides more efficient access. See BODY-F for information on the corner labelling convention.
GTIZYX(name,IZ location,name of array,nydim,nxdim) This GROUND subroutine is used to retrieve body-fitted-grid geometry from EARTH at the IZ slab specified (not necessarily the current slab). The information retrieved is that for one slab, and is put into a 2D array which has been named and dimensioned by the user to nxdim, nydim whose values must be at least NX, NY. The names (and reference numbers) for all geometrical quantities are tabulated below. The quantities refer to cells, rather than the cell corners referred to in GETPT. It should be noted that the variable names are defined in the file GRDBFC, which should be included into any subroutine making use of them.
WARNING: The following table is for use with body-fitted grids only.
| Ref.No. | Name | Geometrical quantity description. |
| 1 | APROJE | east-face area projected normal to u (*) |
| 2 | APROJN | north-face area projected normal to v (*) |
| 3 | APROJH | high-face area projected normal to w (*) |
| 4 | VOLUME | cell volume (*) |
| The above are multiplied automatically by the appropriate porosities in EARTH. The following geometric quantities are NOT so multiplied. | ||
| 5 | ASURFE | surface area of east face |
| 6 | ASURFW | surface area of west face |
| 7 | ASURFN | surface area of north face |
| 8 | ASURFS | surface area of south face |
| 9 | ASURFH | surface area of high face |
| 10 | ASURFL | surface area of low face |
| 11 | DHXPE | twice normal distance of node from E & W walls of cell |
| 12 | XCORNR | cartesian-x coordinate of low-south-west cell corner |
| 13 | DHYPN | twice normal distance of node from N & S walls of cell |
| 14 | YCORNR | cartesian-y coordinate of low-south-west cell corner |
| 15 | DHZPH | twice normal distance of node from H & L walls of cell |
| 16 | ZCORNR | cartesian-z coordinate of low-south-west cell corner |
| 17 | MFAEPE | coefficient for u in formula for east cell face mass flux |
| 18 | MFAEPN | coeff. for v in formula for east cell face mass flux |
| 19 | MFAEPH | coeff. for w in formula for east cell face mass flux |
| 20 | MFANPE | coeff. for u in formula for north cell face mass flux |
| 21 | MFANPN | coeff. for v in formula for north cell face mass flux |
| 22 | MFANPH | coeff. for w in formula for north cell face mass flux |
| 23 | MFAHPE | coeff. for u in formula for high cell face mass flux |
| 24 | MFAHPN | coeff. for v in formula for high cell face mass flux |
| 25 | MFAHPH | coeff. for w in formula for high cell face mass flux |
| 26 | DXGPE | distance from cell node to east-neighbour node (*) |
| 27 | DYGPN | distance from cell node to north-neighbour node (*) |
| 28 | DZGPH | distance from cell node to high-neighbour node (*) |
| 29 | UDIVZE | Earth Coefficients for the U Velocity divergence terms (XZ) |
| 30 | UDIVZN | |
| 31 | UDIVZH | |
| 32 | VDIVZE | Earth Coefficients for the V-velocity divergence terms (YZ) |
| 33 | VDIVZN | |
| 34 | VDIVZH | |
| 35 | WCRVYE | Earth Coefficients for the W-velocity curvature terms |
| 36 | WCRVYN | |
| 37 | WCRVYH | |
| 38 | XCORNO | old cartesian-x coordinate of l-s-w cell corner |
| 39 | YCORNO | old cartesian-y coordinate of l-s-w cell corner |
| 40 | ZCORNO | old cartesian-z coordinate of l-s-w cell corner |
| 41 | UCRVYE | Earth Coefficients for the U-velocity curvature terms |
| 42 | UCRVYN | |
| 43 | UCRVYH | |
| 44 | UNIVHE | XC resolute of unit vector aligned local direction of w |
| 45 | UNIVHN | YC resolute of unit vector aligned local direction of w |
| 46 | UNIVHH | ZC resolute of unit vector aligned local direction of w |
| 47 | VDIVXE | Earth Coefficients for the V-velocity divergence terms (XY) |
| 48 | VDIVXN | |
| 49 | VDIVXH | |
| 50 | WDIVXE | Earth Coefficients for the W-velocity divergence terms (XZ) |
| 51 | WDIVXN | |
| 52 | WDIVXH | |
| 53 | UDIVYE | Earth Coefficients for the U-velocity divergence terms (YX) |
| 54 | UDIVYN | |
| 55 | UDIVYH | |
| 56 | VCRVXE | Earth Coefficients for the V-velocity curvature terms |
| 57 | VCRVXN | |
| 58 | VCRVXH | |
| 59 | MFANPL | Earth Coefficients |
| 60 | MFAHPS | " |
| 61 | MFAHPW | " |
| 62 | WDIVYE | Earth Coefficients for the W-velocity divergence |
| 63 | WDIVYN | terms (YZ) |
| 64 | WDIVYH | |
| 65 | MFAEPS | Earth coefficients |
| 66 | MFAEPL | |
| 67 | MFANPW | |
| 68 | XP | XC coordinate of cell centre |
| 69 | YP | YC coordinate of cell centre |
| 70 | ZP | ZC coordinate of cell centre |
| 71 | UCARTE | coeff for u in eqn for x-directed cartesian resolute |
| 72 | UCARTN | coeff for v in eqn for x-directed cartesian resolute |
| 73 | UCARTH | coeff for w in eqn for x-directed cartesian resolute |
| 74 | VCARTE | coeff for u in eqn for y-directed cartesian resolute |
| 75 | VCARTN | coeff for v in eqn for y-directed cartesian resolute |
| 76 | VCARTH | coeff for w in eqn for y-directed cartesian resolute |
| 77 | WCARTE | coeff for u in eqn for z-directed cartesian resolute |
| 78 | WCARTN | coeff for v in eqn for z-directed cartesian resolute |
| 79 | WCARTH | coeff for u in eqn for z-directed cartesian resolute |
| 80 | MFAEEN | Earth coefficients |
| 81 | MFAEES | |
| 82 | MFAEEH | |
| 83 | MFAEEL | |
| 84 | MFANNE | |
| 85 | MFANNW | |
| 86 | MFANNH | |
| 87 | MFANNL | |
| 88 | MFAHHE | |
| 89 | MFAHHW | |
| 90 | MFAHHN | |
| 91 | MFAHHS | |
| 92 | UNIVEE | XC resolute of unit vector aligned local direction of u |
| 93 | UNIVEN | YC resolute of unit vector aligned local direction of u |
| 94 | UNIVEH | ZC resolute of unit vector aligned local direction of u |
| 95 | UNIVNE | XC resolute of unit vector aligned local direction of v |
| 96 | UNIVNN | YC resolute of unit vector aligned local direction of v |
| 97 | UNIVNH | ZC resolute of unit vector aligned local direction of v |
| 98 | CYCUXE | Earth coeff's associated with Cyclic boundaries (X) |
| 99 | CYCUXN | |
| 100 | CYCUXH | |
| 101 | CYCUYE | Earth coeff's associated with Cyclic boundaries (Y) |
| 102 | CYCUYN | |
| 103 | CYCUYH | |
| 104 | CYCUZE | Earth coeff's associated with Cyclic boundaries (Z) |
| 105 | CYCUZN | |
| 106 | CYCUZH | |
| 107 | CC2FCN | Earth coefficients |
| 108 | CC2FCS | |
| 109 | CC2FCE | |
| 110 | CC2FCW | |
| 111 | CC2FCH | |
| 112 | CC2FCL |
The reference numbers marked with an asterisk (*), refer to quantities that can also be obtained at the current IZ-slab by the use of GETYX with GRDLOC indices AEAST, ANORTH, AHIGH, VOL, DXG2D, DYG2D and DZG2D respectively.
The following 11 subroutines are used in GROUND coding to modify BFC grids during the course of a calculation.
This subroutine MUST be called whenever the grid corner coordinates have been modified before resuming the flow calculation in order to re-calculate all the grid-dependent geometric quantities used in the calculation. Two calls are required, the first with n=1, the second with n=2.
This subroutine writes the complete set of corner coordinates stored in the F-array to a disc file on logical unit 15.
This subroutine dumps the phi (fields) and grid (xyz) files to the appropriate logical units when the filenames CHPHI and CHXYZ (character constants or variables) are not blank.
ICORN(I,J,K) ... integer function to be added to L0B(corner) to yield the F-array location if the corner coordinate for cell I,J,K.
This subroutine prints corner coordinates from the F-array for the slab with index K, where XYZCRN is the corner coordinate index as given in the table below.
| XYZCRN | corner coordinates |
| 12 | XC |
| 14 | YC |
| 16 | ZC |
| 38 | old XC |
| 39 | old YC |
| 40 | old ZC |
This subroutine reads the cell-corner x,y,or z coordinates for the slab with index K from the XYZ file and inserts them into the appropriate segments of the F-array.
This subroutine sets the coordinate in the F-array to 'value', for the corner with indices I,J, and K.
This subroutine writes to the XYZDA file the cell-corner coordinates stored in the F-array for the slab with index K. by calling BGEOM. XYZDA must be set to T in the PREFIX file.
XC(I,J,K) ... real function giving the low-south-west-corner coordinate in the Cartesian-x direction, for cell I,J,K.
XCO(I,J,K) ... real function giving the "old" low-south-west-corner coordinate in the Cartesian-x direction, for cell I,J,K.
YC(I,J,K) ... real function giving the low-south-west-corner coordinate in the Cartesian-y direction, for cell I,J,K.
YCO(I,J,K) ... real function giving the "old" low-south-west-corner coordinate in the Cartesian-y direction, for cell I,J,K.
ZC(I,J,K) ... real function giving the low-south-west-corner coordinate in the Cartesian-z direction, for cell I,J,K.
ZCO(I,J,K) ... real function giving the "old" low-south-west-corner coordinate in the Cartesian-z direction, for cell I,J,K.
The subroutines that follow are for use primarily in body-fitted grid problems, for the generation of grid-related geometrical quantities by means of vector techniques, but can be used for non-BFC cases and in the SATELLITE.
This subroutine calculates the unit vector V normal to a plane passing through the points A, B, C and D, and additionally calculates the AREA of the quadrilateral ABCD. A, B, C, D and V are arrays of dimension 3, the three elements of each array being the resolutes of the vector. See NORMAL for what happens when A, B, C and D are not coplanar.
This subroutine adds vector B to vector A to give vector C, where A, B and C are arrays of dimension 3.
This subroutine averages vectors A and B to give vector C, where A, B and C are arrays of dimension 3.
This subroutine copies vector A in to vector B, where A and B are arrays of dimension 3.
This subroutine sets vector C to the vector product of A and B, where A, B and C are arrays of dimension 3.
This is a FUNCTION subroutine which sets DOT to the scalar product of the vectors A and B, where A and B are arrays of dimension 3.
This is a subroutine which sets MAG to the magnitude (ie. length) of the vector A, and DIR to the unit vector in the same direction as V. A and DIR are arrays of dimension 3.
This is a subroutine which multiplies vector A by scalar S to give vector B, where A and B are arrays of dimension 3.
This is a subroutine which normalizes vector A,where A is an array of dimension 3.
This is a subroutine which calculates the unit vector normal to a plane through points defined by A, B, C and D. If the latter are not coplanar, an intermediate plane is defined by averaging. A, B, C, D and V are arrays of dimension 3.
This subroutine changes vector A to -A, where A is an array of dimension 3.
This subroutine returns the components P, Q and R of vector V, ie. P=V(1), Q=V(2), R=V(3), where V is an array of dimension 3.
This subroutine subtracts vector B from A to give C, where A, B and C are arrays of dimension 3.
This is a FUNCTION subroutine which sets TRIPLE to the triple product of vectors A, B and C ie. A . ( B x C ), where A, B and C are arrays of dimension 3.
This is a FUNCTION subroutine which sets VECM to the magnitude of vector A, where A is an array of dimension 3.
This subroutine defines the vector V from the given scalar cartesian components P, Q and R, where A is an array of dimension 3.
A set of 16 Fortran subroutines effect substitution in a compact way. They have the form:-
SUBnt (n pairs of arguments of the same type) ,where n=1,2,3 or 4
and t is absent for integer arguments,
t is R for real arguments,
t is L for logical arguments and
t is H for character string arguments. Thus: CALL SUB3(I1,4,I3,K+1,J,I2/2) is equivalent
to the 3 statements:
I1=4
I3=K+1
J=I2/2
and
CALL SUB4R(A,3.2,B,C**2,D,FLOAT(I1),E,SQRT(F)) is equivalent to the statements:
A=3.2
B=C**2
D=FLOAT(I1)
E=SQRT(F)
**WARNING** : Do NOT use in even-numbered argument locations variables which are set in
odd-numbered argument locations.Thus: CALL SUB2(I1,2,I2,I1**2) cannot be relied upon to be
the equivalent of,
I1=2
I2=I1**2
A further set of 12 subroutines, facilitating print-out of integer, reals and logicals, has the form:
WRITnt( n pairs of arguments )
where n = 1,2,3 or 4
and t is I for integers,
t is R for reals and
t is L for logicals. The first of each argument pair is always an 8-character string.
Thus:
CALL WRIT2I('ISWEEP ',ISWEEP,'ISTEP ',ISTEP)
will print:
ISWEEP = integer ISTEP = integer
Finally,
CALL WRITBL
will print a blank line;
CALL WRITST
will print a line of asterisks; and
CALL WRIT4('4 characters')
CALL WRIT8('8 characters') and
CALL WRIT40('40 characters')
will print the inserted character strings.
wbs