A shell-and-tube heat-exchanger is used to exemplify the essential ideas of the HEXAGON model devised by D.B Spalding and S.V. Zhubrin.
This may be the first to have shown how the thermo-hydraulics of the shell-side and tube-side fluids can be computed simultaneously.
Moreover the displacements and thermal stresses in tubes and shell are also computed as part of a full SFT (i.e. Solid-Fluid-Thermal) analysis .
The heat exchanger considered is an idealised two-dimensional one, having two baffles within the shell and U-bend tubes connected with a two-part header.
The geometry is as shown.
The tube fluid is water and the shell fluid is air. Their flow patterns are indicated by the calculated velocities, where the white vectors and coloured vectors represent shell- and tube-side fluid respectively.
The MUSES method is employed, the whole volume occupied by the heat exchanger being represented three times, namely for:
The following figures show calculated displacements and
the Y-component of normal thermal stresses.
The model is set up wholly by way of a standard PHOENICS Q1 file, which contains PLANT statements (listed below). These:
Therafter, when the satellite is run, a corresponding GROUND file is written and compiled, without further user intervention. This file is approximately 550 lines in length, and error-free! Obviously, the advantages to the user are very great.
The user puts into the Q1 file in Group 9:
RHO1=GRND
{PRPT01} DEN1=1000.
REGION() 1
{PRPT02} DEN1=1000.
REGION() 2
{PRPT03} DEN1=1.2
REGION() 3
The above three statements, followed by the pointer RHO1=GRND and parameterized REGION qualifiers, instruct PLANT to make the density distributions as the distribution of in-cell marker values dictates.
ENUL=GRND
{PRPT05} VISL=1.*SQRT(U1**2+V1**2)*WDIS
REGION() 1
{PRPT06} VISL=0.01
REGION() 2
{PRPT07} VISL=1.*SQRT(U1**2+V1**2)*WDIS
REGION() 3
The above three statements do the same for viscosities as the previous three have done for densities. Note that the viscosities in the domains marked 1 and 3 are made proportional to the products of local velocity magnitudes and distances to the nearest wall.
The user puts into the Q1 file in Group 13:
PATCH(SS002U,PHASEM,1,10,1,NY,1,NZ,1,1)
{SORC01} CO=.2*(U1**2+V1**2)**0.15
COVAL(SS002U,U1,GRND,0.0)
{SORC02} CO=.2*(U1**2+V1**2)**0.15
COVAL(SS002U,V1,GRND,0.0)
Momentum sinks are introduced by above formulae over all cells having marker value appearing in the number of PATCH name, 002.
PATCH(SS003H,PHASEM,1,10,1,NY,1,NZ,1,1)
{SORC03} CO=2.2*(U1**2+V1**2)**0.25
COVAL(SS003H,U1,GRND,0.0)
{SORC04} CO=2.2*(U1**2+V1**2)**0.25
COVAL(SS003H,V1,GRND,0.0)
Momentum sinks are introduced by above formulae over all cells having marker value appearing in the number of PATCH name, 003.
** Tube fluid heat transfer coefficient
STORE(ALF2);FIINIT(ALF2)=0.0
{SC0601} ALF2=1.+1.*SQRT(U1**2+V1**2+TINY)
REGION() 2
** Shell fluid heat transfer coefficient
STORE(ALF3);FIINIT(ALF3)=0.0
{SC0602} ALF3=1.+3.*SQRT(U1**2+V1**2+TINY)
REGION() 3
Tube and shell fluid heat transfer coefficients, ALF2 and ALF3, are made by aboves dependent on local velocity magnitudes over the REGIONs marked 2 and 3 correspondingly.
** Overall heat transfer coefficient
STORE(HTC);FIINIT(HTC)=0.0
{SC0603} HTC=1./(1/ALF2+1/ALF3[,-IG(1),])
REGION() 2
{SC0604} HTC=1./(1/ALF3 +1/ALF2[,+IG(1),])
REGION() 3
Overall heat transfer coefficient, HTC, distribution are calculated by reference to appropriate local heat transfer coefficients in REGIONs 2 and 3.
** Heat-exchange with shell-fluid, throughout the shell.
PATCH(SS002T,PHASEM,1,NX,1,NY,1,NZ,1,1)
{SORC05} CO =HTC
{SORC05} VAL=TEMP[,-IG(1),]
COVAL(SS002T,TEMP,GRND,GRND)
** Heat-exchange with tube-fluid, throughout the shell.
PATCH(SS003S,PHASEM,1,NX,1,NY,1,NZ,1,1)
{SORC06} CO =HTC
{SORC06} VAL=TEMP[,+IG(1),]
COVAL(SS003S,TEMP,GRND,GRND)
In the above, the PATCH names indicate the sub-domain cell markers , 2 and 3, over which the heat-exchange sources are applied. The indicial operations for TEMP are arranged in appropriate manner.
** Transfer shell fluid temperatures
PATCH(SS005T,CELL,1,NX,1,NY,1,NZ,1,1)
{SORC10} CO=1.e10
{SORC10} VAL=TEMP[-13,-6,]
COVAL(SS005T,TEM1,GRND,GRND)
The temperatures of the shell fluid, TEMP, are transferred into the stress-analysis sub-domain, MARK=5, to be used as TEM1.
** Tube-wall temperature
PATCH(SS100T,CELL,1,NX,1,NY,1,NZ,1,1)
{SORC11} CO=1.e10
{SORC11} VAL=(ALF2[-13,+6,]*TEMP[-13,+6,]+$
ALF3[-13,-6,]*TEMP[-13,-6,])$
/(ALF2[-13,+6,]+ALF3[-13,-6,]+TINY)
COVAL(SS100T,TEM1,GRND,GRND)
Here the tube wall temperatures, TEM1, are calculated in the sub-domain indicated by MARK=100 as PATCH name number specifies. TEM1s are computed via shell and tube fluid temperatures and heat transfer coefficients transferred from corresponding sub-domains as indicial numbers show.
** Transfer the header temperatures
PATCH(SS004T,CELL,1,NX,1,NY,1,NZ,1,1)
{SORC12} CO=1.e10
{SORC12} VAL=TEMP[-13,+3,]
COVAL(SS004T,TEM1,GRND,GRND)
The temperatures of the header-tube fluid, TEMP, are transferred into the stress-analysis sub-domain, MARK=4, to be used as TEM1.
The user puts into the Q1 file in Group 19:
** Tube-fluid velocities transfer
STORE(UU1,VV1)
{SC0605} UU1=U1[-13,+3,]
REGION(14,23,22,29,1,NZ)
{SC0606} VV1=V1[-13,+3,]
REGION(14,23,22,29,1,NZ)
{SC0607} UU1=U1[-13,+6,]
REGION(14,23,7,18,1,NZ)
{SC0608} VV1=V1[-13,+6,]
REGION(14,23,7,18,1,NZ)
Tube fluid velocities are transferred from where they have been calculated for easy visualisation.
** Shell fluid velocities transfer
STORE(U2,V2)
{SC0609} U2=U1[-13,-6,]
REGION(14,23,7,18,1,NZ)
{SC0610} V2=V1[-13,-6,]
REGION(14,23,7,18,1,NZ)
Shell-fluid velocities are transferred from where they have been calculated for easy visualisation.
C Property name: PRPT01
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFMARK=L0F(INAME('MARK'))
LFDEN1 =L0F(AUX(DEN1 ))
DO 90101 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 90101 IY=IYF ,IYL
I=IY+IADD
L0DEN1=LFDEN1+I
L0MARK=LFMARK+I
INMARK=NINT(F(L0MARK))
IF(INMARK.EQ.1 ) THEN
F(L0DEN1 )=1000.
ENDIF
90101 CONTINUE
ENDIF
ENDIF
C Property name: PRPT02
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFMARK=L0F(INAME('MARK'))
LFDEN1 =L0F(AUX(DEN1 ))
DO 90102 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 90102 IY=IYF ,IYL
I=IY+IADD
L0DEN1=LFDEN1+I
L0MARK=LFMARK+I
INMARK=NINT(F(L0MARK))
IF(INMARK.EQ.2 ) THEN
F(L0DEN1 )=1000.
ENDIF
90102 CONTINUE
ENDIF
ENDIF
C Property name: PRPT03
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFMARK=L0F(INAME('MARK'))
LFDEN1 =L0F(AUX(DEN1 ))
DO 90103 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 90103 IY=IYF ,IYL
I=IY+IADD
L0DEN1=LFDEN1+I
L0MARK=LFMARK+I
INMARK=NINT(F(L0MARK))
IF(INMARK.EQ.3 ) THEN
F(L0DEN1 )=1.2
ENDIF
90103 CONTINUE
ENDIF
ENDIF
C Property name: PRPT05
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFMARK=L0F(INAME('MARK'))
LFVISL =L0F(AUX(VISL ))
LFU1 =L0F(U1 )
LFV1 =L0F(V1 )
LFWDIS=L0F(INAME('WDIS '))
DO 90605 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 90605 IY=IYF ,IYL
I=IY+IADD
L0VISL=LFVISL+I
L0MARK=LFMARK+I
INMARK=NINT(F(L0MARK))
IF(INMARK.EQ.1 ) THEN
L0U1 =LFU1 +I
L0V1 =LFV1 +I
L0WDIS=LFWDIS+I
F(L0VISL )=1.*SQRT(F(L0U1)**2+F(L0V1)**2)*F(L0WDIS)
ENDIF
90605 CONTINUE
ENDIF
ENDIF
C Property name: PRPT06
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFMARK=L0F(INAME('MARK'))
LFVISL =L0F(AUX(VISL ))
DO 90606 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 90606 IY=IYF ,IYL
I=IY+IADD
L0VISL=LFVISL+I
L0MARK=LFMARK+I
INMARK=NINT(F(L0MARK))
IF(INMARK.EQ.2 ) THEN
F(L0VISL )=0.01
ENDIF
90606 CONTINUE
ENDIF
ENDIF
C Property name: PRPT07
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFMARK=L0F(INAME('MARK'))
LFVISL =L0F(AUX(VISL ))
LFU1 =L0F(U1 )
LFV1 =L0F(V1 )
LFWDIS=L0F(INAME('WDIS '))
DO 90607 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 90607 IY=IYF ,IYL
I=IY+IADD
L0VISL=LFVISL+I
L0MARK=LFMARK+I
INMARK=NINT(F(L0MARK))
IF(INMARK.EQ.3 ) THEN
L0U1 =LFU1 +I
L0V1 =LFV1 +I
L0WDIS=LFWDIS+I
F(L0VISL )=1.*SQRT(F(L0U1)**2+F(L0V1)**2)*F(L0WDIS)
ENDIF
90607 CONTINUE
ENDIF
ENDIF
C Source name: SORC01
IF(INDVAR.EQ.INAME('U1 ').AND.NPATCH.EQ.'SS002U ') THEN
LFCO =L0F(CO )
LFMARK=L0F(INAME('MARK'))
LFU1 =L0F(U1 )
LFV1 =L0F(V1 )
DO 13701 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13701 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.002) THEN
L0U1 =LFU1 +I
L0V1 =LFV1 +I
F(LFCO +I)=.2*(F(L0U1)**2+F(L0V1)**2)**0.15
ELSE
F(LFCO +I)=0.
ENDIF
13701 CONTINUE
ENDIF
C Source name: SORC02
IF(INDVAR.EQ.INAME('V1 ').AND.NPATCH.EQ.'SS002U ') THEN
LFCO =L0F(CO )
LFMARK=L0F(INAME('MARK'))
LFU1 =L0F(U1 )
LFV1 =L0F(V1 )
DO 13702 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13702 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.002) THEN
L0U1 =LFU1 +I
L0V1 =LFV1 +I
F(LFCO +I)=.2*(F(L0U1)**2+F(L0V1)**2)**0.15
ELSE
F(LFCO +I)=0.
ENDIF
13702 CONTINUE
ENDIF
C Source name: SORC03
IF(INDVAR.EQ.INAME('U1 ').AND.NPATCH.EQ.'SS003H ') THEN
LFCO =L0F(CO )
LFMARK=L0F(INAME('MARK'))
LFU1 =L0F(U1 )
LFV1 =L0F(V1 )
DO 13703 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13703 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.003) THEN
L0U1 =LFU1 +I
L0V1 =LFV1 +I
F(LFCO +I)=2.2*(F(L0U1)**2+F(L0V1)**2)**0.25
ELSE
F(LFCO +I)=0.
ENDIF
13703 CONTINUE
ENDIF
C Source name: SORC04
IF(INDVAR.EQ.INAME('V1 ').AND.NPATCH.EQ.'SS003H ') THEN
LFCO =L0F(CO )
LFMARK=L0F(INAME('MARK'))
LFU1 =L0F(U1 )
LFV1 =L0F(V1 )
DO 13704 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13704 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.003) THEN
L0U1 =LFU1 +I
L0V1 =LFV1 +I
F(LFCO +I)=2.2*(F(L0U1)**2+F(L0V1)**2)**0.25
ELSE
F(LFCO +I)=0.
ENDIF
13704 CONTINUE
ENDIF
C Source name: SORC05
IF(INDVAR.EQ.INAME('TEMP ').AND.NPATCH.EQ.'SS002T ') THEN
LFCO =L0F(CO )
LFMARK=L0F(INAME('MARK'))
LFHTC =L0F(INAME('HTC '))
DO 13705 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13705 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.002) THEN
L0HTC =LFHTC +I
F(LFCO +I)=F(L0HTC)
ELSE
F(LFCO +I)=0.
ENDIF
13705 CONTINUE
ENDIF
C Source name: SORC06
IF(INDVAR.EQ.INAME('TEMP ').AND.NPATCH.EQ.'SS003S ') THEN
LFCO =L0F(CO )
LFMARK=L0F(INAME('MARK'))
LFHTC =L0F(INAME('HTC '))
DO 13706 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13706 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.003) THEN
L0HTC =LFHTC +I
F(LFCO +I)=F(L0HTC)
ELSE
F(LFCO +I)=0.
ENDIF
13706 CONTINUE
ENDIF
C Source name: SORC10
IF(INDVAR.EQ.INAME('TEM1 ').AND.NPATCH.EQ.'SS005T ') THEN
LFCO =L0F(CO )
LFMARK=L0F(INAME('MARK'))
DO 13710 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13710 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.005) THEN
F(LFCO +I)=1.E10
ELSE
F(LFCO +I)=0.
ENDIF
13710 CONTINUE
ENDIF
C Source name: SORC11
IF(INDVAR.EQ.INAME('TEM1 ').AND.NPATCH.EQ.'SS100T ') THEN
LFCO =L0F(CO )
LFMARK=L0F(INAME('MARK'))
DO 13711 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13711 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.100) THEN
F(LFCO +I)=1.E10
ELSE
F(LFCO +I)=0.
ENDIF
13711 CONTINUE
ENDIF
C Source name: SORC12
IF(INDVAR.EQ.INAME('TEM1 ').AND.NPATCH.EQ.'SS004T ') THEN
LFCO =L0F(CO )
LFMARK=L0F(INAME('MARK'))
DO 13712 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13712 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.004) THEN
F(LFCO +I)=1.E10
ELSE
F(LFCO +I)=0.
ENDIF
13712 CONTINUE
ENDIF
C Source name: SORC05
IF(INDVAR.EQ.INAME('TEMP ').AND.NPATCH.EQ.'SS002T ') THEN
LFCO =L0F(CO)
LFVAL =L0F(VAL)
LFMARK=L0F(INAME('MARK'))
LFTEMP=L0F(INAME('TEMP '))
DO 13805 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13805 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.002) THEN
L0TEMP=LFTEMP+I
F(LFVAL+I)=F(L0TEMP-IG(1))
ELSE
F(LFVAL+I)=0.
F(LFCO +I)=0.
ENDIF
13805 CONTINUE
ENDIF
C Source name: SORC06
IF(INDVAR.EQ.INAME('TEMP ').AND.NPATCH.EQ.'SS003S ') THEN
LFCO =L0F(CO)
LFVAL =L0F(VAL)
LFMARK=L0F(INAME('MARK'))
LFTEMP=L0F(INAME('TEMP '))
DO 13806 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13806 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.003) THEN
L0TEMP=LFTEMP+I
F(LFVAL+I)=F(L0TEMP+IG(1))
ELSE
F(LFVAL+I)=0.
F(LFCO +I)=0.
ENDIF
13806 CONTINUE
ENDIF
C Source name: SORC10
IF(INDVAR.EQ.INAME('TEM1 ').AND.NPATCH.EQ.'SS005T ') THEN
LFCO =L0F(CO)
LFVAL =L0F(VAL)
LFMARK=L0F(INAME('MARK'))
LFTEMP=L0F(INAME('TEMP '))
DO 13810 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13810 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.005) THEN
L0TEMP=LFTEMP+I
F(LFVAL+I)=F(L0TEMP-NY*13-6)
ELSE
F(LFVAL+I)=0.
F(LFCO +I)=0.
ENDIF
13810 CONTINUE
ENDIF
C Source name: SORC11
IF(INDVAR.EQ.INAME('TEM1 ').AND.NPATCH.EQ.'SS100T ') THEN
LFCO =L0F(CO)
LFVAL =L0F(VAL)
LFMARK=L0F(INAME('MARK'))
LFALF2=L0F(INAME('ALF2 '))
LFTEMP=L0F(INAME('TEMP '))
LFALF3=L0F(INAME('ALF3 '))
DO 13811 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13811 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.100) THEN
L0ALF2=LFALF2+I
L0TEMP=LFTEMP+I
L0ALF3=LFALF3+I
F(LFVAL+I)=(F(L0ALF2-NY*13+6)*F(L0TEMP-NY*13+6)+
1F(L0ALF3-NY*13-6)*F(L0TEMP-NY*13-6))/(F(L0ALF2-NY*13+6)+
1F(L0ALF3-NY*13-6)+TINY)
ELSE
F(LFVAL+I)=0.
F(LFCO +I)=0.
ENDIF
13811 CONTINUE
ENDIF
C Source name: SORC12
IF(INDVAR.EQ.INAME('TEM1 ').AND.NPATCH.EQ.'SS004T ') THEN
LFCO =L0F(CO)
LFVAL =L0F(VAL)
LFMARK=L0F(INAME('MARK'))
LFTEMP=L0F(INAME('TEMP '))
DO 13812 IX=IXF ,IXL
IADD=NY*(IX-1)
DO 13812 IY=IYF ,IYL
I=IY+IADD
INMARK=NINT(F(LFMARK+I))
IF(INMARK.EQ.004) THEN
L0TEMP=LFTEMP+I
F(LFVAL+I)=F(L0TEMP-NY*13+3)
ELSE
F(LFVAL+I)=0.
F(LFCO +I)=0.
ENDIF
13812 CONTINUE
ENDIF
C Special calls name: SC0601
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFMARK=L0F(INAME('MARK'))
LFALF2 =L0F(INAME('ALF2 '))
LFU1 =L0F(U1 )
LFV1 =L0F(V1 )
DO 19601 IX=1 ,NX
IADD=NY*(IX-1)
DO 19601 IY=1 ,NY
I=IY+IADD
L0ALF2=LFALF2+I
L0MARK=LFMARK+I
INMARK=NINT(F(L0MARK))
IF(INMARK.EQ.2 ) THEN
L0U1 =LFU1 +I
L0V1 =LFV1 +I
F(L0ALF2 )=1.+1.*SQRT(F(L0U1)**2+F(L0V1)**2+TINY)
ENDIF
19601 CONTINUE
ENDIF
ENDIF
C Special calls name: SC0602
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFMARK=L0F(INAME('MARK'))
LFALF3 =L0F(INAME('ALF3 '))
LFU1 =L0F(U1 )
LFV1 =L0F(V1 )
DO 19602 IX=1 ,NX
IADD=NY*(IX-1)
DO 19602 IY=1 ,NY
I=IY+IADD
L0ALF3=LFALF3+I
L0MARK=LFMARK+I
INMARK=NINT(F(L0MARK))
IF(INMARK.EQ.3 ) THEN
L0U1 =LFU1 +I
L0V1 =LFV1 +I
F(L0ALF3 )=1.+3.*SQRT(F(L0U1)**2+F(L0V1)**2+TINY)
ENDIF
19602 CONTINUE
ENDIF
ENDIF
C Special calls name: SC0603
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFMARK=L0F(INAME('MARK'))
LFHTC =L0F(INAME('HTC '))
LFALF2=L0F(INAME('ALF2 '))
LFALF3=L0F(INAME('ALF3 '))
DO 19603 IX=1 ,NX
IADD=NY*(IX-1)
DO 19603 IY=1 ,NY
I=IY+IADD
L0HTC =LFHTC +I
L0MARK=LFMARK+I
INMARK=NINT(F(L0MARK))
IF(INMARK.EQ.2 ) THEN
L0ALF2=LFALF2+I
L0ALF3=LFALF3+I
F(L0HTC )=1./(1/F(L0ALF2)+1/F(L0ALF3-IG(1)))
ENDIF
19603 CONTINUE
ENDIF
ENDIF
C Special calls name: SC0604
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFMARK=L0F(INAME('MARK'))
LFHTC =L0F(INAME('HTC '))
LFALF3=L0F(INAME('ALF3 '))
LFALF2=L0F(INAME('ALF2 '))
DO 19604 IX=1 ,NX
IADD=NY*(IX-1)
DO 19604 IY=1 ,NY
I=IY+IADD
L0HTC =LFHTC +I
L0MARK=LFMARK+I
INMARK=NINT(F(L0MARK))
IF(INMARK.EQ.3 ) THEN
L0ALF3=LFALF3+I
L0ALF2=LFALF2+I
F(L0HTC )=1./(1/F(L0ALF3)+1/F(L0ALF2+IG(1)))
ENDIF
19604 CONTINUE
ENDIF
ENDIF
C Special calls name: SC0605
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFUU1 =L0F(INAME('UU1 '))
LFU1 =L0F(U1 )
DO 19605 IX=14 ,23
IADD=NY*(IX-1)
DO 19605 IY=22 ,29
I=IY+IADD
L0UU1 =LFUU1 +I
L0U1 =LFU1 +I
19605 F(L0UU1 )=F(L0U1-NY*13+3)
ENDIF
ENDIF
C Special calls name: SC0606
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFVV1 =L0F(INAME('VV1 '))
LFV1 =L0F(V1 )
DO 19606 IX=14 ,23
IADD=NY*(IX-1)
DO 19606 IY=22 ,29
I=IY+IADD
L0VV1 =LFVV1 +I
L0V1 =LFV1 +I
19606 F(L0VV1 )=F(L0V1-NY*13+3)
ENDIF
ENDIF
C Special calls name: SC0607
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFUU1 =L0F(INAME('UU1 '))
LFU1 =L0F(U1 )
DO 19607 IX=14 ,23
IADD=NY*(IX-1)
DO 19607 IY=7 ,18
I=IY+IADD
L0UU1 =LFUU1 +I
L0U1 =LFU1 +I
19607 F(L0UU1 )=F(L0U1-NY*13+6)
ENDIF
ENDIF
C Special calls name: SC0608
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFVV1 =L0F(INAME('VV1 '))
LFV1 =L0F(V1 )
DO 19608 IX=14 ,23
IADD=NY*(IX-1)
DO 19608 IY=7 ,18
I=IY+IADD
L0VV1 =LFVV1 +I
L0V1 =LFV1 +I
19608 F(L0VV1 )=F(L0V1-NY*13+6)
ENDIF
ENDIF
C Special calls name: SC0609
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFU2 =L0F(U2 )
LFU1 =L0F(U1 )
DO 19609 IX=14 ,23
IADD=NY*(IX-1)
DO 19609 IY=7 ,18
I=IY+IADD
L0U2 =LFU2 +I
L0U1 =LFU1 +I
19609 F(L0U2 )=F(L0U1-NY*13-6)
ENDIF
ENDIF
C Special calls name: SC0610
IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN
IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN
LFV2 =L0F(V2 )
LFV1 =L0F(V1 )
DO 19610 IX=14 ,23
IADD=NY*(IX-1)
DO 19610 IY=7 ,18
I=IY+IADD
L0V2 =LFV2 +I
L0V1 =LFV1 +I
19610 F(L0V2 )=F(L0V1-NY*13-6)
ENDIF
ENDIF