************************************************************  
 TALK=T;RUN( 1, 1)
 ************************************************************ 
  PQ1 Description
 ************************************************************   
  Title: Heatex
  Authors: EOP/ DBS/ JRS
  Start Date: 10/02/12
  Last editing date: 04/05/13
  Notes by dbs:
  1. The fact that exact values for parallel and counterflow
     are printed in inforout is not mentioned in descr_en.htm
  2. Iy would be useful to activate the metal-temperature feature, 
     at least for transient       
  ***************************************************************
  PQ1 Part 1: Declarations and Settings

  XML-Group General
char(flocon) ! Flow Configuration
flocon=cross! parallel; counter; oblique; 2baffles; 4baffles; 2leakybs
real(timflo) ! timflo=vol*density1/(time_intvl*Mass_flow1)
timflo=0.0

  XML-Group Geometry
real(dxcont) ! Length of Contact Space in x-direction, m
dxcont = 1.0
   Note by dbs 01.05.13 
   The ! command renew lines are all unnecessary, although they
   do no harm.
   The command has been created so as to indicate which settings
   must be ignored by the satellite if it is asked to present the 
   currenly-chosen scenario on the screen. This is something 
   which does not arise in the case of HeatEx which requires no 
   visual display of problem setup.
! command renew
real(dycont) ! Length of Contact Space in y-direction, m
dycont = 1.0
! command renew
real(dzcont) ! Length of Contact Space in z-direction, m
dzcont = 1.0
! command renew
real(ardvol) ! Heat-Transfer Area per unit Volume, 1/m
ardvol=100

  Note by dbs 01.05.13 The following lines are deactivated becsuae,
  although they have relevance to more reaistic heat-exchanger
  simulations, the play no part in HeatEx
  
  real(mevlfr) ! Metal Volume Fraction
  mevlfr = 0.05
  real(tuoudi) ! Tube Outside Diameter, m
  tuoudi = 0.025
  ! command renew
  real(tubthc) ! Tube Wall Thickness, m
  tubthc = 0.002
  ! command renew
  real(lngpit) ! Along-Stream Pitch, m
  lngpit = 0.063
  ! command renew
  real(crspit) ! Cross-Stream Pitch, m
  crspit = 0.058  
  ! command renew
  real(fincoe) ! Fin Coefficient  
  fincoe = 16.0
  real(finpit) ! Fin Pitch
  finpit = 0.0032
  ! command renew
  real(finthc) ! Fin Thickness
  finthc = 0.0008
  ! command renew
  real(finhig) ! Fin Height
  finhig = 0.014  
  ! command renew
  real(conres) ! Contact Resistence in Watts per m**2 deg C
  conres  = 6.e-3 
  boolean(stggrd) ! Staggered Tube Arrangement
  stggrd=t! f  
  ! command renew
  
  XML-Group Variables Solved  
  
  
  XML-Group Material Properties
real(dnsm) ! Metal Density, kg/m**3  
dnsm=7800 
integer(itest) ! Simplify Fluid Properties, 0=no, 1=yes
itest=0 ! 1
real(sphm) ! metal average sp. ht. J/kg. deg C
sphm = 473
real(sph1) ! 1st fluid average sp. ht. J/kg. deg C
sph1 = 4181.8
real(dns1) ! 1st fluid Density, kg/m**3
dns1 = 998.23
real(sph2) ! 2nd fluid average sp. ht. J/kg. deg C
sph2 = 1005.0
real(dns2) ! 2nd fluid Density, kg/m**3
dns2 = 1.189
integer(iprops) ! Constant or Variable 2nd fluid sp. ht., 1=constat, 2=variable
iprops=1! 2
real(tfrB) ! Temperature Fraction of Point B (if variable sp. ht.)
tfrB = 0.45
real(tfrC) ! Temperature Fraction of Point C (if variable sp. ht.)
tfrC = 0.55
real(hfrB) ! Enthalpy Fraction of Point B (if variable sp. ht.)
hfrB = 0.1
real(hfrC) ! Enthalpy Fraction of Point C (if variable sp. ht.)
hfrC = 0.9 

  XML-Group Models
  
  
  XML-Group Initial Conditions
  
  
real(t1ini) ! Temperature of Fluid 1, deg C
t1ini = 20.0
real(t2ini) ! Temperature of Fluid 2, deg C
t2ini = 20.0  
real(tmini) ! Metal Temperature, deg C
tmini = 20.0 

  XML-Group Boundary Conditions
real(min1) ! 1st fluid Mass Flow, kg/s
min1=1.0
real(tin1) ! 1st Fluid Inlet Temperature, deg C
tin1=20.0
real(tin2) ! 2nd Fluid Inlet Temperature, deg C
tin2 = 80.0
real(msp2d1) ! Mass*Specific Heat Ratio (Fluid 2/Fluid 1)
msp2d1=1.0
real(ntumin) ! Number of Transfer units  
ntumin = 1.0  
real(dhtcdt) ! const in U=U0*(1+const*(T-T1)/T2-T1))
dhtcdt = 0.0
real(dhtcdv) ! const in U=U0*(1+const*V2/averageV2)
dhtcdv = 0.0   
  Note by dbs 01.05.13 Whatever roles htrsf1 and 2 may once have
  had, they have none now. Deactivate here and remove fron scene.
  xml
  real(htrsf1) ! fl1 conv. ht.tr resistance / total, dimensionless
  htrsf1 = 0.49
  real(htrsf2) ! fl2 conv. ht.tr resistance / total, dimensionless
  htrsf2 = 0.49   
  
  Note by dbs 01.05.13 No Nusselt or Euler formulae are supplied 
  with or relevant to HeatEx therefore deactivate here and remove 
  from scene.xml
  boolean(nuseul) ! Use Formulae for Nusselt and Euler no.s
  nuseul=f! t  

  XML-Group Output
char(gramon) ! Graphical Monitor Display Curves of
gramon = maxabs! spot; maxmin
boolean(endp) ! Pause at end of Monitor Run
endp = f! t
  integer(nprint) ! Frequency of Print Out
nprint = 1000
boolean(preff) ! Print Effectiveness in inforout?
preff = t! f
boolean(show) ! List Parameter Settings in satlog.dat
show = f! t
boolean(colfil) ! Colour-fill PHOTON Contours?
colfil = t! f 

  XML-Group Numerical
integer(nxcont) ! No. of x-grid intervals in contact space
nxcont = 40
! command renew
integer(nycont) ! No. of y-grid intervals in contact space
nycont = 40
! command renew
  integer(LSTEP) ! No. of Time Steps if not Steady  
LSTEP = 10
  integer(LSWEEP) ! No. of outer iterations (sweeps)
LSWEEP = 200
real(hrlx) ! Linear Under-relaxation Factor
hrlx = 0.10  
  
  XML-End
  
  ************************************************************  
  PQ1 Part 2: Include frommenu.htm
  
incl(frommenu.htm)
  ************************************************************
  PQ1 Part 3: Write Run Script

WRITE(>execPH.bat,@echo on)
WRITE(>>execPH.bat,start /wait %PHOENICS%\d_earth\d_windf\earexe.exe)  
   
  ************************************************************
  PQ1 Part 4: Consequential and Additional Settings
  Note that the declarations below need not be in standard pq1ed 
  form because they are not required to be changed via the PHOENICS
  -Direct menu. Therefore the convenientce of being able to declare
  seveara variable of the same type in a single line is exploited.
  
  Data relating to idealised heat exchangers
real(hin1,hin2,tou1,tou2,min2) 
real(dxhead,area,ohtc)                                     
real(Eff,mindmax)
integer(flotyp,isormeth,igrmon)
integer(nxhead)
boolean(hedh,nusseu)    !
char(parallel,counter,cross,oblique,2baffles,2leakybs)
char(4baffles,formula)       

  Some meanings:
  Flow pattern                     
  flotyp=1 ! parallel-flow  ; flocon=parallel
  flotyp=2 ! counter-flow   ; flocon=counter
  flotyp=3 ! cross-flow     ; flocon=cross
  flotyp=4 ! porous-medium oblique flow  of fluid 2   
                            ; flocon=oblique
  flotyp=5 ! porous-medium 2-baffled flow of fluid 2  
                            ; flocon=2baffles
  flotyp=6 ! porous-medium 2-baffled flow of fluid 2 with leakage
                            ; flocon=2leakybs
  flotyp=7 ! porous-medium 4-baffled flow of fluid 2  
                            ; flocon=4baffles
  
  Fluid properties           
  sph1=  J/kg.degC,specific heat of fluid 1 (water)
  sph2=  J/kg.degC, specific heat of fluid 2 (air)
  dns1=  kg/m**3  as for 20 degC water
  dns2=  kg/m**3  as for 20 degC air at 1 atm pressure
  
  
  tin1= degC, inlet temperature of fluid 1
  tin2= degC, inlet temperature of fluid 1
  min1=mass flow rate of fluid 1 
  min2=mass flow rate of fluid 2
  min2 is later re-set via msp2d1
  ntumin=number of transfer units based smaller mass*specific-heat
  ohtc in  watts/m**2.degC, overall heat-transfer coefficient
  is set via number of transfer units, ntumin
  dhtcdt is the rate of change with dimensionless temperature rise 
         of the overall heat-transfer coefficient.                                              
  dhtcdv is the rate of change with dimensionless velocity rise 
         of the overall heat-transfer coefficient.                                              
  Note: Fluid 1 is thought of as the in-tube fluid which is 
        represented in left-hand grid.
        It is given the properties of water
        
        Fluid 2 is thought of as the in-shell fluid which is  
        represented in right-hand grid.
        For iprops=1, it is given the properties of air; but
        for iprops=2, it is given the properties of a fictitious
        fluid of which the enthalpies at tin1 and tin2 are the 
        same as those of air but between those extremes the 
        enthalpy versus temperature relation is 'three-part-
        linear' as sketched below so as qualitatively to 
        represemt the effect of phase change and 'latent heat'.
        The specific volume of this fluid is that of air in the 
        higher-enthalpy region and of water in the lower-enthalpy
        region. In the intermediate enthalpy region temperature
        varies linearly with enthalpy.                                                                   
  
  min1= kg/s, mass flow rate of fluid 1
  
  dxcont= x-direction length of contact region 
  dycont= y-direction length of contact region 
  dzcont= z-direction length of contact region 
  
  
  data relating to realistic shell&tube heat exchangers
real(dist,mcnd,srfa,st2ti,ti2to,to2fo,fo2ai)

  Note by dbs 01.05.13
  The following appear to have been copied from some more realstic
  heat-echnaner model, for example aircon.
  They are not used here.
  
  tuoudi, the outside diameter of the tubes in the bundle,
  tubthc, the tube-wall thickness
  lngpit, the along-flow pitch of the tubes
  crspit, the cross flow pitch of the tubes            
  finodi, the outside diameter of the fins
  fincoe, the ratio of fin-surface area to tube outside area,
  finpit, the fin pitch itube axis direction
  finhig, the fin height
  stggrd, tube arrangement staggered when true, else in-line 
  conres, tube-to-fin contact resistance tube-to-fin
  nuseul, use Nusselt Euler number formulae
  
  ************************************************************
  Part 2.1 Settings not yet adjustable via menu
  ************************************************************

isormeth=1  ! use the COVAL formulation for enthalpy sources
hedh=F      ! setting to use 2d grid for flotyp=1 and 2 so that
            ! images can be easily compared with other flotyp/
  Note by dbs 01.05.13 If isormeth =1 is the preferred method
  and no provision for changing it via menu id provided, why
  include any other.
  I shall deactivate isormeth-2 below

  ************************************************************

             Note that, although the results will be expressed in 
             dimensionless terms, as being most relevant to the
             instructional purpose, default values are dimensionAL
             so as to impart a sense of realism. Their actual 
             values do not affect the dimensionless effectiveness.

mesgm(  The postulated shape  x horizontal
mesg(      header regions     y vertical  
mesg(     /               \     
mesg(   I--I-------------I--I           ^
mesg(   I  i            >i \I<          .
mesg(   I  i   contact   i  \           .
mesg(   I  i   region    i  I\       dconty 
mesg(   I  i<- dxcont -->i  I dxhead    .
mesg(   I  i             i  I           .
mesg(   I--I-------------I--I           V
if(show) then
show
mesg(before frommenu
flotyp
flocon
dxcont
dycont
dzcont
ardvol
iprops
sph1 
dhtcdt
dhtcdv
dns1 
sph2 
dns2 
tfrB 
tfrC 
hfrB 
hfrC 
min1 
tin1 
tin2 
msp2d1
ntumin
nxcont
nYcont
LSWEEP
hrlx 
igrmon
endp 
preff
endif  
  
mesga(values derived from defaults
area=ardvol*dxcont*dycont*dzcont ! inter-fluid heat-transfer area
  ************************************************************
  Part 3 Interactive re-setting of some variables
  ************************************************************



  FLOTYP
flocon
if(:flocon:.eq.PARALLEL) then
flotyp=1
endif        
if(:flocon:.eq.COUNTER) then  
flotyp=2
endif  
if(:flocon:.eq.CROSS) THen    
flotyp=3
endif
if(:flocon:.eq.OBLIQUE) then  
flotyp=4
endif        
if(:flocon:.eq.2BAFFLES) then 
flotyp=5
endif        
if(:flocon:.eq.2LEAKYBS) then 
flotyp=6
endif        
if(:flocon:.eq.4BAFFLES) then 
flotyp=7
endif        
if(:gramon:.eq.spot) then
igrmon=1
endif
if(:gramon:.eq.maxmin) then
igrmon=2
endif
if(:gramon:.eq.maxabs) then
igrmon=3
endif
gramon
igrmon
if(timflo.gt.0.0) then ! to aid initial testing of transients
itest=1
mesg(itest=:itest:
endif
if(itest.eq.1) then  ! to simplify checking of results
sph1=1
sph2=1
dns1=1
dns2=1
t1ini=(t1ini-tin1)/(tin2-tin1)
t2ini=(t2ini-tin1)/(tin2-tin1)
tin1=0.0
tin2=1.0
tin1
tin2
t1ini
t2ini
endif
min1
sph1
sph2
min2=msp2d1*min1*sph1/sph2     ! effect of msp2d1
min2

hin1=tin1*sph1
hin2=tin2*sph2
min1
hin1
hin2

if(msp2d1.gt.1.0) then
ohtc=ntumin*min1*sph1/area      ! effect of ntumin
mesg(ohtc2=:ohtc:
mindmax=1.0/msp2d1

else
ohtc=ntumin*min2*sph2/area      ! effect of ntumin
mesg(ohtc3=:ohtc:
mindmax=msp2d1
endif

mindmax
mesg(ohtc1=:ohtc:

  
if(show) then
mesg(after frommenu
  flotyp
flocon
dxcont
dycont
dzcont
ardvol
iprops
sph1 
dns1 
dhtcdt
dhtcdv
sph2 
dns2 
tfrB 
tfrC 
hfrB 
hfrC 
min1 
tin1 
tin2 
msp2d1
ntumin
nxcont
nYcont
LSWEEP
hrlx 
igrmon
endp 
show
preff
endif  
  ************************************************************
  Part 4 Consequential and other settings
  ************************************************************
mesga(consequential and other settings

  derived quantities
mesga(derived and other data in left/right terms  
mesg(    header regions      header regions    
mesg(   /               \   /               \  
mesg( I--I-------------I--I--I-------------I--I
mesg( I  i left        i  I  i right       i  I
mesg( I  i contact     i  I  i contact     i  I
mesg( I  i region      i  I  i region      i  I
mesg( I  i             i  I  i             i  I
mesg( I  i             i  I  i             i  I
mesgb( I--I-------------I--I--I-------------I--I
  d=density; c=specific heat; m=mass flow rate
  settings
  ************************************************************
  Part 5 Conventional 25-group Q1
  ************************************************************
  group 1. Run Title and Number
  ************************************************************
TEXT(:irunn: flotyp=:flotyp: ntumin=:ntumin: msp2d1=:msp2d1:

  ************************************************************
  group 2. Time dependence
  !!! erroneous entries
  STEADY = T
  lstep=1
if(timflo.gt.0.0) then
steady=f
timflo
lstep
tlast=lstep*dns1*0.5*xulast*yvlast*zwlast/(timflo*min1)
tlast
grdpwr(t,lstep,tlast,1.0)
idispa=1
itest=1
endif
  ************************************************************
  group 3. X-Direction Grid Spacing
CARTES = T
  mesga(grid settings a
nxhead=1  ! Comment: allowance for /=1 not made fully below 
  if(timflo..eq.0.0) then ! dbs 14.09.16 long-present .. error corrected
if(timflo.eq.0.0) then
if(flotyp.lt.3.and..not.hedh) then
nycont=1 ! economy but currently needed for autoplot
endif
endif
nycont
  steady
  nxhead
  nxcont 
  nycont
  mesg(
  if(flotyp.gt.2) then ! make nxcont divisible by 3 to ease
  nxcont=(nxcont)/3  ! placing of baffles
  nxcont
  nxcont=nxcont*3
  nxcont
  nycont=nxcont
  endif
nregx=6
dxhead=0.1*dxcont ! z-direction length of headers
iregx=1;grdpwr(x,nxhead,dxhead,1.0) ! uniform grid
iregx=2;grdpwr(x,nxcont,dxcont,1.0) ! uniform grid
iregx=3;grdpwr(x,nxhead,dxhead,1.0) ! uniform grid
iregx=4;grdpwr(x,nxhead,dxhead,1.0) ! uniform grid
iregx=5;grdpwr(x,nxcont,dxcont,1.0) ! uniform grid
iregx=6;grdpwr(x,nxhead,dxhead,1.0) ! uniform grid
nregx
nx
xulast

  ************************************************************
  group 4. Y-Direction Grid Spacing
nregy=1
iregy=1;grdpwr(y,nycont,dycont,1.0) ! uniform grid
ny
yvlast
  ************************************************************
  group 5. Z-Direction Grid Spacing
PARAB = F
NZ = 1
ZWLAST =dzcont 
nz
zwlast
  ************************************************************
  group 6. Body-Fitted Coordnsates
  ************************************************************
  group 7. Variables: STOREd,SOLVEd,NAMEd
ONEPHS = T
   * Y in SOLUTN argument list denotes:
   * 1-stored 2-solved 3-whole-field
   * 4-point-by-point 5-explicit 6-harmonic averaging 

solve(p1,u1,v1,h1)
store(hflx,den1,spht, tlr)
if(.not.steady) then
store(tmet)
endif

  dbs 01.05.13 deaqctivated
  if(nuseul) then     ! Define Reynolds numbers
  incl(..\input\reynolds.htm)
  endif
   ! define temperature and density for tube (left-hand) region
(stored var tlr is h1/:sph1: with if(xg.le.xulast/2))
(property den1 is :dns1: with if(xg.le.xulast/2))  
(stored var spht is :sph1: with if(xg.lt.xulast/2))
   ! define teperature and density for shell (right-hand) region
if(iprops.eq.1) then  ! uniform specific heats
(stored var tlr is h1/:sph2: with if(xg.gt.xulast/2))
(property den1 is :dns1: with if(xg.gt.xulast/2))  
(stored var spht is :sph2: with if(xg.gt.xulast/2))
endif

if(iprops.eq.2) then  ! non-uniform  specific heats of fluid 2, 
                      ! suggestive ofphase change
  ! sph2 now represents average specific heat over the range tin1 
  ! to tin2; but the h~t curve is made up of 3 straight lines, not 1
  ! as shown below                   
char(condtn)                      
real(hrA,hrB,hrC,hrD,trA,trB,trC,trD,spAB,spBC,spCD)
  ! r refers to right, shell, and fluid2          
            I
            I               D
            I            *
            I         C
            I        
          hrI        *
            I
            I       *
            I
            I      *
            I
            I     B            
            I  *  
            A____________________________
                        
                        tr
trA=tin1
trD=tin2

trB=trA+tfrb*(trD-trA)  ! 0.45 * possible temp rise
trC=trA+tfrC*(trD-trA)  ! 0.55 * possible temp rise

hrA=tin1*sph2  ! i.e. of shell fluid if reduced to tube inlet temp
hrD=hrA+sph2*(trD-trA)
mesg(hrD=:hrD: hin2=:hin2:
hrB=hrA+hfrB*(hrD-hrA)
hrC=hrA+hfrC*(hrD-hrA)

spAB=(hrB-hrA)/(trB-trA)
spBC=(hrC-hrB)/(trC-trB)
spCD=(hrD-hrC)/(trD-trC)
  density and specific volume
real(dnsB,dnsC,spvlb,spvldh)
dnsB=dns1 ! i.e. water 
dnsC=dns2 ! i.e. air
spvlB=1.0/dnsB
spvldh=(1.0/dnsC-1.0/dnsB)/(hrC-hrB)
if(show)then
hra
hrb
hrc
hrd
tra
trb
trc
trd
spab
spbc
spcd
spvlB
spvldh
endif

patch(shell,volume,nx/2+1,nx,1,ny,1,1,1,lstep)

condtn=h1.lt.:hrb:
(stored var tlr at shell is :trA:+(h1-:hrA:)/:spAB:$ 
 with if(:condtn:))
(property den1 at shell is :dnsB: with if(:condtn:)) 
(stored var spht at shell is :spAB: with if(:condtn:)) 
                                                           
condtn=h1.ge.:hrB:.and.h1.lt.:hrC:                           
(stored var tlr at shell is :trB:+(h1-:hrB:)/:spBC:$
 with if(:condtn:))
(property den1 at shell is 1./(:spvlB:+:spvldh:*(h1-:hrB:))$
 with if(:condtn:)) 
(stored var spht at shell is :spBC: with if(:condtn:)) 
                                                           
condtn=h1.ge.:hrC:                                         
(stored var tlr at shell is :trC:+(h1-:hrC:)/:spCD:$ 
 with if(:condtn:))                  
(property den1 at shell is 1.0/:dns2: with if(:condtn:)) 
(stored var spht at shell is :spCD: with if(:condtn:)) 

endif





store(epor)
mesga(stored and solved variables are:-
stored
  ************************************************************
  group 8. Terms & Devices
   * Y in TERMS argument list denotes:
   * 1-built-in source 2-convection 3-diffusion 4-transient
   * 5-first phase variable 6-interphase transport         
TERMS(h1,n,p,n,p,p,p) ! viscous-dissipation and diffusion inactive
  ************************************************************
  group 9. Properties used if PRPS is not
           stored, and where PRPS = -1.0 if it is!

  **********************************************************
  group 10.Inter-Phase Transfer Processes
  ************************************************************
  group 11.Initial field variables (PHIs)
fiinit(hflx)=0.0
integer(ix1,ix2,iy1,iy2,iz1,iz2,it1,it2)
    Comment; Probably only the epor initial value is needed for
             STEADY=T
    left header of left part
ix1=1;ix2=nxhead;iy1=1;iy2=ny;iz1=1;iz2=nz;it1=1;it2=1
patch(lheadl,inival,ix1,ix2,iy1,iy2,iz1,iz2,it1,it2)
ix1
ix2
(initial of tlr at lheadl is :tin1:)
(initial of h1 at lheadl is :hin1:)
(property den1 at lheadl is :dns1:)  


    contact region of left part
ix1=ix2+1;ix2=ix1+nxcont-1
patch(contl,inival,ix1,ix2,iy1,iy2,iz1,iz2,it1,it2)
  ix1
  ix2
(property den1 at contl is :dns1:)  
if(steady) then
(initial of tlr at contl is :tin1:)
(initial of h1 at contl is :hin1:)
else
(initial of tlr at contl is :t1ini:)
(initial of h1 at contl is :t1ini:*:sph1:)

endif

    right header of left part
ix1=ix2+1;ix2=ix1+nxhead-1
patch(rheadl,inival,ix1,ix2,iy1,iy2,iz1,iz2,it1,it2)
  ix1
  ix2
(property den1 at rheadl is :dns2:)  
(initial of tlr at rheadl is :tin2:)
(initial of h1 at rheadl is :hin2:)
    
    barrier
ix1=nx/2;ix2=ix1    
patch(barrier,inival,ix1,ix2,iy1,iy2,iz1,iz2,it1,it2)
coval(barrier,epor,0.0,0.0)
    Other?
    Note that, although it would be reasonable to set epor
    to eporleft (<1) in the left-hand contact region and to
    1-eporleft in the right-hand contact regions, this has not 
    been done. The u velocit is therefore the 'superficial'
    velocity, as is the v velocity.
    
    left header of right part
ix1=ix2+1;ix2=ix1+nxhead-1
patch(lheadr,inival,ix1,ix2,iy1,iy2,iz1,iz2,it1,it2)
  ix1
  ix2
(initial of tlr at lheadr is :tin2:)

(initial of h1 at lheadr is :hin2:)
(property den1 at lheadr is :dns2:)  

    contact region of right part
ix1=ix2+1;ix2=ix1+nxcont-1
patch(contr,inival,ix1,ix2,iy1,iy2,iz1,iz2,it1,it2)
mesg(contact right
ix1
ix2
(initial of den1 at contr is :dns2:)
if(steady) then
(initial of tlr at contr is :tin2:)
(initial of h1 at contr is :hin2:)
else
(initial of tlr at contr is :t2ini:)
(initial of h1 at contr is :t2ini:*:sph2:)
endif
  
  for flotyp=3 or greater fiinit(u1)=0.0 will suffice

  (initial of h1 at contr is :hin2:)
(property den1 at contr is :dns2:)  

    right header of right part
ix1=ix2+1;ix2=ix1+nxhead-1
patch(rheadr,inival,ix1,ix2,iy1,iy2,iz1,iz2,it1,it2)
  ix1
  ix2
(initial of tlr at rheadr is :tin2:)


(initial of h1 at rheadr is :hin2:)
(property den1 at rheadr is :dns2:)  



  ************************************************************
   group 13. Boundary & Special Sources
XCYCLE = F
   
   13.1 Fix velocities
   (a) all flotyp
   fix v to zero on left
patch(partl,cell,1,nx/2,1,ny,1,1,1,lstep)
coval(partl,v1,fixval,0.0)
   
   (b) flotyp=1 or 2
   fix v to zero on right
if(flotyp.eq.1.or.flotyp.eq.2) then
patch(partr,cell,NX/2+1,nx,1,ny,1,1,1,lstep)
coval(partr,v1,fixval,0.0)
endif   
   (c) flotyp 3
   fix u to zero on right
if(flotyp.eq.3) then
patch(partr,cell,NX/2+1,nx,1,ny,1,1,1,lstep)
coval(partr,u1,fixval,0.0)
endif   
   
   13.2 Fix mass in- and out-flows
   (a) all flotyp
   west in-flow on left; east out on right
patch(inl,west,1,1,1,ny,1,1,1,lstep)
coval(inl,p1,fixflu,min1/dycont)  
coval(inl,h1,onlyms,hin1) 
coval(inl,u1,onlyms,min1/(dns1*dycont)) 
   
patch(oul,west,nx/2,nx/2,1,ny,1,1,1,lstep)
coval(oul,p1,1.e6,0.0)  
   
   (b) flotyp=1 
   west in-flow on right; east out on right
if(flotyp.eq.1) then
patch(inr,west,nx/2+1,nx/2+1,1,ny,1,1,1,lstep)
coval(inr,p1,fixflu,min2/dycont)  
coval(inr,h1,onlyms,hin2) 
coval(inr,u1,onlyms,min2/(dns2*dycont)) 
   
patch(our,west,nx,nx,1,ny,1,1,1,lstep)
coval(our,p1,1.e6,0.0)  
endif   
   (c) flotyp=2
   east inflow on right; west outflow on right
if(flotyp.eq.2) then
patch(inr,west,nx,nx,1,ny,1,1,1,lstep)
coval(inr,p1,fixflu,min2/dycont)  
coval(inr,h1,onlyms,hin2) 
coval(inr,u1,onlyms,-min2/(dns2*dycont)) 
   
patch(our,west,nx/2+1,nx/2+1,1,ny,1,1,1,lstep)
coval(our,p1,1.e6,0.0)  
endif   
   (d) flotyp 3
   south inflow on right; north outflow on right
if(flotyp.eq.3) then
patch(inr,south,nx/2+2,nx-1,1,1,1,1,1,lstep)
coval(inr,p1,fixflu,min2/dxcont)
coval(inr,v1,onlyms,min2/(dns2*dxcont))
coval(inr,h1,onlyms,hin2)  
mesg(hin2 at inr for flotyp 3 is :hin2:                            

patch(our,cell,nx/2+2,nx-1,ny,ny,1,1,1,lstep)
coval(our,p1,1.e6,0.0)  

   extend barrier one cell to right
ix1=nx/2;ix2=ix1+1    
patch(barrier,inival,ix1,ix2,iy1,iy2,iz1,iz2,it1,it2)
   block off right-right header
patch(barrier2,inival,nx-1,nx-1,1,ny,1,1,1,lstep)
coval(barrier2,epor,1.0,0.0)
patch(rrheader,cell,nx,nx,1,ny,1,1,1,lstep)
coval(rrheader,h1,fixval,hin2)
coval(rrheader,p1,fixval,0.0)
patch(lrheader,cell,nx/2+1,nx/2+1,1,ny,1,1,1,lstep)
coval(lrheader,h1,fixval,hin2)
coval(lrheader,p1,fixval,0.0)
endif   
   (e) flotyp 4, 5 and 6
   North-east nozzled inflow on right. north-west nezzled out on right
   zero epor at nx/2+1 and nx-2
if(flotyp.gt.3) then
integer(nxnoz)
nxnoz=nxcont/6
real(dxnoz)
dxnoz=dzcont/6
ix2=nx-1;ix1=ix2-nxnoz+1 
patch(inr,south,ix1,ix2,1,1,1,1,1,lstep)
coval(inr,p1,fixflu,min2/dxnoz)
coval(inr,v1,onlyms,min2/(dns2*dxnoz))
coval(inr,h1,onlyms,hin2)  

ix1=nx/2+2;ix2=ix1+nxnoz-1                             
patch(our,cell,ix1,ix2,ny,ny,1,1,1,lstep)
coval(our,p1,1.e6,0.0)  

   extend barrier one cell to right
ix1=nx/2;ix2=ix1+1    
patch(barrier,inival,ix1,ix2,iy1,iy2,iz1,iz2,it1,it2)
   block off right-right header
patch(barrier2,inival,nx-1,nx-1,1,ny,1,1,1,lstep)
coval(barrier2,epor,1.0,0.0)
patch(rrheader,cell,nx,nx,1,ny,1,1,1,lstep)
coval(rrheader,h1,fixval,hin2)
coval(rrheader,p1,fixval,0.0)
patch(lrheader,cell,nx/2+1,nx/2+1,1,ny,1,1,1,lstep)
coval(lrheader,h1,fixval,hin2)
coval(lrheader,p1,fixval,0.0)
endif   
   
   (f) flotyp 5
   zero epors at two baffle positions, nxcont/3 apart
if(flotyp.gt.4) then
ix1=nx/2+nxcont/3+1;ix2=ix1
patch(baf1,inival,ix1,ix2,ny/6,ny,1,1,1,lstep)
coval(baf1,epor,1.0,0.0)
ix1=nx/2+2*nxcont/3+1;ix2=ix1
patch(baf2,inival,ix1,ix2,1,ny-ny/6+1,1,1,1,lstep)
coval(baf2,epor,1.0,0.0)
endif   
   
   (g) flotyp 6
   Finite epor (>0) at two baffle positions represents leakage
if(flotyp.eq.6) then
coval(baf1,epor,1.0,0.05)
coval(baf2,epor,1.0,0.05)
endif   
   
   (h) flotyp 7
   Set epor=0 at four baffle positions
if(flotyp.eq.7) then
ix1=nx/2+1+nxcont/5;ix2=ix1
nxcont
ix1
patch(baf1,inival,ix1,ix2,ny/6,ny,1,1,1,lstep)
coval(baf1,epor,1.0,0.0)

ix1=nx/2+1+2*nxcont/5;ix2=ix1
ix1
patch(baf2,inival,ix1,ix2,1,ny-ny/6+1,1,1,1,lstep)
coval(baf2,epor,1.0,0.0)

ix1=nx/2+1+3*nxcont/5;ix2=ix1
ix1
patch(baf3,inival,ix1,ix2,ny/6,ny,1,1,1,lstep)
coval(baf3,epor,1.0,0.0)

ix1=nx/2+1+4*nxcont/5;ix2=ix1
ix1
patch(baf4,inival,ix1,ix2,1,ny-ny/6+1,1,1,1,lstep)
coval(baf4,epor,1.0,0.0)
endif   

   13.3 Momentum fluxes
   for all flow types
  if(.not.nusseu) then   
   porous-medium resistance for u and v on right
real(xdirres,ydirres)
xdirres=1.e2;ydirres=xdirres*10   
patch(rightres,volume,nx/2+2,nx-1,1,ny,1,1,1,lstep)
coval(rightres,u1,xdirres,0.0)   
coval(rightres,v1,ydirres,0.0)   

  else  ! call compute resistances in saparate file

  incl(..\input\euler.htm)
  endif   
   13.4 Heat fluxes
   for all flotyp
   location of heat-transfer patches
ix1=nxhead+1;ix2=nxhead+nxcont
patch(left,volume,ix1,ix2,1,ny,1,1,1,lstep)
patch(leftlin,volume,ix1,ix2,1,ny,1,1,1,lstep)
  
ix1=ix1+nx/2;ix2=ix2+nx/2
patch(right,volume,ix1,ix2,1,ny,1,1,1,lstep)
patch(rightlin,volume,ix1,ix2,1,ny,1,1,1,lstep)
  
  temperature difference
(stored var tdif at left is tlr[+:nx/2:,,]-tlr) 
(stored var tdif at right is tlr[-:nx/2:,,]-tlr) 
  
                                           !  heat flux 
if(.not.nusseu) then
if(dhtcdt.eq.0.0.and.dhtcdv.eq.0) then     ! uniform heat-transfer
       
real(vhtc)
vhtc=ohtc*ardvol
       ! vhtc is a constant
(stored var hflx at left is :vhtc:*tdif)   ! coefficient
(stored var hflx at right is :vhtc:*tdif) 

else               
                       ! non-uniform heat-transfer coefficient
 real(dummy)
 if(dhtcdt.ne.0.0) then ! overall heat-transfer coefficient varies 
                       ! with typical dimensionless temperature rise 
                       ! defined as  (.5*(tl+tr)-tin1)/(tin2-tn1)              

 dummy=0.5/(tin2-tin1)                   
 (stored var ttyp at left is :dummy:*(tlr[+:nx/2:,,]+tlr))                   
 (stored var ttyp at right is ttyp[-:nx/2:,,]) ! i.e. that already                  
 (stored var tfac is 1+:dhtcdt:*ttyp)          ! computed for left         
 (stored var vhtc is :vhtc:*tfac)                   

 else
                 ! overall heat-transfer coefficient varies 
                 ! with velocitu rise above cross-flow average 
                 ! defined as  (.5*(tl+tr)-tin1)/(tin2-tn1)              
 store(vabs)
 real(vave)
 vave=min2/(dns2*dycont)
 dummy=dhtcdv/vave
 vave
 dummy
 (stored var vfac at right is 1+dummy*(vabs-:vave:))
 (stored var vfac at left is vfac[+:nx/2:,,])
 (stored var vhtc is :vhtc:*vfac)                   
 endif
       ! vhtc is a calculated 3d variable
(stored var hflx at left is vhtc*tdif) 
(stored var hflx at right is vhtc*tdif) 

endif
else ! if(nusseu) calculate hflux in separate file
incl(..\input\nusselt.htm)
endif
  
  sources
if(isormeth.eq.1.or.iprops.eq.2) then ! This is the preferred method 
  ! distinguished by having one patch with a fixed-flux source and a
  ! second patch with a linearised source, the VAL being the phistar.
  ! Although the patches could be combined with (presumably) a 
  ! computer-time advantage, the VAL would be less easy to interpret.
(source of h1 at left is hflx with fixf) 
(source of h1 at right is hflx with fixf) 
(source of h1 at leftlin is coval(:vhtc:/spht,h1)) 
(source of h1 at rightlin is coval(:vhtc:/spht,h1)) 
endif

  dbs 01.05.13 deactivated next
  if(isormeth.eq.2.and.iprops.eq.1) then  ! use only for iprops=1, 
               ! for which h=sp*t . Otherwise covals are incorrect
  char(coe,valu)          
  coe=:vhtc:/:sph1:
  valu=h1[+:nx/2:,,]*:sph1:/:sph2:  
  (source of h1 at left is coval(:coe:,:valu:))
  coe=:vhtc:/:sph2:
  valu=h1[-:nx/2:,,]*:sph2:/:sph1:  
  (source of h1 at right is coval(:coe:,:valu:))  
  endif

   ************************************************************
   group 14. Downstream Pressure For PARAB
   ************************************************************
   group 15. Terminate Sweeps
 LITHYD = 1 ;LITFLX = 1 ;LITC = 1 ;ITHC1 = 1
 SELREF = f
isg21=lsweep ! to ensure that full sweep number is executed
lsweep=lsweep/lstep
  nprint=lsweep/2  ! to enable convergence to be checked
  ************************************************************
   group 16. Terminate Iterations
  ************************************************************
   group 17. Relaxation
relax(h1,linrlx,hrlx)
conwiz=t  
  ************************************************************
   group 18. Limits
  ************************************************************
   group 19. Data transmitted to GROUND
  ************************************************************
   group 20. Preliminary Printout
  ************************************************************
   group 21. Print-out of Variables
     * Y in OUTPUT argument list denotes:
     * 1-field 2-correction-eq. monitor 3-selective dumping      
     * 4-whole-field residual 5-spot-value table 6-residual table

real(maxhflx)
if(msp2d1.gt.1.0) then
maxhflx=(tin2-tin1)*min1*sph1
maxhflx  
else
maxhflx=(tin2-tin1)*min2*sph2
endif

  
(make effect is 0.0)
(store1 effect at left is nets(h1,left)/:maxhflx:)
if(preff) then
(print of effectiveness is effect)
endif
if(timflo.gt.0.) then
(print of time_flow_/_m1 is :timflo:)
endif
(print of ntumin is :ntumin:)
  (print of flotyp is :flotyp:)
(print of nxcont is :nxcont:)
(print of mspmin_/_mspmax is :mindmax:)
  (print of msp2d1 is :msp2d1:)
  (make fluxr is 0.0)
  (store1 fluxr at contr is nets(h1,right))
  (print of right_heat_flux is fluxr)
  (make fluxl is 0.0)
  (store1 fluxl at contl is nets(h1,left))
  (print of left_heat_flux is fluxl)
  (print of max_poss_heat_flux is :maxhflx:)
  (print of overall_heat_tr_coeff is :ohtc:)
  (print of contact_area is :area:)
  (print of itest is :itest:)
(print of dhtcdt is :dhtcdt:)
(print of iprops is :iprops:)
  (print of isormeth is :isormeth:)
(print of h1_relax_factor is :hrlx:)

real(exactef,expr1)          ! exact expressions for effectiveness
if(flotyp.lt.3) then         ! for flotyp = 1 and 2 only
if(flotyp.eq.2) then         ! for counterflow
 if(mindmax.eq.1.0) then
+ exactef=ntumin/(1.0+ntumin)
 else
+ expr1=exp(-ntumin*(1-mindmax))
+ expr1
+ exactef=(1-expr1)/(1-mindmax*expr1)
 endif 
endif

if(flotyp.eq.1) then     ! for parallel flow
 expr1=exp(-ntumin*(1.0+mindmax))
expr1
 exactef=(1.0-expr1)/(1+mindmax)
endif
(print of exact_effct_ness is :exactef:)
(print of effectvns_ratio is effect/:exactef:)
endif
    create table for plotting changes of nett sum during sweeps
  (table in heatflux is GET(NETS(h1,left)) with HEAD(heatflux)!sweep)

  ************************************************************
   group 22. Monitor Print-Out
tstswp=-1   
SPEDAT(SET,GXMONI,TRANSIENT,L,F)  

  GROUP 22.1. Inforin print-out
write(>inforin,SimScene is HeatEx of April 2013)
if(itest.eq.1) then
write(>>inforin,constant properties are used)
endif
  write(>>inforin,winset is :winset:)
  write(>>inforin,fluid is :flname:)
  write(>>inforin,gravty is :gravty:)  
  write(>>inforin,turbulence model is :turb:)    
write(>>inforin,flow configuration is :flocon:)  
  write(>>inforin,print_general is :prgen:)  
  write(>>inforin,print_boun_conds is :prboco:)  
  write(>>inforin,print_dimensionless is :prdino:)  
  write(>>inforin,priny_grid&geom is :prgrge:)  
  write(>>inforin,print_ref_props is :prrepr:)  
  write(>>inforin,print_numerical is :prnum:)  
  write(>>inforin,print_all is :prntall:)  
  write(>>inforin,   )  
  ************************************************************
   group 23.Field Print-Out & Plot Control
output(p1,n,n,n,n,n,n)   
output(u1,n,n,n,n,n,n)   
output(v1,n,n,n,n,n,n)   
orsiz=0.2
patch(plot,profil,1,nx,1,1,1,1,1,lstep)
coval(plot,tlr,0.0,0.0)
coval(plot,den1,0.0,0.0)
coval(plot,h1,0.0,0.0)
  ************************************************************
   group 24. Dumps For Restarts

if(igrmon.eq.2) then
#maxmin
endif
if(igrmon.eq.3) then
#maxabs
endif
if(endp) then 
 #endpause
endif
  other settings may be in q2
  
  write display macro
integer(nn)
real(vref)
vref=10.0*min2/(dxcont/dns2)
nn=nxcont ! to reduce line length below 
if(steady) then            ! macro for steady
write(>u,p)                                                                               
write(>>u,phi ;;;;)            
if(flotyp.lt.3.and..not.hedh) then
write(>>u, )
write(>>u,AUTOPLOT )
write(>>u, )
write(>>u,file; phi 5 )
write(>>u, )
write(>>u,cl; da 1; tlr; colf 1 )
write(>>u,  msg temperatures: fluid 1 on left fluid 2 on right)
write(>>u,  msg )
write(>>u,  msg                    ignore join in middle)
write(>>u,  pause)
write(>>u,  dump; tempcrve)
write(>>u,cl; da 1;  hflx; colf 1)
write(>>u,  msg heat fluxes: positive on left, negative on right)
write(>>u,  msg )
write(>>u,  msg ignore join in middle and implied zeroes at ends)
write(>>u,  pause)
write(>>u,  dump; hflxcrve)
write(>>u,enduse)
goto next
endif                                                                   
write(>>u, )
write(>>u,window 0 1 .3 .8 )
write(>>u,gr z 1 col 1)
write(>>u,dump grid)

if(colfil) then   ! colour fill on
write(>>u,pause;gr off;cl;gr ou z 1)
write(>>u,con tlr iz 1 ix 2 :nn+1: iy 1 m fi;.001)                                                              
write(>>u,pause)

write(>>u,con tlr iz 1 ix :nn+4: :nx-1: iy 1 m fi;.001)
write(>>u,dump temprtur)
write(>>u,pause)   
           
write(>>u,con off; cl)
write(>>u,con hflx z 1 ix 2 :nn+1: iy 1 m fi;.001)
write(>>u,pause)

write(>>u,con hflx z 1 ix :nn+4: :nx-1: iy 1 m fi;.001)
write(>>u,pause)
write(>>u,dump heatflux)

else  ! colour fill off

write(>>u,pause;gr off;cl;gr ou z 1)
write(>>u,con tlr iz 1 ix 2 :nn+1: iy 1 m  col 1;int 25)                                                              
write(>>u,con tlr iz 1 ix :nn+4: :nx-1: iy 1 m col 1 ;int 25)
write(>>u,pause)   
write(>>u,con off; cl)
write(>>u,con tlr iz 1 ix 2 :nn+1: iy 1 m  col 0;int 25)                                                              
write(>>u,con tlr iz 1 ix :nn+4: :nx-1: iy 1 m col 0 ;int 25)
write(>>u,dump temprtur)
           
write(>>u,con off; cl)
write(>>u,con hflx z 1 ix 2 :nn+1: iy 1 m col 1 ;int 25)
write(>>u,con hflx z 1 ix :nn+4: :nx-1: iy 1 m col 1;int 25)
write(>>u,pause)
write(>>u,con off; cl)
write(>>u,con hflx z 1 ix 2 :nn+1: iy 1 m col 0 ;int 25)
write(>>u,con hflx z 1 ix :nn+4: :nx-1: iy 1 m col 0;int 25)
write(>>u,dump heatflux)
endif                    ! end of colour fill option

write(>>u,pause)
write(>>u,con tdif z 1 ix 3 :nn: iy 1 m fi;.001)
write(>>u,dump tdif)
write(>>u,pause)
write(>>u,con off; cl)
if(dhtcdt.gt.0.0) then
write(>>u,con tfac z 1 ix 3 :nn: iy 1 m fi;.001)
write(>>u,dump tfac)
write(>>u,pause)
endif

if(flotyp.gt.3) then
write(>>u,con off;cl)
write(>>u,window 0.5 0.8 .3 .8 )
write(>>u,set vec ref :vref:)
write(>>u,vec z 1 col 1)
write(>>u,dump vecs)
write(>>u,pause)
write(>>u,vec off;cl)
write(>>u,window 0 1 .3 .8 )
write(>>u,con u1 z 1 fi;.001)
write(>>u,dump u1)
write(>>u,pause)
write(>>u,con off;cl)
write(>>u,con v1 z 1 fi;.001)
write(>>u,dump v1)
endif
label next 
else   ! macro for steady=f uses parada
write(>u,p)                                                                               
write(>>u,parada)     

  if(flotyp.gt.2) then                 ! flotyp>2
integer(zscale)
zscale=2*timflo       
write(>>u,1 1 :zscale:)            
write(>>u, )            
write(>>u,vi -1 1 1)                                                                   
write(>>u,window 0 .85 0.15 0.9)                                                             
write(>>u,do iz=1,:lstep:,1)                                                                    
write(>>u,con tlr z iz ix 2 :nn+1: fi;.001)                                                   
write(>>u,con tlr z iz ix :nn+4: :nx-1: fi;.001)                                                   
write(>>u,enddo) 
  else                                 ! flotyp=1 0r 2
  write(>>u, )            
  write(>>u,window 0 1 .3 .8 )
  write(>>u,con tlr y 1 ix 2 :nn+1: fi;.001)                                                   
  write(>>u,con tlr y 1 ix :nn+4: :nx-1: fi;.001)                                                   
  endif        
                                                                  
write(>>u,dump unstead:flotyp:)     
write(>>u,msg( press return to view step-by-step)                                                              
write(>>u,pause)
write(>>u,clear; red; vi z)
write(>>u,msg(temperature contours of fluids 1 & 2 at successive times)
write(>>u,do iz=1,m)
write(>>u,con tlr z iz ix 2 :nn+1: fi;.001)
write(>>u,con tlr z iz ix :nn+4: :nx-1: fi;.001)
write(>>u,pause)
write(>>u,msg(press return)
write(>>u,enddo)
endif
inifld=t
STOP