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