TEXT( Law-Of-The-Wall Calculator 
REAL(UPLUS,KAY,KUPLUS,YPLUS,EWAL,NUPLUS,LRNT,RE,TEMPP,FACTOR)
REAL(LRNMFAC)
KAY=0.41;EWAL=8.6;enuta=0.05;enutb=4.0
nowipe=t
  DISPLAY
 
  This Q1 computes the various quantities which are deducible from
  the supposition that the "universal law of the wall" can be
  expressed by the single formula:
 
  yplus = uplus + (1/E)*[exp(K*uplus) - 1 - K*uplus/2 -
                                (K*uplus)**2/6 - (K*uplus)**3/24]
 
  (Ref. Spalding DB, Trans ASME, J Appl Mech vol 28 pp 455-458,
                                                             1961)
 
  The quantities in question are:-
  nuplus, the turbulent contribution to the viscosity divided by
          the laminar viscosity;
  re,     the Reynolds Number based on distance from the wall,
          = uplus * yplus;
  lrnt,   the local Reynolds number of turbulence, defined as
          K*yplus, which is the value which nuplus takes at high
          Reynolds number, and which can be regarded as the nominal
          value of nuplus;
  factor, the factor by which the nominal value of nuplus should be
          multiplied to yield the actual value.
 
#$13)
 
  Also computed is LRNMFAC, the local Reynolds Number Multiplication
  Factor computed by PHOENICS, when IENUT is set equal to 6 .
 
  The values used here for K and E are 0.4 and 9.0 respectively;
  and, to calculate the LRNM factor, the values used for ENUTA and
  ENUTB are 0.05 and 4.0 respectively.
 
 
  This Q1 file illustrates how PIL can be used as a general utility
  having nothing directly to do with PHOENICS.
 
 
  Please press RETURN to continue
  readvdu(nx,int,nx)
#$13)
LABEL TOP
 
 
  What uplus do you want, please?
READVDU(UPLUS,REAL,0.0)
KUPLUS=KAY*UPLUS
TEMPP=EXP(KUPLUS)-1 -KUPLUS*(1 + 0.5*KUPLUS + 0.16667*KUPLUS**2)
YPLUS=UPLUS + (TEMPP-KUPLUS**4/24)/EWAL; NUPLUS=KAY*(TEMPP)/EWAL
LRNT=KAY*YPLUS; FACTOR=NUPLUS/LRNT
IF(FACTOR.GT.1) THEN
 FACTOR=1.0+.00001
ENDIF
RE=UPLUS*YPLUS
LRNMFAC=(ENUTA*LRNT)**ENUTB
IF(LRNMFAC.GT.1) THEN
 LRNMFAC=1.0+.00001
ENDIF
UPLUS=UPLUS+.00001
mesg(uplus         yplus         re            nuplus
mesgb(:uplus:    :yplus:    :re:    :nuplus:
mesg(lrnt          factor        lrnmfac
mesgb(:lrnt:    :factor:    :lrnmfac:
  more? (y/n)
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y)THEN
 GOTO TOP
ENDIF
ABORT
  ENDDIS