cgxaslp.for
cFILENAME GXASLP.FOR 291118
C***********************************************************************
C The Algebraic Slip Model is activated in the Q1 file by:
C
C (1) setting ASLP=T
C
C (2) defining volume fraction variable PT* using
C NAME(C1)=PT0 (continuous phase)
C NAME(C2)=PT1, etc.
C and solving using
C SOLVE(PT0)
C SOLVE(PT1), etc.
C
C (3) setting GALA=T
C
C (4) setting SOLVE(VFOL) and TERMS(VFOL,N,N,N,N,P,P) if inflows
C exist (necessary because GALA is used).
C
C (5) setting RHO1=GRND
C RHO2=continuous phase density
C
C (6) setting (if required) ENUL=GRND and
C PRNDTL(PT0)=laminar kinematic viscosity of continuous phase
C
C (7) defining data relating to each dispersed phase:
C PHINT (PT*) = density
C CINT (PT*) = particle, bubble or droplet diameter
C PRNDTL(PT*) = laminar kinematic viscosity (if ENUL=GRND)
C
C (8) using ASM as the first three characters of inlet and outlet
C patch names and specifying incoming VFOL values as the
C reciprocal of the incoming density (compatible with incoming
C PT* values)
C
C (9) setting parameters relating to cell and slab iterations of
C the volume fraction equations using LITER, ENDIT and LINRLX
C values: PT1 is used for the slab iterations, PT0 for the cell
C ones. So for example,
C LITER(PT0)=5; LITER(PT1)=10
C ENDIT(PT0)=1.E-8; ENDIT(PT1)=1.E-6
C RELAX(PT0,LINRLX,0.5); RELAX(PT1,LINRLX,0.3)
C will ensure that the inner (cell) iteration is carried out
C 5 times (or until the sum of the component volume fraction
C changes is less than 1.e-8 in an iteration) using a linear
C relaxation of 0.5, and that the outer (slab) iteration is
C carried out 10 times (or until the sum of the component
C imbalances across the slab is less than 1.e-6) using a linear
C relaxation of 0.3.
C
C Note that because all the volume fractions are solved
C together, residual values are available only for PT0,
C corresponding to the sum of all the component imbalances.
C
C (10) specifying LASLPA and LASLPB:
C LASLPB=T reduces computation time by allowing slip velocity
C only in directions in which gravitational or
C rotational forces are activated using BUOY or ROTA
C patches.
C LASLPA=T uses a generally quicker but less stable iterative
C scheme.
C (11) STORE(RSPT,RCPT) will store the slab and cell erros
C STORE(NIT1,NIT2) will store the cell and slab iteration counters
C
C GXASLP is called from within EARTH at times appropriate to
C GROUND Group 1 Section 2
C GROUND Group 1 Section 3
C GROUND Group 9 Section 1
C GROUND Group 1 Section 6
C GROUND Group 19 Section 6
C
C***********************************************************************
c
SUBROUTINE GXASLP
INCLUDE 'patnos'
INCLUDE 'patcmn'
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'grdear'
INCLUDE 'satgrd'
C
COMMON/ICOVL/M04,I0PHI
COMMON/INDAUX/L0ISL,L0IST,L0SL,L0ST,NSTO,NSOL,L0SLRS,L0TTRS,IFL(3)
COMMON/GENI/NXNY,IGNI1(8),NFM,IGNI11(50)
COMMON/NAMFN/NAMFUN,NAMSUB
COMMON/CONVL/IPT0,IPTL,NPT,L0NORD
COMMON/INTDMN/IDMN,NUMDMN,NXMAX,NYMAX,NZMAX,NXYMAX,NXYZMX,NCLTOT,
1 NDOMMX
COMMON/VRMCMN/L0MSKZ
COMMON /FACES/L0FACE,L0FACZ
LOGICAL NEZ,SLIPU,SLIPV,SLIPW,BDYFX,BDYFY,BDYFZ
COMMON/ASML/LBUOYX,LBUOYY,LBUOYZ,FEXIST(7),BEXIST
LOGICAL LBUOYX,LBUOYY,LBUOYZ
LOGICAL FEXIST,BEXIST,SLD
LOGICAL LVRCFD,MASKPT,L2DVRO,LINVRO,dbnow
REAL NORML(3)
CHARACTER*6 NAMFUN,NAMSUB,CVAR*68
COMMON /TIMSTO/ TIMIN(150),TIMOU(150),PSORIN(2),PSOROU(2),
1 TIMINT(2),TIMOUT(2),VSORIN(2),VSOROU(2)
DOUBLE PRECISION TIMIN,TIMOU,PSORIN,PSOROU,TIMFLIN,TIMFLOU,TIMINT,
1 TIMOUT,VSORIN,VSOROU
C DENS, DIAM and VISC store phase densities, diameters and laminar
c viscosities.
C
SAVE NIT1,NIT2,CRT1,CRT2,L0BDYX,L0BDYY,L0BDYZ,L0BDFL,L0V0X,L0V0Y,
1 L0V0Z,L0V0ZL,L0CHNG,IPT0M1,IPTF,
1 L0SLPS,CD0,SLIPU,SLIPV,SLIPW,BDYFX,BDYFY,BDYFZ,IPTSUM,L0AE,
1 L0AN,L0AH,L0AL,L0MIN,L0MOU,LBVFOL,L0RACC,L0FACL,L0OSLB,
1 L0OCEL,L0OSWP,L0SLP1,L0SLP2,L0SLP3,L0SLP4,L0SLP5,L0SLP6,L0INFL,
1 TERM1,TERM2,TERM3,RECRH2,IASMP1,IASMP2,L0CV4,L0COP,L0COV,
1 L0INFPT,L0OUTPT,L0FPT,L0FPTO,M0L0F,M0L0FO,M0CV4,L0PATN,L0PAT,
1 LBPTEC,LBPTES,LBNIT1,LBNIT2
C
c.... arithmetic-statement functions
DENS(I)=PHINT(I)
DIAM(I)=CINT(I)
VISC(I)=PRNDTL(I)
c
cc call accutm('gxaslp',1)
c write(14,*)'varbl,isw,igr,isc,iz ',
c 1 name(indvar),isweep,igr,isc,iz
NAMSUB='GXASLP'
C>>>>>>>>>>>>>>>>>>>>> Group 1 ------- Section 3 <<<<<<<<<<<<<<<<<<<<<
C DENS (I) = density of ith phase
C DIAM (I) = diameter of ith phase
C VISC (I) = viscosity of ith phase
C L0BDYX = force vector resolved in local xc direction
C L0BDYY = force vector resolved in local yc direction
C L0BDYZ = force vector resolved in local zc direction
C.......................................................................
dbnow=.false.
if(dbnow) write(14,*) 'in gxaslp with igr, isc =',
1 igr,isc
IF(IGR.EQ.1.AND.ISC.EQ.3) THEN
C Set drag coefficient for high slip Reynolds number
CD0 = 0.42
C
C Provide storage for GRound-earth SPare arrays,
IZLAST=-1; JSWEEP=-1
CALL GXMAKE0(L0AL ,NXNY,'AREL')
CALL GXMAKE0(IPTSUM,NXNY,'ISUM')
CALL GXMAKE0(L0RACC ,NXNY,'RACC ')
IF(NZ.GT.1) CALL GXMAKE0(L0FACL,NXNY,'FACL')
IF(NX.GT.1.OR.SOLVE(U1).AND..NOT.CARTES)
1 SLIPU=.TRUE.
IF(NY.GT.1) SLIPV=.TRUE.
IF(NZ.GT.1) SLIPW=.TRUE.
IF(SLIPU) CALL GXMAKE0(L0V0X,NXNY,'I_U0')
IF(SLIPV) CALL GXMAKE0(L0V0Y,NXNY,'I_V0')
IF(SLIPW) CALL GXMAKE0(L0V0Z,NXNY,'I_W0')
IF(SLIPW) CALL GXMAKE0(L0V0ZL,NXNY,'ILW0')
C
C Determine the number of the particle equations
C solved, then fill arrays rhop and DIAM with densities
C and diameters of various groups of particles.
C
C NPT counts the total number of particles, iptf is the indvar
C of the first group of particles, and iptl is the indvar of
C the last group of particles.
C
CALL SUB2(NPT,0,IPTF,0)
DO 10000 IPT = 16 , NPHI
IF(SOLVE(IPT).AND.NAME(IPT)(1:2).EQ.'PT') THEN
IF(IPTF.EQ.0) THEN
IPT0=IPT
IPTF=IPT+1
ENDIF
IPTL = IPT
NPT = NPT+1
ENDIF
10000 CONTINUE
IF(.NOT.NULLPR) THEN
CALL WRITBL
CALL WRYT40('ASM of 29/11/18 has been called ')
CALL WRIT40('ASM of 29/11/18 has been called ')
CALL WRITBL
IF(NUMDMN.GT.1) THEN
CALL WRIT40('ASM not compatible with multi-domain!! ')
CALL WRYT40('ASM not compatible with multi-domain!! ')
CALL SET_ERR(506,
* 'Error. ASM not compatible with multi-domain.',2)
ENDIF
WRITE(14,'('' Number of phases, including continuous = '',
1 I3)') NPT
ENDIF
IPT0M1 = IPT0-1
C
CALL GXMAKE0(L0OSLB,NPT*NXNY,'OSLB')
CALL GXMAKE0(L0OSWP,NPT*NXNY,'OSWP')
CALL GXMAKE0(L0OCEL,NPT,'OCEL')
CALL GXMAKE0(L0CHNG,NPT,'CHNG')
CALL GXMAKE0(L0INFL,NPT,'INFL')
CALL GXMAKE0(L0FPT,NPT,'L0FP')
IF(.NOT.STEADY) CALL GXMAKE0(L0FPTO,NPT,'L0FO')
L0CHNG=L0CHNG+1-IPT0
CALL GXMAKE0(L0SLPS,NPT*7,'SLPS')
CALL SUB4(L0SLP1,L0SLPS, L0SLP2,L0SLPS+NPT, L0SLP3,L0SLPS+2*NPT,
1 L0SLP4, L0SLPS+3*NPT)
CALL SUB3(L0SLP5,L0SLPS+4*NPT, L0SLP6,L0SLPS+5*NPT, L0NORD,
1 L0SLPS+6*NPT)
IF(SOLVE(LBNAME('VFOL'))) LBVFOL = LBNAME('VFOL')
LBPTEC=LBNAME('RCPT'); LBPTES=LBNAME('RSPT')
LBNIT1=LBNAME('NIT1'); LBNIT2=LBNAME('NIT2')
C
Come here to insert the values of the particle densities
C and diameters.
C
NIT1 = LITER(IPT0); NIT2 = LITER(IPT0+1) ! cell & slab iterations
CRT1 = ENDIT(IPT0); CRT2 = ENDIT(IPT0+1) ! cell & slab termination criteria
RLX1 = ABS(DTFALS(IPT0)); RLX2 = ABS(DTFALS(IPT0+1)) ! cell and slab relaxations
C
DO IPT = IPT0,IPTL
c.... switch off solution for ipt
ISLN(IPT) = ISLN(IPT)/3
F(L0ISL+IPT) = - F(L0ISL+IPT)
F(L0TTRS+IABS(ISL(IPT))) = 0.0
ENDDO
IF(.NOT.NULLPR) THEN
CALL WRITBL
CALL WRIT40('Particle Properties')
DO IPT = IPTF,IPTL
CALL WRIT1I('Particle no ',ipt-ipt0)
CALL WRIT3R('density ',dens(ipt),';diameter ',diam(ipt),
1 'viscosity ',visc(ipt))
ENDDO
CALL WRITBL
CALL WRIT40('Carrier Properties')
CALL WRIT2R('density ',RHO2,';kinematic viscosity',
1 PRNDTL(IPT0))
CALL WRITBL
CALL WRIT40('Particle Numerical Settings')
CALL WRIT1I('Cell iterations (LITER(PT0)) ', NIT1)
CALL WRIT1I('Slab iterations (LITER(PT1)) ', NIT2)
CALL WRIT1R('Cell relaxation (DTFALS(PT0)) ', RLX1)
CALL WRIT1R('Slab relaxation (DTFALS(PT1)) ', RLX2)
CALL WRIT1R('Cell cut-off (ENDIT(PT0)) ', CRT1)
CALL WRIT1R('Cell cut-off (ENDIT(PT0)) ', CRT2)
CALL WRITST; CALL WRITBL
ENDIF
TERM1=CD0 / (432. * VISC(IPT0) ** 2)
TERM2=1/CD0
TERM3=1.0 / (18.0 * VISC(IPT0))
RECRH2=1.0/RHO2
BEXIST =.FALSE.
CALL SUB3L(BDYFX,.FALSE.,BDYFY,.FALSE.,BDYFZ,.FALSE.)
CALL SUB2(IASMP1,0, IASMP2,0)
DO 10010 I = 1,NUMPAT
C... Check if patch has SPEDAT(patch_name,ASLP,C,ASM) ie from VR object
CVAR=' '
CALL GETSDC(NAMPAT(I),'ASLP',CVAR)
IF(NAMPAT(I)(1:3).EQ.'ASM'.OR.CVAR.EQ.'ASM') THEN
BEXIST = .TRUE.
IF(IASMP1.EQ.0) IASMP1=I
IASMP2=I
ELSEIF(NAMPAT(I)(1:4).EQ.'BUOY' .OR.
1 NAMPAT(I)(1:4).EQ.'ROTA') THEN
CALL GETCOV(NAMPAT(I),U1,U1CO,U1VA)
CALL GETCOV(NAMPAT(I),V1,V1CO,V1VA)
CALL GETCOV(NAMPAT(I),W1,W1CO,W1VA)
IF(NEZ(U1VA)) BDYFX = .TRUE.
IF(NEZ(V1VA)) BDYFY = .TRUE.
IF(NEZ(W1VA)) BDYFZ = .TRUE.
IF(NAMPAT(I)(1:4).EQ.'BUOY') THEN
LBUOYX=GRNDNO(2,U1VA)
LBUOYY=GRNDNO(2,V1VA)
LBUOYZ=GRNDNO(2,W1VA)
ENDIF
ENDIF
10010 CONTINUE
IF(.NOT.LASLPB) THEN
BDYFX = SLIPU; BDYFY = SLIPV; BDYFZ = SLIPW
ENDIF
IF(BDYFX) CALL GXMAKE0(L0BDYX,NXNY,'GRVX')
IF(BDYFY) CALL GXMAKE0(L0BDYY,NXNY,'GRVY')
IF(BDYFZ) CALL GXMAKE0(L0BDYZ,NXNY,'GRVZ')
IF(BDYFZ) CALL GXMAKE0(L0BDFL,NXNY,'GRVL')
FEXIST(7) =.NOT.STEADY
C>>>>>>>>>>>>>>>>>>>>> Group 1 ------- Section 2 <<<<<<<<<<<<<<<<<<<<<
C Obtain F-array pointers needed for later calculations.
C.......................................................................
ELSEIF(IGR.EQ.1.AND.ISC.EQ.2) THEN
L0MIN0 = L0F(LMIN1 ); L0MOU0 = L0F(LMOUT1)
IF(BDYFX) CALL FN1(-L0BDYX,0.0)
IF(BDYFY) CALL FN1(-L0BDYY,0.0)
IF(BDYFZ) CALL FN1(-L0BDYZ,0.0)
IF(.NOT.NULLPR) THEN
CALL WRITBL
CALL WRIT3L('bdyfx ',bdyfx,',bdyfy ',bdyfy,
1 ',bdyfz ',bdyfz)
CALL WRIT2I('iasmp1 ',iasmp1,',iasmp2 ',iasmp2)
CALL WRITBL
ENDIF
IF(BEXIST) THEN
C... Allocate storage for inflow-outflow nett source calculations
CALL GXMAKE0(L0PATN,NXNY*NZ,'PATN') ! patch no for each cell
CALL GXMAKE0(L0CV4,IASMP2*NPT,'L0CV') ! net source location for patc
CALL GXMAKE0(L0COP,IASMP2*NPT,'L0CO') ! VAL for particle
CALL GXMAKE0(L0COV,IASMP2,'L0CV') ! VAL for VFOL
CALL GXMAKE0(L0INFPT,NXNY*NZ*NPT,'L0IN') ! inflow for each cell & part
CALL GXMAKE0(L0OUTPT,NXNY*NZ*NPT,'L0OU') ! outflow for each cell & par
ENDIF
C
ELSEIF(IGR.EQ.9.AND.ISC.EQ.1) THEN
C>>>>>>>>>>>>>>>>>>>>> Group 9 ------- Section 1 <<<<<<<<<<<<<<<<<<<<<
C Calculate the mixture density rhom if RHO1 = GRND from
C n
C rhom = (1 - sum PTi)*rhol + sum (PTi*rhoi) (8)
C i=1
C
C wherein the summation is over i from 1 to n, n being the
C number of groups of particles.
C.......................................................................
CALL FN1(DEN1,RHO2)
DO 91010 IPT = IPTF,IPTL
CALL FN34(DEN1,IPT,DENS(IPT)-RHO2)
91010 CONTINUE
C
ELSEIF(IGR.EQ.9.AND.ISC.EQ.6) THEN
C>>>>>>>>>>>>>>>>>>>>> Group 9 ------- Section 6 <<<<<<<<<<<<<<<<<<<<<
C Calculate the mixture viscosity nu if ENUL = GRND from
C
C n
C nu = (1 - sum PTi)*null + sum (PTi*nuli) (9)
C i=1
C
C wherein the summation is over i from 1 to n, n being the
C number of groups of particles.
C.......................................................................
CALL FN1(VISL,VISC(IPT0))
DO 96010 IPT = IPTF,IPTL
CALL FN34(VISL,IPT,VISC(IPT)-VISC(IPT0))
96010 CONTINUE
C
ELSEIF(IGR.EQ.19.AND.ISC.EQ.2) THEN
C>>>>>>>>>>>>>>>>>>>>> Group 19 ------- Section 2 <<<<<<<<<<<<<<<<<<<<<
C
C Initialise variables needed to calculate the nett inflow/outflow
C sources for each particle phase. This is done on the first sweep
C of each timestep to allow for the possibility that patches will
C be turned on and off with time.
C.......................................................................
IF(BEXIST.AND.ISWEEP.EQ.FSWEEP) THEN
LBPATNO=LBNAME('PATN') ! Optional 3D STORE(PATN) for patch numbers
IF(LBPATNO.NE.0) L0PATT=L0F(LBPATNO)
CALL ZERNUM(L0PATN,NXNY*NZ) ! set all patch numbers to zero
CALL ZERNUM(L0INFPT,NXNY*NZ*NPT) ! set all inflows to zero
CALL ZERNUM(L0OUTPT,NXNY*NZ*NPT) ! set all outflows to zero
C... Loop over all ASM patches,
IPATT=0
DO IPAT = IASMP1,IASMP2
CALL GETPAT(IPAT,IDUM,TY,I1,I2,J1,J2,K1,K2,L1,L2)
IF(ISTEP.GE.L1.AND.ISTEP.LE.L2) THEN
IPATT=IPATT+1
LVRCFD= ISVROB .AND. MASKPT(IPAT)
C... Loop over cells within patch, storing the patch number
DO IIZ=K1,K2
IF(LVRCFD) THEN
L0MSKZ=L0PVAR(22,IPAT,IIZ)
L0PAT=L0PATNO(IDMN)+(IIZ-1)*NXNY ! index for patch-number store
IPW=0; IBLK=0; L0FACZ=L0FACE
ENDIF
DO IIX=I1,I2
DO IIY=J1,J2
IJK=(IIX-1)*NY+IIY+(IIZ-1)*NXNY
IJKD=(IIX-1)*NY+IIY+(IIZ-1)*NFM
IF(LVRCFD) THEN
C.... For facetted objects, only store cell number if in active cell
IPW=IPW+1
C... Check for 2D facetted object
IF(L2DVRO(IPW)) THEN
F(L0PATN+IJK)=IPAT ! mark cell inside object
C... Check for 3D facetted object - ANGLED_IN/OUT. Check for finite surface area
ELSEIF(LINVRO(IPW).AND..NOT.SLD(IJK)) THEN
CALL GET_SURFACE(IIX,IIY,IIZ,IBLK,L0PAT,CAREA,
1 NORML,IPLUS)
IF(CAREA.GT.0.0) F(L0PATN+IJK)=IPAT ! mark cell with finite area
ENDIF
ELSE
C... otherwise mark all cells within bounding box
F(L0PATN+IJK)=IPAT
ENDIF
C... Optional STORE(PATN) for PHI/RESULT
IF(LBPATNO.NE.0) F(L0PATT+IJKD)=F(L0PATN+IJK)
ENDDO
ENDDO
ENDDO
C... Get VAL for VFOL and store
CALL GETCV(IPAT,LBVFOL,GCO,GVAL)
F(L0COV+IPAT)=GVAL
C... Now loop over all particle variables to store VAL for the particle
C... and the location of the nett source store
DO IPT=IPT0,IPTL
JPT = IPT - IPT0M1
CALL GETCV(IPAT,IPT,GCO,GVAL)
IF(IPT.EQ.IPT0.AND.IPATT.EQ.1) M0CV4=I0PHI+4
F(L0CV4+(IPAT-1)*NPT+JPT)=I0PHI+4-M0CV4
F(L0COP+(IPAT-1)*NPT+JPT)=GVAL
ENDDO
ENDIF ! end of IF(ISTEP
ENDDO ! end of patch loop
ENDIF
ELSEIF(IGR.EQ.19.AND.ISC.EQ. 6) THEN
C>>>>>>>>>>>>>>>>>>>>> Group 19 ------- Section 6 <<<<<<<<<<<<<<<<<<<<<
C
C Solve for cell volume fractions using inner cell iterations
C (linear relaxation RLX1, max number of iterations NIT1)
C and outer slab iterations (linear relaxation RLX2, max number
C of iterations NIT2).
C.......................................................................
IF(ISWEEP.EQ.1.OR.ISWEEP.EQ.LSWEEP.OR.(INDVAR/=9.AND.INDVAR/=1))
1 GO TO 25
C
RESREF(IPT0) = RESREF(P1)
RLX1 = ABS(DTFALS(IPT0))
RLX2 = ABS(DTFALS(IPT0+1))
C
IF(NX.GT.1) THEN
L0U1 = L0F(U1); L0DXG=L0F(DXG2D)
ENDIF
IF(NY.GT.1) THEN
L0V1 = L0F(V1); L0DYG=L0F(DYG2D)
ENDIF
IF(NZ.GT.1) THEN
L0W1 = L0F(W1); L0WL = L0F(LOW(W1))
L0R0L = L0F( LOW(IPT0)); L0R0H = L0F(HIGH(IPT0))
ENDIF
IF(.NOT.STEADY) L0VOL = L0F(VOL)
L0DEN = L0F(DEN1)
L0P1 = L0F(P1)
IF(SLIPU) L0AE = L0F(AEAST )
IF(SLIPV) L0AN = L0F(ANORTH)
IF(SLIPW) L0AH = L0F(AHIGH )
L0MIN = L0F(LMIN1 ); L0MOU = L0F(LMOUT1)
C
IF(IZ.EQ.1) F(L0TTRS+IABS(ISL(IPT0))) = 0.0
FEXIST(6) = IZ.GT. 1; FEXIST(5) = IZ.LT.NZ
c
c.... store beginning-of-sweep values for use in under-relaxation
c
L0OSPT=L0OSWP - NXNY
IF(LBPTEC>0) L0PTEC=L0F(LBPTEC)
IF(LBPTES>0) L0PTES=L0F(LBPTES)
IF(LBNIT1>0) L0NIT1=L0F(LBNIT1)
IF(LBNIT2>0) L0NIT2=L0F(LBNIT2)
DO IPT=IPT0,IPTL
L0IPT=L0F(IPT)
L0OSPT=L0OSPT + NXNY
DO I=1,NXNY
F(L0OSPT+I) = F(L0IPT+I)
ENDDO
C... Store indices for particles
JPT=IPT-IPT0M1
IF(IPT.EQ.IPT0) M0L0F=L0IPT
F(L0FPT+JPT)=L0IPT-M0L0F
IF(.NOT.STEADY) THEN
L0IPTO=L0F(OLD(IPT))
IF(IPT.EQ.IPT0) M0L0FO=L0IPTO
F(L0FPTO+JPT)=L0IPTO-M0L0FO
IF(IZ==1) THEN; TIMIN(IPT)=0.0; TIMOU(IPT)=0.0; ENDIF
ENDIF
ENDDO
C
C The force-and-acceleration loop begins here.
c
DO 907 IX=IXF,IXL
IADD=(IX-1)*NY
DO 907 IY=IYF,IYL
I = IADD+IY
c.... Compute the reciprocal of the absolute acceleration as the
c square root of the sum of the sguares of the pressure gradient
IF(SLD(I)) THEN
F(L0RACC+I) = 0.0
PGRDSQ=0.0
ELSE
PGRDSQ = TINY
C Calculate local pressure-gradient divided by density
c store on l0bdyx, etc
IF(BDYFX) THEN
IF(IX.LT.NX.OR.XCYCLE) THEN
IF(IX.EQ.NX) THEN
IE = IY
ELSE
IE = I + NY
ENDIF
IF(.NOT.SLD(IE)) THEN
F(L0BDYX+I) = ( F(L0P1+IE) - F(L0P1+I) )
1 /F(L0DXG+I)
IF(LBUOYX) F(L0BDYX+I) = F(L0BDYX+I) + BUOYA*BUOYD
F(L0BDYX+I) = 2.0*F(L0BDYX+I)
1 /(F(L0DEN+IE)+F(L0DEN+I))
PGRDSQ = PGRDSQ + F(L0BDYX+I)**2
ENDIF
ENDIF
ENDIF
IF(BDYFY) THEN
IF(IY.LT.NY) THEN
IN = I + 1
IF(.NOT.SLD(IN)) THEN
F(L0BDYY+I) = ( F(L0P1+IN) - F(L0P1+I) )
1 /F(L0DYG+I)
IF(LBUOYY) F(L0BDYY+I) = F(L0BDYY+I) + BUOYB*BUOYD
F(L0BDYY+I) = 2.0*F(L0BDYY+I)
1 /(F(L0DEN+IN)+F(L0DEN+I))
PGRDSQ = PGRDSQ + F(L0BDYY+I)**2
ENDIF
ENDIF
ENDIF
IF(BDYFZ) THEN
IF(IZ.LT.NZ) THEN
IF(.NOT.SLD(I+NXNY)) THEN
F(L0BDYZ+I) = ( F(L0P1+NFM+I) - F(L0P1+I) )
1 /DZG
IF(LBUOYZ) F(L0BDYZ+I) = F(L0BDYZ+I) + BUOYC*BUOYD
F(L0BDYZ+I) = 2.0*F(L0BDYZ+I)
1 /(F(L0DEN+NFM+I)+F(L0DEN+I))
PGRDSQ = PGRDSQ + F(L0BDYZ+I)**2
ENDIF
ENDIF
ENDIF
F(L0RACC+I) = SQRT (1.0/pgrdSQ)
ENDIF
907 CONTINUE
C
C.... The loops for computing particle concentrations start here, viz:
c 1. an outer loop controlling nit2 iterations,
c 2. a middle loop causing all cells in the slab to be visited in
c turn, and
c 3. an inner loop controlling iterations at a single cell.
C
RECDT=1.0/DT
DO 910 IT2 = 1, NIT2
SLBERR = 0.0
F(L0SLRS+IABS(ISL(IPT0))) = 0.0
DO 920 IX = IXF,IXL
IADD = (IX-1)*NY
FEXIST(4) = XCYCLE.OR.IX.GT. 1
FEXIST(3) = XCYCLE.OR.IX.LT.NX
C
C Initialise OLDSLB using values prior to the slab loop
C
DO 920 IY = IYF,IYL
I = IY+IADD
IF(SLD(I)) GO TO 920
DO 925 IPT=IPT0,IPTL
JPT = IPT - IPT0M1
L0PT=NINT(F(L0FPT+JPT))+M0L0F
F(L0OSLB+NPT*(I-1)+JPT) = F(L0PT+I)
925 CONTINUE
FEXIST(2) = IY.GT. 1
FEXIST(1) = IY.LT.NY
C Set neighbouring cell indices for dependent variables in a slab
IN = I+ 1
IS = I- 1
IW = I-NY
IF(IX.EQ. 1) IW = I+(NX-1)*NY
IE = I+NY
IF(IX.EQ.NX) IE = IY
C
C Put slip velocities in array. First initialise to zero.
C
c.... this initialization is necessary
CALL ZERNUM(L0SLPS,6*NPT)
DO 935 IPT=IPTF,IPTL
JPT = IPT - IPT0M1
DENRAT = (DENS(IPT) - RHO2) * RECRH2
DIAM1= DIAM(IPT)
DIAM2=DIAM1**2
c
c.... CRRACC = critical value of the reciprocal of the nagnitude of the
c acceleration for change of drag coefficient formula
c = diam**3 * denrat * CD0 / ( 432 * visc**2 )
c
CRRACC = ABS(DENRAT) * TERM1 * DIAM1 * DIAM2
c
c.... SLDPGx = slip velocity divided by pressure gradient
c x = P for the north, east and high faces;
c = S, W and L for the south, west and low faces
c
IF(F(L0RACC+I).GE.CRRACC) THEN
c.... low Re formula, with slip velocity / pressure gradient , ie SLDPG,
c = constant times density ratio times diameter squared / viscosity
c SLDPGC is common to north, east and high cells, so is calculated
c once only.
SLDPGC = DENRAT * DIAM2 * TERM3
ELSE
FAC1 = SIGN(1.0,DENRAT)
TERM = DENRAT * F(L0RACC+I) * TERM2 * DIAM1
SLDPGC = FAC1 * SQRT(ABS(TERM)*1.333333)
ENDIF
IF(BDYFY) THEN
IF(FEXIST(1)) THEN
F(L0SLPS+1) = 0.0
F(L0SLPS+JPT) = SLDPGC*F(L0BDYY+I)
ENDIF
IF(FEXIST(2)) THEN
IF(F(L0RACC+IS).GE.CRRACC) THEN
SLDPG = DENRAT * DIAM2 * TERM3
ELSE
FAC1 = SIGN(1.0,DENRAT)
TERM = DENRAT * F(L0RACC+IS) * TERM2 * DIAM1
SLDPG = FAC1 * SQRT(ABS(TERM)*1.333333)
ENDIF
L0SLP=L0SLP2
F(L0SLP+1) = 0.0
F(L0SLP+JPT) = SLDPG*F(L0BDYY+IS)
ENDIF
ENDIF
IF(BDYFX) THEN
IF(FEXIST(3)) THEN
L0SLP=L0SLP3
F(L0SLP+1) = 0.0
F(L0SLP+JPT) = SLDPGC*F(L0BDYX+I)
ENDIF
IF(FEXIST(4)) THEN
IF(F(L0RACC+IW).GE.CRRACC) THEN
SLDPG = DENRAT * DIAM2 * TERM3
ELSE
FAC1 = SIGN(1.0,DENRAT)
TERM = DENRAT * F(L0RACC+IW) * TERM2 * DIAM1
SLDPG = FAC1 * SQRT(ABS(TERM)*1.333333)
ENDIF
L0SLP=L0SLP4
F(L0SLP+1) = 0.0
F(L0SLP+JPT) = SLDPG*F(L0BDYX+IW)
ENDIF
ENDIF
IF(BDYFZ) THEN
IF(FEXIST(5)) THEN
L0SLP=L0SLP5
F(L0SLP+1) = 0.0
F(L0SLP+JPT) = SLDPGC*F(L0BDYZ+I)
ENDIF
IF(FEXIST(6)) THEN
IF(F(L0FACL+I).GE.CRRACC) THEN
SLDPG = DENRAT * DIAM2 * TERM3
ELSE
FAC1 = SIGN(1.0,DENRAT)
TERM = DENRAT * F(L0FACL+I) * TERM2 * DIAM1
SLDPG = FAC1 * SQRT(ABS(TERM)*1.333333)
ENDIF
L0SLP=L0SLP6
F(L0SLP+1) = 0.0
F(L0SLP+JPT) = SLDPG*F(L0BDFL+I)
ENDIF
ENDIF
935 CONTINUE
C
C Start cell iteration
C
DO 930 IT1 = 1,NIT1
F(IPTSUM+I) = 0.0
C
C Initialise OLDCEL to value at start of current cell iteration.
C
DO 933 IPT=IPT0,IPTL
JPT = IPT - IPT0M1
L0PT=NINT(F(L0FPT+JPT))+M0L0F
F(L0OCEL+JPT) = F(L0PT+I)
933 CONTINUE
C
C Calculate continuous phase velocity for each cell face in turn
C
C N-face
IF(FEXIST(1)) THEN
IF(BDYFY) THEN
CALL CONVEL(L0SLP1,I,IN,F(L0V1+I),V0)
F(L0V0Y+I) = V0
ELSE
F(L0V0Y+I) = F(L0V1+I)
ENDIF
ENDIF
C S-face
IF(FEXIST(2)) THEN
IF(BDYFY) THEN
CALL CONVEL(L0SLP2,IS,I,F(L0V1+IS),V0)
F(L0V0Y+IS) = V0
ELSE
F(L0V0Y+IS) = F(L0V1+IS)
ENDIF
ENDIF
C E-face
IF(FEXIST(3)) THEN
IF(BDYFX) THEN
CALL CONVEL(L0SLP3,I,IE,F(L0U1+I),V0)
F(L0V0X+I) = V0
ELSE
F(L0V0X+I) = F(L0U1+I)
ENDIF
ENDIF
C W-face
IF(FEXIST(4)) THEN
IF(BDYFX) THEN
CALL CONVEL(L0SLP4,IW,I,F(L0U1+IW),V0)
F(L0V0X+IW) = V0
ELSE
F(L0V0X+IW) = F(L0U1+IW)
ENDIF
ENDIF
C H-face
IF(FEXIST(5)) THEN
IF(BDYFZ) THEN
CALL CONVEL(L0SLP5,I,I+NFM,F(L0W1+I),V0)
F(L0V0Z+I) = V0
ELSE
F(L0V0Z+I) = F(L0W1+I)
ENDIF
ENDIF
C L-face
IF(FEXIST(6)) THEN
IF(BDYFZ) THEN
CALL CONVEL(L0SLP6,I-NFM,I,F(L0WL+I),V0)
F(L0V0ZL+I) = V0
ELSE
F(L0V0ZL+I) = F(L0WL+I)
ENDIF
ENDIF
C
C Calculate imbalance for each component (and its derivative)
C by summing contributions from each cell face, allowing for
C sources/sinks and transience.
C
DBDRM = 0.0
BALTOT = 0.0
DO 940 IPT = IPTL,IPT0,-1
JPT=IPT-IPT0M1
L0SLPT=L0SLPS+JPT
BALANI = 0.0
DBALDR = 0.0
L0RIN = NINT(F(L0FPT+JPT))+M0L0F
RIP = F(L0RIN+I)
C N-face
IF(FEXIST(1)) THEN
ARN = F(L0AN+I)
VPLS= F(L0V0Y+I) + F(L0SLPT)
BALANI = BALANI + ARN*FACBAL(RIP,F(L0RIN+IN),-VPLS)
DBALDR = DBALDR + DFBLDR(-VPLS,ARN)
ENDIF
C S-face
IF(FEXIST(2)) THEN
IF(IS.GT.0) THEN
ARS = F(L0AN +IS)
ELSE
ARS = F(L0AN +1) ! strictly, i2dsb should be used
ENDIF
VPLS= F(L0V0Y+IS) + F(L0SLPT + NPT)
BALANI = BALANI + ARS*FACBAL(RIP,F(L0RIN+IS),VPLS)
DBALDR = DBALDR + DFBLDR(VPLS,ARS)
ENDIF
C E-face
IF(FEXIST(3)) THEN
ARE = F(L0AE +I)
VPLS= F(L0V0X+I) + F(L0SLPT+2*NPT)
BALANI = BALANI + ARE*FACBAL(RIP,F(L0RIN+IE),-VPLS)
DBALDR = DBALDR + DFBLDR(-VPLS,ARE)
ENDIF
C W-face
IF(FEXIST(4)) THEN
IF(IW.GT.0) THEN
ARW = F(L0AE+ IW)
ELSE
ARW = F(L0AE+ I) ! strictly, i2dwb should be used
ENDIF
VPLS= F(L0V0X+IW) + F(L0SLPT+3*NPT)
BALANI = BALANI + ARW*FACBAL(RIP,F(L0RIN+IW),VPLS)
DBALDR = DBALDR + DFBLDR(VPLS,ARW)
ENDIF
c
IF(NZ.GT.1) THEN
C H-face
IF(FEXIST(5)) THEN
L0RIH = NINT(F(L0FPT+JPT))+NFM+M0L0F
ARH = F(L0AH +I)
VPLS= F(L0V0Z+I) + F(L0SLPT+4*NPT)
BALANI = BALANI + ARH*FACBAL(RIP,F(L0RIH+I),
1 -VPLS)
DBALDR = DBALDR + DFBLDR(-VPLS,ARH)
ENDIF
C L-face
IF(FEXIST(6)) THEN
L0RIL = NINT(F(L0FPT+JPT))-NFM+M0L0F
ARL = F(L0AL +I)
VPLS= F(L0V0ZL+I) + F(L0SLPT+5*NPT)
BALANI = BALANI + ARL*FACBAL(RIP,F(L0RIL+I),VPLS)
DBALDR = DBALDR + DFBLDR(VPLS,ARL)
ENDIF
ENDIF
C Inflows
IF(BEXIST) THEN
IF(F(L0MIN+I).GT.0.0) THEN
IF(IT1.EQ.1) THEN
c.... note that it is implied that f(l0min+i) is the volumetric inflow
c to the cell; and that, below, f(l0mou+i) is the volumetric outflow
C... Get patch number for this cell
IJK=(IX-1)*NY+IY+(IZ-1)*NXNY
IPAT=NINT(F(L0PATN+IJK))
IF(IPAT.GT.0) THEN
C... recover VAL for particle and VFOL
RIIN=F(L0COP+(IPAT-1)*NPT+JPT)
RV=F(L0COV+IPAT)
FLOWIN = RIIN*F(L0MIN+I)*RV
BALANI = BALANI + FLOWIN
F(L0INFL+JPT)=FLOWIN
C... store mass inflow to this cell
IJKP=IJK+(JPT-1)*NXNY*NZ
F(L0INFPT+IJKP)=RIIN*F(L0MIN+I)
ENDIF
ELSE
BALANI = BALANI + F(L0INFL+JPT)
ENDIF
ENDIF
C Outflows
IF(F(L0MOU+I).GT.0.0) THEN
BALANI = BALANI - F(L0MOU+I) * RIP
DBALDR = DBALDR - F(L0MOU+I)
C... store mass outflow for this cell
IJK=(IX-1)*NY+IY+(IZ-1)*NXNY
IPAT=NINT(F(L0PATN+IJK))
IF(IPAT.GT.0) THEN
IJKP=IJK+(JPT-1)*NXNY*NZ
F(L0OUTPT+IJKP)=F(L0MOU+I)*RIP*F(L0DEN+I)
ENDIF
ENDIF
ENDIF
C T-face
IF(FEXIST(7)) THEN
L0RIO = NINT(F(L0FPTO+JPT))+M0L0FO
TIMFLO = F(L0VOL+I) * RECDT
BALANI = BALANI + ( F(L0RIO+I) - RIP ) * TIMFLO
DBALDR = DBALDR - TIMFLO
ENDIF
Change
F(L0CHNG+IPT) = BALANI
IF(LASLPA) THEN
DBDRM = DBDRM + ABS(BALANI)*DBALDR
ELSE
DBDRM = MIN(DBALDR,DBDRM)
ENDIF
BALTOT = BALTOT + ABS(BALANI)
940 CONTINUE
ERROR1 = 0.0
IF(LASLPA) DBDRM = DBDRM/(BALTOT+TINY)
C... The following particle-loop uses Newton-Raphson to solve point
C for the volume fractions
F(IPTSUM+I) = 0.0
DO 950 IPT = IPT0,IPTL
JPT=IPT-IPT0M1
L0RIN = NINT(F(L0FPT+JPT))+M0L0F
RNEW = F(L0RIN+I) - RLX1*F(L0CHNG+IPT) / (DBDRM-TINY)
IF(RNEW.LT.0.0) THEN
F(L0RIN+I) = 0.0
GO TO 950
ELSEIF(RNEW.GT.1.0) THEN
F(L0RIN+I) = 1.0
ELSE
F(L0RIN+I) = AMAX1(VARMIN(IPT),AMIN1(RNEW,VARMAX(IPT)))
ENDIF
F(IPTSUM+I) = F(IPTSUM+I) + F(L0RIN+I)
950 CONTINUE
C... The following particle loop seems to adjust volume fractions so they add to 1.0
DO 951 IPT = IPT0,IPTL
JPT = IPT - IPT0M1
L0RIN = NINT(F(L0FPT+JPT))+M0L0F
F(L0RIN+I) = F(L0RIN+I) / F(IPTSUM+I)
ERROR1 = ERROR1 + ABS(F(L0RIN+I)-F(L0OCEL+JPT))
951 CONTINUE
cc call writ3r('error1 ',error1,'crt1 ',crt1,'diff ',
cc 1 error1-crt1)
IF(ERROR1.LT.CRT1) GO TO 960 ! skip out of cell loop
930 CONTINUE ! End of CELL iteration loop
960 CONTINUE
IF(LBPTEC>0) F(L0PTEC+I)=ERROR1
IF(LBNIT1>0) F(L0NIT1+I)=MIN(IT1,NIT1)
F(L0SLRS+IABS(ISL(IPT0))) = SLBRES(IPT0) + BALTOT
DO 965 IPT=IPT0,IPTL
JPT = IPT - IPT0M1
L0RIN = NINT(F(L0FPT+JPT))+M0L0F
SLBERR=SLBERR+ABS(F(L0RIN+I)-F(L0OSLB+NPT*(I-1)+JPT))
965 CONTINUE
920 CONTINUE ! end of IX and IY loops
IF(SLBERR.LT.CRT2) GO TO 970 ! skip out of slab loop
910 CONTINUE ! end of SLAB iteration loop
970 CONTINUE
IF(LBPTES>0.OR.LBNIT2>0) THEN
DO II=1,NXNY
IF(LBPTES>0) F(L0PTES+II)=SLBERR
IF(LBNIT2>0) F(L0NIT2+II)=MIN(IT2,NIT2)
ENDDO
ENDIF
IF(FEXIST(7)) THEN
DO IPT=IPT0,IPTL
JPT=IPT-IPT0M1
L0RIN = NINT(F(L0FPT+JPT))+M0L0F
L0RIO = NINT(F(L0FPTO+JPT))+M0L0FO
! write(14,'(''iz, ipt'',2i3)') iz,ipt
DO I=1,NXNY
IF(SLD(I)) CYCLE
RIP = F(L0RIN+I); RIPO = F(L0RIO+I)
TIMFLO = F(L0VOL+I) * RECDT
! write(14,'(''i='',i3,'' iz='',i3,'' ipt'',i3)')i,iz,ipt
IF(IPT>IPT0) THEN
PDEN=DENS(IPT)
ELSE
PDEN=RHO2
ENDIF
! write(14,'(''pden, rio,rip '',1p3e13.5)') pden,ripo,rip
TIMFLIN=RIPO* TIMFLO*PDEN
TIMFLOU=RIP* TIMFLO*PDEN
TIMIN(IPT)=TIMIN(IPT)+TIMFLIN
TIMOU(IPT)=TIMOU(IPT)-TIMFLOU
ENDDO
! write(14,'(''timin='',1pe13.5,'' timou='',1pe13.5)')
! 1 timin(ipt),timou(ipt)
ENDDO
ENDIF
C
C... loop over all ASM patches to sum up inflows and outflows into nett source
IF(IZ.EQ.NZ) THEN
DO IPAT = IASMP1,IASMP2
C... Check if patch has SPEDAT(patch_name,ASLP,C,ASM)
CVAR=' '
CALL GETSDC(NAMPAT(IPAT),'ASLP',CVAR)
IF(NAMPAT(IPAT)(1:3).EQ.'ASM'.OR.CVAR.EQ.'ASM') THEN
C... Get patch limits
CALL GETPAT(IPAT,IDUM,TY,I1,I2,J1,J2,K1,K2,L1,L2)
IF(ISTEP.GE.L1.AND.ISTEP.LE.L2) THEN
C... Loop over all particles
DO IPT=IPT0,IPTL
JPT = IPT - IPT0M1
C... recover index for nett source, and zero it
I0PHI4=NINT(F(L0CV4+(IPAT-1)*NPT+JPT))+M0CV4
F(I0PHI4)=0.0
C... loop over all cells in patch summing inflow & outflow
DO IIZ=K1,K2
DO IIX=I1,I2
DO IIY=J1,J2
IJK=(IIX-1)*NY+IIY+(IIZ-1)*NXNY
IJKP=IJK+(JPT-1)*NXNY*NZ
IF(NINT(F(L0PATN+IJK)).EQ.IPAT) THEN ! cell has been marked
FLOIN = F(L0INFPT+IJKP)
FLOOUT= F(L0OUTPT+IJKP)
F(I0PHI4)=F(I0PHI4)+FLOIN-FLOOUT
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO ! end of particle loop
ENDIF
ENDIF
ENDDO ! end of patch loop
ENDIF
C
F(L0TTRS+IABS(ISL(IPT0))) = TOTRES(IPT0) + SLBRES(IPT0)/
1 RESREF(IPT0)
C
C RLX2 is a relaxation factor after slab solution
C
L0OSPT=L0OSWP - NXNY
DO IPT=IPT0,IPTL
JPT=IPT-IPT0M1
L0IPT=NINT(F(L0FPT+JPT))+M0L0F
L0OSPT=L0OSPT + NXNY
DO I=1,NXNY
F(L0IPT+I)=F(L0OSPT+I) + RLX2*( F(L0IPT+I)-F(L0OSPT+I))
ENDDO
ENDDO
C
Come here to set the H-face body force of the current
C slab to the L-face body force of the next slab.
C Also define L0FACL for the next slab
C
IF(NZ.GT.1) THEN
CALL FN0(-L0AL,-L0AH)
CALL FN0(-L0FACL,-L0RACC)
IF(BDYFZ) CALL FN0(-L0BDFL,-L0BDYZ)
ENDIF
ENDIF
25 CONTINUE
cc call accutm('gxaslp',2)
cc if(isweep.eq.lsweep.and.igr.eq.19.and.isc.eq.6.and.iz.eq.nz) then
cc write(14,*)'varbl,isw,igr,isc,iz ',
cc 1 name(indvar),isweep,igr,isc,iz
cc call accutm('gxaslp',3)
cc endif
NAMSUB='gxaslp'
END
C***********************************************************************
FUNCTION FACBAL(RIP,RIN,VPLS)
IF(VPLS.GT.0.0) THEN
FACBAL = RIN*VPLS
ELSE
FACBAL = RIP*VPLS
ENDIF
END
C***********************************************************************
FUNCTION DFBLDR(VPLS,A)
DFBLDR = -A*MAX(0.0,-VPLS)
END
C***********************************************************************
c
SUBROUTINE CONVEL(L0SLP,I1,I2,BAL0,V0)
INCLUDE 'farray'
INCLUDE 'grdloc'
CHARACTER*4 CSG1,CSG2,CSG3
COMMON/CONVL/IPT0,IPTL,NPT,L0NORD
COMMON/GENI/IGNI1(9),NFM,IGNI11(50)
c
C The following sequence puts the slip velocities in decreasing order
F(L0NORD+1)=1
DO 10 I=2,NPT
SLP=F(L0SLP+I)
DO 20 J=1,I-1
NORDJ=NINT(F(L0NORD+J))
IF(SLP.GT.F(L0SLP+NORDJ)) THEN
DO 30 K=I,J+1,-1
30 F(L0NORD+K)=NINT(F(L0NORD+K-1))
F(L0NORD+J) = I
GO TO 10
ENDIF
IF(J.EQ.I-1) F(L0NORD+I) = I
20 CONTINUE
10 CONTINUE
c
BALOLD = 0.0
DO 100 JPT=1,NPT
BAL = 0.0
NORDJ=NINT(F(L0NORD+JPT))
V0 = - F(L0SLP + NORDJ)
DO 200 IPT=IPT0,IPTL
L0RPT=L0F(IPT)
JJ = IPT - IPT0 + 1
VPLS = V0 + F(L0SLP+JJ)
IF(VPLS.GT.0.0) THEN
BAL = BAL + F(L0RPT+I1) * VPLS
ELSE
BAL = BAL + F(L0RPT+I2) * VPLS
ENDIF
200 CONTINUE
IF(JPT.EQ.1) THEN
IF(BAL.GE.BAL0) GO TO 101
ELSEIF(JPT.EQ.NPT) THEN
IF(BAL.LE.BAL0) GO TO 101
ENDIF
IF(BAL.GE.BAL0) THEN
IF(BALOLD.LT.BAL0) GO TO 102
ENDIF
BALOLD = BAL
V0OLD = V0
GO TO 100
101 V0 = V0 - (BAL-BAL0)
RETURN
102 V0 = V0OLD + (V0-V0OLD) * (BAL0-BALOLD) / (BAL-BALOLD)
RETURN
100 CONTINUE
CALL WRIT40('failed continuous-phase velocity ')
CALL WRIT40('calculation in subroutine gxaslp ')
CALL SET_ERR(242,'Error. See result file.',2)
C CALL EAROUT(2)
END
c