cgxaslp.for

c
cFILENAME GXASLP.FOR                       120711
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
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 '/phoenics/d_includ/patnos'
      INCLUDE '/phoenics/d_includ/patcmn'
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/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
      LOGICAL LBUOYX,LBUOYY,LBUOYZ
      LOGICAL FEXIST(7),BEXIST,SLD
      LOGICAL LVRCFD,MASKPT,L2DVRO,LINVRO,dbnow
      INTEGER NORML(3)
      CHARACTER*6 NAMFUN,NAMSUB,CVAR*68
c
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,FEXIST,BEXIST,
     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
      SAVE LBUOYX,LBUOYY,LBUOYZ
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,
        CALL SUB2(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 12/07/11 has been called         ')
          CALL WRIT40('ASM of 12/07/11 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')
C
Come  here to insert the values of the particle densities
C     and diameters.
C
        IF(.NOT.NULLPR) CALL WRITBL
        DO 10005 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
          IF(IPT.EQ.IPT0) GO TO 10005
          IF(.NOT.NULLPR) THEN
            CALL WRIT1I('prtcl no',ipt-ipt0)
            CALL WRIT3R('density ',dens(ipt),'diameter',diam(ipt),
     1                  'viscsity',visc(ipt))
          ENDIF
10005   CONTINUE
        IF(.NOT.NULLPR) CALL WRITBL
C
        NIT1 = LITER(IPT0)
        NIT2 = LITER(IPT0+1)
        CRT1 = ENDIT(IPT0)
        CRT2 = ENDIT(IPT0+1)
        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,NPT,'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,T,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) 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
        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
          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)
                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) = RNEW
                  ENDIF
                  F(IPTSUM+I) = F(IPTSUM+I) + F(L0RIN+I)
950             CONTINUE
                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
930           CONTINUE
960           CONTINUE
              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
            IF(SLBERR.LT.CRT2) GO TO 970
910       CONTINUE
970     CONTINUE
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,T,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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/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