TALK=T;RUN(1,1)
  ************************************
  PQ1 description
  ************************************
  Title: Tubeflow test
  Author:  Elena Pankova
  Start Date:  22.02.2013
  Editor: Brian Spalding
    for the purposes of:
    (1) removing split of tube into five sections
    (2) increasing the number of fluids
    (3) reducing the number of turbulence models
    (4) removing out-of-date comments and deactivations
    (5) introducing explanatory comments
    (6) removal of ydimrat
    (7) providing print options for inforout only
    (8) utilising new inforin capability
    (9) removing the parabolic optio entirel, in the grounds that
        * it is not usable only for incompressile fluids, 
        * even for then it provides only an economy of computer 
          memory  which users will not notice; and
        * its nature would have to be explained if it were retained. 
    (10) place flow-configuration option in 'general' rather them
        'other numerical'  
    (11) toend and howmany introduced    
  Last editing date: 23.03.13
  Purpose: To test the 10-item structure
  Brief description:
  Flow is through a straight pipe of circular cross-section with thick 
  wall allowing heat exchange with an external environment.
  
  Calculation modes are:
  1. fully-developed, with uniform longitudinal temperature gradient
  2. fully-developed, with fixed wall temperature
  3. elliptic, i.e whole tube, but not parabolic
  
  ************************************
  PQ1 Part 1 Declarations and settings
  ************************************
vredit=f ! command preventing satellite from double reading
         ! and re-writing of Q1. It renders unnecessary the
         ! introduction of saveXbegin and saveXend around In-Form
         ! statements when VR objects are present.   
  Utilities           
char(formula)
real(dummy)  
  
  xml-group general
  insert declaration and settings below
CHAR(floform) ! flow formulation
floform=developing_flow! full_dev_uniform-Tgrad; full_dev_uniform-T$
wall


  
  xml-group geometry
  insert declaration and settings below

real(xang) ! circumferential extent, radians
xang=3.14159 

real(diam) ! tube inner diameter, m 
diam=0.05
! command renew
  ! dbs 13.03.13 An explanatory note on this is needed

real(thck) ! wall thickness, m 
thck=0.002
! command renew

real(lng1) ! tube length, m
lng1=5.0
  ! dbs 13.03.12 one length segment only
  
  xml-group variables solved
  insert declaration and settings below

  solution
boolean(solp1) ! solve p1
solp1=T ! F

boolean(solu1) ! solve u1
solu1=T ! F

boolean(solv1) ! solve v1
solv1=T ! F

boolean(solw1) ! solve w1
solw1=T ! F

boolean(solt1) ! solve tem1
solt1=T  ! F

                                 
  xml-group material properties
  insert declaration and settings below

char(flname) ! name of fluid in case 089
flname=use_specified_properties ! saturated_water;$
SAE_5W-30_engine_oil;ethylene_glycol;glycerin;air;$
ammonia;carbon_dioxide;carbon_monoxide;hydrogen;$
nitrogen;oxygen;superheated_water_vapour
  
  
real(prno) ! Prandtl number
prno=7.0
  fluid values of water at 20 degC as in PROPS file
real(dnsf) ! fluid density, kg/m**3
dnsf=998.23
  dnsf=1.0 !****
           !**** indicates values sometimes introduced via manual
           !     editing during testing

real(sphf) ! fluid specific heat, J/(kg*degC) 
sphf=4181.8 
  sphf=1.0 !****

real(nulf) ! fluid kinematic viscosity, m**2/s
nulf=1.006E-6
  nulf=1.e-5 !****

real(cndf) ! fluid thermal conductivity, W/(m*degC)
cndf=sphf*nulf*dnsf/prno

   real(prno) ! Prandtl number sphf*dnsf*nulf/cndf
   prno=7.0


   cndf=0.597
   cndf=0.025 !****

real(texf) ! fluid thermal expansient coeff. degC**-1
texf=1.18e-4

  dnsf=1
  sphf=1
  cn
           ! tube-wall values for copper at 27 degC
real(dnst) ! tube wall density, kg/m**3
dnst=8954.0 

real(spht) ! tube wall specific heat, J/(kg*degC) 
spht=383.0

real(cndt) ! tube wall thermal conductivity, W/(m*degC)
cndt=381.0


  xml-group models
  insert declaration and settings below

char(turb) ! turbulence model
  ! dbs 13,02.13 de-activate kl and ko models
turb=lvel!none;kemodl
  turb=none
  turb=kemodl 
boolean(gravty) ! gravity
gravty=F ! T
  gravty=t !****

real(gravax) ! xcart-dir gravitational acceleration, m/s**2
gravax=9.81

real(gravay) ! ycart-dir gravitational acceleration, m/s**2
gravay=0.0

real(gravaz) ! zcart-dir gravitational acceleration, m/s**2
gravaz=0.0


  
  xml-group initial conditions
  insert declaration and settings below
  
  xml-group boundary conditions
  insert declaration and settings below

char(winset) ! inlet velocity setting
winset=from_Re!from_set_value

real(reno)  ! Reynolds number
reno=5000

char(wfin) ! inlet velocity, m/s
wfin=1.0 

real(tfin) ! inlet temperature, degC
tfin=80.0
  tfin=100.0 !***
  tfin=21.0 !***

real(pout) ! outlet gauge pressure, bar
  pout=1.e5
pout=0.0

real(tenv1) ! external temperature section 1, degC
tenv1=20.0
  ! dbs 13.03.13 de-activate 5-length-segment provision

real(envco1) ! external heat-tr, coeff. sect.1, W/(m**2*degC)
envco1=1.e+4
  ! dbs 13.03.13 de-activate 5-length-segment provision

real(foures) ! fouling resistance m**2 degC/W
foures=0.0

  envco1=1.e+7  !*****
  
  ! note that the next sequence has the format which enables pq1ed 
  ! to create the appropriate list lines in tubeflow.xml
  ! The ! following tgrad therefore does NOT signify that a comment 
  ! follows. Instead alternative settings follow, separated by 
  ! semi-colons
  
  xml-group output
  insert declaration and settings below
  
  Print to inforout
integer(howmany) ! how many output sets?
howmany=1
  ! dbs 22.03.13 print-out features relating to RESULT file 
  ! deactivated because these are unlikely to interest SimScene
  ! users  
  
boolean(prgen) ! print general 
prgen=T ! F
  
boolean(prgrge) ! grid and geometry 
prgrge=T ! F

boolean(prrepr) ! ref props 
prrepr=T ! F

boolean(prdino) ! dimensionless nos. 
prdino=T ! F

boolean(prboco) ! boundary conds. 
prboco=T ! F

boolean(prtumo) ! turbulence models 
prtumo=T  ! F

boolean(prnum) ! numerical 
prnum=T ! F

boolean(prntall) ! print all? 
prntall=F ! T     
  
  xml-group computational grid
  insert declaration and settings below

boolean(polar) ! indicator of polar grid
polar=T ! F
     polar=f   !***
integer(nxx) ! circumferential interval no.
nxx=36   
!!! reduced to 1 when absence of gravity render flow 2D
   nxx=1   !******

integer(nrgy) !  radial regions
nrgy=2

integer(nrgy1) ! interval no. in radial region 1
nrgy1=20

integer(nrgy2) ! interval no. in radial region 2
nrgy2=2

  ! dbs 13.03.13 This option now withdrawn for reason below
real(ydimrat) ! y-cell diminution ratio 
  ydimrat=1.0  
            ! Note that 0.9 was at first made the default; but when 
            ! nrgy1 was made large, this led to errors resulting 
            ! from excessively small grid intervals. It seems to be 
            ! a dangerous grid option for that reason.
            
integer(nrgz)  ! number of longitudinal regions 
  ! dbs 13.03.13 now one region only
nrgz=1

integer(nrgz1) ! longitudinal interval no 1
nrgz1=50

  
  xml-group  numerical
  insert declaration and settings below

integer(numits) ! number of iterations
numits=100
 
real(relxw1) ! relaxation factor for longitudinaal velocity
relxw1=0.5

real(relxt1) ! relaxation factor for temperature
relxt1=0.5

boolean(toend) ! continue run to end
toend=T!F

  xml-end
 
  
  ************************************
  PQ1 Part 2 Include Frommenu
  ************************************
INCL(frommenu.htm)
mesg(after frommenu
floform
  ************************************
  PQ1 Part 3 Consequential settings
  ************************************
    GROUP 1. Run title and other preliminaries
    
Text(tubeflow)   
title
floform
BOOLEAN(fdf,tgrad,twall)
if(:floform:.eq.full_dev_uniform_tgrad) then       !-----------
fdf=t                                                         !
tgrad=t                                                       !
  numits=50*numits                                              !
nrgy1=nrgy1*4  ! this increased grid fineness seems to be     !
               ! necessary for accurate Nusselt no,           !
mesg(numits and nrgy1 increased by satellite                  !
endif                                              !-----------

if(:floform:.eq.full_dev_uniform_twall) then       !-----------
fdf=t                                                         !
tgrad=t                                                       !
  numits=50*numits                                              !
nrgy1=nrgy1*4  ! this increased grid fineness seems to be     !
               ! necessary for accurate Nusselt no,           !
mesg(numits and nrgy1 increased by satellite                  !
endif                                              !-----------

if(:floform:.eq.developing_flow) then
fdf=f
endif

fdfsol=fdf


TEXT(fluid is :flname: flow is :floform: grav=:gravty:    
gravty    
if(.not.gravty) then
nxx=1
xang=0.1
  numits=100 !*****
solu1=f
  if(fdf) then
  solp1=f  ! no need to solve for pressure when nx=1
  endif

endif 
mesg(fdf=:fdf:  
if(foures.gt.0.0) then ! add fouling resistance
do ii=1,5
envco:ii:=1/(foures+1.0/envco:ii:)
enddo
endif

 ************************************************************
 
 IRUNN = 1
     
    GROUP 2. Transience; time-step specification
 STEADY = T    
 ************************************************************         
    GROUP 3. X-direction grid specification
 CARTES = .not.polar
   cartes=t !****
   NX = 1
   XULAST =1.
   XFRAC(1)=1.
grdpwr(x,nxx,xang,1.0)  
  nx
xulast 
 ************************************************************    
    GROUP 4. Y-direction grid specification
   NY = 1
   YVLAST =1.
   YFRAC(1)=1.
nregy=nrgy 
  iregy=1;grdpwr(y,nrgy1,-0.5*diam,ydimrat)
  ! dbs 13.03.13 ydimrat replaced by 1.0
iregy=1;grdpwr(y,nrgy1,-0.5*diam,ydimrat)
  iregy=1;grdpwr(y,nrgy1,-0.5*diam,0.95) !*****
  iregy=1;grdpwr(y,nrgy1,-0.5*diam,1.0) !*****
  yvlast  
iregy=2;grdpwr(y,nrgy2,thck,1.0)  
  ny
  yvlast
 ************************************************************     
    GROUP 5. Z-direction grid specification
   NZ = 1
   ZWLAST =1.
   ZFRAC(1)=1.
if(fdf) then         ! fully-developed flow
nrgz1=1;lng1=0.1*diam
grdpwr(z,nrgz1,lng1,1.0)
  ydimrat=1.0
nz

else                      ! elliptic flow 
nregz=nrgz 
iregz=1;grdpwr(z,nrgz1,lng1,1.0)  
  ! dbs 13.03.13     5 z regions reduced to 1
endif                      ! end of first if(fdf)
  nz
zwlast
 ************************************************************           
    GROUP 6. Body-fitted coordinates or grid distortion
 ************************************************************     
    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 
  STORE(A) ! Variable added to show that cn149 changes to cn148   
if(solp1) then
solve(p1)
 if(fdf) then
 solutn(p1,y,n,n,n,n,n) ! store but do not solve
 endif        
endif   

if(solu1) then
solve(u1)
endif   

if(solv1) then
solve(v1)
 if(fdf) then
 solutn(v1,y,n,n,n,n,n) ! store but do not solve
 endif        
endif        

if(solw1) then
solve(w1)
 if(fdf) then
 terms(w1,p,n,p,p,p,p)
 endif        
endif 
   ! Important note: the order in shich the solve and ator staements
   ! appear is very important because thos other than pressures, 
   ! velocities, enhalpies , ke and eps are given variable indices 
   ! of which start large (usually 150) and then diminish; but the 
   ! lower-index variables are calsolated before higher-index ones 
   ! in earth.
   ! it is therefore important that all property variables which 
   ! depend depend upon temeperature should be stored before tem1 
   ! in order that temerature is calculated before they are; since, 
   ! for gaces the enu_expression contains rho1, enul should be 
   ! stored before rho1, so that rho1 is calculated before it. 
   ! Otherwise initial values will be incorrect.

store(enul,cp1,rho1) 
store(cond)
if(solt1) then
solve(tem1)
terms(tem1,n,p,p,p,p,p) ! de-activate kinetic heating
endif
store(shrz)       ! needed for prrad etc, and for nearwall profil

  if(prskin) then    ! various storage settings for print-out purposes
store(skin)
  endif 
store(skin)       ! needed until e1start.for of 12.06.12 reached   !
                  ! earth delivery                                 !

  if(prstan) then
store(stan)
  endif

  if(prhtco) then
store(htco)
  endif

if(:turb:.ne.none) then
 turmod(:turb:)
 store(enut,len1)
 
 if(:turb:.eq.lvel) then
 solve(ltls)
 store(wdis,wgap)
 endif
 
 if(:turb:.eq.kemodl) then
 if(fdf) then
 terms(ke,p,n,p,p,p,p)
 terms(ep,p,n,p,p,p,p)
 endif
 endif
endif
  ! dbs 13.03.13 kl and ko models removed
store(prps)
 ************************************************************        
  mesg(group 8
    GROUP 8. Terms (in differential equations) & devices 
    * Y in TERMS argument list denotes:
    * 1-built-in source 2-convection 3-diffusion 4-transient
    * 5-first phase variable 6-interphase transport         
 DIFCUT =0.5 ;ZDIFAC =1.
 GALA = F ;ADDDIF = F
 ISOLX = -1 ;ISOLY = -1 ;ISOLZ = -1
 ************************************************************     
    GROUP 9. Properties of the medium (or media) 
press0=1.e5 ! additive to p1 to create absolute pressure
temp0=273   ! additive to p1 to create absolute pressure
            ! note that the following declarations appear in
            ! always called-at-the-start library case 014.htm
  CHAR(ABT,ABP)   ! absolute temperature and pressure for 089
  ABT=(TEM1+TEMP0); ABP=(P1+PRESS0)          ! default setting
 
mesg(flname is :flname:

condtn=if(iy.le.:nrgy1:)
condtn
(property enul is 0.0!if(prps.gt.0)) 
if(:flname:.eq.use_specified_properties) then

mesg(use specified properties
  ! the following are provided to enable the boundary-condition 
  ! settings to be the same whether or not constant properties 
  ! are being used
  ! declarations
char(rho_expression)
char(cp_expression)
char(enu_expression)
char(cond_expression)
  ! settings
rho_expression=dnsf 
cp_expression=sphf  
enu_expression=nulf 
cond_expression=cndf


(property rho1 is :dnsf:!:condtn:)
(property enul is :nulf:!:condtn:) 
(property cp1 is :sphf:!:condtn:) 
(property prndtl(tem1) is -:cndf:!:condtn:) 
(property prndtl(tem1) is -:cndt:!if(iy.gt.:nrgy1:)) 
(stored var cond is :cndf:!:condtn:)          ! These 2 lines affect
(stored var cond is :cndt:!if(iy.gt.:nrgy1:)) ! print-out only
dvo1dt=texf       

else ! use properties according to fluid name

fluid_name=:flname: ! set name for use in next-loaded library case       
#$089
                    ! where the following character strings are 
                    !  declared and set
                    !  rho_expression                                 
                    !  cp_expression                                  
                    !  enu_expression                                 
                    !  emu_expression                                 
                    !  cond_expression                                
if(:flname:.eq.saturated_water) then  
dvo1dt=1.18e-4
endif  

if(:flname:.eq.SAE_5W-30_engine_oil) then
dvo1dt=7.0e-4
endif
   
if(:flname:.eq.Ethylene_Glycol) then 
dvo1dt=6.5e-4
endif

if(:flname:.eq.Glycerin) then
mesg(Glycerin 
endif

endif ! end of if(.not.use_specified properties)

 ************************************************************      
    Group 9.1 make reference properties
(make refdns is 0.0)
(store1 refdns is rho1[,1,]!IF(iz.eq.1))

(make refkvis is 0.0)
(store1 refkvis is enul[,1,]!IF(iz.eq.1))
  
(make refsph is 0.0)
(store1 refsph is cp1[,1,]!IF(iz.eq.1))
                                                !
(make refcnd is 0.0)
(store1 refcnd is cond[,1,]!IF(iz.eq.1))
 ************************************************************      
    GROUP 11. Initialization of variable or porosity fields 
fiinit(tem1)=tfin
fiinit(p1)=pout

patch(inisoli,inival,1,nx,nrgy1+1,ny,1,nz,1,1)
coval(inisoli,prps,0.0,103)   ! copper
coval(inisoli,rho1,0.0,dnst)  ! wall density
coval(inisoli,tem1,0.0,tfin)  ! initial wall temperature
coval(inisoli,cond,0.0,cndt)  ! wall conductivity

  mesg(rho_expression =:rho_expression: 
  mesg(enu_expression =:enu_expression: 
  mesg(cond_expression =:cond_expression: 

if(:winset:.eq.from_Re) then
wfin=:reno*nulf/diam:
endif
mesg(wfin =:wfin: 

(initial of w1 is :wfin:)
(initial of rho1 is :rho_expression:!:condtn:)
(initial of enul is :enu_expression:!:condtn:)
(initial of cond is :cond_expression:!:condtn:)
 
 ************************************************************      
    GROUP 13. Boundary conditions and special sources 
 ************************************************************      
    Group 13.1 all flow configurations wall-related settings. Needed?
 EGWF = T
   WALLCO = GRND2
   WALL(innNer,north,1,nx,nrgy1,nrgy1,1,nz,1,lstep) !****
   coval(inNner,w1,grnd2,0.0)          ! which produce inival
   coval(inNner,tem1,grnd2,0.0)        ! patch
   Note: The above lines were introduce for early, not-longer-used
         formulation. They have been left as a reminder that their
         unexpected production of an inival patch should be 
         investigated when time permits.
 ************************************************************
    Group 13.2 Gravity momentum source
  mesg(group 13.2
  gravty
if(gravty) then            ! Gravity source
patch(intube,phasem,1,nx,1,nrgy1,1,nz,1,1)
dummy=nx
dummy=xang/dummy  
if(gravax.gt.0.0) then
formula=:gravax*dvo1dt:*(tem1-:tenv1:)*sin((2*ix+1)*:dummy/2:)  
(source of u1 at intube is :formula:)  
formula=:gravax*dvo1dt:*(tem1-:tenv1:)*cos(ix*:dummy:)  
(source of v1 at intube is :formula:)  
  endif
  if(gravay.gt.0.0) then
  endif
  if(gravaz.gt.0.0) then
  endif
endif                      ! End of gravity source
endif                      ! end of if(gravty)
mesg(after gravity
fdf
  
  (print wfin is wfin)
wfin
 ************************************************************
if(.not.fdf) then          
    
    Group 13.3 Elliptic ! Inlet boundary  
  mesg(at inlet density=:dnsf: velocity=:wfin:
patch(inlet,low,1,nx,1,nrgy1,1,1,1,lstep)
  
  ! the next two give almost identical results  
  (source of w1 at inlet is coval(0,:wfin:)!onlyms)  ! res 2
(source of w1 at inlet is coval(onlyms,:wfin:))  ! res 3
  
  ! either of the next two will serve, results being almost identical
  (source of p1 at inlet is :wfin:*rho1)   ! res1
(source of p1 at inlet is coval(fixflu,:wfin:*rho1)) ! res2

coval(inlet,tem1,onlyms,tfin)                                      !
 
 ************************************************************
    Group 13.4 Elliptic flow ! Heat loss boundary  
integer(izf,izl)                                                   !
izf=1;izl=nrgz1                                                    !
patch(tubewal1,north,1,nx,ny,ny,izf,izl,1,lstep)                   !
coval(tubewal1,tem1,envco1,tenv1)                                  !
 ************************************************************
    Group 13.5 Elliptic flow ! Outlet boundary  
patch(outlet,high,1,nx,1,ny,nz,nz,1,lstep)                        !!
coval(outlet,p1,1.e10,pout)                                        !!
coval(outlet,tem1,onlyms,same)                                    !!
  coval(outlet,w1,onlyms,0.0)                                       !!
mesg(end of 13.5 fdf=:fdf                                                                  !
else             ! end of if(.not.fdf) ----------------------------
 ************************************************************
    Group 13.6 fully-developed-flow momentum-source settings
mesg(13.6

  spmsfl=:wfin:*rho1*:0.5*(diam*0.5)**2*xulast:
  mesg(spmsfl is :spmsfl:                          !
  spmsfl
  dnsf
  (print spmsfl is spmsfl)


(make camsfl)     ! calculated tube-fluid mass flow rate           !
(store1 camsfl is sum(w1*ahigh*rho1))                            !
  (print camsfl is camsfl)

(make spmsfl)     ! specified fluid mass flow rate                 !
(store1 spmsfl is camsfl!if(isweep.eq.1))                            !
  (print spmsfl is spmsfl)

                                                                   !
(make florat is 1.0) ! calculated flow/specified flow              !
(store1 florat is camsfl/spmsfl)                                 !
(print florat is florat)

  eop    

(make pgrmul is 1.0)  ! prgrad multiplier to adjust flow rate      !
  (store1 pgrmul is pgrmul*(0.01+1/florat)/1.01)                   !
(store1 pgrmul is pgrmul*(1+1/florat)/2)                           !
  (store1 pgrmul is pgrmul*(90+1/florat)/91)                       !
  (print presgrad_mult is pgrmul)                                  !
                                                                   !
(make pgrad is 0.0)                                                !
real(volum)                                                        !
volum=0.5*0.25*diam**2*xulast*zwlast                               !
(store1 pgrad is shrz[,:nrgy1:,]*4.0/:diam:) ! first estimate      !
  (print press_grad is pgrad)                                      !
                                                                   !

  dbs
  (make pgrad is 2.0)   
    (make pgrad is 20.0)   
  real(volum)                                                        !
  volum=0.5*0.25*diam**2*xulast*zwlast                               !
    (store1 pgrad is 0.9*shrz[,:nrgy1:,]*4.0/:diam:!if(isweep.eq.2)) 

  (store1 pgrad is (pgrad+1.0-florat)!if(isweep.gt.1))
  ! This is the most sucessful in procuring convergence so far.
  ! I am suspicious of it because it is not dimensionless

(print pgrad is pgrad)

(make frccoe is 0.0)                                               !
dummy=2.0/(:wfin:**2)                                                !
(store1 frccoe is shrz[,:nrgy1:,]*:dummy:/refdns)                  !
  (print frdcoe is frccoe)  

  (store1 fcre is frccoe*reyno)                                      !
                                                                   !
mesg(whole                                                                   !
patch(whole,volume,1,nx,1,nrgy1,1,nz,1,1)                          !
(source of w1 at whole is pgrad*pgrmul)                            !
  (source of w1 at whole is pgrad)                            !
                                                                   !
                                                                   !
 ************************************************************      !
    Group 13.7 fully-developed-flow heat-source settings           !
  There are several ways of computing the heat flow to the external!
  fluid, namely                                                    !
  1. as the net source of the tubewall source                      !
patch(tubewal1,north,1,nx,ny,ny,1,1,1,lstep)                       !
coval(tubewal1,tem1,envco1,tenv1)                                  !
  deactivated                                                      !
  (make heafl1 is 0.0)  ! heat flux to environment                 !
  (store1 heafl1 is (-nets(tem1,tubewal1)))                        !
  (print heafl1 is heafl1)                                         !
                                                                   !
  2. as the temperature difference times the ap at iy=ny           !
  deactivatd                                                       !
  (stored var aptm is apco(tem1))       ! essential                !
  (make heafl2 is 0.0)  ! heat flux to environment                 !
  (store1 heafl2 is (tem1-:tenv1:)*aptm with if(iy.eq.:ny:))       !
  (print heafl2 is heafl2)                                         !
                                                                   !
  3. as the heatflux within the metal                              !
(stored var antm is anco(tem1))                                    !
  (make heafl3 is 0.0)  ! heat flux through wall                   !
  (store1 heafl3 is (tem1-tem1[,+1,])*antm with if(iy.eq.:ny-1:))  !
  (print heafl3 is heafl3)
  
  4. at tube inner surface                                         !
  this one is used                                                 !
(make heafl4 is 0.0)  ! heat flux to environment                   !
(store1 heafl4 is (tem1-tem1[,+1,])*antm with if(iy.eq.:ny-2:))    !
                                                                   !
patch(fluid,high,1,nx,1,nrgy1,1,nz,1,1) ! tube cross-section patch !
                                                                   !
(make entflo is 0.0)                                               !
(store1 entflo at fluid is sum((tem1-:tenv1:)*cp1*rho1*w1*ahigh))  !                                              
  real(entflin) ! enthalpy flowing in                              !
  entflin=spmsfl*sphf*(tfin-tenv1)                                 !
(make entfli is 0.0)                                               !
(store1 entfli at fluid is sum((:tfin:-:tenv1:)*cp1*rho1*w1*ahigh))!                                                
(make heafac is 1.0)  ! should become unity after enough sweeps    !
(store1 heafac is entfli/entflo)                                   !
(make mult is 1.0)   ! multiplies the dtdz                         !
(store1 mult is mult*(10.0+heafac)/11.0)                           !
  (store1 mult is mult*(1.0+heafac)/2.0)                           !
                                                                   !
real(dtdz) ! This if, strictly speaking, MINUS the  temperature    !
           ! gradient                                              !
  ! first estimate the Nusselt number                              !
real(reyn,pran,nulam,nutur,nusse)                                  !
reyn=diam*:wfin:*dnsf/nulf ! Reynolds using specified properties     !
pran=sphf*dnsf*nulf/cndf                                           !
reyn                                                               !
pran                                                               !
nulam=4.5                                                          !
nutur=0.023*reyn**0.8*pran**0.3333                                 !
if(nutur.gt.nulam) then                                            !
nusse=nutur                                                        !
else                                                               !
nusse=nulam                                                        !
endif                                                              !
                                                                   !
  ! Now follows the source calculation for fully developed flow    !
dtdz=(tfin-tenv1)*(nusse*cndf)*diam/(:wfin:*dnsf*sphf)               !
mesg(const tgrad=:tgrad: twall=:twall:                             !
  ! inital guess not affecting result                              !
if(tgrad) then                                                     !
tgrad                                                              !
(source of tem1 at fluid is cp1*rho1*w1*(:dtdz:)*mult)               !
endif                                                              !
   The above is the difference between the enthalpy flow in from   !
   upstream and the ethalpy flow out to downstream.

 ************************************************************      !
    Group 13.4 fully-developed-flow momentum-source settings         !
                                                                   !
  (make effect is 0.0)                                               !
  (print heat_flux is heafl4)
  (print entfli is entfli)
  (store1 effect is heafl4/entfli)                                   !
  (print effectiveness is effect)                                                                   !
if(twall) then                                                     !
twall                                                              !
(source of tem1 at fluid is cp1*rho1*w1*(tem1-:tenv1:)*$           !
mult*:dtdz:/(tem1[1,1,1]-:tenv1:))                                 !
endif                                                              !
                                                                   !
 ************************************************************      !
  Group 13.5 fully-developed-flow temperature-source settings      !
(make maxtem is 0.0)                                               !
(store1 maxtem is tem1[1,1,1])                                     !
                                                                   !
(make maxvel is 0.0)                                               !
(store1 maxvel is w1[1,1,1])                                       !
  (print inlet_vel is :wfin:)                                      !
                                                                   !
(make ohtc is 0.0)                                                 !
(store1 ohtc is heafl4/(anorth[,:nrgy1:,]*(:tfin-tenv1:)))         !
 ************************************************************      !
endif
  ************************************************************                           
    GROUP 15. Termination of sweeps 
 LSWEEP = numits ;ISWC1 = 1
 if(toend) then
 isg21=lsweep
 else
 isg21=1
 endif
 
 numits
 isg21
 
 LITHYD = 1 ;LITFLX = 1 ;LITC = 1 ;ITHC1 = 1
 SELREF = T
   RESFAC =1.0E-03
 ************************************************************      
    GROUP 16. Termination of iterations 
resref(1)=0.0   
selref=f        ! best not used now that residual calculations have
                ! been changed without exploration of consequences
 ************************************************************      

    GROUP 17. Under-relaxation devices 
 OVRRLX =0.
 EXPERT = F ;NNORSL = F
  conwiz=f     ! automatic convergence promotion
  spedat(rlxfac,refvel,r,0.1*:wfin:)
  spedat(rlxfac,reflen,r,zwlast)
 
 relax(w1,linrlx,relxw1)
 relax(tem1,linrlx,relxt1)
 if(fdf) then
 relax(tem1,linrlx,0.2*relxt1)
 endif

    relax(enul,linrlx,0.01)
   relax(enul,linrlx,1.)
if(:turb:.ne.none) then
 relax(enut,linrlx,0.5)
 

endif
 if(:turb:.eq.kemodl) then
 relax(ke,linrlx,0.1)
 relax(ep,linrlx,0.1)
 varmax(ep)= 1.e15
 endif
 ************************************************************      
  mesg(group 18
    GROUP 18. Limits on variables or increments to them 
    how justify these settings?
    
 varmax(tem1)=tenv1+(tfin-tenv1)*10 !****
 varmin(tem1)=tenv1
 varmin(w1)=0.0
 varmax(w1)=2.5*:wfin:
 **********************************************************      
    GROUP 19. Data communicated by satellite to GROUND 
 PARSOL = F
 ************************************************************      
    GROUP 20. Preliminary print-out 
  ! note that the following is rather clumsy, because the inforout-
  ! printing coding allows only reals to be printed, A better method
  ! must be devised 
 ************************************************************    
    GROUP 21. Print-out of variables 
 SUBWGR = F
nprint=lsweep/howmany
    * 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
  if(.not.prprps) then
  output(prps,n,n,y,n,n,n)
  endif    
if(:turb:.ne.none) then  -------------------------------------------
                                                                   !
  if(.not.prenut) then                                               !
  output(enut,n,n,y,n,n,n)                                           
  endif 

  if(.not.prlen1) then                                                !
output(len1,n,n,y,n,n,n)                                            !
  endif                                                              !
                                                                   !
if(turmo.eq.2.) then  !-----------------------------------------                                             !
                                                               !   !
  if(.not.prltls) then                                           !   !
output(ltls,n,n,y,n,n,n)                                       !   !
  endif                                                          !   !
                                                               !   !
  if(.not.prwgap) then                !                          !   !
output(wgap,n,n,y,n,n,n)            !                          !   !
  endif                               !                          !   !
endif    ! end eq 2  -------------------------------------------   !                                                       !
endif    ! end ifturb.ne. none -------------------------------------
 ************************************************************      
    GROUP 22. Spot-value print-out 
 IXMON = 1 ;IYMON = 1 ;IZMON = 1
 NPRMON = 100000 ;NPRMNT = 1 ;TSTSWP = -1
 UWATCH = T ;USTEER = T
 HIGHLO = F
char(moni) 
moni=maxabs
moni=maxmin
#:moni:
 
 ************************************************************      
    GROUP 22.1. Inforin print-out
write(>inforin,SimScene is TubeFlow of March 2013)
write(>>inforin,winset is :winset:)
write(>>inforin,fluid is :flname:)
write(>>inforin,gravty is :gravty:)  
write(>>inforin,turbulence model is :turb:)    
write(>>inforin,flow_form is :floform:)  
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,   )  
  #endpause
  (print wsource is nets(w1,inisoli))
real(xsect)
xsect=0.125*xang*diam**2
  (make force is 0.0)
  (store1 force is (p1[1,1,1]-:pout:)*:xsect:)
  (print force is force)
 
 ************************************************************      
    GROUP 22.2. Inforout print-out
if(prgrge.or.prntall) then ! grid- and geometry-related 
(print used_nrgy1 is :nrgy1:)
(print used_nrgz1 is :nrgz1:)
(print tube_length is :lng1:)
(print tube_inner_diam is :diam:)
if(.not.fdf) then                                  !
(print length/diameter is :zwlast/diam:)           !
endif       ! end not fdf                          !
(print tube_wall_thickness is :thck:) 
endif
if(prrepr.or.prntall) then
(print reference_density is refdns)
(print reference_kin_visc is refkvis)
(print reference_spc_heat is refsph)               !
(print reference_conductivity is refcnd)           !
endif

if(prdino.or.prntall) then  ! dimensionless nos. in inforout
(make reyno is 0.0)
dummy=:wfin:
(store1 reyno is min(1.e10,:dummy*diam:/refkvis)) ! min needed to
(print Reynolds_No is reyno)                     ! prevent overflow
(make prnno is 0.0)
(store1 prnno is refsph*refkvis*refdns/refcnd)
(print Prandtl_No is prnno)
   (print effectiveness is effect)                              
(make nusno is 0.0)                                 
if(fdf) then
(store1 nusno is ohtc*:diam:/refcnd)                
(print Nusselt_No is nusno)                         
(make nuMca is 0.0)
(store1 nuMca is 0.023*reyno^0.8*prnno^0.3333)
(print Nu_from_McAdams is nuMca) 
(print 4*friction_coeff is 4*frccoe)               
(make stanno is 0.0)
(store1 stanno is ohtc/(refdns*refsph*:wfin:))
(print st*pr**2/3 is stanno*prnno^0.667)           
endif
endif

if(prboco.or.prntall) then ! boundary conds.. in inforout
(print inlet_velocity is :wfin:)                    
(print inlet_temperature is :tfin:)                    
(print outlet_pressure is :pout:)                    
endif

if(prntall) then                                                   !
  (print friction_coeff is frccoe)                                   !
 if(fdf) then
 (print spec_mas_flo is :spmsfl:)                                 !
 (print calc_mas_flo is camsfl)                                   !
 (print calc_flow/spec_flow is florat)                              !
 (print presgrad_mult is pgrmul)                                  !
 (print press_grad is pgrad)                                      !
 endif
(print max_temp_ratio is (maxtem-:tenv1:)/:tfin-tenv1:)            !
(print heafl4 is heafl4)                                           !
(print enthalpy_outflow is entflo)      ! print-out not required   !
  (print total_heat_out is entflo+heafl4) ! print-out not required !
(print total_heat_out is entflo) ! print-out not required          !
(print spec_enth_inflow is entfli)   ! print-out not required      !
(print heat_inflow/outflo is heafac)      ! should become unity    !
(print multiplier is mult)              ! print-out not required   !
(print max_vel_ratio is maxvel/:wfin:)                             !
(print 4*friction_coeff is frccoe)                               !
(print frccoe*reyno is fcre)
(print used_lsweep is :lsweep:)
endif                                                              !
 ************************************************************      
    GROUP 23. Field print-out and plot control 
 IDISPA=1;IDISPB=2;IDISPC=NZ
 NUMCLS = 5
 IPLTF = 1 ;IPLTL = -1 ;NPLT = -1
 ISWPRF = 1 ;ISWPRL = 100000
 ITABL = 3 ;IPROF = 1
 ABSIZ =0.5 ;ORSIZ =0.4
 NTZPRF = 1 ;NCOLPF = 50
 ICHR = 2 ;NCOLCO = 45 ;NROWCO = 20
 yzpr=t

iprof=3
patch(crossstr,profil,1,1,1,ny,nz/2,nz/2,1,1)
coval(crossstr,w1,0.0,0.0)   
coval(crossstr,tem1,0.0,0.0)   
if(:turb:.ne.none) Then
coval(crossstr,enut,0.0,0.0)
endif
if(:turb:.eq.lvel) then
coval(crossstr,wdis,0.0,0.0)
endif
if(:turb:.eq.kemodl) then
coval(crossstr,ke,0.0,0.0)
coval(crossstr,len1,0.0,0.0)
endif   

if(nz.gt.1) then   
  patch(axial,profil,1,1,1,1,1,nz,1,1)
  coval(axial,w1,0.0,0.0)   
  coval(axial,tem1,0.0,0.0)   
     
  patch(outersrf,profil,1,1,ny,ny,2,nz,1,1)
  coval(outersrf,#hfl,0.0,0.0)   

  patch(nearwall,profil,1,1,nrgy1,nrgy1,2,nz,1,1) 
  if(fdf) then
  coval(nearwall,shrz,0.0,0.0)                    
  endif
  coval(nearwall,htco,0.0,0.0)                    
     
endif
   inforout
   Use properties at outlet temperature in Reynolds etc numbers
   because these change least

  mesg(inforout  
  Note for nip. It is desirable to be able print integers, logicals 
  and character strings to inforout



                                                    
  (print Euler_No is frccoe*:zwlast/diam:) ! do not print as of no
                                           ! interest for tubes

if(.not.fdf) then                                                  !
real(deltin) ! inlet temperature difference                        !
deltin=tfin-tenv1                                                  !
(make mflo is 0.0)                                                 !
(store1 mflo is nets(r1,inlet))                                    !
(print mass_flow_rate is mflo)                                     !
                                                                   !
(make ohtc is 0.0)                                                 !
real(area)                                                         !
area=xulast*0.5*diam*zwlast                                        !
(store1 ohtc is hflx/(lmtd*:area:))                                !
if(prntall) then                                                   !
(print overall_ht_tr/co is ohtc)                                   !
endif                                                              !
                                                                   !
(make velrat is 0.0)                                               !
(store1 velrat is w1[1,1,9*NZ/10]/:wfin:) ! 9/10 needed ?          !
if(prntall) then                                                   !
(print max/inlet_velocity is velrat)                               !
endif                                                              !
(make deltou is 0.0)                                               !
(store1 deltou is tfou-:tenv1:)                                    !
(print deltou is deltou)
(print deltin is :deltin:)
                                                                   !
(make lmtd is 0.0)                                                 !
(store1 lmtd is (:deltin:-deltou)/loge(:deltin:/deltou))           !
if(prntall) then                                                   !
(print log_mean_t_diff is lmtd)                                    !
endif                                                              !

(make hflx is 0.0)
(store1 hflx is -nets(tem1,tubewal1))
if(prntall) then                                                   !
(print total_heat_flux is hflx)                                    !
                                                                   !
(print inlet_temp. is :tfin:)                                      !
(print external_temp. is :tenv1:)                                  !
(print outlet_temp. is tfou)                                       !
(print effectiveness is (:tfin:-tfou)/:tfin-tenv1:)                !
endif                                                              !
endif     ! end of not fdf -----------------------------------------                                                         !
if(.not.fdf) then !-------------------------------------- 
(make frccoe is 0.0)                                               !
dummy=diam/(zwlast*:wfin:**2)                                        !
formula=:dummy:*(p1[1,1,1]-p1[1,1,nz])/refdns                      !
(store1 frccoe is :formula:)                                       !
(store1 nusno is ohtc*:diam:/refcnd)                
(print Nusselt_No is nusno)                         
(make nuMca is 0.0)
(store1 nuMca is 0.023*reyno^0.8*prnno^0.3333)
(print Nu_from_McAdams is nuMca) 
(print 4*friction_coeff is 4*frccoe)               
(make stanno is 0.0)
(store1 stanno is ohtc/(refdns*refsph*:wfin:))
(print st*pr**2/3 is stanno*prnno^0.667)        
if(prntall) then                                                   !
(print friction_coeff is frccoe)                                   !
endif                                                              !
(print 4*friction_coeff is 4*frccoe)                                   !
                                                                   !
(make tfou is 0.0)                                                 !
(store1 tfou is :tfin:-hflx/(mflo*cp1[,1,])!IF(IZ.EQ.1))           !
endif !  ???? end of not. fdf)  ----------------------------
                                                                   
  izprf=24
  nz=30
  nzprin=1
  izdb1=26
  debug=t
  dbout=t
  flag=t
  fdf
  tstswp=1  
   PHOTON macro
mesg(gravity = :gravty: 
real(vfact)
vfact=zwlast/diam
write(>u,p)                                                                               
write(>>u, )
if(.not.gravty) then   ! 2D only, when gravity is absent
if(fdf) then
write(>>u, )
write(>>u, )
write(>>u,AUTOPLOT )
write(>>u, )
write(>>u, )
write(>>u,file; phi 5 )
write(>>u, )
write(>>u,cl; da 1; w1; colf 1 )
write(>>u,  msg W-velocity profile)
write(>>u,  pause)
write(>>u,  dump; velocrve)
write(>>u,cl; da 1;  tem1; colf 1)
write(>>u,  msg temperature profile)
write(>>u,  pause)
write(>>u,  dump; tempcrve)
write(>>u,cl; da 1;  enul; colf 1)
write(>>u,  msg laminar kinematic viscosity profile)
write(>>u,  pause)
write(>>u,  dump; enulcrve)
if(:turb:.ne.none) then
write(>>u,cl; da 1;  enut; colf 1)
write(>>u,  msg turbulent kinematic viscosity profile)
write(>>u,  pause)
write(>>u,  dump; enutcrve)
endif
write(>>u,msg(                  enter E then RETURN to quit)
write(>>u,enduse)

endif
write(>>u,:vfact: 1)
write(>>u, )
write(>>u,vi -x) 
write(>>u,gr ou x 1 col 1)  
write(>>u,window 0.0 0.95 0.4 0.80)
write(>>u,gr ou x 1 col 1)                                                          
write(>>u,con tem1 x 1 fi;.001)                                                            
write(>>u,pause; dump tem1; con off; red)
if(:flname:.ne.use_specified_properties) then
write(>>u,con enul x 1 y 1 :nrgy1: fi;.001)                                                            
write(>>u,pause; dump enul; con off; red)
endif
write(>>u,con w1 x 1 fi ;.001)                                                           
write(>>u,pause; dump w1; con off; red)
write(>>u,con p1 x 1 fi ;.001)                                                           
write(>>u,pause; dump p1; con off; red)
if(:turb:.ne.none) then
write(>>u,con enut x 1 fi ;.001)                                                           
write(>>u,pause; dump enut; con off; red)
endif
  write(>>u,msg(                  enter E then RETURN to quit)
  write(>>u,enduse)
else  ! 3D when gravity is active
  mesg(about to write u
write(>>u,;;;)
write(>>u, )
write(>>u,vi z)
write(>>u,gr z 1)
write(>>u,pause)
integer(izcon)
if(fdf) then
izcon=1
else
izcon=nz/2
endif
write(>>u,con tem1 z :izcon: fi;.0001)
write(>>u,pause; dump tem1; con off; red)
write(>>u,gr of;red;gr ou z 1)
write(>>u,con u1 z :izcon: fi;.0001)
write(>>u,pause; dump u1; con off; red)
write(>>u,con v1 z :izcon: fi;.0001)
write(>>u,pause; dump v1; con off; red)
write(>>u,con w1 z :izcon: fi;.0001)
write(>>u,pause; dump w1; con off; red)
dummy=:wfin:
write(>>u,set ref vec :0.5*dummy:)
write(>>u,vec z :izcon: sh)
write(>>u,pause; dump vec; con off; red)
write(>>u,msg(                  enter E then RETURN to quit)
write(>>,enduse
endif
 ************************************************************                                                              
    GROUP 24. Dumps for restarts
 SAVE = T ;NOWIPE = F
 NSAVE =CHAM
   ************************************************************                                                              
    GROUP 25. VR-related settings
   ************************************************************                                                              
    GROUP 26. Graphical-display macros
   ************************************************************       
endit(tem1)=0.0
resref(tem1)=0.0
nyprin=1
iyprf=ny-10
STOP