GROUP 1. Run title and other preliminaries
BOOLEAN(TURBL)
INTEGER(NXF,NXL,NYF,NYL,NZF,NZL,ISUM)
REAL(VIN,VFRAC,TIN,CIN)
VIN=10.;TIN=100.;CIN=-1.
IF (EGWF) THEN
+  CASE :CTURB: OF
+  WHEN TURB,4
+     TEXT( TURBULENT TEST OF EGWF - EGWF VARIANT)
+  WHEN LAM,3
+     TEXT( LAMINAR TEST OF EGWF - EGWF VARIANT)
+  ORELSE
+     TEXT( :CTURB: TEST OF EGWF - EGWF VARIANT)
+     TURBL = T
+  ENDCASE
ELSE
+  CASE :CTURB: OF
+  WHEN TURB,4
+     TEXT( TURBULENT TEST OF EGWF - SOWF VARIANT)
+  WHEN LAM,3
+     TEXT( LAMINAR TEST OF EGWF - SOWF VARIANT)
+  ORELSE
+     TEXT( :CTURB: TEST OF EGWF - SOWF VARIANT)
+     TURBL = T
+  ENDCASE
ENDIF
IF (NX .LE. 0) THEN
+  NX=1
ENDIF
IF (NY .LE. 0) THEN
+  NY=1
ENDIF
IF (NZ .LE. 0) THEN
+  NZ=1
ENDIF
IF (N1 .LE. 0) THEN
+  N1=1
ENDIF
 
    GROUP 3. X-direction grid specification
GRDPWR(X,-NX,1.0,2.0)
NXF=1+N1;NXL=NX-N1
    GROUP 4. Y-direction grid specification
GRDPWR(Y,NY,10.0,1.0)
    GROUP 5. Z-direction grid specification
GRDPWR(Z,-NZ,1.0,2.0)
NZF=1+N1;NZL=NZ-N1
    GROUP 7. Variables stored, solved & named
SOLVE(P1);SOLUTN(P1,Y,Y,Y,P,P,Y)
SOLVE(TEM1);SOLUTN(TEM1,Y,Y,Y,P,P,Y)
SOLUTN(C1,Y,Y,Y,P,P,Y)
ISUM=0
IF (NX .GT. 1) THEN
+  SOLVE(U1);ISUM=ISUM+1
+  SOLUTN(U1,P,P,Y,P,P,Y)
ENDIF
IF (NY .GT. 1) THEN
+  SOLVE(V1);ISUM=ISUM+1
+  SOLUTN(V1,P,P,Y,P,P,Y)
ENDIF
IF (NZ .GT. 1) THEN
+  SOLVE(W1);ISUM=ISUM+1
+  SOLUTN(W1,P,P,Y,P,P,Y)
ENDIF
VFRAC=ISUM**0.5
STORE(PRPS)
IF (ISOLBK .NE. 0) THEN
+  STORE(BLOK)
+  IVARBK=-1
ENDIF
IF (TURBL) THEN
+  TURMOD(:CTURB:)
+  SOLUTN(KE,P,P,p,P,P,Y)
+  SOLUTN(EP,P,P,p,P,P,Y)
   ELSE IF (:CTURB: .EQ. TURB) THEN
   +  STORE(LEN1,VIST)
ENDIF
    GROUP 8. Terms (in differential equations) & devices
   *** Get rid of potentially confusing heat sources that are
   *** irrelevant to the tests being performed
TERMS(TEM1,N,P,P,P,P,P)
    GROUP 9. Properties of the medium (or media)
   *** Set up the turbulence models for the turbulent cases
   *** NB. at present limited to an artificially simple mixing
   *** length model
IF (:CTURB: .EQ. TURB .OR. :CTURB: .EQ. KLMODL) THEN
+  IF (NX.GT.1) THEN
+     EL1=LINEARX
+  ELSE
+     EL1=LINEARY
+  ENDIF
+  EL1A=(XULAST**2+YVLAST**2+ZWLAST**2)**0.5;EL1B=0.0
ENDIF
IF (:CTURB: .EQ. TURB) THEN
+  ENUT=PROPLEN;ENUTA=0.0;ENUTB=0.1
ENDIF
    GROUP 11. Initialization of variable or porosity fields
FIINIT(U1)=0.0;FIINIT(V1)=VIN;FIINIT(W1)=0.0
FIINIT(PRPS)=11;PRNDTL(TEM1)=CONDFILE;ENUL=FILE;CP1=FILE
FIINIT(TEM1)=10.;FIINIT(C1)=0.0
 
IF (.NOT. EGWF) THEN
   *** An unsuccesful attempt to bamboozle CONPOR
+  IF (NX .GT. 1) THEN
+     CONPOR(BLCKACP,-1,CELL,1,-(NXF-1),1,NY,1,NZ)
+     CONPOR(BLCKBCP,-1,CELL,-(NXL+1),NX,1,NY,1,NZ)
+  ENDIF
+  IF (NZ .GT. 1) THEN
+     CONPOR(BLCKCCP,-1,CELL,1,NX,1,NY,1,-(NZF-1))
+     CONPOR(BLCKDCP,-1,CELL,1,NX,1,NY,-(NZL+1),NZ)
+  ENDIF
ENDIF
 
CASE :CTURB: OF
WHEN KEMODL,6
+   FIINIT(KE) = 0.01*VIN*VIN
+   FIINIT(EP) = FIINIT(KE)**1.5
WHEN KLMODL,6
+   FIINIT(KE) = 0.01*VIN*VIN
ENDCASE
 
INIADD=F
IF (NX .GT. 1) THEN
+  PATCH(BLCKA ,CELL,1,NXF-1,1,NY,1,NZ,1,1)
+  PATCH(BLCKAI,INIVAL,1,NXF-1,1,NY,1,NZ,1,1)
+  PATCH(BLCKB ,CELL,NXL+1,NX,1,NY,1,NZ,1,1)
+  PATCH(BLCKBI,INIVAL,NXL+1,NX,1,NY,1,NZ,1,1)
+  INIT (BLCKAI,PRPS,0.0,103.);INIT (BLCKBI,PRPS,0.0,103.)
ENDIF
IF (NZ .GT. 1) THEN
+  PATCH(BLCKC ,CELL,1,NX,1,NY,1,NZF-1,1,1)
+  PATCH(BLCKCI,INIVAL,1,NX,1,NY,1,NZF-1,1,1)
+  PATCH(BLCKD ,CELL,1,NX,1,NY,NZL+1,NZ,1,1)
+  PATCH(BLCKDI,INIVAL,1,NX,1,NY,NZL+1,NZ,1,1)
+  INIT (BLCKCI,PRPS,0.0,103.);INIT (BLCKDI,PRPS,0.0,103.)
ENDIF
 
IF (ISOLBK .NE. 0) THEN
+  INTEGER(IBL);IBL=1
+  FIINIT(BLOK)=IBL
+  IF (NX .GT. 1) THEN
+     IBL=IBL+1;INIT(BLCKAI,BLOK,0.0,IBL)
+     IBL=IBL+1;INIT(BLCKBI,BLOK,0.0,IBL)
+  ENDIF
+  IF (NZ .GT. 1) THEN
+     IBL=IBL+1;INIT(BLCKCI,BLOK,0.0,IBL)
+     IBL=IBL+1;INIT(BLCKDI,BLOK,0.0,IBL)
+  ENDIF
ENDIF
    GROUP 13. Boundary conditions and special sources
   *** Set the C1 and TEM1 sbources in the blocks: fixed
   *** flux for heat(??), and fixed values of C1 (heated
   *** rock-salt dissolving into water??)
IF (NX .GT. 1) THEN
+  COVAL(BLCKA,C1,FIXVAL,-CIN)
+  COVAL(BLCKB,C1,FIXVAL,-CIN)
+  COVAL(BLCKA,TEM1,FIXFLU,20.)
+  COVAL(BLCKB,TEM1,FIXFLU,20.)
ENDIF
IF (NZ .GT. 1) THEN
+  COVAL(BLCKC,C1,FIXVAL,-CIN)
+  COVAL(BLCKD,C1,FIXVAL,-CIN)
+  COVAL(BLCKC,TEM1,FIXFLU,20.)
+  COVAL(BLCKD,TEM1,FIXFLU,20.)
ENDIF
   *** Setup inlet conditions
IF (NY .GT. 1) THEN
+  INLET(YINL,SOUTH,NXF,NXL,1,1,nzf,nzl,1,LSTEP)
+  VALUE(YINL,P1,RHO1*VIN)
+  VALUE(YINL,V1,VIN)
+  VALUE(YINL,TEM1,TIN);VALUE(YINL,C1,CIN)
+  IF (FIINIT(KE) .GT. 1.E-10) THEN
+     VALUE(YINL,KE,FIINIT(KE))
+  ENDIF
+  IF (FIINIT(EP) .GT. 1.E-10) THEN
+     VALUE(YINL,EP,FIINIT(EP))
+  ENDIF
+  OUTLET(YOUTL,NORTH,NXF,NXL,NY,NY,nzf,nzl,1,LSTEP)
+  VALUE(YOUTL,P1,0.0)
ENDIF
   *** For the SOWF method set up the WALL PATCHes
IF (.NOT. EGWF) THEN
+  IF (NX .GT. 1) THEN
      WALL(WALLA,WEST,NXF,NXF,1,NY,1,NZ,1,LSTEP)
      COVAL(WALLA,V1,WALLCO,OPPVAL)
      COVAL(WALLA,TEM1,WALLCO,OPPVAL);COVAL(WALLA,C1,WALLCO,OPPVAL)
      WALL(WALLB,EAST,NXL,NXL,1,NY,1,NZ,1,LSTEP)
      COVAL(WALLB,V1,WALLCO,OPPVAL)
      COVAL(WALLB,TEM1,WALLCO,OPPVAL);COVAL(WALLB,C1,WALLCO,OPPVAL)
      INIT(BLCKA,EPOR,0.0,1.0);INIT(BLCKB,EPOR,0.0,1.0)
+  ENDIF
+  IF (NZ .GT. 1) THEN
      WALL(WALLC,LOW ,1,NX,1,NY,NZF,NZF,1,LSTEP)
      COVAL(WALLC,V1,WALLCO,OPPVAL)
      COVAL(WALLC,TEM1,WALLCO,OPPVAL);COVAL(WALLC,C1,WALLCO,OPPVAL)
      WALL(WALLD,HIGH,1,NX,1,NY,NZL,NZL,1,LSTEP)
      COVAL(WALLD,V1,WALLCO,OPPVAL)
      COVAL(WALLD,TEM1,WALLCO,OPPVAL);COVAL(WALLD,C1,WALLCO,OPPVAL)
      INIT(BLCKC,EPOR,0.0,1.0);INIT(BLCKD,EPOR,0.0,1.0)
+  ENDIF
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP=50
    GROUP 16. Termination of iterations
LITER(TEM1)=20;LITER(C1)=20;LITER(U1)=10;LITER(V1)=10
    GROUP 17. Under-relaxation devices
RELAX(U1,FALSDT,5.*XULAST/(NX*VIN))
RELAX(V1,FALSDT,5.*YVLAST/(NY*VIN))
RELAX(W1,FALSDT,5.*ZWLAST/(NZ*VIN))
RELAX(KE,FALSDT,5.*YVLAST/(NY*VIN))
RELAX(EP,FALSDT,5.*YVLAST/(NY*VIN))
RELAX(TEM1,FALSDT,1000.)
RELAX(C1,  FALSDT,1000.)
    GROUP 21. Print-out of variables
OUTPUT(P1,Y,P,P,Y,Y,Y);OUTPUT(U1,Y,P,P,Y,Y,Y)
OUTPUT(V1,Y,P,P,Y,Y,Y);OUTPUT(W1,Y,P,P,Y,Y,Y)
OUTPUT(TEM1,Y,P,P,Y,Y,Y);OUTPUT(C1,Y,P,P,Y,Y,Y)
IF (:CTURB: .EQ. KEMODL .OR. :CTURB: .EQ. KLMODL) THEN
+  OUTPUT(KE,Y,P,P,Y,Y,Y)
ENDIF
    GROUP 22. Spot-value print-out
IXMON=NXF;IYMON=NY/2+1;IZMON=NZ/2+1
TSTSWP=-1
    GROUP 23. Field print-out and plot control
NXPRIN=1;NYPRIN=1;NZPRIN=1;IYPRF=NY-2;IYPRL=NY
IF (NX .GT. 1 .AND. NY .GT. 1) THEN
+  PATCH(CONZ,CONTUR,1,NX,1,NY,NZ/2+1,NZ/2+1,1,1)
+  PLOT (CONZ,P1,0.0,20.)
+  PLOT (CONZ,TEM1,0.0,20.)
+  PLOT (CONZ,C1,0.0,20.)
ENDIF
IF (NY .GT. 1 .AND. NZ .GT. 1) THEN
+  PATCH(CONX,CONTUR,NX/2+1,NX/2+1,1,NY,1,NZ,1,1)
+  PLOT (CONX,P1,0.0,20.)
+  PLOT (CONX,TEM1,0.0,20.)
+  PLOT (CONX,C1,0.0,20.)
ENDIF
IF (NZ .GT. 1 .AND. NX.GT. 1) THEN
+  PATCH(CONY,CONTUR,1,NX,NY/2+1,NY/2+1,1,NZ,1,1)
+  PLOT (CONY,P1,0.0,20.)
+  PLOT (CONY,TEM1,0.0,20.)
+  PLOT (CONY,C1,0.0,20.)
ENDIF
IF (NX .GT. 1) THEN
+  PATCH(PRFX,PROFIL,1,NX,NY/2+1,NY/2+1,NZ/2+1,NZ/2+1,1,1)
+  PLOT (PRFX,TEM1,100.0,1000.)
+  PLOT (PRFX,C1,-1.,1.)
ENDIF
IF (NY .GT. 1) THEN
+  PATCH(PRFY,PROFIL,NX/2+1,NX/2+1,1,NY,NZ/2+1,NZ/2+1,1,1)
+  PLOT (PRFY,TEM1,100.0,1000.)
+  PLOT (PRFY,C1,-1.,1.)
ENDIF
IF (NZ .GT. 1) THEN
+  PATCH(PRFZ,PROFIL,NX/2+1,NX/2+1,NY/2+1,NY/2+1,1,NZ,1,1)
+  PLOT (PRFZ,TEM1,100.0,1000.)
+  PLOT (PRFZ,C1,-1.,1.)
ENDIF