c C name estrgr.for................................... 190219 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
      USE intercept_data  !<<------------- Include module for gridline intercepts
      INCLUDE 'objnam'
      INCLUDE 'farray'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE 'lbnamer'
      INCLUDE 'd_earth/clstr' ! sparsol common block
C
C... Include some EARTH commons to get no. of variables and
C    storage information
      COMMON /VDAT/VERDAT,VER
      CHARACTER VERDAT*27,VER*5
      COMMON/F0/KXDX,KXXG,KXXU,IF01(2),KYDY,KYYG,IF02(6),KYYV,IF03(2),
     1          KZDZ,KZZG,IF04(6),KZZW,IF05(280)
      COMMON /GENL/ LDUM1(20), INCORJ,INCORM, LDUM2(38)
      LOGICAL LDUM1,LDUM2,INCORJ,INCORM
      COMMON /NAMFN/ NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB,LINE*256
      CHARACTER*68 FNAME,BUFF*256
C
      COMMON /GEODMN0/ I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
     1                I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
      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
     1       /FACES/L0FACE,L0FACZ
      COMMON/DRHODP/ITEMP,IDEN/DVMOD/IDVCGR
      COMMON /STOGEO/STORGEO
      LOGICAL STORGEO
      COMMON /WORDI1/ NWDS,NCHARS(20),NSEMI,H,NLINES
      INTEGER H
      COMMON /WORDL1/ ERROR,MORWDS
      LOGICAL ERROR,MORWDS
      COMMON /WORDC1/ WD(20),INLINE
      CHARACTER*20 WD, INLINE*120
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,HIFREQ
      LOGICAL GDBG,SLD,LSOLID,VAC,LVAC,POR,LPOR,GRN
      LOGICAL QEQ,QNE,NEZ,LABOVE,LDOHI,NEWREAD,ISRLZZ
      REAL JZCATH
 
      CHARACTER*256 FORM
C
C----Following array dimensions are NYM,NXM
C--- except  GZWNZ(NZM)
      REAL, ALLOCATABLE :: GYV(:,:),GBX(:,:),GBY(:,:),GBZ(:,:),
     1          GVAL(:,:),GJCATH(:,:),GZWNZ(:),HANODE(:,:)
C
      INTEGER, ALLOCATABLE :: IANODE(:,:),IFREEZ(:,:)
 
      REAL, ALLOCATABLE :: XCD(:),YCD(:),ZCD(:)
      REAL, ALLOCATABLE :: BXD(:,:,:),BYD(:,:,:),BZD(:,:,:)
      REAL, ALLOCATABLE :: BXC(:,:,:),BYC(:,:,:),BZC(:,:,:)
      REAL, ALLOCATABLE :: BXDY(:,:,:),BYDY(:,:,:),BZDY(:,:,:),
     1                     BXDX(:,:,:),BYDX(:,:,:),BZDX(:,:,:)
      REAL, ALLOCATABLE :: JZC(:,:),JZCY(:,:)
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 (RHOCOL,RG(33)),(CONCOL,RG(34)),
     1            (RHOROD,RG(35)),(CONROD,RG(36)),
     1            (RHOBUS,RG(37)),(CONBUS,RG(38)),
     1            (RHOSTB,RG(71)),(CONSTB,RG(72)),
     1            (RHOCOP,RG(73)),(CONCOP,RG(74)),
     1            (RHCST1,RG(75)),(CNCST1,RG(76)),
     1            (RHCST2,RG(77)),(CNCST2,RG(78))
      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 IV194 /0/
      DATA IREF,IREF2 /0, 0/
      DATA IERR1,IERR2,IERR3,IERR4,IERR5,IERR6,IERR7,IERR8 /8*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==13) GO TO 13
      IF(IGR==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)/=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   :   '//VER)
      CALL WRYT40('GROUND file is ESTRGR.HTM of:    190219 ')
      CALL WRIT40('PHOENICS version number is   :   '//VER)
      CALL WRIT40('GROUND file is ESTRGR.HTM of:    190219 ')
C
      IF(.NOT.ALLOCATED(GYV)) THEN
        ALLOCATE(GYV(NY,NX),GBX(NY,NX),GBY(NY,NX),GBZ(NY,NX),
     1   GVAL(NY,NX),GJCATH(NY,NX),GZWNZ(NZ),STAT=IERR)
        ALLOCATE(IANODE(NY,NX),HANODE(NY,NX),IFREEZ(NY,NX),STAT=IERR)
        IANODE=0; HANODE=0.0; IFREEZ=0
      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
C--- Ensure that UDIFF is set when BEMF required or two phase
      IF(NEZ(BEMF).OR..NOT.ONEPHS) UDIFF=.TRUE.
C
C--- Stores for conductivity
      CALL GXMAKE(L0CND,NX*NY,'COND')
      CALL GXMAKE(L0CNDH,NX*NY,'CNDH')
      IF(UDIFF) THEN
        CALL GXMAKE(L0CONDE,NX*NY*NZ,'CONDE')
        CALL GXMAKE(L0CONDN,NX*NY*NZ,'CONDN')
        CALL GXMAKE(L0CONDH,NX*NY*NZ,'CONDH')
      ENDIF
C
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
C
C... activate 2D stores for Z distances
      AXDZ=GRND
C
C--- Default location for air surface
      IF(IZ3==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'); IHINT=INAME('HINT')
      IF(EPOT>0) THEN
        LBEPOT=EPOT; ISEPOT=.TRUE.
      ENDIF
C
C--- Set atmospheric conditions for free surface adjustment
      PAIR=0.0
C
C--- Ensure specific heat is set
      IF(IPRP==0.AND..NOT.GRN(CP1)) CP1=GRND
 
      FNAME='NO'; CALL GETSDC('ESTER','HI-FILE',FNAME)
      LDOHI=FNAME/='NO'
      IF(LDOHI) THEN
        HIFREQ=0; CALL GETSDI('ESTER','HI-FREQ',HIFREQ)
        IF(HIFREQ==0) HIFREQ=LSWEEP
      ENDIF
      FNAME='NO'; CALL GETSDC('ESTER','CBARJ-FILE',FNAME)
      DOCBJ=FNAME/='NO'
      FNAME='NO'; CALL GETSDC('ESTER','RODSJ-FILE',FNAME)
      DORODJ=FNAME/='NO'
      FNAME='NO'; CALL GETSDC('ESTER','RODSV-FILE',FNAME)
      DORODV=FNAME/='NO'
      RETURN
C
 1002 CONTINUE
      NAMSUB='G1.2  '
C
      RETURN
C*****************************************************************
C
C--- GROUP 5. Z-direction grid specification
C
    5 CONTINUE
C--- Recompute LXYDZ,LXYDZG,LXYDZL
      IF(IZ==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=IZ0.AND.IZ<=IZ2) THEN
C.... Adjust cut-link stores of Z-distances for Sparsol
        L0XYDZ=L0F(LXYDZ)
        L0PRP=L0F(IPRPS); L0PRPH=L0F(HIGH(IPRPS))
        L0HI=L0F(HI);     L0HIH=L0F(HIGH(HI))
        L0DX=L0F(LXYDX);  L0DY=L0F(LXYDY)
        DO ICLSORT=1,NCUTLNK
          ICL=ICUTLIST(ICLSORT) ! in order to increase IZ
          IZZ=CLIZ(ICLSORT)     ! iz index of current cut-link
          IF(IZZ/=IZ) CYCLE
          ICLTYP=ICL_type(ICL,.TRUE.,0) ! get type of cut-link
          IF(ICLTYP==1.OR.ICLTYP==5) CYCLE ! cut-link through solids o,
           ! or link is cut by fluid objects, link is not considered
          IJ=ICL_ij(ICL,.TRUE.,0) ! get index ij
          IDIR=ICL_idir(ICL,.TRUE.,0) ! Get direction of current cut-link
          LABOVE=ICLTYP/=3 ! object above intersection
          IF(IDIR==5) THEN ! High direction
            DISTW=0.5*F(L0XYDZ+IJ)               ! half local slab depth
!            CGtoI=RCLdstvl(ICL,.FALSE.,DISTW)    ! distance between intersection and
                                                 ! cell node of off-line velocity component
            IF(LABOVE) THEN ! solids, cutting link, above intersection
              DISTF=DISTW ! half current slab depth
              V=RCLdstdw(ICL,.FALSE.,DISTF)      ! distance in fluid between interface
                                                 ! and current cell center
              CONDS=CNDCT(IJ,L0PRP+NFM,L0PRPH+NFM) ! conductivity of solid
              DISTS=0.5*(F(L0HIH+IJ)-F(L0HI+IJ)) ! distance in solid
            ELSE ! solids, cutting link, below intersection
              DISTF=0.5*(F(L0HIH+IJ)-F(L0HI+IJ)) ! half upper slab depth
              V=RCLdstup(ICL,.FALSE.,DISTF)      ! distance in fluid between interface
                                                 ! and neighboring cell center
              CONDS=CNDCT(IJ,L0PRP,L0PRPH)       ! conductivity of solid
              DISTS=DISTW                        ! distance in solid
            ENDIF
            CGtoI=RCLdstvl(ICL,.FALSE.,DISTF)    ! distance between intersection and
!          ELSEIF(IDIR==1) THEN ! North
!            IF(LABOVE) THEN
!              CONDS=CNDCT(IJ+1,L0PRP,L0PRPH)
!              DISTS=0.5*F(L0DY+IJ+1)
!            ELSE
!              CONDS=CNDCT(IJ,L0PRP,L0PRPH)
!              DISTS=0.5*F(L0DY+IJ)
!            ENDIF
!          ELSEIF(IDIR==3) THEN ! East
!            IF(LABOVE) THEN
!              CONDS=CNDCT(IJ+NY,L0PRP,L0PRPH)
!              DISTS=0.5*F(L0DX+IJ+NY)
!            ELSE
!              CONDS=CNDCT(IJ,L0PRP,L0PRPH)
!              DISTS=0.5*F(L0DX+IJ)
!            ENDIF
            RESIST=DISTS/CONDS
            V=RCLtres(ICL,.FALSE.,RESIST) ! save resistance of solid
          ENDIF
        ENDDO
      ENDIF
      IF(STORGEO) THEN
        CALL FN026(VOL,AHIGH,LXYDZ) ! Vol = Ahigh * DZ
      ENDIF
      RETURN
C*****************************************************************
C
C--- GROUP 8. Terms (in differential equations) & devices
C
    8 IF(ISC==5.OR.ISC==6) THEN
C   * -----------GROUP 8  SECTION 5 & 6 --- ADDITIONAL VELOCITY
C--- for W1AD<=GRND--- phase 1 additional velocity (VELAD).
C--- for W2AD<=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==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==R1.OR.INDVAR==R2) THEN
C.... Cut diffusive links in R1/R2 equations in metal
          IF(IZ<=IZ1) THEN
            IF(NDIREC==1) THEN
              LVAL=LAN
            ELSEIF(NDIREC==3) THEN
              LVAL=LAE
            ELSEIF(NDIREC==5) THEN
              LVAL=LD11
            ENDIF
            CALL FN1(LVAL,0.0)
          ENDIF
        ELSEIF(INDVAR==EPOT.AND.(ISWEEP>1.OR.ISTEP>1)) THEN
C... Recompute diffusive links for EPOT, without volume fraction multiplier
C    Also recalculate areas to avoid Sparsol modifications
          IF(IZ>IZ3.OR.IZ0) l0ar=l0f(lbar)
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
            if(ndirec==5) then
              AREA=F(L0D1+I)*F(L0D2+I) ! DY*DZ or DY*DZ or DX*DY
            else
              AREA=F(L0D1+I)*0.5*(F(L0D2+I)+F(L0D2+I+IADD)) ! DY*DZ or DY*DZ or DX*DY
            endif
            if(l0ar>0) f(l0ar+i)=area
            F(L0DIF+I)=F(L0COND+I+(IZ-1)*NX*NY)*AREA/F(L0DEL+I)
            IF(IADD<=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
      IF(QEQ(CP1,GRND)) CALL FN1(ISPH1,1.0)
      RETURN
   91 CONTINUE
      NAMSUB='G9.1  '
C   * -----------GROUP 9  SECTION  1 ---------------------------
C--- for RHO1<=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
      IF(QEQ(RHO1,GRND)) THEN
        L0DEN=L0F(DEN1)
        L0PRP=L0F(IPRPS)
        DO 910 I=1,NX*NY
          IPRP=NINT(F(L0PRP+I))
          IF(IPRP==51.OR.IPRP==-1) THEN
C--- Metal
          F(L0DEN+I)=RHOMET
        ELSEIF(IPRP==52) THEN
C--- Bath
          F(L0DEN+I)=RHOELC
        ELSEIF(IPRP==151.OR.IPRP==104) THEN
C--- Carbon block
          F(L0DEN+I)=RHOCAT
        ELSEIF(IPRP==100) THEN
C--- Anode
          F(L0DEN+I)=RHOANO
        ELSEIF(IPRP==198.OR.IPRP==199) THEN
          F(L0DEN+I)=0.0
        ENDIF
  910   CONTINUE
      ENDIF
      RETURN
   92 CONTINUE
      NAMSUB='G9.2  '
C   * -----------GROUP 9  SECTION  2 ---------------------------
C--- for DRH1DP<=GRND--- D(LN(DEN))/DP for phase 1 (D1DP).
      RETURN
C
   95 CONTINUE
      NAMSUB='G9.5  '
C   * -----------GROUP 9  SECTION  5 ---------------------------
C    For ENUT<=GRND--- reference turbulent kinematic viscosity
C                                                  Index VIST
      IF(QEQ(ENUT,GRND)) THEN
C... Different constant turbulent viscosity in metal and bath
        IF(IZ<=IZ1)THEN
C ... Metal
          CALL FN1(VIST,ENUTA)
        ELSE IF (IZ<=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<=GRND--- reference laminar kinematic viscosity.
C
      IF(IZ<=IZ0) THEN
C--- Cathode - carbon block
       CALL FN1(VISL,0.0)
C--- Metal
      ELSEIF(IZ<=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( )<=GRND--- laminar PRANDTL nos., or diffusivity.
!      IF(QEQ(PRNDTL(EPOT),-GRND).AND.INDVAR==EPOT) THEN
      IF(INDVAR==EPOT) THEN
        IF(QEQ(PRNDTL(EPOT),-GRND)) THEN
C--- Get conductivities for current Z-slab
          L0PRPS=L0F(IPRPS)
          IF(IZ=MIN(BX,BY,BZ).AND.INDVAR<=MAX(BX,BY,BZ)) THEN
C--- Initialise Bx, By and Bz
        IF(MAGF==0) THEN
          IF(IMAGF==0) THEN
C--- Calculate Bx By & Bz from algebraic expressions (PDR 20 case)
            L0XU=L0F(XU2D); L0YV=L0F(YV2D); L0VAL=L0F(VAL)
            IF(INDVAR==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==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==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==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==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==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
        ELSEIF(MAGF>0) THEN
C--- Read Bx By & Bz from disc file
          IF(INDVAR==MIN(BX,BY,BZ)) THEN ! first B field variable
            IF(IZ==1) THEN
              FNAME=CG(1); CALL GETSDC('ESTER','MAGF',FNAME)
              LL=LENGZZ(FNAME); CALL LOWCSE(FNAME,-LL) ! make lowercase
 1102         CONTINUE
              CALL DSCFLS(MAGF,FNAME,12,0)
              IF(FNAME(1:2)=='q1') THEN
                DO WHILE(LINE/='  *MAGF')
                  READ(MAGF,'(A)',ERR=1104,END=1104) LINE
                ENDDO
                GO TO 1105
 1104           CONTINUE
                IF(FNAME=='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)
              ENDIF
 1105         CONTINUE ! established file name
              LINE=' '
              DO WHILE (LINE==' ')
                READ(MAGF,'(A)',ERR=1108,IOSTAT=IOS) LINE
                ICOM=INDEX(LINE,'!')       ! check for comments
                IF(ICOM>0) LINE(ICOM:)=' ' ! delete any comment
              ENDDO
              LL=LENGZZ(LINE)
              CALL SPLTZZZ (LINE,WD,NWDS,NCHARS,LL,' ,;='//CHAR(9),20) ! split into words
              IF(WD(1)=='TITLE') THEN ! New style
                NEWREAD=.TRUE.; IREAD=0
!                DO WHILE(WD(1)/='VARIABLES')
!                  READ(MAGF,'(A)',ERR=1108,IOSTAT=IOS) LINE
!                  ICOM=INDEX(LINE,'!')       ! check for comments
!                  IF(ICOM>0) LINE(ICOM:)=' ' ! delete any comment
!                  LL=LENGZZ(LINE)
!                  CALL SPLTZZZ (LINE,WD,NWDS,NCHARS,LL,' ,;='//CHAR(9),20) ! split into words
!                ENDDO
C... Get grid size of data from ZONE
                NXD=0; NYD=0; NZD=0
                DO WHILE(WD(1)/='ZONE')
                  READ(MAGF,'(A)',ERR=1108,IOSTAT=IOS) LINE
                  ICOM=INDEX(LINE,'!'); IF(ICOM>0) LINE(ICOM:)=' '
                  IF(LINE(1:1)=='#') LINE=' '
                  LL=LENGZZ(LINE)
                  CALL SPLTZZZ (LINE,WD,NWDS,NCHARS,LL,' ,;='//CHAR(9),
     1                                                               20) ! split into words
                ENDDO
                IF(WD(1)=='ZONE') THEN
                  DO IWD=2,NWDS
                    IF(WD(IWD)=='I') THEN
                      NXD=IRDZZZ(IWD+1); IF(ERROR) GO TO 1108
                    ENDIF
                    IF(WD(IWD)=='J') THEN
                      NYD=IRDZZZ(IWD+1); IF(ERROR) GO TO 1108
                    ENDIF
                    IF(WD(IWD)=='K') THEN
                      NZD=IRDZZZ(IWD+1); IF(ERROR) GO TO 1108
                    ENDIF
                  ENDDO
                ENDIF
                IF(NXD>0.AND.NYD>0.AND.NZD>0) THEN
C... Allocate storage for data
                  ALLOCATE(XCD(NXD),YCD(NYD),ZCD(NZD),BXD(NXD,NYD,NZD),
     1                       BYD(NXD,NYD,NZD),BZD(NXD,NYD,NZD),STAT=IOS)
                  IF(IOS/=0) GO TO 1108
C... Read the data file into XCD(1D), YCD(1D), ZCD(1D), BXD(3D), BYD(3D) & BZD(3D)
C... Allow for blank lines and comments after !. Ignore lines starting #
                  LINE=' '
                  DO IZD=1,NZD; DO IYD=1,NYD; DO IXD=1,NXD
                    DO WHILE (LINE==' ')
                      READ(MAGF,'(A)',ERR=1108,IOSTAT=IOS) LINE
                      ICOM=INDEX(LINE,'!'); IF(ICOM>0) LINE(ICOM:)=' '
                      IF(LINE(1:1)=='#') LINE=' '
                    ENDDO
                    READ(LINE,*,ERR=1108,IOSTAT=IOS)
     1                      XCD(IXD),YCD(IYD),ZCD(IZD),BXD(IXD,IYD,IZD),
     1                              BYD(IXD,IYD,IZD),BZD(IXD,IYD,IZD),BM
                    LINE=' '
                  ENDDO; ENDDO; ENDDO
                  XOFF=0.5*XULAST; YOFF=0.5*YVLAST ! recentre on domain
                  XCD=XCD+XOFF; YCD=YCD+YOFF
                  IF(IZ0/=0) THEN ! offset by height of cathode block
                    ZOFF = F(KZZW+IZ0) ! Height of cathode carbon block
                    ZCD=ZCD+ZOFF
                  ENDIF
C... Interpolate in Y in BXD,BYD,BZD into BXDY,BYDY,BZDY
                  ALLOCATE(BXDY(NXD,NY,NZD),BYDY(NXD,NY,NZD),
     1                                        BZDY(NXD,NY,NZD),STAT=IOS)
                  IF(IOS/=0) GO TO 1108
                  DO IZD=1,NZD
                    DO IXD=1,NXD ! loop over all X in data
                      DO IY=1,NY ! loop over actual Y in solution
                        YLOC=F(KYYG+IY) ! actual Y cell-centre
                        IF(YLOC<=YCD(1)) THEN ! below start of data, take first point
                          BXDY(IXD,IY,IZD)=BXD(IXD,1,IZD)
                          BYDY(IXD,IY,IZD)=BYD(IXD,1,IZD)
                          BZDY(IXD,IY,IZD)=BZD(IXD,1,IZD)
                        ELSEIF(YLOC>=YCD(NYD)) THEN ! above end of data, take last point
                          BXDY(IXD,IY,IZD)=BXDY(IXD,NYD,IZD)
                          BYDY(IXD,IY,IZD)=BYDY(IXD,NYD,IZD)
                          BZDY(IXD,IY,IZD)=BZDY(IXD,NYD,IZD)
                        ELSE ! interpolate
                          DO IYD=2,NYD ! loop over all Y in data
                            IF(YLOC<=YCD(IYD)) THEN ! found interval in data
                              FC=(YLOC-YCD(IYD-1))/(YCD(IYD)-YCD(IYD-1))
                              BXDY(IXD,IY,IZD)=BXD(IXD,IYD-1,IZD)+
     1                          (BXD(IXD,IYD,IZD)-BXD(IXD,IYD-1,IZD))*FC
                              BYDY(IXD,IY,IZD)=BYD(IXD,IYD-1,IZD)+
     1                          (BYD(IXD,IYD,IZD)-BYD(IXD,IYD-1,IZD))*FC
                              BZDY(IXD,IY,IZD)=BZD(IXD,IYD-1,IZD)+
     1                          (BZD(IXD,IYD,IZD)-BZD(IXD,IYD-1,IZD))*FC
                             EXIT
                           ENDIF
                          ENDDO
                        ENDIF
                      ENDDO ! end IY loop
                    ENDDO   ! end IXD loop
                  ENDDO     ! end IZD loop
C... interpolate in X in BXDY,BYDY,BZDY into BXDX,BYDY,BZDZ
                  ALLOCATE(BXDX(NX,NY,NZD),BYDX(NX,NY,NZD),
     1                                        BZDX(NX,NY,NZD),STAT=IOS)
                  IF(IOS/=0) GO TO 1108
                  DO IZD=1,NZD
                    DO IY=1,NY ! loop over actual Y
                      DO IX=1,NX ! loop over actual X
                        XLOC=F(KXXG+IX) ! actual X cell centre
                        IF(XLOC<=XCD(1)) THEN ! below start of data, take first
                          BXDX(IX,IY,IZD)=BXDY(1,IY,IZD)
                          BYDX(IX,IY,IZD)=BYDY(1,IY,IZD)
                          BZDX(IX,IY,IZD)=BZDY(1,IY,IZD)
                        ELSEIF(XLOC>=XCD(NXD)) THEN ! above end of data , take last
                          BXDX(IX,IY,IZD)=BXDY(NXD,IY,IZD)
                          BYDX(IX,IY,IZD)=BYDY(NXD,IY,IZD)
                          BZDX(IX,IY,IZD)=BZDY(NXD,IY,IZD)
                        ELSE ! interpolate
                          DO IXD=2,NXD ! loop over X in data
                            IF(XLOC<=XCD(IXD)) THEN ! found interval
                              FC=(XLOC-XCD(IXD-1))/(XCD(IXD)-XCD(IXD-1))
                              BXDX(IX,IY,IZD)=BXDY(IXD-1,IY,IZD)+
     1                          (BXDY(IXD,IY,IZD)-BXDY(IXD-1,IY,IZD))*FC
                              BYDX(IX,IY,IZD)=BYDY(IXD-1,IY,IZD)+
     1                          (BYDY(IXD,IY,IZD)-BYDY(IXD-1,IY,IZD))*FC
                              BZDX(IX,IY,IZD)=BZDY(IXD-1,IY,IZD)+
     1                          (BZDY(IXD,IY,IZD)-BZDY(IXD-1,IY,IZD))*FC
                              EXIT
                            ENDIF
                          ENDDO
                        ENDIF
                      ENDDO ! end IX loop
                    ENDDO   ! end IY loop
                  ENDDO     ! end IZD loop
C... Finally interpolate in Z in BXDX,BYDX,BZDX into BXC,BYC,BZC
                  ALLOCATE(BXC(NX,NY,NZ),BYC(NX,NY,NZ),
     1                                        BZC(NX,NY,NZ),STAT=IOS)
                  IF(IOS/=0) GO TO 1108
                  IZSTART=ITWO(1,IZ0+1,IZ0==0) ! skip slabs in cathode
                  BXC=0.0; BYC=0.0; BZC=0.0    ! initialise all slabs to zero
                  DO IZC=IZSTART,NZ ! loop over actual Z
                    ZLOC=F(KZZG+IZC) ! actual Z cell centre
                    DO IY=1,NY ! loop over actual Y
                      DO IX=1,NX ! loop over actual X
                        IF(ZLOC<=ZCD(1)) THEN ! below start of data, take first
                          BXC(IX,IY,IZC)=BXDX(1,IY,IZD)
                          BYC(IX,IY,IZC)=BYDX(1,IY,IZD)
                          BZC(IX,IY,IZC)=BZDX(1,IY,IZD)
                        ELSEIF(ZLOC>=ZCD(NZD)) THEN ! above end of data , take last
                          BXC(IX,IY,IZC)=BXDX(NXD,IY,IZD)
                          BYC(IX,IY,IZC)=BYDX(NXD,IY,IZD)
                          BZC(IX,IY,IZC)=BZDX(NXD,IY,IZD)
                        ELSE ! interpolate
                          DO IZD=2,NZD ! loop over X in data
                            IF(ZLOC<=ZCD(IZD)) THEN ! found interval
                              FC=(ZLOC-ZCD(IZD-1))/(ZCD(IZD)-ZCD(IZD-1))
                              BXC(IX,IY,IZC)=BXDX(IX,IY,IZD-1)+
     1                            (BXDX(IX,IY,IZD)-BXDX(IX,IY,IZD-1))*FC
                              BYC(IX,IY,IZC)=BYDX(IX,IY,IZD-1)+
     1                            (BYDX(IX,IY,IZD)-BYDX(IX,IY,IZD-1))*FC
                              BZC(IX,IY,IZC)=BZDX(IX,IY,IZD-1)+
     1                            (BZDX(IX,IY,IZD)-BZDX(IX,IY,IZD-1))*FC
                              EXIT
                            ENDIF
                          ENDDO
                        ENDIF
                      ENDDO ! end IX loop
                    ENDDO   ! end IY loop
                  ENDDO     ! end IZ loop
                  DEALLOCATE(XCD,YCD,ZCD,BXD,BYD,BZD,BXDY,BYDY,BZDY,
     1                                          BXDX,BYDX,BZDX,STAT=IOS)
                  IF(IOS/=0) GO TO 1108
                  WRITE(LUPR1,
     1       '(/'' MAGNETIC FIELDS READ AND INTERPOLATED FROM FILE '',
     2                A/)') FNAME(1:LENGZZ(FNAME))
                ELSE ! didn't find grid sizes
                  IREAD=1; GO TO 1108
                ENDIF
              ELSE ! old-style file, read in old way
                BACKSPACE(MAGF); NEWREAD=.FALSE.
              ENDIF
            ENDIF ! end IZ==1
          ENDIF ! end of actions for first B field variable
          IF(.NOT.NEWREAD) THEN ! Read B field old way
            IF(INDVAR==MIN(BX,BY,BZ)) THEN
              IREAD=0
              DO ICOMPON = 1,3
                DO IY =NY, 1, -1
                  IF(ICOMPON == 1)THEN
                    READ(MAGF,*,END=1110,ERR=1108,IOSTAT=IOS)
     1                                              (GBX(IY,IX),IX=1,NX)
                  ELSE IF(ICOMPON == 2)THEN
                    READ(MAGF,*,END=1110,ERR=1108,IOSTAT=IOS)
     1                                              (GBY(IY,IX),IX=1,NX)
                  ELSE
                    READ(MAGF,*,END=1110,ERR=1108,IOSTAT=IOS)
     1                                              (GBZ(IY,IX),IX=1,NX)
                  END IF
                  IREAD=IREAD+NX
                ENDDO
              ENDDO
              IF(MOD(IREAD,NXNY)==0) THEN
                WRITE(LUPR1,9980) FNAME(1:LENGZZ(FNAME)),IZ
 9980           FORMAT(/' MAGNETIC FIELDS READ FROM FILE ',A,' AT IZ = ',
     1                                                              I3/)
              ENDIF
            ENDIF
 1110       CONTINUE ! jump here if end-of-file found, i.e. no data for this slab
            IF(INDVAR==BX) CALL SETYX(VAL,GBX,NY,NX)
            IF(INDVAR==BY) CALL SETYX(VAL,GBY,NY,NX)
            IF(INDVAR==BZ) CALL SETYX(VAL,GBZ,NY,NX)
            IF(IZ==NZ.AND.INDVAR==MAX(BX,BY,BZ)) THEN
              CALL DSCFLS(MAGF,FNAME,14,0)
              DEALLOCATE(GBX,GBY,GBZ,STAT=IOS)
              IF(IOS/=0) GOT O 1108
            ENDIF
          ELSE ! Set B fields for new style file
            L0VAL=L0F(VAL) ! copy from current slab BXC, BYC or BZC into VAL
            IF(INDVAR==BX) THEN
              DO IX=1,NX; DO IY=1,NY
                I=(IX-1)*NY+IY; F(L0VAL+I)=BXC(IX,IY,IZ)
              ENDDO; ENDDO
            ELSEIF(INDVAR==BY) THEN
              DO IX=1,NX; DO IY=1,NY
                I=(IX-1)*NY+IY; F(L0VAL+I)=BYC(IX,IY,IZ)
              ENDDO; ENDDO
            ELSEIF(INDVAR==BZ) THEN
              DO IX=1,NX; DO IY=1,NY
                I=(IX-1)*NY+IY; F(L0VAL+I)=BZC(IX,IY,IZ)
              ENDDO; ENDDO
            ENDIF
            IF(IZ==NZ.AND.INDVAR==MAX(BX,NY,BZ)) THEN ! clear up scratch
              CALL DSCFLS(MAGF,FNAME,14,0)  ! and close file
              DEALLOCATE(BXC,BYC,BZC,STAT=IOS); IF(IOS/=0) GO TO 1108
            ENDIF
          ENDIF !
        ENDIF ! end of MAGF setting or reading
C--- Convert from GAUSS to TESLA
        IF(IMAGU==1) THEN
          CALL FN2(VAL,VAL,0.0,1.E-4)
C--- Convert from milliTesla to TESLA
        ELSEIF(IMAGU==2) THEN
          CALL FN2(VAL,VAL,0.0,1.E-3)
        ENDIF
      ENDIF
C
      IF(INDVAR==HI) THEN
C--- Initialise HI store : cell face heights
       IF(IZ==1) THEN
C--- Initilise interface, anode and free surface height
        CALL GETZ(ZWNZ,GZWNZ,NZ)
C--- Store interface height in GRSP2
        IF(NPATCH(1:5)/='INIF1') THEN
          IF(HUMP==0) THEN
C---  Flat interface
            CALL FN1(GRSP2,GZWNZ(IZ1))
          ELSE
C---  Read interface height from 'HUMP' file
            FNAME=CG(6)
            CALL GETSDC('ESTER','HUMP',FNAME)
            LL=LENGZZ(FNAME); CALL LOWCSE(FNAME,-LL) ! make lowercase
 1142       CONTINUE
            CALL DSCFLS(HUMP,FNAME,12,0)
            IF(FNAME(1:2)=='q1') THEN
              DO WHILE(LINE/='  *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=='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)
 1146       CONTINUE
!            CALL SETYX(GRSP2,GVAL,NYM,NXM)
            CALL SETYX(GRSP2,GVAL,NY,NX)
          ENDIF
        ELSE
C
C--- Sequence for setting inclined interface for transient test case
C
         GZ0=GZWNZ(IZ1)
         IF(NPATCH(6:6)==' '.OR.NPATCH(6:6)=='Y') THEN
C--- PATCH name is INIF1 or INIF1Y. Slope in Y direction
           CALL GETYX(YG2D,GYV,NY,NX)
           CALL GETONE(YV2D,GLEN,NY,NX)
         ELSE
C--- Slope in X direction
           CALL GETYX(XG2D,GYV,NY,NX)
           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,NY,NX)
        ENDIF
C
C--- Store anode height in GRSP3
        IF(HANO/=0) THEN
C---  Read anode underside heights from disc file
          FNAME=CG(5); IREAD=0
          CALL GETSDC('ESTER','HANO',FNAME)
          LL=LENGZZ(FNAME); CALL LOWCSE(FNAME,-LL) ! make lowercase
 1111     CONTINUE
          CALL DSCFLS(HANO,FNAME,12,0)
          IF(FNAME(1:2)=='q1') THEN
            DO WHILE (LINE/='  *HANO')
              READ(HANO,'(A)',ERR=1112,END=1112,IOSTAT=IOS) LINE
            ENDDO
            GO TO 1113
 1112       CONTINUE
            IF(FNAME=='q1') THEN
              CALL DSCFLS(MAGF,FNAME,14,0)
              FNAME='q1ear'
              GO TO 1111
            ENDIF
            GO TO 1010
          ENDIF
 1113     CONTINUE
          LUNITR=HANO; IREAD=0; LINE=' '
          DO WHILE (LINE==' ')
            READ(LUNITR,'(A)',ERR=1010,END=1010,IOSTAT=IOS) LINE ! read first line
            ICOM=INDEX(LINE,'!'); IF(ICOM>0) LINE(ICOM:)=' '
            IF(LINE(1:1)=='#') LINE=' '
          ENDDO
          LL=LENGZZ(LINE)
          CALL SPLTZZZ (LINE,WD,NWDS,NCHARS,LL,' ,;='//CHAR(9),20) ! split into words
          IF(FNAME(1:2)=='q1') THEN
            NEWREAD=.FALSE.
          ELSEIF(WD(1)=='TITLE') THEN
            NEWREAD=.TRUE.
C... Get grid size of data from ZONE
            NXD=0; NYD=0; NZD=0
            DO WHILE(WD(1)/='ZONE')
              READ(LUNITR,'(A)',ERR=1010,IOSTAT=IOS) LINE
              ICOM=INDEX(LINE,'!'); IF(ICOM>0) LINE(ICOM:)=' '
              IF(LINE(1:1)=='#') LINE=' '
              LL=LENGZZ(LINE)
              CALL SPLTZZZ (LINE,WD,NWDS,NCHARS,LL,' ,;='//CHAR(9),20) ! split into words
            ENDDO
            IF(WD(1)=='ZONE') THEN
              DO IWD=2,NWDS
                IF(WD(IWD)=='I') THEN
                  NXD=IRDZZZ(IWD+1); IF(ERROR) GO TO 1010
                ENDIF
                IF(WD(IWD)=='J') THEN
                  NYD=IRDZZZ(IWD+1); IF(ERROR) GO TO 1010
                ENDIF
              ENDDO
            ENDIF
          ELSEIF(NWDS==2) THEN ! if only 2 words, new style file with NX, NY on 1st line
            NEWREAD=.TRUE.
            NXD=IRDZZZ(1); IF(ERROR) GO TO 1010 ! get NX,NY for data
            NYD=IRDZZZ(2); IF(ERROR) GO TO 1010
          ELSE
            NEWREAD=.FALSE.
          ENDIF
          IF(NEWREAD) THEN
            ALLOCATE(XCD(NXD),YCD(NYD),JZC(NXD,NYD), STAT=IOS) ! allocate storage
            IF(IOS/=0) GO TO 1010 ! error allocating memory
            LINE=' '
C... Read the file into XCD(1D), YCD(1D) & JZC(2D)
C... Allow for blank lines and comments after !. Ignore lines starting #
            DO IYD=1,NYD; DO IXD=1,NXD
              DO WHILE(LINE==' ')
                READ(LUNITR,'(A)',END=1010,ERR=1010,IOSTAT=IOS) LINE
                ICOM=INDEX(LINE,'!'); IF(ICOM>0) LINE(ICOM:)=' '
                IF(LINE(1:1)=='#') LINE=' '
              ENDDO
              READ(LINE,*,END=1010,ERR=1010,IOSTAT=IOS) XCD(IXD),
     1                                          YCD(IYD),JZC(IXD,IYD)
              LINE=' '
            ENDDO; ENDDO
            XOFF=0.5*XULAST; YOFF=0.5*YVLAST ! recentre on domain
            XCD=XCD+XOFF; YCD=YCD+YOFF ! add offset
            ALLOCATE(JZCY(NXD,NY), STAT=IOS) ! allocate scratch store for interpolation
            IF(IOS/=0) GO TO 1010 ! error allocating memory
C... Interpolate in Y in JZC into JZCY
            DO IXD=1,NXD ! loop over all X in data
              DO IY=1,NY ! loop over actual Y in solution
                YLOC=F(KYYG+IY) ! actual Y cell-centre
                IF(YLOC<=YCD(1)) THEN ! below start of data, take first point
                  JZCY(IXD,IY)=JZC(IXD,1)
                ELSEIF(YLOC>=YCD(NYD)) THEN ! above end of data, take last point
                  JZCY(IXD,IY)=JZC(IXD,NYD)
                ELSE ! interpolate
                  DO IYD=2,NYD ! loop over all Y in data
                    IF(YLOC<=YCD(IYD)) THEN ! found interval in data
                      JZCY(IXD,IY)=JZC(IXD,IYD-1)+
     1                  (JZC(IXD,IYD)-JZC(IXD,IYD-1))*
     1                  (YLOC-YCD(IYD-1))/(YCD(IYD)-YCD(IYD-1))
                      EXIT
                    ENDIF
                  ENDDO
                ENDIF
              ENDDO ! end IY loop
            ENDDO   ! end IXD loop
C... interpolate in X in JZCY
            DO IY=1,NY ! loop over actual Y
              DO IX=1,NX ! loop over actual X
                XLOC=F(KXXG+IX) ! actual X cell centre
                IF(XLOC<=XCD(1)) THEN ! below start of data, take first
                  HANODE(IY,IX)=JZCY(1,1)
                ELSEIF(XLOC>=XCD(NXD)) THEN ! above end of data , take last
                  HANODE(IY,IX)=JZCY(NXD,1)
                ELSE ! interpolate
                  DO IXD=2,NXD ! loop over X in data
                    IF(XLOC<=XCD(IXD)) THEN ! found interval
                      HANODE(IY,IX)=JZCY(IXD-1,IY)+
     1                  (JZCY(IXD,IY)-JZCY(IXD-1,IY))*
     1                  (XLOC-XCD(IXD-1))/(XCD(IXD)-XCD(IXD-1))
                      EXIT
                    ENDIF
                  ENDDO
                ENDIF
              ENDDO ! end IX loop
            ENDDO   ! end IY loop
            WRITE(LUPR1,
     1 '(/,'' ANODE HEIGHTS  READ AND INTERPOLATED FROM FILE '',A/)')
     1                                                             FNAME
            CALL DSCFLS(LUNITR,FNAME,14,0) ! close the file
            DEALLOCATE(XCD,YCD,JZC,JZCY, STAT=IERR) ! clear memory
          ELSE
            BACKSPACE(HANO)
            DO IY=NY,1,-1
             READ(HANO,*,ERR=1114,END=1114,IOSTAT=IOS) (HANODE(IY,IX),
     1                                                          IX=1,NX)
             IREAD=IREAD+NX
           ENDDO
           IF(MOD(IREAD,NXNY)==0) THEN
             CALL DSCFLS(HANO,FNAME,14,0)
           WRITE(LUPR1,'(/'' ANODE HEIGHTS READ FROM FILE '',A/)') FNAME
             GO TO 11132
           ENDIF
 1114      CONTINUE
           WRITE(LUPR1,'('' ERROR IN READING HANO FILE '',A)') FNAME
           IF(IOS/=0) THEN
             CALL IOEMZZ(IOS,BUFF)
             WRITE(LUPR1,'('' Error is: '',A)') BUFF
           ENDIF
           IF(MOD(IREAD,NXNY)/=0) THEN
             WRITE(LUPR1,'('' File may be for different mesh size'')')
           ENDIF
           CALL SET_ERR(319, 'ERROR IN READING HANO FILE.',1)
         ENDIF
11132    CONTINUE
         CALL SETYX(GRSP3,HANODE,NY,NX)
         IF(IZ1>0) THEN
           IF(IZ0>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<=IZ0) THEN
        CALL FN1(VAL,GZWNZ(IZ))
       ELSEIF(IZ<=IZ1) THEN
        CALL FN2(VAL,GRSP2,0.0,GZWNZ(IZ)/GZWNZ(IZ1))
       ELSE IF(IZ<=IZ2) THEN
        ALFA=(GZWNZ(IZ)-GZWNZ(IZ1))/(GZWNZ(IZ2)-GZWNZ(IZ1))
        CALL FN10(VAL,GRSP3,GRSP2,0.0,ALFA,1.-ALFA)
       ELSE IF(IZ<=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==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==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==FZ) THEN
!C--- Fz
!       IF(IZIZ0) GDEN=51.
        IF(IZ>IZ1) GDEN=52.
        IF(IZ>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>IZ0.AND.IZ<=IFREEZ(IY,IX)+IZ0) F(L0VAL+I)=198.
            IF(IZ>IZ2 .AND. IANODE(IY,IX)==1) F(L0VAL+I)=100.
          ENDDO
        ENDDO
      ENDIF
C
      IF(INDVAR==JZ) THEN
C--- Set Jz to cathode distribution initially - GJCATH read in G1.1
!        CALL SETYX(VAL,GJCATH,NY,NX)
        CALL FN1(VAL,0.0)
      ENDIF
C
      RETURN
C*****************************************************************
C
C--- GROUP 13. Boundary conditions and special sources
C
   13 CONTINUE
      IF(ISC/=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)=='ANOPTNT'.OR.OBJNAM(IPAT)=='ANOPTNTL').AND.
     1                                INDVAR==EPOT.AND.IZ==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=='SURFACE'.OR.OBJNAM(IPAT)=='SURFACE').AND.
     1                                              INDVAR==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/=12) RETURN
C---------------------------------------- value = GRND
      NAMSUB='G13V  '
C
      IF(NPATCH(1:7)=='ANOPTNT'.OR.OBJNAM(IPAT)=='ANOPTNTL') THEN
        IF (INDVAR == 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)
          CALL POTENT (GVAL,NY,NX)
C--- Set into EARTH
!          CALL SETYX(VAL,GVAL,NYM,NXM)
          CALL SETYX(VAL,GVAL,NY,NX)
        ENDIF
      ENDIF
C
      IF(NPATCH=='LORENTZ'.OR.OBJNAM(IPAT)=='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>1) L0JZL=L0F(LOW(JZ))
        L0BX=L0F(BX); L0BY=L0F(BY); L0BZ=L0F(BZ)
        L0VAL=L0F(VAL)
C
        IF(INDVAR==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==NX) THEN
C--- Inside anode or freeze - Fx = 0.0
              F(L0VAL+I)=0.0
            ELSE
C--- Average JY to cell face
              IF(IY>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>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==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==NY) THEN
C--- Inside anode or freeze - Jy = 0.0
              F(L0VAL+I)=0.0
            ELSE
C--- Average JX to cell face
              IF(IX>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>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==W1) THEN
C--- Fz = Jx * By - Jy * Bx
          IF(IZ1) 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>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==W2.AND.NPATCH=='GRAVW2') THEN
        CALL FN2(VAL,R2,0.0,AGRAVZ*(RHO2-RHOELC))
      ENDIF
C
C--- Gravity sources for U2
      IF(INDVAR==U2.AND.NPATCH=='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==V2.AND.NPATCH=='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=='INDUCE'.OR.OBJNAM(IPAT)=='INDUCE').AND.
     1                                           INDVAR==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>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>1.AND..NOT.SLD(I-NY))
     1                         GSU=GSU-F(L0AE+I-NY)*F(L0JIX+I-NY)
              IF(IY>1.AND..NOT.SLD(I-1))
     1                         GSU=GSU-F(L0AN+I-1)*F(L0JIY+I-1)
              IF(IZ>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==P2.AND.NPATCH=='GAS') THEN
        L0PRP=L0F(HIGH(IPRPS))
        L0VAL=L0F(VAL)
        IF(GMDOT>0.0) THEN
C... Constant release rate set in Q1
          DO 13115 I=1,NX*NY
            IF(NINT(F(L0PRP+I))==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))==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=='SURFACE'.OR.OBJNAM(IPAT)=='SURFACE').AND.
     1                                                INDVAR==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==W2.AND.NPATCH=='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==EPOT.AND.(NPATCH(1:7)=='CATHCUR'.OR.OBJNAM(IPAT)
     1                                         =='CATHCURR')) THEN
C--- GJCATH set in Group 19.1. The cathode current should be zero
C    inside the freeze, or else very large potentials will arise.
        CALL SETYX(VAL,GJCATH,NY,NX)
        lbj=lbname('JCAT')
        if(lbj>0.and.isweep==fsweep) then
          call fn0(lbj,val)
          l0jc=l0f(lbj); l0al=i2dalb
          totj=0.0; tota=0.; totjf=0.0; totfree=0.
          do i=1,nx*ny
            totj=totj+f(l0jc+i)*f(l0al+i); tota=tota+f(l0al+i)
            if(.not.sld(i)) then
              totjf=totjf+f(l0jc+i)*f(l0al+i); totfree=totfree+f(l0al+i)
            endif
          enddo
          write(14,'(''total cathode current '',1pe13.6)') totj
          write(14,'(''total cathode area    '',1pe13.6)') tota
          write(14,'(''exposed cathode area  '',1pe13.6)') totfree
          write(14,'(''current thru free area'',1pe13.6)') totjf
        endif
      ENDIF
C
C--- Calculate W at interface
      IF((NPATCH=='INTFACE'.OR.OBJNAM(IPAT)=='INTFACE').AND.
     1                              (INDVAR==W1.OR.INDVAR==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,197,198),ISC
      GO TO (191,20,193,194,20,196,197,198,20,20,20),ISC
  191 CONTINUE
      IF(ISTEP==FSTEP) THEN
        IF(JZCATH<=0.0) THEN
          GJCELL=0.0; CALL GETSDR('ESTER','JCELL',GJCELL) ! total cell current through collector bar
          IF(GJCELL/=0.0) JZCATH=-ABS(GJCELL/(XULAST*YVLAST)) ! must be -ve
C--- Set cathode current to JZCATH directly
          GJCATH=JZCATH
        ELSE
C--- Read cathode current distribution from disc file
          LUNITR=IFIX(JZCATH); IREAD=0
          FNAME=CG(2); CALL GETSDC('ESTER','JZCATH',FNAME)
          LL=LENGZZ(FNAME); CALL LOWCSE(FNAME,-LL) ! make lowercase
 1005     CONTINUE
          CALL DSCFLS(LUNITR,FNAME,12,0) ! open the file
          IF(FNAME(1:2)=='q1') THEN ! if Q1, scan for marker
            DO WHILE (LINE/='  *JZC')
              READ(LUNITR,'(A)',ERR=1010,END=1006,IOSTAT=IOS) LINE
            ENDDO
            GO TO 1007
 1006       CONTINUE
            IF(FNAME=='q1') THEN
              CALL DSCFLS(MAGF,FNAME,14,0)
              FNAME='q1ear'
              GO TO 1005
            ENDIF
            GO TO 1010
          ENDIF
 1007     CONTINUE ! established file name
          LINE=' '; TOTJC=GRND
          DO WHILE (LINE==' ')
            READ(LUNITR,'(A)',ERR=1010,END=1010,IOSTAT=IOS) LINE ! read first line
            ICOM=INDEX(LINE,'!'); IF(ICOM>0) LINE(ICOM:)=' '
            IF(LINE(1:1)=='#') LINE=' '
          ENDDO
          LL=LENGZZ(LINE)
          CALL SPLTZZZ (LINE,WD,NWDS,NCHARS,LL,' ,;='//CHAR(9),20) ! split into words
          IF(FNAME(1:2)=='q1') THEN
            NEWREAD=.FALSE.
          ELSEIF(WD(1)=='TITLE') THEN
            NEWREAD=.TRUE.
            IF(NWDS==2.AND.ISRLZZ((WD(2)))) THEN
              RRR=RRDZZZ(2); IF(ERROR) GO TO 1010; TOTJC=RRR
            ENDIF
C... Get grid size of data from ZONE
            NXD=0; NYD=0; NZD=0
            DO WHILE(WD(1)/='ZONE')
              READ(LUNITR,'(A)',ERR=1010,IOSTAT=IOS) LINE
              ICOM=INDEX(LINE,'!'); IF(ICOM>0) LINE(ICOM:)=' '
              IF(LINE(1:1)=='#') LINE=' '
              LL=LENGZZ(LINE)
              CALL SPLTZZZ (LINE,WD,NWDS,NCHARS,LL,' ,;='//CHAR(9),20) ! split into words
            ENDDO
            IF(WD(1)=='ZONE') THEN
              DO IWD=2,NWDS
                IF(WD(IWD)=='I') THEN
                  NXD=IRDZZZ(IWD+1); IF(ERROR) GO TO 1010
                ENDIF
                IF(WD(IWD)=='J') THEN
                  NYD=IRDZZZ(IWD+1); IF(ERROR) GO TO 1010
                ENDIF
              ENDDO
              IF(TOTJC==GRND) THEN
                DO IWD=2,NWDS
                  IF(WD(IWD)=='T') THEN
                    DO IC=1,NCHARS(IWD+1)
                      IF(WD(IWD+1)(IC:IC)=='"') WD(IWD+1)(IC:IC)=' '
                    ENDDO
                    CALL COMPR2(WD(IWD+1),NCHARS(IWD+1))
                    RRR=RRDZZZ(IWD+1)
                    IF(ERROR) GO TO 1010; TOTJC=RRR
                    EXIT
                  ENDIF
                ENDDO
              ENDIF
            ENDIF
          ELSEIF(NWDS<=3) THEN ! if only 2 or 3 words, new style file with NX, NY on 1st line
            NEWREAD=.TRUE.
              NXD=IRDZZZ(1); IF(ERROR) GO TO 1010 ! get NX,NY for data
              NYD=IRDZZZ(2); IF(ERROR) GO TO 1010
            IF(NWDS==3) THEN ! found 3 items, total current is third
              TOTJC=RRDZZZ(3); IF(ERROR) GO TO 1010
            ELSE             ! found 2 items, total current may be on next line
              LINE=' '
              DO WHILE (LINE==' ')
                READ(LUNITR,'(A)',ERR=1010,END=1010,IOSTAT=IOS) LINE ! read first line
                ICOM=INDEX(LINE,'!'); IF(ICOM>0) LINE(ICOM:)=' '
                IF(LINE(1:1)=='#') LINE=' '
              ENDDO
              LL=LENGZZ(LINE); CALL SPLTZZ (LINE,WD,NWDS,NCHARS,LL) ! split into words
              IF(NWDS==1) THEN ! found single item, get total current
                TOTJC=RRDZZZ(1); IF(ERROR) GO TO 1010
              ELSE             ! total current was not set
                BACKSPACE(LUNITR)
              ENDIF
            ENDIF
          ELSE
            NEWREAD=.FALSE.
          ENDIF
          IF(NEWREAD) THEN
            ALLOCATE(XCD(NXD),YCD(NYD),JZC(NXD,NYD), STAT=IOS) ! allocate storage
            IF(IOS/=0) GO TO 1010 ! error allocating memory
            LINE=' '
C... Read the file into XCD(1D), YCD(1D) & JZC(2D)
C... Allow for blank lines and comments after !. Ignore lines starting #
            DO IYD=1,NYD; DO IXD=1,NXD
              DO WHILE(LINE==' ')
                READ(LUNITR,'(A)',END=1010,ERR=1010,IOSTAT=IOS) LINE
                ICOM=INDEX(LINE,'!'); IF(ICOM>0) LINE(ICOM:)=' '
                IF(LINE(1:1)=='#') LINE=' '
              ENDDO
              READ(LINE,*,END=1010,ERR=1010,IOSTAT=IOS) XCD(IXD),
     1                                          YCD(IYD),JZC(IXD,IYD)
              LINE=' '
            ENDDO; ENDDO
            XOFF=0.5*XULAST; YOFF=0.5*YVLAST ! recentre on domain
            XCD=XCD+XOFF; YCD=YCD+YOFF ! add offset
C... sum up total current from data file
            TOTARD=0.0; TOTJD=0.0; TOTARD0=0; TOTJD0=0.
            DO IXD=1,NXD
              IF(IXD==1) THEN
                DXX=(XCD(IXD+1)-XCD(IXD))
              ELSEIF(IXD=YCD(NYD)) THEN ! above end of data, take last point
                  JZCY(IXD,IY)=JZC(IXD,NYD)
                ELSE ! interpolate
                  DO IYD=2,NYD ! loop over all Y in data
                    IF(YLOC<=YCD(IYD)) THEN ! found interval in data
                      JZCY(IXD,IY)=JZC(IXD,IYD-1)+
     1                  (JZC(IXD,IYD)-JZC(IXD,IYD-1))*
     1                  (YLOC-YCD(IYD-1))/(YCD(IYD)-YCD(IYD-1))
                      EXIT
                    ENDIF
                  ENDDO
                ENDIF
              ENDDO ! end IY loop
            ENDDO   ! end IXD loop
C... interpolate in X in JZCY
            DO IY=1,NY ! loop over actual Y
              DO IX=1,NX ! loop over actual X
                XLOC=F(KXXG+IX) ! actual X cell centre
                IF(XLOC<=XCD(1)) THEN ! below start of data, take first
                  GJCATH(IY,IX)=JZCY(1,1)
                ELSEIF(XLOC>=XCD(NXD)) THEN ! above end of data , take last
                  GJCATH(IY,IX)=JZCY(NXD,1)
                ELSE ! interpolate
                  DO IXD=2,NXD ! loop over X in data
                    IF(XLOC<=XCD(IXD)) THEN ! found interval
                      GJCATH(IY,IX)=JZCY(IXD-1,IY)+
     1                  (JZCY(IXD,IY)-JZCY(IXD-1,IY))*
     1                  (XLOC-XCD(IXD-1))/(XCD(IXD)-XCD(IXD-1))
                      EXIT
                    ENDIF
                  ENDDO
                ENDIF
              ENDDO ! end IX loop
            ENDDO   ! end IY loop
            WRITE(LUPR1,
     1 '(/,'' CATHODE CURRENTS READ AND INTERPOLATED FROM FILE '',A/)')
     1                                                             FNAME
            CALL DSCFLS(LUNITR,FNAME,14,0) ! close the file
            if(.true.) then
              open(lunitr,file='jcath.plt',status='unknown')
              close(lunitr,status='delete',iostat=ios)
              open(lunitr,file='jcath.plt',status='new')
              write(lunitr,'(''TITLE=Intepolated cathode current'')')
              write(lunitr,'(''VARIABLES="X","Y","Jb"'')')
              write(lunitr,'(''ZONE T="B1", I='',I6,'', J='',I6,
     1                                           '', F=POINT'')') nx,ny
              do ix=1,nx
                do iy=1,ny
                  xloc=f(kxxg+ix)-0.5*xulast; yloc=f(kyyg+iy)-0.5*yvlast
                  write(lunitr,'(1p3e13.5)') xloc,yloc,gjcath(iy,ix)
                enddo
              enddo
              close(lunitr,status='keep',iostat=ios)
            endif
            DEALLOCATE(XCD,YCD,JZC,JZCY, STAT=IERR) ! clear memory
            L0AL=I2DALB; TOTJ=0.0; TOTA=0.; TOTFREEA=0.; TOTFREEJ=0.0
            L0FACZ0=L0FACZ; L0FACZ=L0FACE
            DO IX=1,NX; DO IY=1,NY
              I=(IX-1)*NY+IY
              TOTJ=TOTJ+GJCATH(IY,IX)*F(L0AL+I); TOTA=TOTA+F(L0AL+I) ! sum total area & current
              IF(.NOT.VAC(I).AND..NOT.POR(I)) THEN
                TOTFREEA=TOTFREEA+F(L0AL+I) ! sum free area and current
                TOTFREEJ=TOTFREEJ+GJCATH(IY,IX)*F(L0AL+I)
              ENDIF
            ENDDO; ENDDO
            TOTJADJ=0.0
            IF(TOTJC==GRND) THEN
              TOTJ=TOTJD
            ELSE
              TOTJ=TOTJC
            ENDIF
            DO IX=1,NX; DO IY=1,NY
              I=(IX-1)*NY+IY
              GJCATH(IY,IX)=GJCATH(IY,IX)*TOTJ/TOTFREEJ
              IF(.NOT.VAC(I).AND..NOT.POR(I)) THEN
                TOTJADJ=TOTJADJ+GJCATH(IY,IX)*F(L0AL+I)
              else
                gjcath(iy,ix)=0.0
              ENDIF
            ENDDO; ENDDO
            if(.true.) then
              open(lunitr,file='jcath_adj.plt',status='unknown')
              close(lunitr,status='delete',iostat=ios)
              open(lunitr,file='jcath_adj.plt',status='new')
              write(lunitr,
     1       '(''TITLE=Adjusted interpolated cathode current'')')
              write(lunitr,'(''VARIABLES="X","Y","Jb"'')')
              write(lunitr,'(''ZONE T="B1", I='',I6,'', J='',I6,
     1                                           '', F=POINT'')') nx,ny
              do ix=1,nx
                do iy=1,ny
                  xloc=f(kxxg+ix)-0.5*xulast; yloc=f(kyyg+iy)-0.5*yvlast
                  write(lunitr,'(1p3e13.5)') xloc,yloc,gjcath(iy,ix)
                enddo
              enddo
              close(lunitr,status='keep',iostat=ios)
            endif
            L0FACZ=L0FACZ0
            IF(TOTJC==GRND) THEN
              WRITE(14,'('' Estimated Total cathode current     '',
     1                                           1PE13.6,'' A'')') TOTJD
              WRITE(14,'('' Estimated cathode area              '',
     1                                          1PE13.6,'' m^2'')') TOTA
              WRITE(14,'('' Total interpolated cathode current  '',
     1                                            1PE13.6,'' A'')') TOTJ
            ELSE
              WRITE(14,'('' Total cathode current from data     '',
     1                                           1PE13.6,'' A'')') TOTJC
            ENDIF
            WRITE(14,'('' Exposed cathode area                '',1PE13.6
     1                                             ,'' m^2'')') TOTFREEA
            WRITE(14,'('' Actual interpolated cathode current '',1PE13.6
     1                                               ,'' A'')') TOTFREEJ
            WRITE(14,'('' After adjustment                    '',1PE13.6
     1                                               ,'' A'')') TOTJADJ
            CALL WRITST
            GO TO 1015
          ELSE ! old-style file
            BACKSPACE(LUNITR) ! go back one line
            DO IY=NY,1,-1 ! read the data
              READ(LUNITR,*,ERR=1010,END=1010,IOSTAT=IOS)(GJCATH(IY,IX),
     1                                                          IX=1,NX)
              IREAD=IREAD+NX
            ENDDO
            IF(MOD(IREAD,NXNY)==0) THEN
              CALL DSCFLS(LUNITR,FNAME,14,0) ! close the file
              WRITE(LUPR1,
     1             '(/,'' CATHODE CURRENTS READ FROM FILE '',A/)') FNAME
              GO TO 1015
            ENDIF
          ENDIF ! end JZCATH reading
 1015     CONTINUE
        ENDIF ! end
        IF(ISTEP==1.AND.JZ>0) THEN ! copy cathode current to JZ at all slabs
          DO IZZ=1,NZ
            L0JZ=L0F(ANYZ(JZ,IZZ))
            IF(FX>0) THEN
              L0FX=L0F(ANYZ(FX,IZZ)); L0FY=L0F(ANYZ(FY,IZZ))
              L0FZ=L0F(ANYZ(FZ,IZZ))
              L0BX=L0F(ANYZ(BX,IZZ)); L0BY=L0F(ANYZ(BY,IZZ))
            ENDIF
            DO IX=1,NX; DO IY=1,NY
              I=(IX-1)*NY+IY
              F(L0JZ+I)=GJCATH(IY,IX)
              IF(FX>0) THEN
C... Fx = Jy*Bz - Jz*By  (Jy = 0)
                F(L0FX+I)=-F(L0JZ+I)*F(L0BY+I)
C... Fy = Jz*Bx - Jx*Bz  (Jx = 0)
                F(L0FY+I)= F(L0JZ+I)*F(L0BX+I)
C... Fz = Jx*By - Jy*Bx  (Jx=0, Jy = 0)
                F(L0FZ+I)= 0.0
              ENDIF
            ENDDO; ENDDO
          ENDDO
        ENDIF
      ENDIF
      RETURN
  193 CONTINUE
C   * ------------------- SECTION 3 ---- START OF SLAB
C
!      IF(INDVAR==EPOT) RETURN
C---  Set debug switch
      GDBG=IZ>=IG(17).AND.IZ<=IG(18).AND.ISWEEP>=IG(19).AND.
     1      ISWEEP<=IG(20)
C
C--- For RESTART only, recover interface,anode & free surface
C     heights on first sweep
      IF(ISWEEP==FSWEEP .AND. ISTEP==FSTEP .AND.
     1   QEQ(FIINIT(HI),READFI) ) THEN
       IF(IZ==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==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))==100) IANODE(IY,IX)=1
        IF(NINT(F(L0PRP+I))==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==1.AND..NOT.STEADY.AND.IZ==1.AND.ISTEP==FSTEP) THEN
       CALL FN34(GRSP2,ANYZ(W1,IZ1),DT)
      ENDIF
C
C--- Adjust interface height
C    -----------------------
C
      IF(MOD(ISWEEP,NIH)==0.AND.IZ==IZ1.AND.ISWEEP>=IHF.and.indvar==r1)
     1                                                              THEN
       IF(IREF==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>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>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
        DHLIM= MIN(GZWNZ(IZ1)-Z0,GZWNZ(IZ2)-GZWNZ(IZ1))/FDH
        HMIN=Z0+(1.-FHLIM)*(GZWNZ(IZ1)-Z0)
        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)
       CALL GETZ(ZWNZ,GZWNZ,NZ)
       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))
       lbc=lbname('DH1')
       if(lbc>0) l0bc=l0f(lbc)
       lbpd=lbname('PDIF')
       if(lbpd>0) l0pd=l0f(lbpd)
       lbfd=lbname('FDIF')
       if(lbfd>0) l0fd=l0f(lbfd)
       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
       if(lbpd>0) f(l0pd+i)=DPREF-(F(L0P1H+I)-F(L0P1+I))
       if(lbfd>0) f(l0fd+i)=DFZREF-DFZ
       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<0.) THEN
        CHANGE=AMAX1(CHANGE,-DHLIM)
       ELSE
        CHANGE=AMIN1(CHANGE, DHLIM)
       ENDIF
       if(lbc>0) f(l0bc+i)=change
C--- Get new height
       HLOCAL=F(L0HI+I)+CHANGE
C--- Ensure height limits not exceeded
       IF(.NOT.LG(4)) THEN ! limit between cathode & anode underside
!         F(L0HI+I)=AMAX1(HMIN,AMIN1(HLOCAL,FHLIM*(F(L0HA+I)-GZWNZ(IZ1))+
!     1                                                    GZWNZ(IZ1)))
         hlocal=AMAX1(HMIN,AMIN1(HLOCAL,
     1                         FHLIM*(F(L0HA+I)-GZWNZ(IZ1))+GZWNZ(IZ1)))
       ELSE  ! ACD constant so only limit to cathode
!         F(L0HI+I)=AMAX1(HMIN,HLOCAL)
         hlocal=AMAX1(HMIN,HLOCAL)
       ENDIF
!       if(lbc>0) f(l0bc+i-nfm)=hlocal-f(l0hi+i)
       f(l0hi+i)=hlocal
       VOLSUM=VOLSUM+F(L0AH+I)*(F(L0HI+I)-Z0)
19306  CONTINUE
C--- Apply block correction to preserve total volume of metal
       volsum0=volsum; iadj=0
19307  continue
       iadj=iadj+1
       DH=(VLOWER-VOLSUM)/GAREAM
       volsum=0.; VOLSUM2=0.
       DO 19308 I=1,NX*NY
       IF(SLD(I)) GO TO 19308
!       F(L0HI+I)=F(L0HI+I)+DH
       HLOCAL=F(L0HI+I)+dh
C--- Ensure correction does not violate height limits
       IF(.NOT.LG(4)) THEN ! limit between cathode and anode underside
!         F(L0HI+I)=AMAX1(HMIN,AMIN1(F(L0HI+I),
!     1                         FHLIM*(F(L0HA+I)-GZWNZ(IZ1))+GZWNZ(IZ1)))
         hlocal=AMAX1(HMIN,AMIN1(hlocal,
     1                         FHLIM*(F(L0HA+I)-GZWNZ(IZ1))+GZWNZ(IZ1)))
       ELSE  ! limit to above cathode
!         F(L0HI+I)=AMAX1(HMIN,F(L0HI+I))
         hlocal=AMAX1(HMIN,hlocal)
       ENDIF
       if(lbc>0) f(l0bc+i+nfm)=hlocal-f(l0hi+i) ! local second adjustment
       f(l0hi+i)=hlocal ! save local height
       VOLSUM2=VOLSUM2+F(L0AH+I)*(F(L0HI+I)-Z0)
       VOLSUM=VOLSUM+F(L0AH+I)*(F(L0HI+I)-Z0)
19308  CONTINUE
!!!       if(iadj<3) go to 19307
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 (constant ACD)
       IF(LG(4)) CALL FN2(GRSP3,GRSP2,GZWNZ(IZ2)-GZWNZ(IZ1),1.0)
       LBSURF=LBNAME('SURF')
       IF(LBSURF>0) THEN ! SURF exists
C... SURF contains 1 in cells full of metal or air, 0 in cells full of electrolyte.
C    cells containing the interface havr SURF interpolated between 0 - 1 based on
C    height. A VRV or PHOTON iso-surface of SURF=0.5 will show the shape of the
C    interface.
         L0SURF=L0F(ANYZ(LBSURF,1)) ! start address at IZ=1
         DO I=1,NXNY
           GHI=F(L0HI+I) ! current interface height
           DO IZZ=1,NZ
             IF(IZZ==1) THEN
               ZWWM=0.0
             ELSE
               ZWWM=GZWNZ(IZZ-1) ! lower face of cell
             ENDIF
             ZWW=GZWNZ(IZZ) ! upper face of cell
             IZADD=(IZZ-1)*NFM
             IF(GHI<=ZWWM) THEN ! interface below lower face - all electrolyte
               F(L0SURF+I+IZADD)=0.0
             ELSEIF(GHI>ZWW) THEN ! interface above upper face - all metal
               F(L0SURF+I+IZADD)=1.0
             ELSE ! interface in current cell - interpolate between 1 - 0
               F(L0SURF+I+IZADD)=(GHI-ZWWM)/(ZWW-ZWWM)
             ENDIF
           ENDDO
         ENDDO
       ENDIF
 
      ENDIF
C.... store values into InForm variables for optional table plot
      IF(IERR1==0) CALL SET_INF('DH',DH,IERR1)
      IF(IERR2==0) CALL SET_INF('VLOWER',VLOWER,IERR2)
      IF(IERR3==0) CALL SET_INF('VOLSUM',VOLSUM,IERR3)
      IF(IERR3a==0) CALL SET_INF('VOLSM0',VOLSUM0,IERR3a)
      IF(IERR4==0) CALL SET_INF('VOLSM2',VOLSUM2,IERR4)
C
C--- Adjust free surface height
C    --------------------------
C
      IF(LG(2)) THEN
       IF(MOD(ISWEEP,NIH)==0.AND.IZ==IZ3.AND.ISWEEP>=IHF.and.indvar==r1)
     1                                                              THEN
        IF(IREF2==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)-GZWNZ(IZ1))
         FREEAH=FREEAH+F(L0AH+I)
19316    CONTINUE
C--- Get variables for later use
         DHLIM1= (GZWNZ(IZ3)-GZWNZ(IZ2))/FDH
         HMIN1=GZWNZ(IZ2)+(1.-FHLIM)*(GZWNZ(IZ3)-GZWNZ(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))
        lbc=lbname('DH2')
        if(lbc>0) l0bc=l0f(lbc)
        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<0.) THEN
           CHANGE=AMAX1(CHANGE,-DHLIM1)
          ELSE
           CHANGE=AMIN1(CHANGE, DHLIM1)
          ENDIF
C--- Get new height
          HLOCAL=F(L0HIH+I)+CHANGE
          if(lbc>0) f(l0bc+i)=change
C--- Ensure height limits not exceeded
          IF(.NOT.LG(4)) THEN
            F(L0HIH+I)=AMAX1(HMIN1,AMIN1(HLOCAL,(1.+FHLIM)*GZWNZ(NZ)))
          ELSE
            F(L0HIH+I)=AMAX1(HMIN1,HLOCAL)
          ENDIF
          VUPSUM=VUPSUM+F(L0AH+I)*(F(L0HIH+I)-F(L0HI+I))
19318   CONTINUE
C--- Apply block correction to preserve total volume of electrolyte.
        DHU=(VUPPER-VUPSUM)/FREEAH
        VUPSM2=0.0
        DO 19320 I=1,NX*NY
          IF(SLD(I)) GO TO 19320
          F(L0HIH+I)=F(L0HIH+I)+DHU
C--- Ensure correction does not violate height limits
          IF(.NOT.LG(4)) THEN
            F(L0HIH+I)=AMAX1(HMIN1,AMIN1(F(L0HIH+I),
     1                                            (1.+FHLIM)*GZWNZ(NZ)))
          ELSE
            F(L0HIH+I)=AMAX1(HMIN1,F(L0HIH+I))
          ENDIF
          VUPSM2=VUPSM2+F(L0AH+I)*(F(L0HIH+I)-F(L0HI+I))
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
        LBSURF=LBNAME('SURF')
        IF(LBSURF>0) THEN ! SURF exists
C... SURF contains 1 in cells full of metal or air, 0 in cells full of electrolyte.
C    cells containing the interface havr SURF interpolated between 0 - 1 based on
C    height. A VRV or PHOTON iso-surface of SURF=0.5 will show the shape of the
C    interface.
          L0SURF=L0F(ANYZ(LBSURF,1)) ! start address at IZ=1
          DO I=1,NXNY
            GHI=F(L0HIH+I) ! current free surface height
            DO IZZ=1,NZ
              IF(IZZ==1) THEN
                ZWWM=0.0
              ELSE
                ZWWM=GZWNZ(IZZ-1) ! lower face of cell
              ENDIF
              ZWW=GZWNZ(IZZ) ! upper face of cell
              IZADD=(IZZ-1)*NFM
              IF(GHI<=ZWWM) THEN ! interface below lower face - all air
                F(L0SURF+I+IZADD)=1.0
              ELSEIF(GHI>ZWW) THEN ! interface above upper face - all electrolyte
!                F(L0SURF+I+IZADD)=1.0 ! do nothing as already set above
              ELSE ! interface in current cell - interpolate between 0 - 1
                F(L0SURF+I+IZADD)=(ZWW-GHI)/(ZWW-ZWWM)
              ENDIF
            ENDDO
          ENDDO
        ENDIF
       ENDIF
      ENDIF
C.... store values into InForm variables for optional table plot
      IF(IERR5==0) CALL SET_INF('DHU',DHU,IERR5)
      IF(IERR6==0) CALL SET_INF('VUPPER',VUPPER,IERR6)
      IF(IERR7==0) CALL SET_INF('VUPSUM',VUPSUM,IERR7)
      IF(IERR8==0) CALL SET_INF('VUPSM2',VUPSM2,IERR8)
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
      CALL GETZ(ZWNZ,GZWNZ,NZ)
      IF(IZ0>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<=IZ0) THEN
C--- Cathode
       CALL FN1(CPOR,1.0)
       CALL FN1(HI,GZWNZ(IZ))
      ELSEIF(IZ<=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<=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<=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
      IF(IZ>IZ0) THEN
        L0EPOR=L0F(EPOR); L0NPOR=L0F(NPOR)
        L0AE=L0F(AEAST); L0AN=L0F(ANORTH)
        DO I=1,NX*NY
          F(L0AE+I)=F(L0AE+I)/(F(L0EPOR+I)+TINY)
          F(L0AN+I)=F(L0AN+I)/(F(L0NPOR+I)+TINY)
        ENDDO
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--- Apply new porosities
        CALL POR_GEOM(EPOR,IZ) ! multiply by new porosity
        CALL POR_GEOM(NPOR,IZ)
      ELSEIF(ISG62>0.AND.IZ<=IZ0.AND.INDVAR/=EPOT) THEN ! adjust North areas in collector bars
        IF(ISWEEP==FSWEEP.AND.ISTEP==FSTEP) THEN ! only need to do this once
          NCB=0
          CBWID=0.19; CALL GETSDR('ESTER','CBARWID',CBWID)
          CPWID=0.05; CALL GETSDR('ESTER','COPRWID',CPWID)
          DO IPAT=1,NUMPAT ! loop over patches looking for collector bars
            IF(OBJNAM(IPAT)(1:4)=='CBAR') THEN
              CALL GETPAT(IPAT,IR,TYP,JX1,JX2,JY1,JY2,JZ1,JZ2,JT1,JT2)
              IF(IZ>JZ2) CYCLE
              NCB=NCB+1; L0PRPS=L0F(IPRPS); L0DEL=L0F(XU2D)
              XMIN=999.; XMAX=-999.
              DO IX=JX1,JX2
                I=(IX-1)*NY+1
                IF(NINT(F(L0PRPS+I))==103) THEN
                  XMIN=MIN(XMIN,F(L0DEL+I-NY));XMAX=MAX(XMAX,F(L0DEL+I))
                ENDIF
              ENDDO
              DELCB=XMAX-XMIN ! modelled wdth of collector bar
              L0NPOR=L0F(NPOR)
              IF(JY2==NY) JY2=NY-1 ! don't adjust porosity at NY as this upsets net source
              DO IX=JX1,JX2; DO IY=JY1,JY2
                I=(IX-1)*NY+IY
                IF(NINT(F(L0PRPS+I))==103) THEN
                  F(L0NPOR+I)=CBWID/DELCB
                ENDIF
              ENDDO; ENDDO
            ENDIF
          ENDDO
          DO IPAT=1,NUMPAT ! loop over patches looking for collector bars
            IF(OBJNAM(IPAT)(1:4)=='COPR') THEN
              CALL GETPAT(IPAT,IR,TYP,JX1,JX2,JY1,JY2,JZ1,JZ2,JT1,JT2)
              IF(IZ>JZ2) CYCLE
              L0PRPS=L0F(IPRPS); L0DEL=L0F(XU2D)
              XMIN=999.; XMAX=-999.
              DO IX=JX1,JX2
                I=(IX-1)*NY+1
                IF(NINT(F(L0PRPS+I))==106) THEN
                  XMIN=MIN(XMIN,F(L0DEL+I-NY));XMAX=MAX(XMAX,F(L0DEL+I))
                ENDIF
              ENDDO
              DELCP=XMAX-XMIN ! modelled width of copper bar
              L0NPOR=L0F(NPOR)
              IF(JY2==NY) JY2=NY-1 ! don't adjust porosity at NY as this upsets net source
              DO IX=JX1,JX2; DO IY=JY1,JY2
                I=(IX-1)*NY+IY
                IF(NINT(F(L0PRPS+I))==106) THEN
                  F(L0NPOR+I)=CPWID/DELCP
                ENDIF
              ENDDO; ENDDO
            ENDIF
          ENDDO
          IF(NCB>0) CALL POR_GEOM(NPOR,IZ) ! adjust geometry if collector bars found
        ENDIF
      ENDIF
!      write(14,'(''isweep, iz, indvar '',3i3)') isweep,iz,indvar
!      call chklb1(l0npor,'npor')
!      call chklb1(l0an,'anorth')
C
      RETURN
  194 CONTINUE
C
C   * ------------------- SECTION 4 ---- START OF ITERATION.
C
C---  Compute currents and induced currents
      IF(LG(5).OR.(ISWEEP==1.AND.ISTEP==1)) RETURN
C---  Only compute currents when solving for potential
      IF(INDVAR/=EPOT) RETURN
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(IZ1) 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)
      IF(UDIFF) THEN
        DO IY=1,NY
          I=(NX-1)*NY+IY
          F(L0CONDE+I+(IZ-1)*NX*NY)=0.0
        ENDDO
      ENDIF
      lbcone=lbname('#CE')
      l0cone=0; if(lbcone>0) l0cone=l0f(lbcone)
      if(lbcone>0) call fn1(lbcone,0.0)
      DO 1940 IX=1,NXM1
      DO 1940 IY=1,NY
      I=(IX-1)*NY+IY
C--- Get harmonic average of 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))==52.AND.NINT(F(L0PRPS+I+NY))
     1                                                     ==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
      IF(UDIFF) F(L0CONDE+I+(IZ-1)*NX*NY)=GCON
      if(lbcone>0) f(l0cone+i)=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>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>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)
      IF(UDIFF) THEN
        DO IX=1,NX
          I=(IX-1)*NY+NY
          F(L0CONDN+I+(IZ-1)*NX*NY)=0.0
        ENDDO
      ENDIF
      lbconn=lbname('#CN')
      l0conn=0; if(lbconn>0) l0conn=l0f(lbconn)
      if(lbconn>0) call fn1(lbconn,0.0)
      DO 1942 IX=1,NX
      DO 1942 IY=1,NYM1
      I=(IX-1)*NY+IY
C--- Get harmonic average of 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))==52.AND.NINT(F(L0PRPS+I+1))
     1                                                     ==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
      IF(UDIFF) F(L0CONDN+I+(IZ-1)*NX*NY)=GCON
      if(lbconn>0) f(l0conn+i)=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>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>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(UDIFF.AND.IZ==NZ) THEN
        DO I=1,NXNY
          F(L0CONDH+I+(IZ-1)*NX*NY)=0.0
        ENDDO
      ENDIF
      lbconh=lbname('#CH');  lbdlg=lbname('#DZG')
      l0conh=0; if(lbconh>0) l0conh=l0f(lbconh)
      l0dlg=0; if(lbdlg>0) l0dlg=l0f(lbdlg)
      if(lbconh>0) call fn1(lbconh,0.0)
      if(lbdlg>0) call fn1(lbdlg,0.0)
      IF(IZ0) f(l0conh+i)=gcon
C--- Get Jz
      F(L0JZ+I)=-GCON*(F(L0EPOTH+I)-F(L0EPOT+I))/F(L0DELG+I)
      if(lbdlg>0) f(l0dlg+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>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>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
C... Save interface height offset to HINT if it is STOREd
      IF(IHINT>0) THEN
        L0HINT=L0F(IHINT); L0HI=L0F(GRSP2); ZMET=GZWNZ(IZ1)
        DO I=1,NXNY
          IF(.NOT.SLD(I)) THEN
            F(L0HINT+I)=F(L0HI+I)-ZMET
          ELSE
            F(L0HINT+I)=0.0
          ENDIF
        ENDDO
      ENDIF
      RETURN
  197 CONTINUE
C   * ------------------- SECTION 7 ---- Finish of sweep.
      IF(STEADY.AND.LDOHI) THEN
        IF(MOD(ISWEEP,HIFREQ)==0.OR.ISWEEP==LSWEEP) THEN
          WRITE(FNAME,'(''HINT'',I4.4,''.CSV'')') ISWEEP
          LL=LENGZZ(FNAME); CALL LOWCSE(FNAME,-LL)
          LU=111
          OPEN(LU,FILE=FNAME,STATUS='NEW',IOSTAT=IOS)
          IF(IOS/=0) THEN
            OPEN(LU,FILE=FNAME,STATUS='UNKNOWN',IOSTAT=IOS)
            CLOSE(LU,STATUS='DELETE',IOSTAT=IOS)
            OPEN(LU,FILE=FNAME,STATUS='NEW',IOSTAT=IOS)
            IF(IOS/=0) RETURN
          ENDIF
          L0XG=L0F(XG2D); L0YG=L0F(YG2D); L0HI=L0F(GRSP2)
          FORM='(999(F8.4,1H,),F8.4)'
          WRITE(FORM(2:4),'(I3)') NX
          CALL GETZ(ZWNZ,GZWNZ,NZ); ZMET=GZWNZ(IZ1)
          XM=0.5*XULAST; YM=0.5*YVLAST
          WRITE(LU,FORM) 0.,((F(L0XG+1+(IX-1)*NY)-XM),IX=1,NX)
          DO IY=NY,1,-1
            WRITE(LU,FORM) (F(L0YG+IY)-YM),
     1                     ((F(L0HI+(IX-1)*NY+IY)-ZMET),IX=1,NX)
          ENDDO
          CLOSE(LU,STATUS='KEEP',IOSTAT=IOS)
        ENDIF
      ENDIF
      RETURN
  198 CONTINUE
C   * ------------------- SECTION 8 ---- Finish of time step.
      IF(.NOT.STEADY.AND.LDOHI) THEN
        IF(MOD(ISTEP,HIFREQ)==0.OR.ISTEP==LSTEP) THEN
          WRITE(FNAME,'(''HINT'',I4.4,''.TXT'')') ISTEP
          LL=LENGZZ(FNAME); CALL LOWCSE(FNAME,-LL)
          LU=111
          OPEN(LU,FILE=FNAME,STATUS='NEW',IOSTAT=IOS)
          IF(IOS/=0) THEN
            OPEN(LU,FILE=FNAME,STATUS='UNKNOWN',IOSTAT=IOS)
            CLOSE(LU,STATUS='DELETE')
            OPEN(LU,FILE=FNAME,STATUS='NEW',IOSTAT=IOS)
            IF(IOS/=0) RETURN
          ENDIF
          L0XG=L0F(XG2D); L0YG=L0F(YG2D); L0HI=L0F(GRSP2)
          FORM='(999(F8.4,1H,),F8.4)'
          WRITE(FORM(2:4),'(I3)') NX
          CALL GETZ(ZWNZ,GZWNZ,NZ); ZMET=GZWNZ(IZ1)
          XM=0.5*XULAST; YM=0.5*YVLAST
          WRITE(LU,FORM) 0.,((F(L0XG+1+(IX-1)*NY)-XM),IX=1,NX)
          DO IY=NY,1,-1
            WRITE(LU,FORM) (F(L0YG+IY)-YM),
     1                     ((F(L0HI+(IX-1)*NY+IY)-ZMET),IX=1,NX)
          ENDDO
          CLOSE(LU,STATUS='KEEP',IOSTAT=IOS)
        ENDIF
      ENDIF
      RETURN
   20 CONTINUE
      RETURN
C... Error returns
 1108 WRITE(LUPR1,9981) FNAME
 9981 FORMAT(' ERROR IN READING MAGNETIC FIELDS FILE ',A)
      IF(IOS/=0) THEN
        CALL IOEMZZ(IOS,BUFF)
        WRITE(LUPR1,'('' Error is: '',A)') BUFF
      ENDIF
      IF(MOD(IREAD,NXNY)/=0) THEN
        WRITE(LUPR1,'('' File may be for different mesh size'')')
      ENDIF
      CALL DSCFLS(MAGF,FNAME,14,0)
      CALL SET_ERR(317, 'Error. See result file.',1)
C
 1010 CONTINUE ! there was an error or the marker was not found
      WRITE(LUPR1,'(/,'' ERROR IN READING CURRENTS FILE '',A)')FNAME
      IF(LINE/=' ') WRITE(LUPR1,'(''Reading from: '',A)') LINE
      IF(IOS/=0) THEN
        CALL IOEMZZ(IOS,BUFF)
        WRITE(LUPR1,'('' Error is: '',A)') BUFF
      ENDIF
      IF(MOD(IREAD,NXNY)/=0) THEN
        WRITE(LUPR1,'('' File may be for different mesh size'')')
      ENDIF
      CALL SET_ERR(317, 'Error. See result file.',1)
 
      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 '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<0) RETURN
      CALL L0F2(K1,K2,I,I2M1,IADD,'FNAVXY')
          IF(IDIR==0) GO TO 3
        DO 2 IX=IXF,IXL
      I=I+IADD
        DO 2 IY=IYF,IYL
      I=I+1
          IF(IY==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==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 'farray'
      COMMON /IGE/ IXF,IXL,IYF,IYL,IGFILL(21)
      COMMON /NAMFN/ NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      SAVE
      NAMSUB='FNSUM2'
          IF(IXL<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 'farray'
      COMMON /IGE/ IXF,IXL,IYF,IYL,IGFILL(21)
      COMMON /NAMFN/ NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      SAVE
      NAMSUB='FNSUM1'
          IF(IXL<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 'farray'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE '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)),
     1            (RHOCOL,RG(33)),(CONCOL,RG(34)),
     1            (RHOROD,RG(35)),(CONROD,RG(36)),
     1            (RHOBUS,RG(37)),(CONBUS,RG(38)),
     1            (RHOBUS,RG(37)),(CONBUS,RG(38)),
     1            (RHOSTB,RG(71)),(CONSTB,RG(72)),
     1            (RHOCOP,RG(73)),(CONCOP,RG(74)),
     1            (RHCST1,RG(75)),(CNCST1,RG(76)),
     1            (RHCST2,RG(77)),(CNCST2,RG(78))
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
          F(L0CON+I)=CNDCT(I,L0PRP,L0PRPH)
        ENDDO
      ENDDO
      RETURN
      END
C---------------------------------------------------------------------
      REAL FUNCTION CNDCT(I,L0PRP,L0PRPH)
      INCLUDE 'farray'
      COMMON /IDATA/ NX,NY,NZ,IDFILL(117)
      COMMON/RDATA/TINY,RDFIL1(84)
      COMMON/RGRND/RG(100)
      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))
      EQUIVALENCE (RHOCOL,RG(33)),(CONCOL,RG(34)),
     1            (RHOROD,RG(35)),(CONROD,RG(36)),
     1            (RHOBUS,RG(37)),(CONBUS,RG(38)),
     1            (RHOSTB,RG(71)),(CONSTB,RG(72)),
     1            (RHOCOP,RG(73)),(CONCOP,RG(74)),
     1            (RHCST1,RG(75)),(CNCST1,RG(76)),
     1            (RHCST2,RG(77)),(CNCST2,RG(78))
      LOGICAL NEZ
      IPROP=NINT(F(L0PRP+I))
      IF(IPROP==51.OR.IPROP==-1) THEN
C--- Metal
        CNDCT=CONMET
      ELSEIF(IPROP==52) THEN
C--- Bath
        CNDCT=CONELC
        IF(NEZ(CONFAC)) THEN
          IY=1+MOD(I-1,NY); IX=1+(I-1)/NY
          IF((IY 1.AND.NINT(F(L0PRP+I- 1))==100).OR.
     1       (IX 1.AND.NINT(F(L0PRP+I-NY))==100).OR.
     1        NINT(F(L0PRPH+I))==100) CNDCT=CNDCT*CONFAC
        ENDIF
      ELSEIF(IPROP==100) THEN
C--- Anode
        CNDCT=CONANO
      ELSEIF(IPROP==101) THEN
C--- Anode rod
        CNDCT=CONROD
      ELSEIF(IPROP==102) THEN
C--- Anode busbar
        CNDCT=CONBUS
      ELSEIF(IPROP==103) THEN
C--- Collector bar
        CNDCT=CONCOL
      ELSEIF(IPROP==151.OR.IPROP==104) THEN
C--- Cathode
        CNDCT=CONCAT
      ELSEIF(IPROP==105) THEN
C--- Anode stub
        CNDCT=CONSTB
      ELSEIF(IPROP==106) THEN
C--- Copper insert
        CNDCT=CONCOP
      ELSEIF(IPROP==107) THEN
C--- Cast iron around collector bar
        CNDCT=CNCST1
      ELSEIF(IPROP==108) THEN
C--- Cast iron around anode stubs
        CNDCT=CNCST2
      ELSEIF(IPROP==105) THEN
C--- Anode stub
        CNDCT=CONSTB
      ELSEIF(IPROP==198.OR.IPROP==199)THEN
C--- Freeze / airspace, treated as insulator
        CNDCT=TINY
      ENDIF
      END
C---------------------------------------------------------------------
      SUBROUTINE POTENT (GPOT,NYY,NXX)
C--- Example subroutine for setting anode potentials
      INCLUDE 'farray'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE '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, BUFF*256
C
      IF(IGO==12) THEN
C... Open existing file, STATUS = 'OLD'
        OPEN(LUNIT, FILE = FNAME, STATUS = 'OLD', ERR=100, IOSTAT=IOS)
        REWIND LUNIT
      ELSE IF(IGO==14) THEN
        CLOSE(LUNIT,STATUS='KEEP')
      ENDIF
      GO TO 999
C
C... Error in opening file
  100 WRITE(BUFF,'(''Error number '',I3,'' in opening file :'',A)')
     1   IOS,FNAME
      WRITE(LUPR1,'(A)') BUFF(1:LENGZZ(BUFF))
      CALL IOEMZZ(IOS,BUFF)
      WRITE(LUPR1,'(A)') BUFF(1:LENGZZ(BUFF))
      CALL SET_ERR(322, 'Error in opening file.',2)
  999 RETURN
      END
c