GREX3.F is a user-accessible Fortran file exemplifying how GROUND -style subroutines can be used to augment or modify the flow- simulating power of EARTH.
This file contains the single subroutine GREX3, which is called by EARTH. It is structured like the empty GROUND, the PHENC entry on which should be read before this one.
The listing, which can be inspected via /phoenics/d_earth/d_core/grex3.for, is extensively commented. The present article describes the contents of the subroutine, group-by-group, in a more general way.
The FORTRAN subroutine GREX3 is called whenever USEGRX equals T (which is the default setting). It has been created for two purposes:
The first pieces of advice to the prospective user of GROUND are therefore:
Because the GREX3 subroutine is long and the features which it exemplifies very numerous, this 'careful inspection' is neither easy nor quick. The following notes have been written mainly so as to reduce its difficulty, not to make it unnecessary; for it is the listing itself which defines precisely what it is that PHOENICS does.
It should be remarked here, to save both user and computer-time, that no-one is obliged to make use of GREX3. If none of its features are to be utilised, computer-time will be saved by setting USEGRX=F in the SATELLITE; and if the user desires to reduce the size of his EARTH module, he can achieve this by 'commenting out' or removing all the coding in GREX8!between the lines:
IF(.NOT.USEGRX) GO TO 25 and: 25 CONTINUE
and thereafter re-compiling and re-linking the EARTH module. Beginners who want to keep their minds and files uncluttered, and veterans who appreciate the associated economy, may alike find this practice advantageous.
Inspection of the listing of GREX3 reveals that, after a page of preliminary material, its whole content is divided up into precisely the same 24 groups as have already been encountered in the Q1 files and in the FORTRAN files SATLIT and GROUND.
Moreover, many of the groups are further divided into sections, the divisions being exactly the same as those which appear in GROUND. GREX3 is indeed nothing but a particular version of GROUND, differing from it in having many of its groups and sections occupied by coding sequences, whereas those of GROUND were all empty.
The GREX3 subroutine is called repeatedly from within EARTH. What governs the part of the subroutine to which control is sent is the value of IGR, an index which is set within EARTH and is the group number required; and which section of a group is then visited depends on the value of another index, ISC, which is also set within EARTH.
GREX3 calls further subroutines, which appear in four groupings, as follows:
The Fortran ciding of these subroutines can be viewed by way of POLIS (PHOENICS-2/Fortran).
This group consists of two sections, the first of which is called at the start of the run after the SATELLITE data have been read in but before the EARTH F-array segment-location indices have been allocated. It is called at this point so that the user can instruct EARTH to allocate storage.
The user is free to introduce into GROUND subroutines storage for arrays which he wishes to use there; to do so, he should provide DIMENSION statements obeying the rules governing the version of Fortran which is accepted by the user's computer.
However, it is possible for the user to activate certain arrays, for which a 'dormant' provision has been made in EARTH.
There are 10 such arrays , with names EASP1, EASP2, EASP3, ... , EASP10 (see GRDLOC in Appendix 2); and they are activated in Group 1, by means of a call to subroutine MAKE. Inspection of Group 1 of GREX3 will show precisely how this is to be done.
That inspection will reveal that, as well as being called for the variables EASP1, etc, MAKE is additionally called with the arguments XG2D, XU2D, etc, the names in question also being among those which appear in GRDLOC.
Calls of this kind do create storage for the quantities in question; but they further ensure that the array is filled by appropriate values. The quantities in question are all geometrical ones which, if the grid is of the cartesian or polar varieties, are normally stored only one-dimensionally, for reasons of economy; however when MAKE is called, the variables are suitably placed in the elements of a two-dimensional array of dimension NX times NY.
For example, the statement:
CALL MAKE(XG2D)
has the effect of putting the x-coordinates of the grid points into a 2D array, which can thereafter be utilised, if it is desired, in a CALL to one of the service subroutines. This utilisation is also exemplified by the following statement:
CALL FN2(LEN1,XG2D,EL1A,EL1B)
which computes the first-phase length scale from x-coordinate values.
As will be seen by inspection of the listings of GREX3 and GXLEN, both the CALL MAKE and the CALL FN2 are preceded by the same condition, namely IF(GRNDNO(1,EL1)). Obviously, the first CALL has been inserted to provide storage for and values for the variable XG2D, which are made use of in the second CALL; and the signal is given in the SATELLITE by the setting of EL1 to GRND1.
Settings of this kind will be found in the Q1 files of the PHOENICS Input Library.
The 2D-store block-location indices that are recognized as arguments of MAKE are; XG2D, XU2D, DXG2D, DXU2D, YG2D, YV2D, DYG2D, DYV2D, RG2D and RV2D.
This group is called at the start of each time step when the user has set TLAST=GRND in PIL, to instruct EARTH to visit this group of GROUND for a setting of the time step DT.
Normally the time-step is fixed by the TFRAC setting, made in PIL, but this option povides the opportunity for DT to be made a function of other calculated quantities, or of some complicated distribution awkward to set in PIL.
It should be noted that the number of time steps is still fixed by TFRAC(1). Thus the PIL commands:
TFRAC(1)=-250.0; TFRAC(2)=1.0;TLAST=GRND
will instruct EARTH to perform 250 time steps, the size of each being determined by the setting of DT in group 2.
When control returns to EARTH, it sets TIM=TIM+DT.
This group is CALLed at the start of the current z-slab at the stage when the geometry is being calculated. It is useful in parabolic calculations (PARAB=T) to expand the x-extent of the grid.
The EARTH CALL to this group is activated when AZXU=GRND. The user's job is to set the variable XRAT which is the ratio of the x-extent at the current z-slab to that of the previous slab.
Thereafter EARTH modifies all geometrical entities (eg. internodal distances, cell-face areas, XULAST etc) to conform with XFRAC.
Typically, it is desired to make the domain width a function of the downstream distance of the current slab (ZW). For example,
XRAT=(RG(1)+RG(2)*ZW+RG(3)*ZW**2)/XULAST
gives a quadratic dependence. The RGs are PIL-set parameters.
More often it might be desired to set XRAT to vary with z so as to accommodate a boundary layer, which would necessitate local determination of the width of the layer in order that edge velocities should be sufficiently close to the free-stream values.
It should be noted that group 3 is called after group 5 (where ZW is set), so that in group 3 ZW is the current-slab z (ZWL stores the z-location of the low face of the current z-slab, ie the ZW of the previous slab).
Group 4 performs an exactly corresponding function for the y- direction domain width, YVLAST, to that performed by group 3 for XULAST. EARTH visits the group when AZYV is set to GRND in PIL.
This group is visited by EARTH when AZDZ has been set .LE. GRND in SATELLITE; it sets the step size DZ in parabolic calculations. Two examples are provided in GREX3 (activated by AZDZ=GRND1 or AZDZ=GRND2).
They make the step size a fixed fraction (RSG10) of the calculation domain widths XULAST and YVLAST respectively. The user would need to provide coding here if the options provided in GREX3 were inadequate for his needs.
This group is called from EARTH when UGEOM=T, at the start of the current IZ slab, just after the current-slab geometry (eg inter- nodal distances, cell face areas etc) has been set. It therefore provides the possibility of amending these geometrical entities.
In non-BFC calculations, geometrical modifications are also often performed, by amendment to the porosity factors.
When BFC equals T, however, the geometrical quantities are calculated and stored whole-field at the start of the calculations; they are multiplied by the porosity factors there, once for all. Thus, this group provides opportunities for amendment of the geometrical quantities once BFC equals F.
GREX3 contains no example of the use of this Group, but a CALL to GXPORA from GROUP 19 has a similar effect.
This Group is not entered
This group contains 15 sections. Provision has been made in Group 8 of GROUND for extensive intervention in the procedures of formulation and solution of the finite-domain equations which are central to the flow-simulation process of PHOENICS.
Whoever wishes to use this will find helpful advice in the comments which appear in group 8 of GREX3, and specifically in sections 8 to 14.
These comments explain what must be done to ensure that the relevant sections of GREX3 (and GROUND, if USEGRD=T) are to be entered: logical variables such as UCONV, UDIFF etc are to be set .TRUE. in the SATELLITE.
These variables allow the user to set, in accordance with whatever rule he prefers: the convective coefficients, the diffusive coefficients, the convection neighbour values, the diffusion neighbour values, the two components of the linearised source terms, any of the coefficients in the variable-correction equation, the solution procedure, or the results of that procedure.
Values of the coefficients, neighbour values, etc will already have been set by EARTH in accordance with the standard prescription at the time at which control passes to GROUND. What these values are can be determined by the activation of the calls to the PRN subroutine which have been, deactivated by C's, in the first column (see the listing of GREX3).
Even if the user decides he does not want to activate these statements, he will find the way in which the various quantities are referred to be instructive about how the functions L0F and L0FZ should be used in arguments of other functions, for example FN0, when values are to be restored to EARTH.
This group originally contained 13 sections, all of which were concerned with setting material properties or other auxiliary quantities used in calculations.
Nowadays these calls are made directly from withn EARTH; so this group is empty.
This group contains four sections, all of which are concerned with setting quantities that determine the intensity of inter- phase transport. Thus section 1 sets the coefficient of inter- phase diffusion of momentum, ie friction, which is also used after multiplication by CINT for diffusion of other variables when their CINTs are not equal to GRND.
Section 2 sets the inter-phase convection, ie the mass-transfer rate. Sections 3 and 4 are used to set phase-to-interface diffusion transfer coefficients, normally for non-velocity variables (an isotropic velocity-transfer coefficient being set in section 1).
The reader is advised at this stage to study GREX3. The listing shows how PIL is used to instruct EARTH to visit each section (eg CMDOT.LE.GRND) for the interphase friction coefficient, and the index (or index function) which permits GROUND to refer to the storage location of the quantity in EARTH.
The coefficients of the diffusive transfer of momentum must be set in accordance with its definition:
(Total interphase friction force for the cell)
(Velocity difference between the two phases).
The coefficient of the diffusive transport of the other PHI variables must be in accordance with the definition:
(Diffusive flux of PHI to interface for the cell)
(PHI difference between the bulk of the phase and its ) ( interface value PHINT ) .
The first two sections are called at the start of the hydrodynamic iteration of the current IZ slab. When visited, EARTH expects GROUND to return an array of values in the F-array segment address locations determined by the indices INTFRC and INTMDT
As an example, consider an inter-phase friction coefficient equal to a constant (CFIP1A) times the in-cell mass of the first-phase times the in-cell volume fraction of the second phase. The following statement inserted in section 1 effects this dependence:
CALL FN21(INTFRC,MASS1,R2,0.0,CFIP1A)
The subroutine FN21(y,x1,x2,a,b) has the following mathematical significance:
y=a+b*x1*x2
For each cell in the current IZ slab, the index MASS1 refers to the EARTH store of the mass of phase one in each cell at the current slab.
An equivalent but more understandable sequence is affected by (see also previous section):
L0FIP=L0F(INTFRC)
L0MAS=L0F(MASS1)
L0R2=L0F(R2)
NXNY=NX*NY
DO 109 I=1,NXNY
109 F(L0FIP+I)=CFIP1A*F(L0MAS+I)*F(L0R2+I)
In the above sequence, the first three statements locate: the zero-locations in the F-array of:- the friction coefficients; the mass of phase 1 in the cell; and the volume fraction of the second phase.
Another technique is to use the subroutine GETYX to get, and store locally, the arrays for MASS1 and R2; to calculate and store the required result in the array GFIP; and to set the EARTH array to the data contained in it. An example now follows
CALL GETYX(MASS1,GM1,NYDIM,NXDIM)
CALL GETYX(R2,GR2,NYDIM,NXDIM)
DO 102 IX=1,NX
DO 102 IY=1,NY
102 GFIP(IY,IX)=CFIP1A*GM1(IY,IX)*GR2(IY,IX)
CALL SETYX(INTFRC,GFIP,NYDIM,NXDIM)
The arrays GFIP, GM1 and GR2 are dimensioned to NYDIM, NXDIM (which must be greater than NY, NX respectively ) at the top of the subroutine. This technique although familiar to users of PHOENICS 81, is not recommended, because the first two methods are more economical.
Sections 3 and 4 are CALLed when the interphase terms are being assembled for variable PHI, when CINT(PHI) is less than equal to GRND. Thus, different formulae can be supplied for different dependent variables. EARTH sets the GROUND variable INDVAR so that the FORTRAN can distinguish one PHI entry from another. The indices CO1I and CO2I permit the F-array zero-locations to be deduced as L0F(CO1I) and L0F(CO2I).
GREX3 supplies options for CMDOT, CFIPS and CINT(INDVAR), which are described under their respective entries in chapter 2.
This group is for setting non-uniform initial conditions for variables that are stored whole-field. This group is visited by EARTH for variables for which the fourth argument of INIT (set in the SATELLITE) is GRND (see INIT for background information).
EARTH visits group 11 for the field values of variable INDVAR over the current patch at the current IZ step. The field values are to be set in EARTH at a segment address located by means of the integer index VAL. The following example will clarify what has to be done.
Suppose that it is desired to initialise the w-velocity field to a parabolic profile over the last half of the domain: the PIL commands
PATCH(LASTHALF,INIVAL,1,NX,1,NY,NZ/2,NZ,1,1)
INIT(LASTHALF,W1,0.0,GRND)
instruct EARTH to visit group 11 of GROUND for an array of values for the field W1 at each slab within the sub-domain indicated by arguments 3 to 8 of PATCH. Prior to calling group 11, EARTH sets NPATCH(1:8)=LASTHALF, INDVAR=W1, IXF=1, IXL=NX, IYF=1, IYL=NY; and IZ contains the current z-slab that EARTH is considering.
The following coding does what is needed in group 11:
IF(NPATCH.EQ.'LASTHALF') THEN
IF(INDVAR.EQ.W1) THEN
CALL FN4(VAL,YG2D,RG(1),RG(2),RG(3))
ENDIF
ENDIF
RETURN
The subroutine FN4(y,x,a,b,c) has the mathematical significance:
y=a+b*x+c*x**2
for each cell in the current PATCH at IZ. The index YG2D refers to the EARTH array of length NX*NY elements that contains the y- coordinate of the cell centres at the slab. It should be noted that RG(1), RG(2) and RG(3) are PIL parameters.
An equivalent but more transparent sequence is effected by first determining the segment address of VAL and YG2D and then providing a DO loop that sets the field directly:
L0VAL=L0F(VAL) L0Y=L0F(YG2D) DO 111 IX=IXF,IXL DO 111 IY=IYF,IYL ICELL=IY+(IX-1)*NY XX=F(L0Y+ICELL)**2 111 F(L0VAL+ICELL)=RG(1)+RG(2)*F(L0Y+ICELL) 1+RG(3)*XX
Yet another technique to achieve the same is to use the subroutine GETYX to get and store locally the array YG2D, to overwrite this local array with the array of values required, and then to set the EARTH store of VAL equal to the data in this array, thus:
CALL GETYX(YG2D,GY,NYDIM,NXDIM)
DO 111 IX=IXF,IXL
DO 111 IY=IYF,IYL
XX=GY(IY,IX)**2
111 GY(IY,IX)=RG(1)+RG(2)*GY(IY,IX)+RG(3)*XX
CALL SETYX(VAL,GY,NYDIM,NXDIM)
where GY is an array dimensioned to NYDIM, NXDIM at the top of subroutine GROUND ( NYDIM, NXDIM must be geater than or equal to NY, NX respectively ).
Group 11 can be visited for any number of variables for a given PATCH for which non-uniform fields are wanted. Any number of PATCHes may be used. The user must use the parameters NPATCH, INDVAR and IZ to distinguish one patch from another, one variable from another, and one slab from another.
Group 13 of GROUND is the place where the user can provide non- linear sources and boundary-condition information for PATCHes of the domain for variables identified by COVAL.
The PIL instructions
PATCH(name,type,ixf,ixl,...,itf,itl)
COVAL(name,PHI,GRND,...)
instruct EARTH to visit group 13 of GROUND for an array of 'coefficients' for each z-slab indicated for the variable PHI. The index L0F(CO) gives the zero location of the appropriate segment of the F-array into which the coefficients must be put.
The PIL instructions
PATCH(name,type,ixf,ixl,...,itf,itl)
COVAL(name,PHI,...,GRND)
correspondingly instruct EARTH to visit group 13 of GROUND for an array of 'values'.
The PIL instructions
PATCH(name,type,ixf,ixl,...,itf,itl)
COVAL(name,PHI,GRND,GRND)
instruct EARTH to visit group 13 once for an array of COefficients and again for an array of VALues.
The PIL instructions
PATCH(name,type,ixf,ixl,...,itf,itl)
COVAL(name,PHI,GRND,GRND)
COVAL(name,PHIA,0,GRND)
COVAL(name,PHIB,GRND,0)
causes EARTH to visit groups 13 for COs and VALs for variable PHI, VALs for variable PHIA and COs for variable PHIB. In this case, there are four CALLs from EARTH for the PATCH in question for each slab IZ in the range of arguments 7 and 8 in PATCH. At each slab the COefficient and VALue arrays need to be set over the extent ixf to ixl, iyf to iyl; ie (ixl-ixf+1)*(iyl-iyf+1) values are to be set.
Before EARTH calls group 13, it sets NPATCH (character*8), IZ and INDVAR to: the current PATCH name, the current IZ considered, and the current variable in question respectively. Reference to these variables in the FORTRAN coding can distinguish between the possibilities selected.
Group 13 of GROUND is subdivided into 22 sections. The first 11 sections are provided for the setting of COefficient-array options and the last 11 are provided for setting VALue-array options, according to the following scheme.
ISC=1 Section 1 COVAL(name,PHI,GRND,...)
ISC=2 Section 2 COVAL(name,PHI,GRND1,...)
ISC=3 Section 3 COVAL(name,PHI,GRND2,...)
ISC=4 Section 4 COVAL(name,PHI,GRND3,...)
ISC=5 Section 5 COVAL(name,PHI,GRND4,...)
ISC=6 Section 6 COVAL(name,PHI,GRND5,...)
ISC=7 Section 7 COVAL(name,PHI,GRND6,...)
ISC=8 Section 8 COVAL(name,PHI,GRND7,...)
ISC=9 Section 9 COVAL(name,PHI,GRND8,...)
ISC=10 Section 10 COVAL(name,PHI,GRND9,...)
ISC=11 Section 11 COVAL(name,PHI,GRND10,...)
ISC=12 Section 12 COVAL(name,PHI,...,GRND) ISC=13 Section 13 COVAL(name,PHI,...,GRND1) ISC=14 Section 14 COVAL(name,PHI,...,GRND2) ISC=15 Section 15 COVAL(name,PHI,...,GRND3) ISC=16 Section 16 COVAL(name,PHI,...,GRND4) ISC=17 Section 17 COVAL(name,PHI,...,GRND5) ISC=18 Section 18 COVAL(name,PHI,...,GRND6) ISC=19 Section 19 COVAL(name,PHI,...,GRND7) ISC=20 Section 20 COVAL(name,PHI,...,GRND8) ISC=21 Section 21 COVAL(name,PHI,...,GRND9) ISC=22 Section 22 COVAL(name,PHI,...,GRND10)
In any particular section, coefficient (or value) arrays can be set for any number of different PHIs for which there are COVALs.
The possibilities offered by this feature are heavily exploited in GREX3, where an additional degree of freedom is added by the recognition of special PATCH names. The beginner to GROUND is advised to use GRND only, and hence to set his coefficients and values in section 1 and section 12 of group 13 respectively.
In what follows, two examples are provided of group 13 coding. The reader can examine group 13 of GREX3 and the GX library subroutines called from there for further examples.
The first example is one in which a known internal heat source per unit volume is present over a restricted portion of the domain, namely at IX=3 to 7, IY=2 to 20 and IZ=3 to 6. Suppose that this heat source per unit volume, q, is known to vary with position as follows:
q = ax + by + cz
This is not a non-linear source for it does not depend upon any solved quantity. It is however non-uniform, and without the possibility of GROUND coding would necessitate (7-3+1)*(20-2+1)*(6-3+1) = 580 separate PATCHs to set the heat flux in the cells covered. The following PIL commands instruct EARTH to look for an array of values of the PATCH instead.
PATCH(HEAT,VOLUME,3,7,2,20,3,6,5,10)
COVAL(HEAT,H1,FIXFLU,GRND)
The last two arguments of PATCH dictate that the source will be applied only during time steps 5 to 10 inclusive. The coding in group 13 to effect the formula is as follows:
IF(NPATCH(1:4).EQ.'HEAT') THEN
IF(INDVAR.EQ.H1) THEN
CALL FN10(VAL,XG2D,YG2D,RG(3)*ZW,RG(1),RG(2))
ENDIF
ENDIF
RETURN
The subroutine FN10(y,x1,x2,a,b,c) has the following mathematical significance:
y=a+b*x1+c*x2
for all cells in the zone IXF to IXL, IYF to IYL at the current slab. RG(1), RG(2) and RG(3) are PIL variables representing p, q, and r. (It should be noted that VAL is not multiplied by 1.0E10 i.e 1.0/FIXFLU: this is done automatically in EARTH.)
An equivalent, but more transparent sequence, is effected by first determining the segment address for the EARTH arrays referred to by VAL, XG2D, and then by setting VAL directly for all cells required:
L0VAL=L0F(VAL)
L0x=L0F(XG2D)
L0Y=L0F(YG2D)
DO 131 IX=IXF,IXL
DO 131 IY=IYF,IYL
ICELL=IY+(IX-1)*NY
131 F(L0VAL+ICELL)=RG(1)*F(L0X+ICELL)
1+RG(2)*F(L0Y+ICELL)+RG(3)*ZW
Yet another technique to achieve the same end is to use the subroutine GETYX to get and store locally the groups referred to by XG2D and YG2D, to fill a local array with required values of the flux and then to fill the F-array segment address indicated by VAL to this local array:
CALL GETYX(XG2D,GX,NYDIM,NXDIM)
CALL GETYX(YG2D,GY,NYDIM,NXDIM)
DO 131 IX=IXF,IXL
DO 131 IY=IYF,IYL
131 GVAL(IY,IX)=RG(1)*GX(IY,IX)+RG(2)*GY(IY,IX)+RG(3)*ZW
CALL SETYX(VAL,GVAL,NYDIM,NXDIM)
The arrays GX, GY and GVAL must be dimensioned to NXDIM,NYDIM at the top of subroutine GROUND ( NYDIM, NXDIM must be greater than or equal to NY, NX respectively ). Although this method will be familiar to users of PHOENICS 81, it is not recommended, for the other two methods are much more economical as they work directly with the EARTH stores, obviating the additional step of GETting and SETting
The second example concerns the representation of a parabolic profile, in the y-direction, of W1 at the inlet plane IZ=1 by means of sources in group 13. The PIL statements required are:
PATCH(INLET,LOW,1,NX,1,NY,1,1,1,1)
COVAL(INLET,P1,FIXFLU,GRND)
COVAL(INLET,W1,ONLYMS,GRND)
The group 13 coding to bring this about is as follows:
IF(NPATCH(1:5).EQ.'INLET') THEN
IF(INDVAR.EQ.P1.OR.INDVAR.EQ.W1) THEN
RECRHO=1.0
IF(INDVAR.EQ.W1) RECRHO=1./RHO1
CALL FN4(VAL,YG2D,RG(1)*RECRHO,RG(2)*RECRHO,RG(3)*RECRHO)
ENDIF
ENDIF
RETURN
SUBROUTINE FN4(y,x,a,b,c) has the mathematical significance:
y=a+b*x+c*x**2.
For each cell in the PATCH at IZ=1, the integer name YG2D refers to the EARTH array of length NX*NY elements that contains the y- coordinate of the cell centres at the slab. When INDVAR=P1, the mass flow rate per unit area is fixed; when INDVAR=W1 this is divided by density to give the w momentum convected in by this mass flow rate.
The third example concerns the representation by a sink of momentum of a resistive force on the U1 velocity (due for example to the presence of stationary baffles). The momentum form of the force is:
-a*u**b.
(It should be noted that, if b=2, the source is quadratic and can be represented from PIL; see COVAL for details). In linearised form, this source can be represented by a coefficient of a*u1**(b-1) and a value of zero. Thus the PIL statements are:
PATCH(BAFFLE,PHASEM,1,NX,1,NY,1,NZ,1,1)
COVAL(BAFFLE,U1,GRND,0.0)
The GROUND coding that represents this source in group 13 is as follows:
IF(NPATCH(1:6).EQ.'BAFFLE') THEN
IF(INDVAR.EQ.U1) THEN
L0U=L0F(U1)
L0CO=L0F(CO)
NXNY=NX*NY
GB=RG(2)-1.0
DO 131 I=1,NXNY
XX=F(L0U+I)**GB
131 F(L0CO+I)=RG(1)*XX
ENDIF
ENDIF
RETURN
The first move is to obtain the zero-locations of the EARTH stores for the coefficient and the U1 velocity (at the current IZ). The final move is to set the EARTH coefficient array directly.
This could also have been done by a CALL to FN8, or by means of local arrays filled by GETYX and returned to EARTH by SETYX. Examples earlier in this group, and in groups 9, 10 and 11, explain how these alternatives work.
This section is used in calculations of unconfined parabolic flows (IPARAB=1) in which the user wishes to impose a finite z- direction pressure gradient, corresponding to that in the free stream.
The section is accessed, if the user has set AZPH=GRND in Q1, at the start of work at each forward step. It is the job of sequences in group 14 to return in PBAR the value of the imposed pressure that applies at the next forward slab, ie at IZ+1.
The boundary conditions on P1 at the edge of the current IZ must be consistent with the prescribed variation of PBAR with z.
For three-dimensional boundary-layer flows, the imposed free-stream pressure may be a function of x (or y), in addition to being a function of z. The current version of PHOENICS does not have this possibility as a built-in option.
However it is possible to handle this case, at least for incompressible flows, as follows:
IF(NPATCH(1:3).EQ.'PXZ') THEN
IF(INDVAR.EQ.'W1') THEN
L0VAL=L0F(VAL)
DO 131 IX=1,NX
DO 131 IY=1,NY
ICELL=(IX-1)*NY+IY
131 F(L0VAL+ICELL)=-DPDZ(IX)
ENDIF
ENDIF
RETURN
where DPDZ(IX) is a user-prescribed one-dimensional array which specifies the imposed pressure gradient as a function of IX. The user would reset this at each forward z-step, in section 3 of group 19 of GROUND, so as to prescribe the z-wise variation of the imposed pressure gradient.
There are 10 sections in this group, called by EARTH at the stages of the solution sequence indicated in the comments on the listing.
Examples of the use of these sections may be found in of GREX3. Of especial interest are:
The call to GXPIST and hence ZMOVE in section 1 (ie at the start of the time step) for the definition of the moving grid.
The modification of the porosities at the start of the of the z- slab as a function of the pressure (to represent free surface and the effect of compliant walls etc) (the in-EARTH multiplication of the nominal areas by the porosity fractions comes after this entry into section 3).
The call to GXGENK for the calculation of the strain rate expression used in the generation of KE sources and for the source of enthalpy for the dissipation of mechanical energy into heat. In this last respect the user should recognise that he has to access this quantity which he may want to amend near walls.
This group is called at the beginning of the run, after the data have been read in and modified by EARTH to ensure consistency. It is from this group in subroutine GREX3 that DATPRN is called to ECHO the SATELLITE data.
wbs