c C name estrgr.for................................... 181007 c

C***************************************************************
C
C  This GROUND subroutine contains the FORTRAN sequences needed
C  to convert EARTH into a simulator of flows in aluminium-
C  reduction cells of the Hall- & Soderberg-type. This GROUND is
C  selected in the input file by the statement NAMGRD=ESTR.
C  The input file contains notes on the process simulated, &
C  gives additional references.
C
C    Group 11 contains sequences for the initialization of the
C  three components of the magnetic field, the metal-electrolyte
C  interface height, the slopes of the anode undersides and the
C  Lorentz forces. Group 13 provides for the wall friction,
C  for the calculation of the Lorentz forces from the prescribed
C  magnetic field and the electric-current density, for the
C  representation of electomagnetic induction, for the sources of
C  under-anode gas evolution, for the current distribution at the
C  cathode, etc. The metal-electrolyte interface height is
C  computed in group 19, along with the current density deduced
C  from the product of the gradient of electric potential and
C  the electrical conductivity of the medium (ie. the
C  electrolyte or the metal according to z location).
C
C  This subroutine contains three examples of user-coded
C  functions, namely FNAVXY , FNSUM1 & FNSUM2
C
      SUBROUTINE ESTRGR
      INCLUDE '/phoenics/d_includ/objnam'
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
C
C... Include some EARTH commons to get no. of variables and
C    storage information
      COMMON /GENL/ LDUM1(20), INCORJ,INCORM, LDUM2(38)
      LOGICAL LDUM1,LDUM2,INCORJ,INCORM
      COMMON /NAMFN/ NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB,line*68
      CHARACTER FNAME*20
C
      COMMON/GENI/  NXNY,IGF1(8),NFM,IF2(5),NM,IF3(15),IPRL,IBTAU,ILTLS,
     1              IGFIL(15),ITEM1,
     1              ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,
     1              IRADZ,IVFOL
      COMMON/DRHODP/ITEMP,IDEN/DVMOD/IDVCGR
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION STARTS:
C
C 1   Set dimensions of data-for-GROUND arrays here. WARNING: the
C     corresponding arrays in the MAIN program of the satellite
C     and EARTH must have the same dimensions.
      COMMON/LGRND/LG(20)/IGRND/IG(20)/RGRND/RG(100)/CGRND/CG(10)
      LOGICAL LG
      CHARACTER*4 CG
C
C 2   User dimensions own arrays here, for example:
C     DIMENSION UUH(10,10),UUC(10,10),UUX(10,10),UUZ(10)
C
      INTEGER MAGF,HI,CPOR,EPOT,HANO,FX,FY,FZ,BX,BY,BZ,HUMP
      LOGICAL GDBG,SETJ,SLD,LSOLID,VAC,LVAC,POR,LPOR,GRN
      LOGICAL QEQ,QNE,NEZ
      REAL JZCATH
C
C--- Set dimensions of user-declared arrays
      PARAMETER (NYM=50, NXM=100, NXYM=NYM * NXM, NZM=25)
C
C----Following array dimensions are NYM,NXM
C--- except  GZWNZ(NZM)
      DIMENSION GYV(NYM,NXM),GBX(NYM,NXM),GBY(NYM,NXM),GBZ(NYM,NXM),
     1          GVAL(NYM,NXM),GJCATH(NYM,NXM),GZWNZ(NZM)
C
      DIMENSION IANODE(NYM,NXM),HANODE(NYM,NXM)
      DIMENSION IFREEZ(NYM,NXM)
C
C--- Equivalence elements of IG and RG arrays to their Q1 names
      EQUIVALENCE (IZ1,IG(1)),(IZ2,IG(2)),(NIH,IG(3)),(IHF,IG(4)),
     1 (MAGF,IG(5)),(IZ3,IG(6)),(IFRZ,IG(9)),(IANO,IG(10)),
     1 (HANO,IG(11)),(HUMP,IG(12)),(IZ0,IG(13)),(IMAGU,IG(14)),
     1 (IMAGF,IG(15))
      EQUIVALENCE (RHOMET,RG(1)), (RHOELC,RG(2)), (RHOANO,RG(3)),
     1            (CONMET,RG(4)), (CONELC,RG(5)), (CONANO,RG(6)),
     1            (AGRAVZ,RG(7)), (GMDOT,RG(8)),  (JZCATH,RG(9)),
     1            (SLOPE,RG(10)), (SLOH,RG(11)),  (FDH,RG(12)),
     1            (FHLIM,RG(13)), (FRCON,RG(14)), (RELF,RG(15)),
     1            (GDZ1,RG(16)),  (GDZ2,RG(17)),  (ENUMET,RG(18)),
     1            (ENUELC,RG(19)),(RHOFRZ,RG(20)),(CONFRZ,RG(21)),
     1            (ANOPOT,RG(22)),(RHOAIR,RG(23)),(CONAIR,RG(24)),
     1            (GPCOEF,RG(25)),(CONST1,RG(26)),(CONST2,RG(27)),
     1            (CONST3,RG(28)),(CONFAC,RG(29)),(BEMF,RG(30)),
     1            (CONCAT,RG(31)),(RHOCAT,RG(32))
      EQUIVALENCE (AA0,RG(40)),(AA1,RG(41)),(AA2,RG(42)),(AA3,RG(43)),
     1            (AA4,RG(44)),(AA5,RG(45)),(AA6,RG(46)),(AA7,RG(47)),
     1            (BB0,RG(50)),(BB1,RG(51)),(BB2,RG(52)),(BB3,RG(53)),
     1            (BB4,RG(54)),(BB5,RG(55)),(BB6,RG(56)),(BB7,RG(57)),
     1            (CC0,RG(60)),(CC1,RG(61)),(CC2,RG(62)),(CC3,RG(63)),
     1            (CC4,RG(64)),(CC5,RG(65)),(CC6,RG(66)),(CC7,RG(67))
C
      SAVE
C
C 3   User places his data statements here, for example:
C     DATA NXDIM,NYDIM/10,10/
C
C--- Initialise working arrays to 0.0 Length of each array is NXM * NYM
      DATA IANODE,IFREEZ  /NXYM*0  , NXYM*0/
      DATA IV194 /0/
      DATA IREF,IREF2 /0, 0/
C
C
C 4   Insert own coding below as desired, guided by GREX examples.
C     Note that the satellite-to-GREX special data in the labelled
C     COMMONs /RSG/, /ISG/, /LSG/ and /CSG/ can be included and
C     used below but the user must check GREX for any conflicting
C     uses. The same comment applies to the EARTH-spare working
C     arrays EASP1, EASP2,....EASP10. If the call to GREX has been
C     deactivated then they can all be used without reservation.
C
      IXL=IABS(IXL)
      IF(IGR.EQ.13) GO TO 13
      IF(IGR.EQ.19) GO TO 19
      GO TO (1,2,2,2,5,2,2,8,9,2,11,2,13,2,2,2,2,2,19,2,2,
     1 2,2,2),IGR
    2 CONTINUE
      RETURN
C*****************************************************************
C
C--- GROUP 1. Run title
C
    1 GO TO (1001,1002,1000),ISC
 1000 IF(.NOT.ONEPHS.AND.MOD(ISLN(R1),13).NE.0) THEN
        ISLN(R1)=ISLN(R1)*13
        ISLN(R2)=ISLN(R1)
      ENDIF
      RETURN
 1001 CONTINUE
      NAMSUB='G1.1  '
C
      CALL WRYT40('PHOENICS version number is   :   2007   ')
      CALL WRYT40('GROUND file is ESTRGR.HTM of:    181007 ')
      CALL WRIT40('PHOENICS version number is   :   2007   ')
      CALL WRIT40('GROUND file is ESTRGR.HTM of:    181007 ')
C
C--- Check dimensioning of local arrays
C
      IF(NX.GT.NXM.OR.NY.GT.NYM.OR.NZ.GT.NZM) THEN
       CALL WRIT40('Local arrays in GROUND are too small!!  ')
       CALL WRIT40('Please reset NXM NYM & NZM in GROUND,   ')
       CALL WRIT40('CONDCT & POTENT, recompile and relink   ')
      CALL SET_ERR(312,
     *   'Error. Local arrays in GROUND are too small.',1)
C       CALL EAROUT(1)
      ENDIF
C
C--- Allocate additional EARTH storage
      IF(LG(9)) THEN
C--- Stores for vertical Lorentz force contribution to interface
C    balance
        CALL GXMAKE(L0PDIF,NX*NY,'PDIF')
        IF(LG(3)) CALL GXMAKE(L0FZDIF,NX*NY,'FZDIF')
      ENDIF
C--- Stores for conductivity
      CALL GXMAKE(L0CND,NX*NY,'COND')
      CALL GXMAKE(L0CNDH,NX*NY,'CNDH')
      CALL GXMAKE(L0CONDE,NX*NY*NZ,'CONDE')
      CALL GXMAKE(L0CONDN,NX*NY*NZ,'CONDN')
      CALL GXMAKE(L0CONDH,NX*NY*NZ,'CONDH')
 
C--- GRSP2 stores interface height
      CALL MAKE(GRSP2)
C--- GRSP3 stores height of anode underside
      CALL MAKE(GRSP3)
C--- GRSP4 stores free surface height
      CALL MAKE(GRSP4)
C
      NXM1=NX-1
      NYM1=NY-1
      AXDZ=GRND
C--- Default location for air surface
      IF(IZ3.EQ.0) IZ3=NZ
C--- Find storage locations for STOREd variables
      EPOT=INAME('EPOT')
      CPOR=INAME('CPOR')
      HI  =INAME('HI  ')
      JX  =INAME('JX  ')
      JY  =INAME('JY  ')
      JZ  =INAME('JZ  ')
      JIX =INAME('JIX ')
      JIY =INAME('JIY ')
      JIZ =INAME('JIZ ')
      BX  =INAME('BX  ')
      BY  =INAME('BY  ')
      BZ  =INAME('BZ  ')
      FX  =INAME('FX  ')
      FY  =INAME('FY  ')
      FZ  =INAME('FZ  ')
      IPRPS=INAME('PRPS')
C
C--- Set atmospheric conditions for free surface adjustment
      PAIR=0.0
C
C--- Initialise JZ and set cathode current density for use in Group 13
C--- MUST be done here, as Group 11 not visited on RESTRT!!!!
      IF(JZCATH.LE.0.0) THEN
C--- Set cathode current to JZCATH directly
        DO 1004 IX=1,NX
        DO 1004 IY=1,NY
 1004     GJCATH(IY,IX)=JZCATH
      ELSE
C--- Read cathode current distribution from disc file
        LUNITR=IFIX(JZCATH)
        FNAME=CG(2)
 1005   CONTINUE
        CALL DSCFLS(LUNITR,FNAME,12,0)
        IF(FNAME(1:2).EQ.'q1') THEN
          DO WHILE (LINE.NE.'  *JZC')
            READ(LUNITR,'(A)',ERR=1010,END=1010) LINE
          ENDDO
        ENDIF
        DO IY=NY,1,-1
          READ(LUNITR,*,ERR=1010,END=1010)(GJCATH(IY,IX),IX=1,NX)
        ENDDO
        CALL DSCFLS(LUNITR,FNAME,14,0)
        WRITE(LUPR1,'(/,'' CATHODE CURRENTS READ FROM FILE '',A/)')
     1                                                             FNAME
        GO TO 1015
 1010   CONTINUE
        IF(FNAME.EQ.'q1') THEN
          CALL DSCFLS(LUNITR,FNAME,14,0)
          FNAME='q1ear'
          GO TO 1005
        ENDIF
        WRITE(LUPR1,'(/,'' ERROR IN READING CURRENTS FILE '',A)') FNAME
      CALL SET_ERR(313, 'Error in reading file',1)
C        CALL EAROUT(1)
 1015   CONTINUE
      ENDIF
C
C--- Ensure that UDIFF is set when BEMF required
      IF(NEZ(BEMF)) UDIFF=.TRUE.
C--- Ensure specific heat is set
      IF(.NOT.GRN(CP1)) CP1=GRND
      RETURN
C
 1002 CONTINUE
      NAMSUB='G1.2  '
C
C... Check if dependent variables stored incore. If on disc, then
C    indexing functions fail, so stop.
      IF(.NOT.INCORM) THEN
       CALL WRIT40('Dependent variables stored on disc.     ')
       CALL WRIT40('Please increase F array as shown above, ')
       CALL WRIT40('then recompile and relink               ')
       CALL SET_ERR(314, 'Error. See result file.',1)
C       CALL EAROUT(1)
      ENDIF
C... Same applies to old-time variables for transient.
      IF(.NOT.STEADY .AND. .NOT.INCORJ) THEN
       CALL WRIT40('Old-time variables stored on disc.      ')
       CALL WRIT40('Please increase F array as shown above, ')
       CALL WRIT40('then recompile and relink               ')
       CALL SET_ERR(315, 'Error. See result file.',1)
C       CALL EAROUT(1)
      ENDIF
      RETURN
C*****************************************************************
C
C--- GROUP 5. Z-direction grid specification
C
    5 CONTINUE
C--- Recompute LXYDZ,LXYDZG,LXYDZL
      IF(IZ.EQ.1) THEN
       CALL FN0(LXYDZ,HI)
       CALL FN2(LXYDZG,HIGH(HI),0.0,0.5)
       CALL FN1(LXYDZL,0.0)
      ELSE
       CALL FN10(LXYDZ,HI,LOW(HI),0.0,1.0,-1.0)
       CALL FN0(LXYDZL,LXYDZG)
       IF(IZ.LT.NZ) CALL FN10(LXYDZG,HIGH(HI),LOW(HI),0.0,0.5,-0.5)
      ENDIF
C
C... Ensure grid size does not go to zero
      CALL FN22(LXYDZG,1.E-10)
      CALL FN22(LXYDZL,1.E-10)
      CALL FN22(LXYDZ ,1.E-10)
C
      RETURN
C*****************************************************************
C
C--- GROUP 8. Terms (in differential equations) & devices
C
    8 IF(ISC.EQ.5.OR.ISC.EQ.6) THEN
C   * -----------GROUP 8  SECTION 5 & 6 --- ADDITIONAL VELOCITY
C--- for W1AD.LE.GRND--- phase 1 additional velocity (VELAD).
C--- for W2AD.LE.GRND--- phase 2 additional velocity (VELAD).
C
C--- Calculate local grid velocity : -Wg = (Hold - Hnew) / DT
        NAMSUB='G85&6 '
        IF(.NOT.STEADY) THEN
          CALL FN10(VELAD,OLD(HI),HI,0.0,1./DT,-1./DT)
          IF(GDBG.AND.LG(15)) THEN
            WRITE(LUPR1,*) 'GROUP 8.5 AT IZ = ',IZ,' ISWEEP = ',ISWEEP
            WRITE(LUPR1,*) 'DT = ',DT
            CALL PRN('HI  ',HI)
            CALL PRN('HIOL',OLD(HI))
            CALL PRN('W1AD',VELAD)
          ENDIF
        ENDIF
      ELSE IF(ISC.EQ.9) THEN
C   * -----------GROUP 8  SECTION 9 --- DIFFUSION COEFFICIENTS
C.... Entered when UDIFF =.TRUE.; block-location indices are:
C     LAE for east, LAN for north, and LD11 for high.
C     The user should provide INDVAR and NDIREC IF's as appropriate.
C     EARTH will apply the DIFCUT and GP12 modifications after the user
        IF(INDVAR.EQ.R1.OR.INDVAR.EQ.R2) THEN
C.... Cut diffusive links in R1/R2 equations in metal
          IF(IZ.LE.IZ1) THEN
            IF(NDIREC.EQ.1) THEN
              LVAL=LAN
            ELSEIF(NDIREC.EQ.3) THEN
              LVAL=LAE
            ELSEIF(NDIREC.EQ.5) THEN
              LVAL=LD11
            ENDIF
            CALL FN1(LVAL,0.0)
          ENDIF
        ELSEIF(INDVAR.EQ.EPOT.AND.ISWEEP.GT.1) THEN
C... Recompute diffusive links for EPOT, without volume fraction multiplier
          IF(NDIREC.EQ.1) THEN
            L0COND = L0CONDN
            L0DIF  = L0F(LAN)
            L0AR   = L0F(ANORTH)
            L0DEL  = L0F(DYG2D)
            IADD   = 1
          ELSEIF(NDIREC.EQ.3) THEN
            L0COND = L0CONDE
            L0DIF  = L0F(LAE)
            L0AR   = L0F(AEAST)
            L0DEL  = L0F(DXG2D)
            IADD   = NY
          ELSE
            L0COND = L0CONDH
            L0DIF  = L0F(LD11)
            L0AR   = L0F(AHIGH)
            L0DEL  = L0F(LXYDZG)
            IADD   = (IZ-1)*NX*NY
          ENDIF
C... Compute diffusive link as sigma*A*/inter-nodal-distance
C... Note that E & N porosities change to allow interface to float,
C    but H porosity is always 1.0. Internodal distance is changed
C    instead. Sigma is stored in Group 19, Section 4 below.
          DO I = 1,NX*NY
            F(L0DIF+I)=F(L0COND+I+(IZ-1)*NX*NY)*F(L0AR+I)/F(L0DEL+I)
            IF(IADD.LE.NX*NY) THEN
              IF(VAC(I).OR. VAC(I+IADD).OR.POR(I).OR. POR(I+IADD))
     1                                                    F(L0DIF+I)=0.0
            ELSE
              IF(VAC(I).OR.LVAC(I+IADD).OR.POR(I).OR.LPOR(I+IADD))
     1                                                    F(L0DIF+I)=0.0
            ENDIF
          ENDDO
C          CALL PRN('DIF ',-L0DIF)
        ENDIF
      ENDIF
      RETURN
C*****************************************************************
C
C--- GROUP 9. Properties of the medium (or media)
C
    9 GO TO (91,92,90,90,95,96,97,90,90,90,90,90,90,904,90),ISC
   90 CONTINUE
      RETURN
  904 CONTINUE
      NAMSUB='G9.104'
C   * ------------------- SECTION 14 ---------------------------
C     ------ phase-1 specific heat
C--- Set to 1.0, so that EPOT solution (in H1 slot) is not affected
C--- by new H1 diffusion treatment
      CALL FN1(ISPH1,1.0)
      RETURN
   91 CONTINUE
      NAMSUB='G9.1  '
C   * -----------GROUP 9  SECTION  1 ---------------------------
C--- for RHO1.LE.GRND--- density for phase 1          Index DEN1.
C
C--- For RHO1=GRND, DEN1 store must be updated,or EARTH will flag
C    an error. Relaxation of 0.0 ensures retention of initial values
      L0DEN=L0F(DEN1)
      L0PRP=L0F(IPRPS)
      DO 910 I=1,NX*NY
      IF(NINT(F(L0PRP+I)).EQ.51) THEN
C--- Metal
        F(L0DEN+I)=RHOMET
      ELSEIF(NINT(F(L0PRP+I)).EQ.52) THEN
C--- Bath
        F(L0DEN+I)=RHOELC
      ELSEIF(NINT(F(L0PRP+I)).EQ.151) THEN
C--- Carbon block
        F(L0DEN+I)=RHOCAT
      ELSEIF(NINT(F(L0PRP+I)).EQ.100) THEN
C--- Anode
        F(L0DEN+I)=RHOANO
      ELSEIF(NINT(F(L0PRP+I)).EQ.198.OR.NINT(F(L0PRP+I)).EQ.199) THEN
        F(L0DEN+I)=0.0
      ENDIF
  910 CONTINUE
      RETURN
   92 CONTINUE
      NAMSUB='G9.2  '
C   * -----------GROUP 9  SECTION  2 ---------------------------
C--- for DRH1DP.LE.GRND--- D(LN(DEN))/DP for phase 1 (D1DP).
      RETURN
C
   95 CONTINUE
      NAMSUB='G9.5  '
C   * -----------GROUP 9  SECTION  5 ---------------------------
C    For ENUT.LE.GRND--- reference turbulent kinematic viscosity
C                                                  Index VIST
      IF(QEQ(ENUT,GRND)) THEN
C... Different constant turbulent viscosity in metal and bath
        IF(IZ.LE.IZ1)THEN
C ... Metal
          CALL FN1(VIST,ENUTA)
        ELSE IF (IZ.LE.IZ3)THEN
C ... Bath
          CALL FN1(VIST,ENUTB)
        ELSE
C... Air space - non participating
          CALL FN1(VIST,0.0)
        ENDIF
      ENDIF
      RETURN
C
   96 CONTINUE
      NAMSUB='G9.6  '
C   * ------------------- SECTION  6 ---------------------------
C    For ENUL.LE.GRND--- reference laminar kinematic viscosity.
C
      IF(IZ.LE.IZ0) THEN
C--- Cathode - carbon block
       CALL FN1(VISL,0.0)
C--- Metal
      ELSEIF(IZ.LE.IZ1) THEN
       CALL FN1(VISL,ENUMET)
      ELSE
C--- Electrolyte
       CALL FN1(VISL,ENUELC)
      ENDIF
      RETURN
   97 CONTINUE
      NAMSUB='G9.7  '
C   * -----------GROUP 9  SECTION  7 ---------------------------
C
C--- for PRNDTL( ).LE.GRND--- laminar PRANDTL nos., or diffusivity.
      IF(QEQ(PRNDTL(EPOT),-GRND).AND.INDVAR.EQ.EPOT) THEN
C--- Get conductivities for current Z-slab
        L0PRPS=L0F(IPRPS)
        IF(IZ.LT.NZ) THEN
          L0PRPH=L0F(HIGH(IPRPS))
        ELSE
          L0PRPH=L0PRPS
        ENDIF
        CALL CONDCT(L0CND,L0PRPS,L0PRPH)
        CALL FN0(LAMPR,-L0CND)
      ENDIF
      RETURN
C*****************************************************************
C
C--- GROUP 11. Initialization of variable or porosity fields
C
   11 CONTINUE
      NAMSUB='G11   '
      IF(INDVAR.GE.MIN(BX,BY,BZ).AND.INDVAR.LE.MAX(BX,BY,BZ)) THEN
C--- Initialise Bx, By and Bz
        IF(MAGF.EQ.0) THEN
          IF(IMAGF.EQ.0) THEN
C--- Calculate Bx By & Bz from algebraic expressions (PDR 20 case)
            L0XU=L0F(XU2D)
            L0YV=L0F(YV2D)
            L0VAL=L0F(VAL)
            IF(INDVAR.EQ.BX) THEN
C--- Bx
              DO I=1,NX*NY
                F(L0VAL+I)=CONST1*(1.-2.*F(L0YV+I)/F(L0YV+NX*NY))
              ENDDO
            ELSE IF(INDVAR.EQ.BY) THEN
C--- By
              DO I=1,NX*NY
                F(L0VAL+I)=CONST2*(1.-2.*F(L0XU+I)/F(L0XU+NX*NY))
              ENDDO
            ELSE IF(INDVAR.EQ.BZ) THEN
C--- Bz
              DO I=1,NX*NY
                F(L0VAL+I)=CONST3*(1.-2.*F(L0XU+I)/F(L0XU+NX*NY))*
     1                            (1.-2.*F(L0YV+I)/F(L0YV+NX*NY))
              ENDDO
            ENDIF
          ELSE
C--- Calculate Bx By & Bz from polynomial algebraic expressions
            L0XU=L0F(XU2D)
            L0YV=L0F(YV2D)
            L0XG=L0F(XG2D)
            L0YG=L0F(YG2D)
            L0VAL=L0F(VAL)
            XMAX=0.5*F(L0XU+NX*NY)
            YMAX=0.5*F(L0YV+NX*NY)
            IF(INDVAR.EQ.BX) THEN
C--- Bx = a0+a1.x+a2.y+a3.x^2+a4.xy+a5.y^2+a6.x^2+a7.xy^2
              DO I=1,NX*NY
                XG=F(L0XG+I)-XMAX
                YG=F(L0YG+I)-YMAX
                F(L0VAL+I)=AA0 + AA1*XG + AA2*YG + AA3*XG**2 + AA4*XG*YG
     1                   + AA5*YG**2 +AA6*(XG**2)*YG + AA7*XG*YG**2
              ENDDO
            ELSE IF(INDVAR.EQ.BY) THEN
C--- By =  b0+b1.x+b2.y+b3.x^2+b4.xy+b5.y^2+b6.x^2+b7.xy^2
              DO I=1,NX*NY
                XG=F(L0XG+I)-XMAX
                YG=F(L0YG+I)-YMAX
                F(L0VAL+I)=BB0 + BB1*XG + BB2*YG + BB3*XG**2 + BB4*XG*YG
     1                   + BB5*YG**2 +BB6*(XG**2)*YG + BB7*XG*YG**2
              ENDDO
            ELSE IF(INDVAR.EQ.BZ) THEN
C--- Bz =  c0+c1.x+c2.y+c3.x^2+c4.xy+c5.y^2+c6.x^2+c7.xy^2
              DO I=1,NX*NY
                XG=F(L0XG+I)-XMAX
                YG=F(L0YG+I)-YMAX
                F(L0VAL+I)=CC0 + CC1*XG + CC2*YG + CC3*XG**2 + CC4*XG*YG
     1                   + CC5*YG**2 +CC6*(XG**2)*YG + CC7*XG*YG**2
              ENDDO
            ENDIF
          ENDIF
        ELSE IF(MAGF.GT.0) THEN
C--- Read Bx By & Bz from disc file
          IF(IZ.EQ.1.AND.INDVAR.EQ.MIN(BX,BY,BZ)) THEN
            FNAME=CG(1)
 1102       CONTINUE
            CALL DSCFLS(MAGF,FNAME,12,0)
            IF(FNAME(1:2).EQ.'q1') THEN
              DO WHILE(LINE.NE.'  *MAGF')
                READ(MAGF,'(A)',ERR=1104,END=1104) LINE
              ENDDO
              GO TO 1105
 1104         CONTINUE
              IF(FNAME.EQ.'q1') THEN
                CALL DSCFLS(MAGF,FNAME,14,0)
                FNAME='q1ear'
                GO TO 1102
              ENDIF
              WRITE(LUPR1,9981) FNAME
              CALL DSCFLS(MAGF,FNAME,14,0)
              CALL SET_ERR(316, 'Error. See result file.',1)
C            CALL EAROUT(1)
            ENDIF
 1105       CONTINUE
          ENDIF
          IF(INDVAR.EQ.MIN(BX,BY,BZ)) THEN
            DO ICOMPON = 1,3
              DO IY =NY, 1, -1
                IF(ICOMPON .EQ. 1)THEN
                  READ(MAGF,*,END=1110,ERR=1108) (GBX(IY,IX),IX=1,NX)
                ELSE IF(ICOMPON .EQ. 2)THEN
                  READ(MAGF,*,END=1110,ERR=1108) (GBY(IY,IX),IX=1,NX)
                ELSE
                  READ(MAGF,*,END=1110,ERR=1108) (GBZ(IY,IX),IX=1,NX)
                END IF
              ENDDO
            ENDDO
            WRITE(LUPR1,9980) FNAME,IZ
 9980       FORMAT(/' MAGNETIC FIELDS READ FROM FILE ',A,' IZ = ',I3/)
            GO TO 1110
 1108       WRITE(LUPR1,9981) FNAME
 9981       FORMAT(' ERROR IN READING FIELDS FILE ',A)
            CALL DSCFLS(MAGF,FNAME,14,0)
            CALL SET_ERR(317, 'Error. See result file.',1)
C            CALL EAROUT(1)
          ENDIF
 1110     IF(INDVAR.EQ.BX) CALL SETYX(VAL,GBX,NYM,NXM)
          IF(INDVAR.EQ.BY) CALL SETYX(VAL,GBY,NYM,NXM)
          IF(INDVAR.EQ.BZ) CALL SETYX(VAL,GBZ,NYM,NXM)
          IF(IZ.EQ.NZ.AND.INDVAR.EQ.MIN(BX,BY,BZ)) THEN
            CALL DSCFLS(MAGF,FNAME,14,0)
          ENDIF
        ENDIF
C--- Convert from GAUSS to TESLA
        IF(IMAGU.EQ.1) THEN
          CALL FN2(VAL,VAL,0.0,1.E-4)
C--- Convert from milliTesla to TESLA
        ELSEIF(IMAGU.EQ.2) THEN
          CALL FN2(VAL,VAL,0.0,1.E-3)
        ENDIF
      ENDIF
C
      IF(INDVAR.EQ.HI) THEN
C--- Initialise HI store : cell face heights
       IF(IZ.EQ.1) THEN
C--- Initilise interface, anode and free surface height
        CALL GETZ(ZWNZ,GZWNZ,NZM)
C--- Store interface height in GRSP2
        IF(NPATCH(1:5).NE.'INIF1') THEN
          IF(HUMP.EQ.0) THEN
C---  Flat interface
            CALL FN1(GRSP2,GZWNZ(IZ1))
          ELSE
C---  Read interface height from 'HUMP' file
            FNAME=CG(6)
 1142       CONTINUE
            CALL DSCFLS(HUMP,FNAME,12,0)
            IF(FNAME(1:2).EQ.'q1') THEN
              DO WHILE(LINE.NE.'  *HUMP')
                READ(HUMP,'(A)',ERR=1144,END=1144) LINE
              ENDDO
            ENDIF
            DO IY =NY, 1, -1
              READ(HUMP,*,ERR=1144,END=1144)(GVAL(IY,IX),IX=1,NX)
            ENDDO
            CALL DSCFLS(HUMP,FNAME,14,0)
            WRITE(LUPR1,'(/'' INITIAL INTERFACE HEIGHT READ FROM FILE ''
     1                                                     ,A)') FNAME
            GO TO 1146
 1144       CONTINUE
            IF(FNAME.EQ.'q1') THEN
              CALL DSCFLS(HUMP,FNAME,14,0)
              FNAME='q1ear'
              GO TO 1142
            ENDIF
            WRITE(LUPR1,'('' ERROR IN READING HUMP FILE '',A/)') FNAME
            CALL SET_ERR(318, 'ERROR IN READING HUMP FILE.',1)
C            CALL EAROUT(1)
 1146       CONTINUE
            CALL SETYX(GRSP2,GVAL,NYM,NXM)
          ENDIF
        ELSE
C
C--- Sequence for setting inclined interface for transient test case
C
         GZ0=GZWNZ(IZ1)
         IF(NPATCH(6:6).EQ.' '.OR.NPATCH(6:6).EQ.'Y') THEN
C--- PATCH name is INIF1 or INIF1Y. Slope in Y direction
           CALL GETYX(YG2D,GYV,NYM,NXM)
           CALL GETONE(YV2D,GLEN,NY,NX)
         ELSE
C--- Slope in X direction
           CALL GETYX(XG2D,GYV,NYM,NXM)
           CALL GETONE(XU2D,GLEN,NY,NX)
         ENDIF
         GDHDY=(GDZ2-GDZ1)/GLEN
         GZ1=GZ0+GDZ1
         DO IX=1,NX
           DO IY=1,NY
             GVAL(IY,IX)=GZ1+GDHDY*GYV(IY,IX)
           ENDDO
         ENDDO
         CALL SETYX(GRSP2,GVAL,NYM,NXM)
        ENDIF
C
C--- Store anode height in GRSP3
        IF(HANO.NE.0) THEN
C---  Read anode underside heights from disc file
         FNAME=CG(5)
 1111    CONTINUE
         CALL DSCFLS(HANO,FNAME,12,0)
         IF(FNAME(1:2).EQ.'q1') THEN
           DO WHILE (LINE.NE.'  *HANO')
             READ(HANO,'(A)',ERR=1113,END=1113) LINE
           ENDDO
         ENDIF
         DO IY=NY,1,-1
           READ(HANO,*,ERR=1113,END=1113) (HANODE(IY,IX),IX=1,NX)
         ENDDO
         CALL DSCFLS(HANO,FNAME,14,0)
         WRITE(LUPR1,'(/'' ANODE HEIGHTS READ FROM FILE '',A/)') FNAME
         GO TO 11132
 1113    CONTINUE
         IF(FNAME.EQ.'q1') THEN
           CALL DSCFLS(HANO,FNAME,14,0)
           FNAME='q1ear'
           GO TO 1111
         ENDIF
         WRITE(LUPR1,'('' ERROR IN READING HANO FILE '',A)') FNAME
         CALL SET_ERR(319, 'ERROR IN READING HANO FILE.',1)
C         CALL EAROUT(1)
11132    CONTINUE
         CALL SETYX(GRSP3,HANODE,NYM,NXM)
         IF(IZ1.GT.0) THEN
           IF(IZ0.GT.0) THEN
             CALL FN2(GRSP3,GRSP3,GZWNZ(IZ0),1.0)
           ELSE
             CALL FN2(GRSP3,GRSP3,0.0,1.0)
           ENDIF
         ENDIF
        ELSE
C--- Flat anode underside
         CALL FN1(GRSP3,GZWNZ(IZ2))
        ENDIF
C
C--- Store free surface height in GRSP4
C---  flat free surface
        CALL FN1(GRSP4,GZWNZ(IZ3))
       ENDIF
C--- Compute HI from GRSP2,GRSP3 AND GRSP4
       IF(IZ.LE.IZ0) THEN
        CALL FN1(VAL,GZWNZ(IZ))
       ELSEIF(IZ.LE.IZ1) THEN
        CALL FN2(VAL,GRSP2,0.0,GZWNZ(IZ)/GZWNZ(IZ1))
       ELSE IF(IZ.LE.IZ2) THEN
        ALFA=(GZWNZ(IZ)-GZWNZ(IZ1))/(GZWNZ(IZ2)-GZWNZ(IZ1))
        CALL FN10(VAL,GRSP3,GRSP2,0.0,ALFA,1.-ALFA)
       ELSE IF(IZ.LE.IZ3) THEN
        ALFA=(GZWNZ(IZ)-GZWNZ(IZ2))/(GZWNZ(IZ3)-GZWNZ(IZ2))
        CALL FN10(VAL,GRSP4,GRSP3,0.0,ALFA,1.-ALFA)
       ELSE
        ALFA=(GZWNZ(IZ)-GZWNZ(IZ3))/(GZWNZ(NZ)-GZWNZ(IZ3))
        CALL FN2(VAL,GRSP4,0.0,ALFA)
       ENDIF
      ENDIF
C
C--- Initialise Lorentz Forces from J*B, ensuring zero values at
C    domain edges
      IF(INDVAR.EQ.FX) THEN
C--- Fx
       CALL FN48(VAL,JY,BZ,JZ,BY,1.0,-1.0)
       IXF=NX
       CALL FN1(VAL,0.0)
       IXF=1
      ELSE IF(INDVAR.EQ.FY) THEN
C--- Fy
       CALL FN48(VAL,JZ,BX,JX,BZ,1.0,-1.0)
       IYF=NY
       CALL FN1(VAL,0.0)
       IYF=1
      ELSE IF(INDVAR.EQ.FZ) THEN
C--- Fz
       IF(IZ.LT.NZ) THEN
        CALL FN48(VAL,JX,BY,JY,BX,1.0,-1.0)
       ELSE
        CALL FN1(VAL,0.0)
       ENDIF
      ENDIF
C
      IF(INDVAR.EQ.IPRPS) THEN
C--- Initialise PRPS fields
        IF(IZ.EQ.1) THEN
C--- Read IANODE and IFRZ arrays once
          IF(IANO.NE.0) THEN
C--- Read IANODE array from disc file
            FNAME=CG(4)
 1114       CONTINUE
            CALL DSCFLS(IANO,FNAME,12,0)
            IF(FNAME(1:2).EQ.'q1') THEN
              DO WHILE (LINE.NE.'  *IANO')
                READ(IANO,'(A)',ERR=1118,END=1118) LINE
              ENDDO
            ENDIF
            DO IY=NY,1,-1
              READ(IANO,*,ERR=1118,END=1118) (IANODE(IY,IX),IX=1,NX)
            ENDDO
            CALL DSCFLS(IANO,FNAME,14,0)
            WRITE(LUPR1,'(/'' ANODE LOCATIONS READ FROM FILE '',A/)')
     1                                                             FNAME
            GO TO 1122
 1118       CONTINUE
            IF(FNAME.EQ.'q1') THEN
              CALL DSCFLS(IANO,FNAME,14,0)
              FNAME='q1ear'
              GO TO 1114
            ENDIF
            WRITE(LUPR1,'('' ERROR IN READING IANO FILE '',A)') FNAME
            CALL SET_ERR(320, 'ERROR IN READING IANO FILE.',1)
C            CALL EAROUT(1)
          ENDIF
 1122     IF(IFRZ.NE.0) THEN
C--- Read IFREEZ array from disc file
            FNAME=CG(3)
 1124       CONTINUE
            CALL DSCFLS(IFRZ,FNAME,12,0)
            IF(FNAME(1:2).EQ.'q1') THEN
              DO WHILE (LINE.NE.'  *IFRZ')
                READ(IFRZ,'(A)',ERR=1130,END=1130) LINE
              ENDDO
            ENDIF
            DO IY=NY,1,-1
              READ(IFRZ,*,ERR=1130,END=1130) (IFREEZ(IY,IX),IX=1,NX)
            ENDDO
            CALL DSCFLS(IFRZ,FNAME,14,0)
            WRITE(LUPR1,'(/,'' FREEZE HEIGHTS READ FROM FILE '',A/)')
     1                                                             FNAME
            GO TO 1134
 1130       CONTINUE
            IF(FNAME.EQ.'q1') THEN
              CALL DSCFLS(IFRZ,FNAME,14,0)
              FNAME='q1ear'
              GO TO 1124
            ENDIF
            WRITE(LUPR1,'('' ERROR IN READING FILE '',A)') FNAME
            CALL SET_ERR(321, 'ERROR IN READING FILE.',1)
C            CALL EAROUT(1)
          ENDIF
 1134     CONTINUE
        ENDIF
C--- Setting PRPS field. PRPS used only to identify liquid/solid, not to
C    get properties. These are set as before in ESTRGR, Group 9
C        51 - metal, 52 - electrolyte, fluids
C       151 - cathode,  participating solid
C       100 - anode, participating solid
C       198 - freeze, non-participating solid with friction
C       199 - air, non-participating solid with no friction
        GDEN=151.
        IF(IZ.GT.IZ0) GDEN=51.
        IF(IZ.GT.IZ1) GDEN=52.
        IF(IZ.GT.IZ3) GDEN=199.
        L0VAL=L0F(VAL)
        DO IX=1,NX
          DO IY=1,NY
            I=(IX-1)*NY+IY
            F(L0VAL+I)=GDEN
            IF(IZ.GT.IZ0.AND.IZ.LE.IFREEZ(IY,IX)+IZ0) F(L0VAL+I)=198.
            IF(IZ.GT.IZ2 .AND. IANODE(IY,IX).EQ.1) F(L0VAL+I)=100.
          ENDDO
        ENDDO
      ENDIF
C
      IF(INDVAR.EQ.JZ) THEN
C--- Set Jz to cathode distribution initially - GJCATH read in G1.1
        CALL SETYX(VAL,GJCATH,NYM,NXM)
      ENDIF
C
      RETURN
C*****************************************************************
C
C--- GROUP 13. Boundary conditions and special sources
C
   13 CONTINUE
      IF(ISC.NE.1) GO TO 1311
      NAMSUB='G13C  '
C--------------------------------- coefficient = GRND
C
      L0CO=L0F(CO)
      L0PRPS=L0F(IPRPS)
C--- Set anode potential : VALue set in Q1 or below.
      IF((NPATCH(1:7).EQ.'ANOPTNT'.OR.OBJNAM(IPAT).EQ.'ANOPTNTL').AND.
     1                                INDVAR.EQ.EPOT.AND.IZ.EQ.NZ) THEN
       L0DZ=L0F(LXYDZ)
       DO I = 1,NX*NY
         IF(QEQ(F(L0PRPS+I),100.)) THEN
           F(L0CO+I)=2.0*CONANO/F(L0DZ+I)
         ELSE
           F(L0CO+I)=0.0
         ENDIF
       ENDDO
       RETURN
      ENDIF
C
C--- Gas outlet
      IF((NPATCH.EQ.'SURFACE'.OR.OBJNAM(IPAT).EQ.'SURFACE').AND.
     1                                              INDVAR.EQ.P2) THEN
       DO I = 1,NX*NY
         IF(QNE(F(L0PRPS+I),100.).AND.QNE(F(L0PRPS+I),198.)) THEN
           F(L0CO+I)=GPCOEF
         ELSE
           F(L0CO+I)=0.0
         ENDIF
       ENDDO
       RETURN
      ENDIF
      RETURN
C
 1311 IF(ISC.NE.12) RETURN
C---------------------------------------- value = GRND
      NAMSUB='G13V  '
C
      IF(NPATCH(1:7).EQ.'ANOPTNT'.OR.OBJNAM(IPAT).EQ.'ANOPTNTL') THEN
        IF (INDVAR .EQ. EPOT) THEN
C--- Call subroutine to obtain potentials at upper face of anodes
C    Version supplied returns constant value set in Q1.
          CALL POTENT (GVAL,NYM,NXM)
C--- Set into EARTH
          CALL SETYX(VAL,GVAL,NYM,NXM)
        ENDIF
      ENDIF
C
      IF(NPATCH.EQ.'LORENTZ'.OR.OBJNAM(IPAT).EQ.'LORENTZ') THEN
C--- Calculate and add Lorentz body-forces from F = J * B
C--- B field stored at cell centres. J values stored at cell faces.
C    F values computed at required cell faces, using interpolated
C    values of B and J. F is set to 0.0 inside and on faces of anodes
C    and freeze.
C
C--- Get variables for future use
        L0JX=L0F(JX)
        L0JY=L0F(JY)
        L0JZ=L0F(JZ)
        IF(IZ.GT.1) L0JZL=L0F(LOW(JZ))
        L0BX=L0F(BX)
        L0BY=L0F(BY)
        L0BZ=L0F(BZ)
        L0VAL=L0F(VAL)
C
        IF(INDVAR.EQ.U1) THEN
C--- Fx = Jy * Bz - Jz * By
          L0DXU=L0F(DXU2D)
          L0DXG=L0F(DXG2D)
          L0FX=L0F(FX)
C
          DO 13111 IX=1,NX
          DO 13111 IY=1,NY
            I=(IX-1)*NY+IY
            IF(SLD(I).OR.SLD(I+NY).OR.IX.EQ.NX) THEN
C--- Inside anode or freeze - Fx = 0.0
              F(L0VAL+I)=0.0
            ELSE
C--- Average JY to cell face
              IF(IY.GT.1) THEN
                GJYB=0.25*(F(L0DXU+I+NY) *(F(L0JY+I)+F(L0JY+I-1))
     1           +F(L0DXU+I)  *(F(L0JY+I+NY)+F(L0JY+I-1+NY)))/F(L0DXG+I)
              ELSE
                GJYB=0.25*(F(L0DXU+I+NY) *F(L0JY+I)
     1           +F(L0DXU+I)  *F(L0JY+I+NY))/F(L0DXG+I)
              ENDIF
C--- Average JZ to cell face
              IF(IZ.GT.1) THEN
                GJZB=0.25*(F(L0DXU+I+NY)*(F(L0JZ+I)+F(L0JZL+I))
     1           +F(L0DXU+I)*(F(L0JZ+I+NY)+F(L0JZL+I+NY)))/F(L0DXG+I)
              ELSE
                GJZB=0.25*(F(L0DXU+I+NY)*(F(L0JZ+I)+GJCATH(IY,IX))
     1           +F(L0DXU+I)*(F(L0JZ+I+NY)+GJCATH(IY,IX+1)))/F(L0DXG+I)
              ENDIF
C--- Average Bz and By to cell face
              GBZB=0.5*(F(L0DXU+I+NY)*F(L0BZ+I)+F(L0DXU+I)*F(L0BZ+I+NY))
     1            /F(L0DXG+I)
              GBYB=0.5*(F(L0DXU+I+NY)*F(L0BY+I)+F(L0DXU+I)*F(L0BY+I+NY))
     1            /F(L0DXG+I)
C--- Compute Fx
              GFX=GJYB*GBZB-GJZB*GBYB
C--- Transfer to VALue store
              F(L0VAL+I)=RELF*GFX + (1.0-RELF)*F(L0FX+I)
            ENDIF
13111     CONTINUE
C--- Store Fx in FX for printout
          CALL FN0(FX, VAL)
C
        ELSE IF(INDVAR.EQ.V1) THEN
C--- Fy = Jz * Bx - Jx * Bz
          L0DYV=L0F(DYV2D)
          L0DYG=L0F(DYG2D)
          L0FY=L0F(FY)
          DO 13112 IX=1,NX
          DO 13112 IY=1,NY
            I=(IX-1)*NY+IY
            IF(SLD(I).OR.SLD(I+1).OR.IY.EQ.NY) THEN
C--- Inside anode or freeze - Jy = 0.0
              F(L0VAL+I)=0.0
            ELSE
C--- Average JX to cell face
              IF(IX.GT.1) THEN
                GJXB=0.25*(F(L0DYV+I+1)*(F(L0JX+I)+F(L0JX+I-NY))
     1           +F(L0DYV+I)*(F(L0JX+I+1)+F(L0JX+I+1-NY)))/F(L0DYG+I)
              ELSE
                GJXB=0.25*(F(L0DYV+I+1)*F(L0JX+I)
     1           +F(L0DYV+I)*F(L0JX+I+1))/F(L0DYG+I)
              ENDIF
C--- Average JZ to cell face
              IF(IZ.GT.1) THEN
                GJZB=0.25*(F(L0DYV+I+1)*(F(L0JZ+I)+F(L0JZL+I))
     1           +F(L0DYV+I)*(F(L0JZ+I+1)+F(L0JZL+I+1)))/F(L0DYG+I)
              ELSE
                GJZB=0.25*(F(L0DYV+I+1)*(F(L0JZ+I)+GJCATH(IY,IX))
     1           +F(L0DYV+I)*(F(L0JZ+I+1)+GJCATH(IY+1,IX)))/F(L0DYG+I)
              ENDIF
C--- Average Bz & Bx to cell faces
              GBZB=0.5*(F(L0DYV+I+1)*F(L0BZ+I)+F(L0DYV+I)*F(L0BZ+I+1))/
     1          F(L0DYG+I)
              GBXB=0.5*(F(L0DYV+I+1)*F(L0BX+I)+F(L0DYV+I)*F(L0BX+I+1))/
     1          F(L0DYG+I)
C--- Compute Fy
              GFY=GJZB*GBXB-GJXB*GBZB
C--- Transfer to VALue store
              F(L0VAL+I)=RELF*GFY + (1.0-RELF)*F(L0FY+I)
            ENDIF
13112     CONTINUE
C--- Store Fy in FY for printout
          CALL FN0(FY,VAL)
C
        ELSE IF(INDVAR.EQ.W1) THEN
C--- Fz = Jx * By - Jy * Bx
          IF(IZ.LT.NZ) THEN
            L0DZW=L0F(LXYDZ)
            L0DZG=L0F(LXYDZG)
            L0FZ=L0F(FZ)
            L0HI=L0F(HI)
            L0HIH=L0F(HIGH(HI))
            L0JXH=L0F(HIGH(JX))
            L0JYH=L0F(HIGH(JY))
            L0BXH=L0F(HIGH(BX))
            L0BYH=L0F(HIGH(BY))
            IADD=(IZ-1)*NX*NY
            DO 13113 IX=1,NX
            DO 13113 IY=1,NY
              I=(IX-1)*NY+IY
              IF(SLD(I).OR.LSOLID(I+IADD)) THEN
C--- Inside anode or freeze - Fz = 0.0
                F(L0VAL+I)=0.0
              ELSE
                GDELH=F(L0HIH+I)-F(L0HI+I)
C--- Average JY to cell face
                IF(IY.GT.1) THEN
                  GJYB=0.25*(GDELH*(F(L0JY+I)+F(L0JY+I-1))
     1            +F(L0DZW+I)*(F(L0JYH+I)+F(L0JYH+I-1)))/F(L0DZG+I)
                ELSE
                  GJYB=0.25*(GDELH*F(L0JY+I)
     1            +F(L0DZW+I)*F(L0JYH+I))/F(L0DZG+I)
                ENDIF
C--- Average JX to cell face
                IF(IX.GT.1) THEN
                  GJXB=0.25*(GDELH*(F(L0JX+I)+F(L0JX+I-NY))
     1            +F(L0DZW+I)*(F(L0JXH+I)+F(L0JXH+I-NY)))/F(L0DZG+I)
                ELSE
                  GJXB=0.25*(GDELH*F(L0JX+I)
     1            +F(L0DZW+I)*F(L0JXH+I))/F(L0DZG+I)
                ENDIF
C--- Average Bx & By to cell face
                GBXB=0.5*(GDELH*F(L0BX+I)+F(L0DZW+I)*F(L0BXH+I))/
     1           F(L0DZG+I)
                GBYB=0.5*(GDELH*F(L0BY+I)+F(L0DZW+I)*F(L0BYH+I))/
     1           F(L0DZG+I)
C--- Compute Fz
                GFZ=GJXB*GBYB-GJYB*GBXB
C--- Transfer to VALue store
                F(L0VAL+I)=RELF*GFZ + (1.0-RELF)*F(L0FZ+I)
              ENDIF
13113       CONTINUE
C--- Store Fz in FZ for printout
            CALL FN0(FZ,VAL)
          ENDIF
        ENDIF
C--- Multiply by volume fraction for 2-Phase
        IF(.NOT.ONEPHS) CALL FN21(VAL,VAL,R1,0.0,1.0)
      ENDIF
C
C--- Gravity sources for W2
      IF(INDVAR.EQ.W2.AND.NPATCH.EQ.'GRAVW2') THEN
        CALL FN2(VAL,R2,0.0,AGRAVZ*(RHO2-RHOELC))
      ENDIF
C
C--- Gravity sources for U2
      IF(INDVAR.EQ.U2.AND.NPATCH.EQ.'GRAVU2') THEN
        GDRH2P=AGRAVZ*(RHO2-RHOELC)
        CALL FNAVXY(VAL,HI,0.0,-GDRH2P,GDRH2P,X)
        CALL FN21(VAL,VAL,R2,0.0,1.0)
      ENDIF
C
C--- Gravity sources for V2
      IF(INDVAR.EQ.V2.AND.NPATCH.EQ.'GRAVV2') THEN
        GDRH2P=AGRAVZ*(RHO2-RHOELC)
        CALL FNAVXY(VAL,HI,0.0,-GDRH2P,GDRH2P,Y)
        CALL FN21(VAL,VAL,R2,0.0,1.0)
      ENDIF
C
C--- Compute Induced-Current source for Electric Potential
C    S = div (Ji)
      IF((NPATCH.EQ.'INDUCE'.OR.OBJNAM(IPAT).EQ.'INDUCE').AND.
     1                                           INDVAR.EQ.EPOT) THEN
        IF(LG(1)) THEN
          L0AE=L0F(AEAST)
          L0AN=L0F(ANORTH)
          L0AH=L0F(AHIGH)
          L0JIX=L0F(JIX)
          L0JIY=L0F(JIY)
          L0JIZ=L0F(JIZ)
          IF(IZ.GT.1) L0JIZL=L0F(LOW(JIZ))
          L0VAL=L0F(VAL)
          IADD=(IZ-1)*NX*NY
          DO 13114 IX=1,NX
          DO 13114 IY=1,NY
            I=(IX-1)*NY+IY
            IF(SLD(I)) THEN
              GSU=0.0
            ELSE
              GSU=F(L0AE+I)*F(L0JIX+I)+F(L0AN+I)*F(L0JIY+I)+F(L0AH+I)
     1                                                      *F(L0JIZ+I)
              IF(IX.GT.1.AND..NOT.SLD(I-NY))
     1                         GSU=GSU-F(L0AE+I-NY)*F(L0JIX+I-NY)
              IF(IY.GT.1.AND..NOT.SLD(I-1))
     1                         GSU=GSU-F(L0AN+I-1)*F(L0JIY+I-1)
              IF(IZ.GT.1.AND..NOT.LSOLID(I-IADD))
     1                         GSU=GSU-F(L0AH+I)*F(L0JIZL+I)
            ENDIF
13114       F(L0VAL+I)=GSU
          IF(.NOT.ONEPHS) CALL FN26(VAL,R1)
        ELSE
          CALL FN1(VAL,0.0)
        ENDIF
        RETURN
      ENDIF
C
C--- Under-anode gas sources
      IF(INDVAR.EQ.P2.AND.NPATCH.EQ.'GAS') THEN
        L0PRP=L0F(HIGH(IPRPS))
        L0VAL=L0F(VAL)
        IF(GMDOT.GT.0.0) THEN
C... Constant release rate set in Q1
          DO 13115 I=1,NX*NY
            IF(NINT(F(L0PRP+I)).EQ.100) THEN
              F(L0VAL+I)=GMDOT
            ELSE
              F(L0VAL+I)=0.0
            ENDIF
13115     CONTINUE
          RETURN
        ELSE
C... Release rate calculated from Vol=2.62E-4 Jz m^3/m^2 , Jz in kA.
          L0JZ=L0F(JZ)
C... multiply constant by 1E-3 to convert current in A to kA
          GCONST=2.62E-4*1.0E-3
          DO I=1,NX*NY
            IF(NINT(F(L0PRP+I)).EQ.100) THEN
              F(L0VAL+I)=RHO2*GCONST*ABS(F(L0JZ+I))
            ELSE
              F(L0VAL+I)=0.0
            ENDIF
          ENDDO
          RETURN
        ENDIF
      ENDIF
C
C--- Gas outlet
      IF((NPATCH.EQ.'SURFACE'.OR.OBJNAM(IPAT).EQ.'SURFACE').AND.
     1                                                INDVAR.EQ.P2) THEN
        L0R2=L0F(R2)
        L0VAL=L0F(VAL)
        DO 13116 I=1,NX*NY
          IF(SLD(I)) THEN
            F(L0VAL+I)=0.0
          ELSE
            F(L0VAL+I)=AGRAVZ*RHO2*F(L0R2+I)/CFIPS
          ENDIF
13116   CONTINUE
        RETURN
      ENDIF
C
C--- Fix gas rise velocity
      IF(INDVAR.EQ.W2.AND.NPATCH.EQ.'RISE') THEN
        L0VAL=L0F(VAL)
        DO 13117 I=1,NX*NY
          IF(SLD(I)) THEN
            F(L0VAL+I)=0.0
          ELSE
            F(L0VAL+I)=ABS(AGRAVZ/CFIPS)
          ENDIF
13117   CONTINUE
        RETURN
      ENDIF
C
C--- Set cathode current density for potential eqn.
      IF(INDVAR.EQ.EPOT.AND.(NPATCH(1:7).EQ.'CATHCUR'.OR.OBJNAM(IPAT)
     1                                         .EQ.'CATHCURR')) THEN
C--- GJCATH set in Group 1. The cathode current should be zero
C    inside the freeze, or else very large potentials will arise.
        CALL SETYX(VAL,GJCATH,NYM,NXM)
      ENDIF
C
C--- Calculate W at interface
      IF((NPATCH.EQ.'INTFACE'.OR.OBJNAM(IPAT).EQ.'INTFACE').AND.
     1                              (INDVAR.EQ.W1.OR.INDVAR.EQ.W2)) THEN
        IF(.NOT.STEADY) THEN
C---  W = ( HI - HIold ) / dt  for transient
          CALL FN10(VAL,HI,OLD(HI),0.0,1./DT,-1./DT)
          IF(LG(11).AND.GDBG) THEN
C---  Debug output
            CALL PRN('CURH',HI)
            CALL PRN('OLDH',OLD(HI))
            CALL PRN('VAL ',VAL)
          ENDIF
        ELSE
C---  W = 0.0  for steady state
          CALL FN1(VAL,0.0)
        ENDIF
      ENDIF
C
C
      RETURN
C***************************************************************
C
C--- GROUP 19. Special calls to GROUND from EARTH
C
   19 NAMSUB='G19   '
      GO TO (191,191,193,194,191,196,191,191),ISC
  191 CONTINUE
      RETURN
  193 CONTINUE
C   * ------------------- SECTION 3 ---- START OF SLAB
C
C---  Set debug switch
      GDBG=IZ.GE.IG(17).AND.IZ.LE.IG(18).AND.ISWEEP.GE.IG(19).AND.
     1      ISWEEP.LE.IG(20)
C
C--- For RESTART only, recover interface,anode & free surface
C     heights on first sweep
      IF(ISWEEP.EQ.FSWEEP .AND. ISTEP.EQ.FSTEP .AND.
     1   QEQ(FIINIT(HI),READFI) ) THEN
       IF(IZ.EQ.1) THEN
C--- Interface - GRSP2
        CALL FN0(GRSP2,ANYZ(HI,IZ1))
        CALL PRN('HINT',GRSP2)
C--- Anode undersurface - GRSP3
        CALL FN0(GRSP3,ANYZ(HI,IZ2))
        CALL PRN('HANO',GRSP3)
C--- Free surface - GRSP4
        CALL FN0(GRSP4,ANYZ(HI,IZ3))
        CALL PRN('HAIR',GRSP4)
       ENDIF
       IF(IZ.EQ.IZ2+1) THEN
C--- Recreate IANODE and IFREEZ arrays in case needed in POTENT
        L0PRP=L0F(IPRPS)
        DO 19301 IX=1,NX
        DO 19301 IY=1,NY
        I=(IX-1)*NY+IY
        IF(NINT(F(L0PRP+I)).EQ.100) IANODE(IY,IX)=1
        IF(NINT(F(L0PRP+I)).EQ.198) IFREEZ(IY,IX)=1
19301   CONTINUE
       ENDIF
      ENDIF
C
C--- Guess new interface position on 1st sweep of new time-step
C     Hnew = Hold + dH/dT * DT   : [ dH/dT = W1 at IZ = IZ1 ]
      IF(ISWEEP.EQ.1.AND..NOT.STEADY.AND.IZ.EQ.1) THEN
       CALL FN34(GRSP2,ANYZ(W1,IZ1),DT)
      ENDIF
C
C--- Adjust interface height
C    -----------------------
C
      IF(MOD(ISWEEP,NIH).EQ.0.AND.IZ.EQ.IZ1.AND.ISWEEP.GE.IHF) THEN
       IF(IREF.EQ.0) THEN
C--- Search for reference pressure point (not in blockage)
        CALL SUB3R(VLOWER,0.,GDRHO,AGRAVZ*(RHOMET-RHOELC),GAREAM,0.)
        DO 19302 I=1,NX*NY
        IF(.NOT.SLD(I)) THEN
          IF(IZ1.GT.1) THEN
            IADD=(IZ-2)*NX*NY
            IF(.NOT.LSOLID(I+IADD)) GOTO 19304
          ELSE
            GO TO 19304
          ENDIF
        ENDIF
19302   CONTINUE
C--- Reference point found, store coordinates
19304   IREF=I
C--- Compute and store total volume below interface and area of metal
        L0AH=L0F(AHIGH)
        L0HI=L0F(GRSP2)
        Z0=0.0
        IF(IZ0.GT.0) Z0=GZWNZ(IZ0)
        DO 19305 I=1,NX*NY
          IF(SLD(I)) GO TO 19305
          VLOWER=VLOWER+(F(L0HI+I)-Z0)*F(L0AH+I)
          GAREAM=GAREAM+F(L0AH+I)
19305   CONTINUE
C--- Get variables for later use
        CALL SUB2R(DHLIM,FDH*GZWNZ(IZ1),HMIN,(1.-FHLIM)*GZWNZ(IZ1))
        L0ZW=L0F(ZWNZ)
       ENDIF
C
C--- Get new interface height from hydrostatic balance. Lorentz
C    forces included in balance if LG(3)=.T.
C
C--- Get required variables
       CALL GETZ(ZWNZ,GZWNZ,NZM)
       L0ZW=L0F(ZWNZ)
       L0YV=L0F(YV2D)
       L0XU=L0F(XU2D)
       GDISTY=F(L0YV+(NX-1)*NY+NY)
       GDISTX=F(L0XU+(NX-1)*NY+NY)
       GYDIST=AMAX1(GDISTY,GDISTX)
       IF(LG(3)) THEN
        L0FZ =L0F(FZ)
        L0FZH=L0F(HIGH(FZ))
        L0FZL=L0F(LOW(FZ))
        L0HI =L0F(HI)
        L0HIH=L0F(HIGH(HI))
        L0HIL=L0F(LOW(HI))
        GZ1=0.5*(F(L0HI+IREF)+F(L0HIL+IREF))
        GZ2=0.5*(F(L0HIH+IREF)+F(L0HI+IREF))
        GFZ1=0.5*(F(L0FZ+IREF)+F(L0FZL+IREF))
        GFZ2=0.5*(F(L0FZH+IREF)+F(L0FZ+IREF))
        DFZREF=GFZ1*GZ1-GFZ2*GZ2
        IF(LG(9).AND.GDBG) THEN
          WRITE(14,'(''z1,z2,f1,f2,dfzref '',1P5E12.4)') GZ1,GZ2,GFZ1,
     1         GFZ2,DFZREF
        ENDIF
       ELSE
        DFZREF=0.0
        DFZ=0.0
       ENDIF
       L0HI=L0F(GRSP2)
       L0HA=L0F(GRSP3)
       L0P1=L0F(P1)
       L0P1H=L0F(HIGH(P1))
       CALL SUB3R(VOLSUM,0.,HREF,F(L0HI+IREF),DPREF,F(L0P1H+IREF)-
     1            F(L0P1+IREF))
C
       IF(LG(9).AND.GDBG) THEN
C--- Debug for adjustment
        WRITE(LUPR1,9983) ISWEEP
 9983   FORMAT(1X,' INTERFACE ADJUSTMENT AT ISWEEP= ',I4)
        WRITE(LUPR1,9984) VLOWER,HREF,DPREF,GDRHO
 9984   FORMAT(1X,' VOL ',1PE13.6,' HREF ',1PE13.6,' DPREF ',
     1    1PE13.6,' GDRHO ',1PE13.6)
        CALL PRN('P1  ',P1)
        CALL PRN('P1H ',HIGH(P1))
        IF(LG(3)) THEN
         CALL PRN('FZH ',HIGH(FZ))
         CALL PRN('FZ  ',FZ)
         CALL PRN('FZL ',LOW(FZ))
        ENDIF
        CALL PRN('HI  ',GRSP2)
       ENDIF
C
C--- Perform adjustment
       FROUDE=RHOMET*(GYDIST/DT)**2/(ABS(GDRHO)*GZWNZ(IZ1))
C       FROUDE=RHOMET*(GYDIST/DT)**2/(ABS(GDRHO)*F(L0ZW+IZ1))
       DO 19306 I=1,NX*NY
       IF(SLD(I)) GO TO 19306
       FAC=GDRHO
       IF(LG(3)) THEN
         FAC=FAC-0.5*(F(L0FZH+I)-F(L0FZL+I))
         GZ3 = 0.5*(F(L0HI+I) +F(L0HIL+I))
         GZ4 = 0.5*(F(L0HIH+I)+F(L0HI+I))
         GFZ3= 0.5*(F(L0FZ+I) +F(L0FZL+I))
         GFZ4= 0.5*(F(L0FZH+I)+F(L0FZ+I))
         DFZ = GFZ3*GZ3-GFZ4*GZ4
         IF(LG(9).AND.GDBG) THEN
           F(L0FZDIF+I)=DFZREF-DFZ
           F(L0PDIF+I)=F(L0P1H+I)-DPREF
         ENDIF
       ELSE
         DFZREF=0.0
         DFZ=0.0
       ENDIF
       GRES=(HREF-F(L0HI+I))*FAC-DPREF+F(L0P1H+I)-
     1               F(L0P1+I)+DFZREF-DFZ
       CHANGE=GRES/(FAC*(1.0+FROUDE*FRCON))
C--- Relax change
       CHANGE=SLOH*CHANGE
       IF(CHANGE.LT.0.) THEN
        CHANGE=AMAX1(CHANGE,-DHLIM)
       ELSE
        CHANGE=AMIN1(CHANGE, DHLIM)
       ENDIF
C--- Get new height
       HLOCAL=F(L0HI+I)+CHANGE
C--- Ensure height limits not exceeded
       F(L0HI+I)=AMAX1(HMIN,AMIN1(HLOCAL,FHLIM*F(L0HA+I)))
       VOLSUM=VOLSUM+F(L0AH+I)*(F(L0HI+I)-Z0)
19306  CONTINUE
C--- Apply block correction to preserve total volume of metal
       DH=(VLOWER-VOLSUM)/GAREAM
       DO 19308 I=1,NX*NY
       IF(SLD(I)) GO TO 19308
       F(L0HI+I)=F(L0HI+I)+DH
C--- Ensure correction does not violate height limits
       F(L0HI+I)=AMAX1(HMIN,AMIN1(F(L0HI+I),FHLIM*F(L0HA+I)))
19308  CONTINUE
C
       IF(LG(9).AND.GDBG) THEN
C--- Debug after adjustment
        WRITE(LUPR1,9986) DH,VOLSUM,FROUDE
 9986   FORMAT(1X,' DH ',1PE13.6,' VOL* ',1PE13.6,' FROUDE ',1PE13.6)
        WRITE(LUPR1,*) ' DHLIM= ',DHLIM,' FHLIM=',FHLIM,' HMIN=',HMIN
        CALL PRN('HI* ',GRSP2)
        if(lg(3))then
          write(lupr1,'(''pdif'',1p5e13.6)') (f(l0pdif+ii),ii=1,NX*NY)
          write(lupr1,'(''fzdf'',1p5e13.6)') (f(l0fzdif+ii),ii=1,NX*NY)
        endif
       ENDIF
C
C--- Adjust anode height to preserve anode-interface gap
       IF(LG(4)) CALL FN2(GRSP3,GRSP2,GZWNZ(IZ2)-GZWNZ(IZ1),1.0)
      ENDIF
C
C--- Adjust free surface height
C    --------------------------
C
      IF(LG(2)) THEN
       IF(MOD(ISWEEP,NIH).EQ.0.AND.IZ.EQ.IZ3.AND.ISWEEP.GT.IHF) THEN
        IF(IREF2.EQ.0) THEN
C--- Search for reference pressure point (not in blockage or anode)
         L0AH=L0F(AHIGH)
         CALL SUB2R(VUPPER,0.,GDRHO1,AGRAVZ*(RHOELC-RHOAIR))
         FREEAH=0.0
         DO 19310 I=1,NX*NY
         IF(.NOT.SLD(I)) GO TO 19312
19310    CONTINUE
19312    IREF2=I
C--- Compute and store total volume of electrolyte in inter-anode
C    and anode-wall gaps, and free surface area.
         L0HIH=L0F(GRSP4)
         L0HI=L0F(GRSP2)
         DO 19316 I=1,NX*NY
         IF(SLD(I)) GO TO 19316
         VUPPER=VUPPER+F(L0AH+I)*(F(L0HIH+I)-F(L0HI+I))
         FREEAH=FREEAH+F(L0AH+I)
19316    CONTINUE
C--- Get variables for later use
         CALL SUB2R(DHLIM1,FDH*GZWNZ(NZ),HMIN1,(2.-FHLIM)*GZWNZ(IZ2))
C         CALL SUB2R(DHLIM1,FDH*F(L0ZW+NZ),HMIN1,(2.-FHLIM)*F(L0ZW+IZ2))
        ENDIF
C
C--- Get new free-surface height from hydrostatic balance.
C
C--- Get required variables
        L0HIH=L0F(GRSP4)
        L0HI=L0F(GRSP2)
        L0P1=L0F(P1)
        CALL SUB3R(VUPSUM,0.,HREF,F(L0HIH+IREF2),DPREF,PAIR-
     1             F(L0P1+IREF2))
C
        IF(LG(9).AND.GDBG) THEN
C--- Debug for adjustment
         WRITE(LUPR1,9987) ISWEEP
 9987    FORMAT(1X,' FREE SURFACE ADJUSTMENT AT ISWEEP= ',I4)
         WRITE(LUPR1,9988) VUPPER,HREF,DPREF,GDRHO1
 9988    FORMAT(1X,' VOL ',1PE13.6,' HREF ',1PE13.6,' DPREF ',
     1    1PE13.6,' GDRHO1 ',1PE13.6)
         CALL PRN('P1  ',P1)
         CALL PRN('HIFS',GRSP4)
        ENDIF
C
C--- Perform adjustment over free surface only
        FROUDE=RHOELC*(GYDIST/DT)**2/(ABS(GDRHO1)*GZWNZ(IZ3))
C        FROUDE=RHOELC*(GYDIST/DT)**2/(ABS(GDRHO1)*F(L0ZW+IZ3))
        DO 19318 I=1,NX*NY
        IF(SLD(I)) GO TO 19318
        FAC=GDRHO1
        GRES=(HREF-F(L0HIH+I))*FAC-DPREF+PAIR-F(L0P1+I)
        CHANGE=GRES/(FAC*(1.0+FROUDE*FRCON))
C--- Relax change
        CHANGE=SLOH*CHANGE
        IF(CHANGE.LT.0.) THEN
         CHANGE=AMAX1(CHANGE,-DHLIM1)
        ELSE
         CHANGE=AMIN1(CHANGE, DHLIM1)
        ENDIF
C--- Get new height
        HLOCAL=F(L0HIH+I)+CHANGE
C--- Enusre height limits not exceeded
        F(L0HIH+I)=AMAX1(HMIN1,AMIN1(HLOCAL,(1.+FHLIM)*GZWNZ(NZ)))
        VUPSUM=VUPSUM+F(L0AH+I)*(F(L0HIH+I)-F(L0HI+I))
19318   CONTINUE
C--- Apply block correction to preserve total volume of electrolyte.
        DH=(VUPPER-VUPSUM)/FREEAH
        DO 19320 I=1,NX*NY
        IF(SLD(I)) GO TO 19320
        F(L0HIH+I)=F(L0HIH+I)+DH
C--- Ensure correction does not violate height limits
        F(L0HIH+I)=AMAX1(HMIN1,AMIN1(F(L0HIH+I),(1.+FHLIM)*GZWNZ(NZ)))
19320   CONTINUE
C
        IF(LG(9).AND.GDBG) THEN
C--- Debug after adjustment
         WRITE(LUPR1,9990) DH,VUPSUM,FROUDE
 9990    FORMAT(1X,' DH ',1PE13.6,' VOL* ',1PE13.6,' FROUDE ',1PE13.6)
         CALL PRN('HFS*',GRSP4)
        ENDIF
       ENDIF
      ENDIF
C
C--- Compute and apply grid distortion factors, stored as porosities
C    Storage locations are :
C     HI - 30 ; CPOR - 31 ; EPOR - 32 ; NPOR - 33
C
      SETPOR=.TRUE.
      CALL GETZ(ZWNZ,GZWNZ,NZM)
      IF(IZ0.GT.0) THEN
        HMET=GZWNZ(IZ1)-GZWNZ(IZ0)
        HCAT=GZWNZ(IZ0)
      ELSE
        HMET=GZWNZ(IZ1)
        HCAT=0.0
      ENDIF
      HELEC=GZWNZ(IZ2)-GZWNZ(IZ1)
      HANOD=GZWNZ(IZ3)-GZWNZ(IZ2)
      IF(IZ.LE.IZ0) THEN
C--- Cathode
       CALL FN1(CPOR,1.0)
       CALL FN1(HI,GZWNZ(IZ))
      ELSEIF(IZ.LE.IZ1) THEN
C--- Metal
       CALL FN2(CPOR,GRSP2,-HCAT/HMET,1.0/HMET)
C       CALL FN2(HI,CPOR,0.0,GZWNZ(IZ))
       CALL FN10(HI,CPOR,CPOR,HCAT,GZWNZ(IZ),-HCAT)
      ELSE IF(IZ.LE.IZ2) THEN
C--- Electrolyte
       CALL FN10(CPOR,GRSP3,GRSP2,0.0,1./HELEC,-1./HELEC)
       HEL=GZWNZ(IZ)-GZWNZ(IZ1)
       CALL FN10(HI,GRSP2,CPOR,0.0,1.0,HEL)
      ELSE IF(IZ.LE.IZ3) THEN
C--- Anodes immersed in electrolyte
       CALL FN10(CPOR,GRSP4,GRSP3,0.0,1./HANOD,-1./HANOD)
       HAN=GZWNZ(IZ)-GZWNZ(IZ2)
       CALL FN10(HI,GRSP3,CPOR,0.0,1.0,HAN)
      ELSE
C--- Anodes in air space
       CALL FN1(CPOR,1.0)
       CALL FN1(HI,GZWNZ(IZ))
      ENDIF
C--- Average cell-centre factors to cell faces
      CALL FNAVXY(EPOR,CPOR,0.0,0.5,0.5,X)
      CALL FNAVXY(NPOR,CPOR,0.0,0.5,0.5,Y)
C
      RETURN
  194 CONTINUE
C
C   * ------------------- SECTION 4 ---- START OF ITERATION.
C
C---  Compute currents and induced currents
      IF(LG(5).OR.ISWEEP.EQ.1) RETURN
      SETJ=ISWEEP.GT.1
C---  Only compute currents when solving for potential
      IV194=IV194 + 1
      IF(IV194.LE.NZ) RETURN
      IF(IV194.EQ.2*NZ) IV194=0
C---  Compute Jx,Jy and Jz from -sigma * grad(phi)
C     also Jix,Jiy and Jiz from sigma * (V x B) if LG(1)=T
C---  Get required variables
      L0EPOT=L0F(EPOT)
      IF(IZ.LT.NZ) L0EPOTH=L0F(HIGH(EPOT))
      L0JX=L0F(JX)
      L0PRPS=L0F(IPRPS)
      L0PRPH=L0F(HIGH(IPRPS))
      CALL CONDCT(L0CND,L0PRPS,L0PRPH)
      IF(IZ.LT.NZ) THEN
        L0PRPS=L0F(HIGH(IPRPS))
        IF(IZ.LE.NZ-1) THEN
          L0PRPH=L0F(ANYZ(IPRPS,IZ+2))
        ELSE
          L0PRPH=L0PRPS
        ENDIF
        CALL CONDCT(L0CNDH,L0PRPS,L0PRPH)
      ENDIF
      L0DEL=L0F(DXU2D)
      L0DELG=L0F(DXG2D)
      IF(LG(1)) THEN
C---  Variables required for induced currents only
        L0JIX=L0F(JIX)
        L0U1=L0F(U1)
        L0V1=L0F(V1)
        L0W1=L0F(W1)
        IF(IZ.GT.1) L0W1L=L0F(LOW(W1))
        L0BX=L0F(BX)
        L0BY=L0F(BY)
        L0BZ=L0F(BZ)
      ENDIF
C
      IF(LG(10).AND.GDBG) THEN
C--- Debug output
       WRITE(LUPR1,9991)  IZ,ISWEEP,IV194,INDVAR
 9991  FORMAT(1X,' IZ= ',I3,' ISWEEP= ',I5,' IV194= ',I4,' INDVAR= ',
     1    I4)
       CALL PRN('EPOT',EPOT)
       CALL PRN('EPTH',HIGH(EPOT))
       CALL PRN('COND',-L0CND)
       CALL PRN('CNDH',-L0CNDH)
       CALL PRN('DXU ',DXU2D)
       CALL PRN('DXG ',DXG2D)
      ENDIF
C
C--- Jx = -sigma * (Ee-Ep)/DXG
C--- Jix = sigma * (V * Bz - W * By)
      DO IY=1,NY
        I=(NX-1)*NY+IY
        F(L0CONDE+I+(IZ-1)*NX*NY)=0.0
      ENDDO
      DO 1940 IX=1,NXM1
      DO 1940 IY=1,NY
      I=(IX-1)*NY+IY
C--- Get conductivity
      GCON=2.*F(L0DELG+I)/(F(L0DEL+I)/F(L0CND+I)+F(L0DEL+I+NY)/
     1 F(L0CND+I+NY))
C--- Introduce BEMF at anode-electrolyte interface
      IF(NEZ(BEMF).AND.NINT(F(L0PRPS+I)).EQ.52.AND.NINT(F(L0PRPS+I+NY))
     1                                                     .EQ.100) THEN
       GCON=GCON*F(L0JX+I)*F(L0DELG+I)/(F(L0JX+I)*F(L0DELG+I)+GCON*BEMF)
      ENDIF
C--- save for use as link in EPOT equation
      F(L0CONDE+I+(IZ-1)*NX*NY)=GCON
C--- Get Jx
      F(L0JX+I)=-GCON*(F(L0EPOT+I+NY)-F(L0EPOT+I))/F(L0DELG+I)
      IF(LG(1)) THEN
C--- Get Jix
       IF(SLD(I)) THEN
        F(L0JIX+I)=0.0
       ELSE
C--- Average V1 to cell face
        IF(IY.GT.1) THEN
         GVB=0.25*(F(L0DEL+I+NY)*(F(L0V1+I)+F(L0V1+I-1))
     1        +F(L0DEL+I)*(F(L0V1+I+NY)+F(L0V1+I-1+NY)))/F(L0DELG+I)
        ELSE
         GVB=0.25*(F(L0DEL+I+NY)*F(L0V1+I)
     1        +F(L0DEL+I)*F(L0V1+I+NY))/F(L0DELG+I)
        ENDIF
C--- Average W1 to cell face
        IF(IZ.GT.1) THEN
         GWB=0.25*(F(L0DEL+I+NY)*(F(L0W1+I)+F(L0W1L+I))
     1         +F(L0DEL+I)*(F(L0W1+I+NY)+F(L0W1L+I+NY)))/F(L0DELG+I)
        ELSE
         GWB=0.25*(F(L0DEL+I+NY)*F(L0W1+I)
     1         +F(L0DEL+I)*F(L0W1+I+NY))/F(L0DELG+I)
        ENDIF
C--- Average Bz and By to cell face
        GBZB=0.5*(F(L0DEL+I+NY)*F(L0BZ+I)+F(L0DEL+I)*F(L0BZ+I+NY))/
     1   F(L0DELG+I)
        GBYB=0.5*(F(L0DEL+I+NY)*F(L0BY+I)+F(L0DEL+I)*F(L0BY+I+NY))/
     1   F(L0DELG+I)
C--- Compute Jix
        F(L0JIX+I)=GCON*(GVB*GBZB-GWB*GBYB)
       ENDIF
C--- Jx = Jx + Jix
       F(L0JX+I)=F(L0JX+I)+F(L0JIX+I)
      ENDIF
 1940 CONTINUE
C
      IF(LG(10).AND.GDBG) THEN
C---  Debug
       CALL PRN('JX  ',JX)
       IF(LG(1)) CALL PRN('JIX ',JIX)
      ENDIF
C
C--- Jy = -sigma * (En-Ep)/DYG
C--- Jiy = sigma * (W * Bx - U * Bz)
      L0DEL=L0F(DYV2D)
      L0DELG=L0F(DYG2D)
      L0JY=L0F(JY)
      IF(LG(1)) L0JIY=L0F(JIY)
      DO IX=1,NX
        I=(IX-1)*NY+NY
        F(L0CONDN+I+(IZ-1)*NX*NY)=0.0
      ENDDO
      DO 1942 IX=1,NX
      DO 1942 IY=1,NYM1
      I=(IX-1)*NY+IY
C--- Get conductivity
      GCON=2.*F(L0DELG+I)/(F(L0DEL+I)/F(L0CND+I)+F(L0DEL+I+1)/
     1 F(L0CND+I+1))
C--- Introduce BEMF at anode-electrolyte interface
      IF(NEZ(BEMF).AND.NINT(F(L0PRPS+I)).EQ.52.AND.NINT(F(L0PRPS+I+1))
     1                                                     .EQ.100) THEN
       GCON=GCON*F(L0JY+I)*F(L0DELG+I)/(F(L0JY+I)*F(L0DELG+I)+GCON*BEMF)
      ENDIF
C--- save for use as link in EPOT equation
      F(L0CONDN+I+(IZ-1)*NX*NY)=GCON
C--- Get Jy
      F(L0JY+I)=-GCON*(F(L0EPOT+I+1)-F(L0EPOT+I))/F(L0DELG+I)
      IF(LG(1)) THEN
C--- Get Jiy
       IF(SLD(I)) THEN
        F(L0JIY+I)=0.0
       ELSE
C--- Average U1 to cell face
        IF(IX.GT.1) THEN
         GUB=0.25*(F(L0DEL+I+1)*(F(L0U1+I)+F(L0U1+I-NY))
     1        +F(L0DEL+I)*(F(L0U1+I+1)+F(L0U1+I+1-NY)))/F(L0DELG+I)
        ELSE
         GUB=0.25*(F(L0DEL+I+1)*F(L0U1+I)
     1        +F(L0DEL+I)*F(L0U1+I+1))/F(L0DELG+I)
        ENDIF
C--- Average W1 to cell face
        IF(IZ.GT.1) THEN
         GWB=0.25*(F(L0DEL+I+1)*(F(L0W1+I)+F(L0W1L+I))
     1         +F(L0DEL+I)*(F(L0W1+I+1)+F(L0W1L+I+1)))/F(L0DELG+I)
        ELSE
         GWB=0.25*(F(L0DEL+I+1)*(F(L0W1+I)
     1         +F(L0DEL+I)*F(L0W1+I+1)))/F(L0DELG+I)
        ENDIF
C--- Average Bz & Bx to cell faces
        GBZB=0.5*(F(L0DEL+I+1)*F(L0BZ+I)+F(L0DEL+I)*F(L0BZ+I+1))/
     1   F(L0DELG+I)
        GBXB=0.5*(F(L0DEL+I+1)*F(L0BX+I)+F(L0DEL+I)*F(L0BX+I+1))/
     1   F(L0DELG+I)
C--- Compute Jiy
        F(L0JIY+I)=GCON*(GWB*GBXB-GUB*GBZB)
       ENDIF
C--- Jy = Jy + Jiy
       F(L0JY+I)=F(L0JY+I)+F(L0JIY+I)
      ENDIF
 1942 CONTINUE
C
      IF(LG(10).AND.GDBG) THEN
C---  Debug
       CALL PRN('DYV ',DYV2D)
       CALL PRN('DYG ',DYG2D)
       CALL PRN('JY  ',JY)
       IF(LG(1)) CALL PRN('JIY ',JIY)
      ENDIF
C
C--- Jz = -sigma * (Eh-Ep)/DZG
C--- Jiz = sigma * (U * By - V * Bx)
      IF(IZ.LT.NZ) THEN
      L0JZ=L0F(JZ)
      L0DEL=L0F(LXYDZ)
      L0DELG=L0F(LXYDZG)
      L0HI=L0F(HI)
      L0HIH=L0F(HIGH(HI))
      IF(LG(1)) THEN
        L0U1H=L0F(HIGH(U1))
        L0V1H=L0F(HIGH(V1))
        L0BXH=L0F(HIGH(BX))
        L0BYH=L0F(HIGH(BY))
        L0JIZ=L0F(JIZ)
      ENDIF
      DO 1944 IX=1,NX
      DO 1944 IY=1,NY
      I=(IX-1)*NY+IY
C--- Get conductivity
      GDELH=AMAX1( F(L0HIH+I)-F(L0HI+I), 1.E-10)
C
      GCON=2.*F(L0DELG+I)/(F(L0DEL+I)/F(L0CND+I)+GDELH/
     1 F(L0CNDH+I))
C--- Introduce BEMF at anode-electrolyte interface
      IF(NEZ(BEMF).AND.NINT(F(L0PRPS+I)).EQ.52.AND.NINT(F(L0PRPH+I))
     1                                                     .EQ.100) THEN
       GCON=GCON*F(L0JZ+I)*F(L0DELG+I)/(F(L0JZ+I)*F(L0DELG+I)+GCON*BEMF)
      ENDIF
C--- save for use as link in EPOT equation
      F(L0CONDH+I+(IZ-1)*NX*NY)=GCON
C--- Get Jz
      F(L0JZ+I)=-GCON*(F(L0EPOTH+I)-F(L0EPOT+I))/F(L0DELG+I)
      IF(LG(1)) THEN
C--- Get Jiz
       IF(SLD(I)) THEN
        F(L0JIZ+I)=0.0
       ELSE
C--- Average V1 to cell face
        IF(IY.GT.1) THEN
         GVB=0.25*(GDELH*(F(L0V1+I)+F(L0V1+I-1))
     1         +F(L0DEL+I)*(F(L0V1H+I)+F(L0V1H+I-1)))/F(L0DELG+I)
        ELSE
         GVB=0.25*(GDELH*F(L0V1+I)
     1         +F(L0DEL+I)*F(L0V1H+I))/F(L0DELG+I)
        ENDIF
C--- Average U1 to cell face
        IF(IX.GT.1) THEN
         GUB=0.25*(GDELH*(F(L0U1+I)+F(L0U1+I-NY))
     1         +F(L0DEL+I)*(F(L0U1H+I)+F(L0U1H+I-NY)))/F(L0DELG+I)
        ELSE
         GUB=0.25*(GDELH*F(L0U1+I)
     1         +F(L0DEL+I)*F(L0U1H+I))/F(L0DELG+I)
        ENDIF
C--- Average Bx & By to cell face
        GBXB=0.5*(GDELH*F(L0BX+I)+F(L0DEL+I)*F(L0BXH+I))/
     1   F(L0DELG+I)
        GBYB=0.5*(GDELH*F(L0BY+I)+F(L0DEL+I)*F(L0BYH+I))/
     1   F(L0DELG+I)
C--- Compute Jiz
        F(L0JIZ+I)=GCON*(GUB*GBYB-GVB*GBXB)
       ENDIF
C--- Jz = Jz + Jiz
       F(L0JZ+I)=F(L0JZ+I)+F(L0JIZ+I)
      ENDIF
 1944 CONTINUE
C
      IF(LG(10).AND.GDBG) THEN
C---  Debug
       CALL PRN('DZW ',LXYDZ)
       CALL PRN('DZG ',LXYDZG)
       CALL PRN('HI  ',HI)
       CALL PRN('HIH ',HIGH(HI))
       CALL PRN('JZ  ',JZ)
       IF(LG(1)) CALL PRN('JIZ ',JIZ)
      ENDIF
C
      ENDIF
      RETURN
C
196   CONTINUE
C
C   * ------------------- SECTION 6 ---- FINISH OF IZ SLAB
C
      RETURN
      END
C---------------------------------------------------------------------
C--- User routines start
C---------------------------------------------------------------------
      SUBROUTINE FNAVXY(K1,K2,A,B1,B2,IDIR)
C... this subroutine averages K2 into K1 in the IDIR direction
C    using  K1 = A + B1*K2p + B2*K2idir
      INCLUDE '/phoenics/d_includ/farray'
      COMMON /IGE/ IXF,IXL,IYF,IYL,IGFILL(21)
      COMMON /IDATA/ NX,NY,NZ, IDFILL(117)
      COMMON /NAMFN/ NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      SAVE
      NAMSUB='FNAVXY'
          IF(IXL.LT.0) RETURN
      CALL L0F2(K1,K2,I,I2M1,IADD,'FNAVXY')
          IF(IDIR.EQ.0) GO TO 3
        DO 2 IX=IXF,IXL
      I=I+IADD
        DO 2 IY=IYF,IYL
      I=I+1
          IF(IY.EQ.NY) GO TO 1
      F(I)=A+B1*F(I2M1+I)+B2*F(I2M1+I+1)
      GO TO 2
    1 F(I)=A+B1*F(I2M1+I)+B2*F(I2M1+I)
    2 CONTINUE
      RETURN
    3 CONTINUE
        DO 5 IX=IXF,IXL
      I=I+IADD
        DO 5 IY=IYF,IYL
      I=I+1
          IF(IX.EQ.NX) GO TO 4
      F(I)=A+B1*F(I2M1+I)+B2*F(I2M1+I+NY)
      GO TO 5
    4 F(I)=A+B1*F(I2M1+I)+B2*F(I2M1+I)
    5 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE FNSUM2(A,K1,K2)
C... This subroutine sums the arrays K1 and K2 into the REAL variable A
      INCLUDE '/phoenics/d_includ/farray'
      COMMON /IGE/ IXF,IXL,IYF,IYL,IGFILL(21)
      COMMON /NAMFN/ NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      SAVE
      NAMSUB='FNSUM2'
          IF(IXL.LT.0) RETURN
      CALL L0F2(K1,K2,I,I2M1,IADD,'FNSUM2')
        DO 2 IX=IXF,IXL
      I=I+IADD
        DO 2 IY=IYF,IYL
      I=I+1
      A=A+F(I)*F(I2M1+I)
    2 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE FNSUM1(A,K1)
C... This subroutine sums the array K1 into the REAL variable A
      INCLUDE '/phoenics/d_includ/farray'
      COMMON /IGE/ IXF,IXL,IYF,IYL,IGFILL(21)
      COMMON /NAMFN/ NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      SAVE
      NAMSUB='FNSUM1'
          IF(IXL.LT.0) RETURN
      CALL L0F1(K1,I,IADD,'FNSUM1')
        DO 2 IX=IXF,IXL
      I=I+IADD
        DO 2 IY=IYF,IYL
      I=I+1
      A=A+F(I)
    2 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE CONDCT(L0CON,L0PRP,L0PRPH)
C--- Subroutine for setting electric conductivities
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
C
C 1   Set dimensions of satellite-to-GROUND data arrays to those
C     of the satellite.
      COMMON/LGRND/LG(20)/IGRND/IG(20)/RGRND/RG(100)/CGRND/CG(10)
      LOGICAL LG,NEZ
      CHARACTER*4 CG
      COMMON /NAMFN/ NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      EQUIVALENCE (RHOMET,RG(1)), (RHOELC,RG(2)), (RHOANO,RG(3)),
     1            (CONMET,RG(4)), (CONELC,RG(5)), (CONANO,RG(6)),
     1            (RHOFRZ,RG(20)),(CONFRZ,RG(21)),(RHOAIR,RG(23)),
     1            (CONAIR,RG(24)),(CONFAC,RG(29)),(BEMF,RG(30)),
     1            (CONCAT,RG(31)),(RHOCAT,RG(32))
C--- Set conductivity according to location - DEN1 used as flag
C
      NAMSUB='CONDCT'
      DO IX=1,NX
        DO IY=1,NY
          I=(IX-1)*NY+IY
          IPROP=NINT(F(L0PRP+I))
          IF(IPROP.EQ.51) THEN
C--- Metal
            F(L0CON+I)=CONMET
          ELSEIF(IPROP.EQ.52) THEN
C--- Bath
            F(L0CON+I)=CONELC
            IF(NEZ(CONFAC)) THEN
              IF((IY.LT.NY.AND.NINT(F(L0PRP+I+ 1)).EQ.100).OR.
     1           (IY.GT. 1.AND.NINT(F(L0PRP+I- 1)).EQ.100).OR.
     1           (IX.LT.NX.AND.NINT(F(L0PRP+I+NY)).EQ.100).OR.
     1           (IX.GT. 1.AND.NINT(F(L0PRP+I-NY)).EQ.100).OR.
     1            NINT(F(L0PRPH+I)).EQ.100) F(L0CON+I)=F(L0CON+I)*CONFAC
            ENDIF
          ELSEIF(IPROP.EQ.151) THEN
C--- Cathode
            F(L0CON+I)=CONCAT
          ELSEIF(IPROP.EQ.100) THEN
C--- Anode
            F(L0CON+I)=CONANO
          ELSEIF(NINT(F(L0PRP+I)).EQ.198.OR.NINT(F(L0PRP+I)).EQ.199)THEN
            F(L0CON+I)=0.0
          ENDIF
        ENDDO
      ENDDO
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE POTENT (GPOT,NYY,NXX)
C--- Example subroutine for setting anode potentials
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
C 1   Set dimensions of satellite-to-GROUND data arrays to those
C     of the satellite.
      COMMON/LGRND/LG(20)/IGRND/IG(20)/RGRND/RG(100)/CGRND/CG(10)
      LOGICAL LG
      LOGICAL QEQ
      CHARACTER*4 CG
      COMMON /NAMFN/ NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      DIMENSION GPOT(NYY,NXX)
      EQUIVALENCE (RHOMET,RG(1)),(RHOELC,RG(2)),(RHOANO,RG(3)),
     1            (CONMET,RG(4)),(CONELC,RG(5)),(CONANO,RG(6)),
     1            (RHOFRZ,RG(20)),(ANOPOT,RG(22))
C
C--- Set potential according to location - DEN1 used as flag
C
      NAMSUB='POTENT'
      L0D=L0F(DEN1)
      DO 100 IX = 1, NX
        DO 100 IY = 1, NY
          I=(IX-1)*NY+IY
          IF (QEQ(F(L0D+I),RHOANO)) THEN
C--- Anode
            GPOT(IY,IX) = ANOPOT
          ELSE
            GPOT(IY,IX)=0.0
          ENDIF
  100 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DSCFLS (LUNIT, FNAME, IGO,ID2)
C
C... Subroutine to OPEN files at run-time. May be machine dependent.
C
      COMMON /IDATA/ IDUM1(3), LUPR1, IDUM2(116)
      CHARACTER *(*) FNAME
C
      IF(IGO.EQ.12) THEN
C... Open existing file, STATUS = 'OLD'
        OPEN(LUNIT, FILE = FNAME, STATUS = 'OLD', ERR=100, IOSTAT=IOS)
        REWIND LUNIT
      ELSE IF(IGO.EQ.14) THEN
        CLOSE(LUNIT,STATUS='KEEP')
      ENDIF
      GO TO 999
C
C... Error in opening file
  100 WRITE(LUPR1,*) 'Error number ',IOS,' in opening file :',FNAME
      CALL SET_ERR(322, 'Error in opening file.',2)
C      CALL EAROUT(2)
  999 RETURN
      END
c