cc

C.... File Name ................ GREX3.FOR ...................... 280623
c
C     This subroutine is called by EARTH within the solution cycle.
C
C     Subroutine GREX3, with the subsidiary subroutines GXxxxx,
C     allows the PHOENICS EARTH program to access a wide range of
C     special physical models, algorithm adjustments and print-out
C     features.
C
C     'EX' stands for 'example', signifying that the coding practices
C     exhibited here may be used as patterns by PHOENICS users who
C     wish to introduce their own models, algorithms, etc.
C
C     The GX subroutines that are called in group 13 are activated
C     by settings in the SATELLITE, mainly by way of PATCH names.
C     Many of them, like GREX3 itself, call on 'utility' subroutines'
C     named FNxxx.
C
C     The FNxxx subroutines perform various arithmetic operations on
C     segments of the F-array corresponding to variables of like kind,
C     indicated by the arguments. The FUNCT entry of the Encyclopaedia
C     provides full explanations.
C
c
      SUBROUTINE GREX3
C--------------------------------------- INCLUDE COMMON BLOCK files
      INCLUDE 'patcmn'
      INCLUDE 'objnam'
      INCLUDE 'farray'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE 'grdbfc'
      INCLUDE 'parear'
      INCLUDE 'lbnamer'
      CHARACTER*48 COMMAND
      COMMON/IMAIN2/IMAIN1(9),ISWSTP,IMAIN2(33)
      COMMON/F0/KF01(70),KAP,KF02(233)
     1      /MRKCMN/LBMARK,L0MARK
      COMMON/LBFC/STORSA(6),STORWD(6),LBFCSP /ARRT/ARRIVT
      LOGICAL STORSA,STORWD,LBFCSP,LGTEMP,SPCLGR,ARRIVT
      COMMON/RDLOG/LOGQ(10)/RDINT/INTQ(10)/RDREAL/REALQ(10)
      COMMON/RAKEEP/RADIAT/DVMOD/IDVCGR
      LOGICAL RADIAT,NEZ,dbgloc
      COMMON /LRNTM3/L0UTAU,NMWALL,L0WALL,L0DSKN,IVPRST
      LOGICAL NOCALC,DRAGON,DRG1ST,PBFC,DFAIL,CALLSL
      COMMON/DRAGCM/NOCALC,DRAGON,DRG1ST,PBFC,DFAIL
      COMMON /LUNITS/LUNIT(60)
     1       /GENL/LFIL1(3),TIMPRN,LGF1(10),debgtz,SWPRIN,LGF2(40),
     1        TURB,LGF3(3)
      LOGICAL SLIPVL,LG,KEBUOY,LGF1,debgtz,TIMPRN,SWPRIN,LGF2,TURB,LGF3,
     1        HOCS,ESCRS,FSQSOR,RADCVD,BLIN,LGXIO,LSUN,LWAVE,SURFT
      SAVE SLIPVL,NUMSPPR,KEBUOY,SPCLGR,ESCRS,FSQSOR,HOCS,PSIVEL,
     1     BLIN,LWAVE,SURFT
      COMMON/GRXFLAGS/LGXIO,LSUN
      COMMON/LBINGREX/LBSLU,LBSLV,LBSLW,LBWAVE,LBGENK,LBGNK2,LBMACH,
     1                LBPTOT,LBVABS,LBVLSQ,LBPTO2,LBVAB2,LBMAC2,LBFSQ,
     1                LBARRT,LBWCRT,LBPSIV,LBVOXY,LBVOYZ,LBVOZX
      COMMON/LRNTM1/L0WDIS,L0FMU,L0FONE,L0FTWO,L0REYN,L0REYT,L0UD1,L0UD2
     1             ,L0UD3,L0UD4
      COMMON/PASQLF/ITPRO,PASQBUOY,BUOSSG,MOLEN
      LOGICAL PASQBUOY,BUOSSG
      COMMON/EXPLOSION/EXPLOS
      LOGICAL EXPLOS,PSIVEL
C.... Provision for the use of GRAPH
      DIMENSION TIME(100),OR(100,10),LABOR(10),IOSC(10)
      CHARACTER*4 LABOR,CG
      PARAMETER (NLG=100, NIG=100, NRG=200, NCG=100)
      COMMON/LGRND/LG(NLG)/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      COMMON/GENI/NXNY,IFIL1(8),NFM,IGFIL1(21),IPRL,IBTAU,ILTLS,
     1 IGFIL(15),ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,
     1 IRADZ,IVFOL
      PARAMETER (MXXLET=20)
      COMMON /LETL0C/ L0C(MXXLET)
      CHARACTER*(MXXLET) CINDEX, WREF*2
      COMMON /LETBLK/ CINDEX
      COMMON/IELAPS/IREC,ITIME1,NUMSOL
      COMMON/GRX/NUMSEC
      COMMON/RHILO/HI3D,RLO3D
      LOGICAL DONE192,GXMON,SLD,DOIT
      COMMON/GXMON2/PLTCLSX; LOGICAL PLTCLSX
      SAVE DONE192
      SAVE TIME,OR,LABOR
      DIMENSION DIGIT(10)
      CHARACTER*2 DIGIT
      DATA DIGIT/'1','2','3','4','5','6','7','8','9','10'/
      DATA DONE192 /.FALSE./
c
      dbgloc=debgtz.and.dbgrnd
      if(dbgloc) then
        call writbl
        call banner(1,'grex  ',120121)
        write(14,*)'entry to  grex for: igr,isc,indvar ',igr,isc,indvar
      endif
      if(test) call fcheck(0,'grex3 ')
c----------------------------------------------------------------------
      NAMSUB='GREX3'
C.... GO TO the group indicated by IGR, then to section ISC
      IF(IGR==19) THEN
        GO TO 19
      ELSEIF(IGR==13) THEN
        GO TO 13
      ELSE
        GO TO (1,2,3,4,5,6,25,8,9,10,11,12,13,14,25,25,25,25,19,
     1         20,21,25,25,24) IGR
      ENDIF
C****************************************************************
C
C---- GROUP 1. Run title and other preliminaries
C
    1 GO TO (1001,1002,1003),ISC
C
C.... group 1 sections are called in the order 1, 2, 3.
C****************************************************************
C
C   * -----------GROUP 1  SECTION  1 ------ Initializations and
C     allocations of storage for variables which it is essential to
C     dump to the PHI (or PHIDA) file for restart purposes
C
 1001 IX=1 ; IY=1
C----------------------------------- Identification of GREX3 to VDU
      IF(.NOT.NULLPR.AND.UWATCH)
     1          CALL WRYT40('GREX3 of 280623 has been called')
C
      IF(READQ1) CALL GXRDQ1 ! Read Q1EAR file when READQ1 = T in Q1
c
c.... Set absolute value below which zeroes will appear in field-value
c     print-out (unless debug=t, which sets print0 to 0.0).
c     Users wishing to change print0 without re-compiling should insert
c       PRT0BEGIN
c       PRINT0 [desired value]
c       PRT0END
c     at the top of their Q1 file
c
      IF(.NOT.DEBUG) PRINT0=1.E-11
      CALL RQ1R('PRT0','PRINT0',PRINT0)
C.... Set up PATCH-wise storage, if necessary, for GXBFC
      IF(BFC) CALL GXBFC
      PSIVEL=LBPSIV/=0
C------------------------ PHI and XYZ files for parabolic cases
C
      IF(PARAB.OR.(.NOT.STEADY.AND.PNAM=='    ')) THEN
        IF(IDISPA>0) CALL GXPARA
      ENDIF
      NUMSEC=0
C.... Initializations dependent on the existence of patch names
      RADIAT=.FALSE. ; KEBUOY=.FALSE. ; HOCS=.FALSE. ; FSQSOR=.FALSE.
      BLIN=.FALSE. ;LGXIO=.FALSE.; LSUN=.FALSE.; LWAVE=.FALSE.
      SURFT=.FALSE. ; PASQBUOY=.FALSE.
      DO 10011 I = 1,NUMPAT
        IF(NAMPAT(I)(1:6)=='KEBUOY') THEN
          KEBUOY=.TRUE.
        ELSEIF(NAMPAT(I)(1:4)=='RADI') THEN
          RADIAT=.TRUE.
          CALL GXRADI(CARTES,TEMP1,DEN1,RHO1) ! Radiation initialization
        ELSEIF(NAMPAT(I)(1:4)=='HOCS') THEN
          HOCS=.TRUE.
        ELSEIF(NAMPAT(I)(1:4)=='STEN') THEN
          SURFT=.TRUE.
        ELSEIF(NAMPAT(I)(1:5)=='FSQSO') THEN
          FSQSOR=.TRUE.
        ELSEIF(NAMPAT(I)(1:4)=='BLIN'.OR.NAMPAT(I)(1:4)=='GXBL')THEN
          BLIN=.TRUE.          ! Boundary-layer inlet profile facility
        ELSEIF(.NOT.LWAVE.AND.NAMPAT(I)(1:4)=='WAVE'.AND.IVFOL>0) THEN
          CALL GETCOV(NAMPAT(I),IVFOL,GCO,GVAL)
          IF(GCO/=-999.) LWAVE=.TRUE.
        ELSEIF(.NOT.BFC) THEN
          COMMAND=' '
          CALL GETSDC(OBJNAM(I),'INLET',COMMAND)
          IF(COMMAND==' ') CALL GETSDC(OBJNAM(I),'OUTLET',COMMAND) ! test for
                                                            ! angled-out object
          IF(COMMAND=='GXI'.OR.COMMAND=='GXO') LGXIO=.TRUE.
        ENDIF
10011 CONTINUE
      LGXIO=LGXIO.OR.(MIMD.AND.NPROC>1)
      COMMAND=' '
      CALL GETSDC('SUNLIGHT','TIME',COMMAND)
      IF(COMMAND/=' ') LSUN=.TRUE.
C---------------------------------- SLIPVL setting
      IF(.NOT.ONEPHS) THEN
        SLIPVL = LBNAME('SLPU')/=0.OR.LBNAME('SLPV')/=0.OR.
     1           LBNAME('SLPW')/=0
        IF(SLIPVL) THEN
          LBSLU=0 ; LBSLV=0 ;LBSLW=0
          IF(SOLVE(3).AND.SOLVE(4)) LBSLU=LBNAME('SLPU')
          IF(SOLVE(5).AND.SOLVE(6)) LBSLV=LBNAME('SLPV')
          IF(SOLVE(7).AND.SOLVE(8)) LBSLW=LBNAME('SLPW')
        ENDIF
        IF(NEZ(CPIP)) CALL GXINTP
        IF(NEZ(CLIFT)) CALL GXLIFT
      ENDIF
      IF(BLIN) CALL GXBLIN   ! Boundary-layer inlet profile facility
      IF(FDFSOL) CALL GXFDFS ! Fully-developed floq
      IF(HOCS) CALL GXHOCS
      IF(NOX1) CALL GXNOX
      IF(NOX2) CALL GXNOX
      IF(CHMKIN) CALL GXCHKI(0,0.0)
!!!      IF(LSG87) CALL GXTCALC(0,0.0)
      IF(LBARRT/=0) CALL GXARRT
      IF((GRNDNO(9,TMP1).OR.GRNDNO(10,TMP1)).AND..NOT.LSG87) THEN
        ESCRS=.TRUE.
        IF(GRNDNO(9,TMP1)) THEN
C Define operating pressure in atmospheres via CHSOC for CHEMKIN
          CHSOC=PRESS0/1.01325E5
          CALL GXCHKI(0,0.0)
        ENDIF
        CALL GXSCRI
      ELSE
        ESCRS=.FALSE.
      ENDIF
      EXPLOS=.FALSE.
C***************************************************************
      CALL UCASZZ(NAMGRD,48)
      SPCLGR= YPLS.OR.WALPRN.OR.S2SR
     1    .OR.NAMGRD=='FURN'.OR.NAMGRD=='ESTR'.OR.NAMGRD=='FLAR'
     1    .OR.NAMGRD=='HTBX'.OR.NAMGRD=='TACT'.OR.NAMGRD=='MICA'
     1    .OR.NAMGRD=='CONV'.OR.NAMGRD=='CVD '.OR.NAMGRD=='F1'
     1    .OR.NAMGRD=='WATS'.OR.NAMGRD=='WAVE'
      LL=LENGZZ(NAMGRD)
      IF(SPCLGR.AND.NAMGRD/='NONE') THEN
        IF(.NOT.NULLPR) CALL WRIT40
     1     ('Special GROUND '//NAMGRD(1:LL)//' is in use           ')
      ELSEIF(NAMGRD/='NONE') THEN
        CALL WRITST
        CALL WRITBL
        CALL WRIT40
     1   ('Stop because '//NAMGRD(1:LL)//' is not known to GREX3 ')
        CALL WRYT40
     1   ('Stop because '//NAMGRD(1:LL)//' is not known to GREX3 ')
        CALL SET_ERR(213,'Error. See result file.',1)
      ENDIF
      IF(MYID==0) CALL LATEST_DUMP(1,0,0,' ',' ',' ')
      GO TO 25
C
C***********************************************************************
C
C   * -----------GROUP 1  SECTION  3 ---------------------------
c
C---- Use Section 3 to create storage via MAKE or GXMAKE which it is
C     not necessary to dump to the PHI(da) file for restarts.
C
 1003 CONTINUE
C
C.... Provision of special EARTH arrays by means of the
C     storage-setting subroutine MAKE, used thus:
C     CALL MAKE(variable name)
C
C---------------------------------------------------------BFC=T
      IF(BFC) THEN
        CALL GXBFC
C.... Segment EASP4 is used for buoyancy sources when BFC=T
C     (see GXBUOY)
        CALL MAKE(EASP4)
C
      ENDIF
      IF(SURFT) CALL GXSURFT
C
C-------------------------------------------------EASP1, EASP2, EASP3
C
C.... F-array segments with zero-locations EASP1, EASP2 & EASP3 are
C     always available, for use as local temporary stores
C     (See for example their uses in GXWALL & GXPORA)
      CALL MAKE(EASP1)
      CALL MAKE(EASP2)
      CALL MAKE(EASP3)
      CALL MAKE(EASP5) ! to store SKINT in gxwall.for
      IF(UCONV.AND.NZ>1) CALL MAKE(EASP20)
C
C-------------------------------------------------------EASP7
C.... Store for k-eps source-term linearisation
      IF(BFC.OR.KELIN==1) CALL MAKE(EASP7)
C
C---------------------------------------------------Letter Mask
      MXL=MXXLET
      CINDEX=' '
      DO 111 IPAT=1,NUMPAT
        KK=INDEX(NAMPAT(IPAT)(:2),'%')
        IF(KK/=0) THEN
          WREF=NAMPAT(IPAT)(KK+1:KK+2)
C*        WREF(:1)  = $  LETTER MASK IN XY PLANE
C*                  = &  LETTER MASK IN ZY PLANE
C*                  = #  LETTER MASK IN XZ PLANE
          IF(WREF(:1)=='$' .OR. WREF(:1)=='&'
     1                    .OR. WREF(:1)=='#') THEN
            K=INDEX(CINDEX,WREF(2:2))
            IF(K==0) CALL LETMSK(CINDEX,WREF(2:2),L0C,MXL)
          ENDIF
        ENDIF
  111 CONTINUE
C---------------------------------------------------------LGEN1
C.... Provision of storage for square-of-velocity-gradient
C     expressions for viscous-dissipation or turbulence-
C     energy source.
C
      IF(GENK .OR. DUDX.OR.DVDX.OR.DWDX.OR.DUDY.OR.DVDY.OR.DWDY.OR.
     1             DUDZ.OR.DVDZ.OR.DWDZ) THEN
        IF(BFC) THEN
          LIMIT1=1
          LIMIT2=18
        ELSE
          LIMIT1=5
          LIMIT2=9
        ENDIF
        DO 1000 I= LIMIT1,LIMIT2
          CALL MAKE(LXYEA(I))
 1000   CONTINUE
      ENDIF
      IF((SOLVE(H1).AND.(GRNDNO(5,TMP1).OR.GRNDNO(6,TMP1)))
     1  .OR.(SOLVE(H2).AND.(GRNDNO(5,TMP2).OR.GRNDNO(6,TMP2)))) NCRT=1
      IF(IENUTA.EQ.14) CALL GX_REALISABLE_KE
      IF(IENUTA.EQ.10.OR.IENUTA.EQ.11) CALL GXKW_WILCOX ! not needed.
      IF(IENUTA.EQ.15.OR.IENUTA.EQ.16) CALL GXKW_WILCOX
      IF(IENUTA.GE.17.AND.IENUTA.LE.20) CALL GXKW_MENTER
      GO TO 25
C
C***********************************************************************
C   * -----------GROUP 1  SECTION  2 Further initializations which
C     require that storage shall aready have been allocated
c
 1002 CONTINUE
c      call fcheck(10,namsub)
      IF(KEBUOY) CALL GXKEGB
c      call fcheck(20,namsub)
      namsub='grex3'
      IF(FSQSOR) CALL GXFSQS
C.... Wall-function initializations
      IF(TURB .AND. NMWALL>0) CALL GXWFUN
      IF(HOCS) CALL GXHOCS
      IF(.NOT.ONEPHS) LBWAVE=LBNAME('WAVE')
C.... Possibly set BFC grid by way of coding in GXBFGR
      IF(SETBFC.OR.MOVBFC) CALL GXBFGR
      IF(LGXIO) CALL GXIO
      IF(LSUN) CALL GXSUN
      IF(BLIN) CALL GXBLIN   ! Boundary-layer inlet profile facility
      CALL GXOUTLET
      CALL WRTPAT
      GO TO 25
c
C*****************************************************************
C
C
C--- GROUP 2. Transience; time-step specification
C
    2 CONTINUE
C   * Set DT if TLAST has been made <=GRND in satellite
C... Set timestep from Courant Number
      IF(TLAST-TFIRST==GRND1) CALL GXTIMESTEP
      GO TO 25
C*****************************************************************
C
C
C--- GROUP 3. X-direction grid specification
C
    3 CONTINUE
C   * Set XRAT if AZXU has been made <=GRND in satellite
      GO TO 25
C*****************************************************************
C
C
C--- GROUP 4. Y-direction grid specification
C
    4 CONTINUE
C   * Set YRAT if AZYV has been made <=GRND in satellite
C
      GO TO 25
C*****************************************************************
C
C
C--- GROUP 5. Z-direction grid specification
C
    5 CONTINUE
C   * Set DZ if AZDZ has been made <=GRND in satellite
      IF(GRNDNO(1,AZDZ)) THEN     ! azdz = grnd1 or propx
        DZ=DZW1*XULAST            ! i.e. proportional to xulast
      ELSEIF(GRNDNO(2,AZDZ)) THEN ! azdz = grnd2 or propy
        DZ=DZW1*YVLAST            ! i.e. proportional to yvlast
      ENDIF                       ! dzw1 = proportionality factor
c   * Insert DZ in the appropriate slab-wise array
      CALL FN1(LXYDZ,DZ)
      GO TO 25
C*****************************************************************
C
C
C--- GROUP 6. Body-fitted coordinates or grid distortion
C    This group is visited when UGEOM is set .TRUE.. The visits
C    occur at the start of each z slab after its geometry has been
C    computed. Hence, this is the right place to access AEAST, AHIGH,..
C    etc.
    6 CONTINUE
C
      GO TO 25
C*****************************************************************
C   * Make changes for this group only in group 19.
C
C--- GROUP 7. Variables stored, solved & named
C*****************************************************************
C
C
C--- GROUP 8. Terms (in differential equations) & devices
C
    8 GO TO (81,82,83,84,85,86,87,88,89,810,811,812,813,814,815,816)
     1         ISC
C   * Add velocities relative to grid for phase 1 and/or phase 2
C   * -----------GROUP 8  SECTION  1 ---------------------------
C--- for U1AD<=GRND--- phase 1 additional velocity (VELAD).
   81 CONTINUE
C--- for U1AD<=GRND--- phase 1 additional velocity (VELAD).
c.... The option U1AD=GRND1 adds to U1 a component equal  -W1*TAN(ALF)
c     where ALF is the east-face grid angle. It is needed for the
c     parabolic solver for grids for which constant-x lines are
c     not precisely at right angles to constant-z lines. The
c     correction is also available for U2, V1 and V2 (see below).
c
      IF(GRNDNO(1,U1AD)) CALL GXPVAD(XULAST,XRAT,XU2D,W1,'X')
      GO TO 25
C   * -----------GROUP 8  SECTION  2 ---------------------------
C--- for U2AD<=GRND--- phase 2 additional velocity (VELAD).
   82 CONTINUE
      IF(GRNDNO(1,U2AD)) CALL GXPVAD(XULAST,XRAT,XU2D,W2,'X')
      GO TO 25
C   * -----------GROUP 8  SECTION  3 ---------------------------
   83 CONTINUE
C--- for V1AD<=GRND--- phase 1 additional velocity (VELAD).
      IF(GRNDNO(1,V1AD)) CALL GXPVAD(YVLAST,YRAT,YV2D,W1,'Y')
      GO TO 25
C   * -----------GROUP 8  SECTION  4 ---------------------------
C--- for V2AD<=GRND--- phase 2 additional velocity (VELAD).
   84 CONTINUE
      IF(GRNDNO(1,V2AD)) CALL GXPVAD(YVLAST,YRAT,YV2D,W2,'Y')
      GO TO 25
C   * -----------GROUP 8  SECTION  5 ---------------------------
C--- for W1AD<=GRND--- phase 1 additional velocity (VELAD).
   85 CONTINUE
      GO TO 25
C   *  ----------GROUP 8  SECTION  6 ---------------------------
C--- for W2AD<=GRND--- phase 2 additional velocity (VELAD).
   86 CONTINUE
      GO TO 25
C   * -----------GROUP 8  SECTION 7 ---- VOLUMETRIC SOURCE FOR GALA
C---- Entered when GALA =.TRUE.; block-location index is LSU or LSU2
   87 CONTINUE
      GO TO 25
C   * -----------GROUP 8  SECTION 8 --- CONVECTION COEFFICIENTS
C--- Entered when UCONV =.TRUE.
   88 CONTINUE
c      call fcheck(881,'grex3 ')
C.... Correction applicable to velocities in compressible flow
      IF(COMPRS) THEN                          ! COMPRS is a PIL variable
        IF(INDVAR>=U1.AND.INDVAR<=W2) THEN  ! do this only for velocities
          IF((INDVAR==U1.OR.INDVAR==U2).AND.NDIREC==3) THEN
            CALL GXCMPR(LD11,EAST(P1),P1,1.4)
            CALL GXCMPR(LD12,P1,EAST(P1),1.4)
          ELSEIF((INDVAR==V1.OR.INDVAR==V2).AND.NDIREC==1) THEN
            CALL GXCMPR(LD11,NORTH(P1),P1,1.4)
            CALL GXCMPR(LD12,P1,NORTH(P1),1.4)
          ELSEIF(INDVAR==W1.OR.INDVAR==W2) THEN
            IF(NDIREC==5) THEN
              CALL GXCMPR(EASP20,HIGH(P1),P1,1.4)
            ELSEIF(NDIREC==6) THEN
              CALL GXCMPR(EASP20,P1,HIGH(P1),1.4)
            ENDIF
          ENDIF
        ENDIF
      ENDIF
c      call fcheck(882,'grex3 ')
C
      GO TO 25
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
C     has made his settings.
C
   89 CONTINUE
      IF(CHMKIN) CALL GXCHKI(0,0.0)
      GO TO 25
C   * -----------GROUP 8  SECTION 10 --- CONVECTION NEIGHBOURS
C---- Entered when UCONNE =.TRUE.; block-location indices are LD7
  810 CONTINUE
      GO TO 25
C   * -----------GROUP 8  SECTION 11 --- DIFFUSION NEIGHBOURS
C---- Entered when UDIFNE =.TRUE.; block-location indices are LD7
C     for south, west, high or low, and LD8 for north or east.
C     User should provide INDVAR and NDIREC IF's as above.
  811 CONTINUE
      IF(CHMKIN) CALL GXCHKI(0,0.0)
      GO TO 25
C   * -----------GROUP 8  SECTION 12 --- LINEARISED SOURCES
C---- Entered when USOURC =.TRUE.
  812 CONTINUE
      IF(FDFSOL) CALL GXFDFS
      IF(CHMKIN) CALL GXCHKI(0,0.0)
      GO TO 25
C   * -----------GROUP 8  SECTION 13 --- CORRECTION COEFFICIENTS
C---- Entered when UCORCO =.TRUE.
  813 CONTINUE
      GO TO 25
c
C   * -----------GROUP 8  SECTION 14 --- USER'S SOLVER
C---- Entered when USOLVE =.TRUE.;
  814 CONTINUE
c........................ SLVR and CSG3 are EQUIVALENCEd in GRDLOC
C-------------------------    Call Gauss-Seidel solver
      IF(SLVR=='GAUS'.AND..NOT.SLBSOL) CALL GXGAUS(OVRRLX,
     1 LITER(INDVAR),IXMON,IYMON,IZMON,ENDIT(INDVAR),XCYCLE,LUPR1,LUPR3)
C
      USER=.FALSE.
C
      CALLSL = INDVAR==1 .OR. INDVAR==3 .OR. INDVAR==5 .OR.
     1         INDVAR==7 .OR.(INDVAR>=ICNGRB .AND. INDVAR<=ICNGRC)
      IF(CALLSL) THEN
        NZZ=NZ
        IF( SLBSOL ) NZZ=1
        IF( SLVR=='CGGR') THEN
          CALL GXCGGR(NX,NY,NZZ)
        ELSEIF( SLVR=='CRGR') THEN
          CALL GXCRGR(NX,NY,NZZ)
        ENDIF
C....  entry to GXCHKI for N-R solvers
      ELSEIF(CHMKIN) then
        CALL GXCHKI(0,0.0)
      ENDIF
C
      GO TO 25
C   * -----------GROUP 8  SECTION 15 --- Change solution result
C---- Entered when UCORR =.TRUE.; block-location indices are as above.
  815 CONTINUE
C.... entry to GXCHKI for N-R solvers
C
      IF(CHMKIN) CALL GXCHKI(0,0.0)
      GO TO 25
  816 CONTINUE
C   * ------------------- SECTION 16 --- Change DVEL/DPs if UDVDP = .TRUE.
      GO TO 25
C   * Make all other group-8 changes in group 19.
C***********************************************************************
C
C
C--- GROUP 9. Properties of the medium (or media)
C
C... EARTH no longer calls Group 9 and 10 of GREX for properties. Instead
C    it directly calls subroutine SLBPRP (in EARTH) which itself calls the
c    the appropriate one of the following:
c    SLBDEN (in GXDENS.HTM) for DENST1
c    SLBDRP                     CMPRS1
c    SLBTHE                     THRME1
c    SLBVST                     VISTRB
c    SLBVSL                     VISCLM
c    SLBPRL                     PRNLAM
c    SLBTMP                     TEMPR1
c    SLBSPH                     SPEHT1
c    SLBLEN                     MIXLN1
c    SLBDEN                     DENST2
c    SLBDRP                     CMPRS2
c    SLBTHE                     THRME2
c    SLBTMP                     TEMPR2
c    SLBSPH                     SPEHT2
c    SLBLEN                     MIXLN2
c    SLBIVL                     INTVL1
c    SLBFRC                     INTFCO
c    SLBIMS                     INTMTR
c    SLBIDF                     INTRC1
c    SLBVRM                     IVRMCO
c
C    property flag is set to one of the GRND values. However, if a property
C    flag is set to GRND, SLBPRP calls the GROUND subroutine, expecting to
C    find there user-provided coding for the property in question.
C
    9 GO TO 25
C***********************************************************************
C
C
C--- GROUP 10. Inter-phase-transfer processes and properties
C
C... See comment to Group 9.
   10 GO TO 25
C***********************************************************************
C
C
C--- GROUP 11. Initialization of variable or porosity fields
C
   11 IF(NPATCH(1:4)=='IBFC') THEN ! Initial field of uniform
        CALL GXBFC                   ! flow in a curvilinear grid
      ELSEIF(NPATCH(1:6)=='INIPOL') THEN
        CALL GXIPOL
      ELSEIF(NPATCH(1:6)=='ICHEMK') THEN
        IF(CHMKIN) CALL GXCHKI(0,0.0)
      ELSEIF(NPATCH(1:4)=='BLIN'.OR.NPATCH(1:4)=='GXBL') THEN  ! Boundary-layer profile
        CALL GXBLIN   ! Boundary-layer inlet profile facility
      ENDIF
      GO TO 25
C*****************************************************************
C
C
C--- GROUP 12. Convection and diffusion adjustments
C
   12 CONTINUE
      GO TO 25
C*****************************************************************
C
C
C--- GROUP 13. Boundary conditions and special sources
C
C   * XCYCLE may be changed in group 19.
C
C     Sections 1  to 11 for CO
C     Sections 12 to 22 for VAL
   13 NPAT= NPATCH(1:4)
c.... "Un-comment" call to accutm('grp 13') in three places below, in
c     order to compute and print cumulative time spent in group 13.
c     Such calls can be placed anywhere; but, since accutm itself takes
c     a not-insignificant amount of computer time, they are best
c     "commented" when not needed
c
c      call accutm('grp 13',1)
c
      IF(ISVROB) THEN ! Momentum sources due to swirling fan object
        IF(NAMGRD/='HTBX'.AND.NAMGRD/='FLAR') CALL SWIRFAN
      ENDIF
      IF(NPATCH(1:4)=='BLIN'.OR.NPATCH(1:4)=='GXBL') THEN  ! Boundary-layer inlet profile
        CALL GXBLIN
      ELSEIF(NPATCH(1:3)=='FIR') THEN
        CALL GXFIRE                                    ! FIRE
      ELSEIF(NPATCH(1:6)=='DENDIF') THEN
        IF(INDVAR>2.AND.INDVAR<9) CALL GXDENDIF  ! density difference
      ELSEIF(NPAT=='KESO') THEN ! Sources of turbulence energy and
                                  !dissipation rate   GXKESO
c.... "Un-comment" call to accutm('ke sor') in three places below, in
c     order to compute and print cumulative time spent in GXKESO
c        call accutm('ke sor',1)
        CALL GXKESO(VIST,LEN1,VISL,FIXVAL)
        IF(.NOT.LSG4.AND.ISG60==1.AND.ISG62==0)
     1    CALL PB2_KESO(VIST,LEN1,VISL,FIXVAL)
c        call accutm('ke sor',2)
c        if(isweep==lsweep-1.and.istep==lstep)
c     1        call accutm('ke sor',1)
      ELSEIF(NPAT=='KEBU') THEN
        CALL GXKEGB
C
      ELSEIF(NPAT=='RNGM') THEN ! Extra source of EP in RNG k-e model
        IF(INDVAR==EP) CALL GXRNGM(VIST,LEN1,FIXFLU)
      ELSEIF(NPAT=='KECH') THEN ! Extra source of EP in Chen-Kim k-e model
        IF(INDVAR==EP.AND.ISC==16) THEN
          L0G=L0F(GEN1); L0K=L0F(KE); L0VAL=L0F(VAL); L0VIST=L0F(VIST)
          L0CO=L0F(CO)
          DOIT=IENUTA==3.OR.IENUTA==4
          DO I=1,NXNY
            IF(SLD(I).OR.F(KAP+I)>FIXVAL) THEN
              F(L0VAL+I)=0.0; F(L0CO+I)=0.0
            ELSE
              F(L0VAL+I)=0.25*F(L0G+I)*F(L0G+I)*F(L0VIST+I)*F(L0VIST+I)
     1                                                  /(F(L0K+I)+TINY)
              IF(DOIT) F(L0VAL+I)=F(L0VAL+I)*F(L0FONE+I)
            ENDIF
          ENDDO
        ENDIF
      ELSEIF(NPAT=='REKE') THEN
        IF(IENUTA==14) CALL GX_REALISABLE_KE
      ELSEIF(NPAT=='KWSO') THEN
        IF(IENUTA==10.OR.IENUTA==11) CALL GXKW_WILCOX
        IF(IENUTA==15.OR.IENUTA==16) CALL GXKW_WILCOX
        IF(IENUTA.GE.17.AND.IENUTA.LE.20) CALL GXKW_MENTER
      ELSEIF(NPAT=='TSKE') THEN
        IF(IENUTA==7) CALL GXTSKE(VIST,LEN1,FIXFLU)
      ELSEIF(NPAT=='EYAP') THEN
        IF(ILTLS>0) CALL GXEYAP
      ELSEIF(NPAT=='BUOY') THEN  ! Sources of momentum caused by buoyancy              GXBUOY
        CALL GXBUOY(BFC,CARTES,PARAB,DEN1,DEN2)
      ELSEIF(NPAT=='STEN') THEN  ! Sources of momentum caused by surface tension       GXSURFT
        CALL GXSURFT
      ELSEIF(NPAT=='CLDA') THEN  ! Conservative low-dispersion algorithm               GXCLDA
        CALL GXCLDA
      ELSEIF(NPAT=='HOCS') THEN  ! Higher-order schemes for convection  GXHOCS
        CALL GXHOCS
C.... Exit boundary condition for Supersonic flow Z-direction.
      ELSEIF(NPAT=='OUTL') THEN
        IF(INDVAR==P1) THEN
          IF(ISC==13) THEN
            INDW=ITWO(W1,LBNAME('WCRT'),.NOT.BFC)
            CALL FN21(VAL,LOW(INDW),LOW(DEN1),0.0,-0.995)
          ENDIF
        ENDIF
C
C.... The following call to                               GXNEPA
C     allows the "values" used in a COVAL command to be chosen as
C     the "neighbour-values" of the cell in specified space, time or
C     "variable-space" directions. The condition for the call is that
C     the PATCH name should begin with the characters NE.
      ELSEIF(NPAT(1:2)=='NE') THEN
        CALL GXNEPA(NPAT)
C
C.... The following call to                               GXPOLR
C     fixes the u- and v-velocities at the circumferential boundary of a
C     cylindrical-polar domain for uniform flow at speed POLRA
      ELSEIF(NPAT(2:4)=='POL') THEN
        IF(.NOT.CARTES) CALL GXPOLR(NPAT(1:2))
C
C.... Sources used to set profiles                         GXPROF
      ELSEIF(NPAT=='PROF') THEN
        CALL GXPROF(DEN1,RHO1,FIINIT(H1))
C
C.... Calculate radiation sources .                       GXRADI
      ELSEIF(NPAT=='RADI') THEN
        CALL GXRADI(CARTES,TEMP1,DEN1,RHO1)
C
C.... Centrifugal & Coriolis forces due to rotation of
C     coordinate system about axis of cylindrical-polar
C     system.                                            GXROTA
      ELSEIF(NPAT=='ROTA' .AND. .NOT.GCV) THEN
        IF(.NOT.CARTES.OR.BFC) CALL GXROTA(BFC,NX)
C
C.... Velocity resolutes at a curved inlet boundary       GXBFC
      ELSEIF(NPAT(1:3)=='BFC') THEN
        CALL GXBFC
C
C.... Time-varying sources                                GXTIM
      ELSEIF(NPAT(1:3)=='TIM') THEN
        IF(NEZ(TIMA)) CALL GXTIM
C
C.... Sources provided for single-slab solution of fully-
C     developed flow                                      GXFDFS
      ELSEIF(NPAT(1:3)=='FDF') THEN
        CALL GXFDFS
C
C.... Chemical-reaction source                            GXCHSO
      ELSEIF(NPAT=='CHSO') THEN
        IF(.NOT.EXPLOS) THEN
          CALL GXCHSO(TEMP1)
        ENDIF
C
C.... Scalar fluctuations source                          GXFSQS
      ELSEIF(NPAT=='FSQS') THEN
        CALL GXFSQS
C.... Chemical-reaction sources, inlet properties and
C     something else                      - from CHEMKIN  GXCHKI
      ELSEIF(NPATCH(1:5)=='CHEMK'.OR.NPATCH(1:4)=='NOCP'.OR.
     1  NPATCH(1:6)=='CHEMID'.OR.NPATCH(1:5)=='CKWAL') THEN
        IF(CHMKIN) CALL GXCHKI(0,0.0)
C.... Finite-rate chemical source terms for extended SCRS model
      ELSEIF(ESCRS.AND.NPAT=='SCRS') THEN
        CALL GXSCRI
!!!      ELSEIF(LSG87) THEN
!!!        CALL GXTCALC(0,0.0)
C....                                      Two-phase-flow patches
      ELSEIF(.NOT.ONEPHS) THEN
C.... Momentum source caused by lateral gravity
C     in 2-phase flow                                     GXLATG
        IF(NPAT=='LATG') THEN
          CALL GXLATG
C.... Tentative length-scale source for two-fluid turbulence model
        ELSEIF(NPAT=='LESO') THEN
          CALL GXLESO('SOUTH')
C.... Shear source of lateral velocity for the 2-fluid
C     turbulence model                                     GXSHSO
        ELSEIF(NPAT=='SHSO') THEN
          CALL GXSHSO(INTFRC,LEN1)
C.... Sources providing for upwinding of the interphase
C     transport terms                                     GXIPST
        ELSEIF(NPAT(1:4)=='IPST') THEN
          CALL GXIPST(INTFRC,CINT(INDVAR),NPATCH)
C.... Sources of momentum for interfacial pressure        GXINTP
        ELSEIF(NPAT=='INTP') THEN
          IF(NEZ(CPIP).AND.ISC==16) CALL GXINTP
C.... Sources of momentum for interfacial lift            GXLIFT
        ELSEIF(NPAT=='LIFT') THEN
          IF(NEZ(CLIFT).AND.ISC==16) CALL GXLIFT
C.... Turbulence modulation due to presence of particles  GXDISP
        ELSEIF(NPAT=='KEDI') THEN
          CALL GXDISP
        ENDIF
      ENDIF
      IF(ISC==13.AND.INDVAR>=U1.AND.INDVAR<=W2) THEN
C.... Deduced inflow velocity at fixed-pressure outlet    GXOUTLET
        CALL GXOUTLET
      ENDIF
C.... Angled inlet and outlet                            GXIO
      IF(LGXIO) CALL GXIO
      IF(LSUN) CALL GXSUN
C
      IF(LSG62) CALL GXNOX
      IF(LSG63) CALL GXNOX
c      call accutm('grp 13',2)
c      if(isweep==lsweep-1.and.istep==lstep) call accutm('grp 13',3)
      GO TO 25
C***************************************************************
C
C
C--- GROUP 14. Downstream pressure for PARAB=.TRUE.
C
   14 CONTINUE
      GO TO 25
C***************************************************************
C   * Make changes for these groups 15, 16, 17 & 18 only in group 19.
C***************************************************************
C
C
C--- GROUP 19. Special calls to GROUND from EARTH
C
   19 GO TO (191,192,193,194,195,196,197,198,199,1910,1911),ISC
C
C   * ----------GROUP 19  SECTION 1 ---- START OF TIME STEP.
  191 CONTINUE
C.... Print or compute new grid and geometry
      IF(MOVBFC) CALL GXBFGR
C.... Call to GXMOBS to modify blockages:
C      * either to introduce complex initial field for steady cases;
C      * or moving obstacles/blockages for unsteady cases; or both.
      IF(IPORIB==1.OR.IPORIB==2) CALL GXMOBS
C.... Multiply DYG2D by EPOR, for thin-film lubrication
      IF(ADJEPR) CALL FN26(DYG2D,EPOR)
C.... Call GXPIST to expand and contract the grid in accordance
C     with the kinematics of crankshaft-connecting-rod mechanisms.
      IF(NEZ(W1AD)) THEN
        IF(NEZ(AZW1))
     1  CALL GXPIST(ISTEP,NZ,TIM,DT,ISWEEP,LSWEEP,NPRINT,NTPRIN)
      ENDIF
      IF(RESET.AND..NOT.EXPLOS) THEN
        CALL GXRSET
      ENDIF
      IF(.NOT.STEADY.AND.LSUN) CALL GXSUN
      IF(.NOT.STEADY.AND.BLIN) CALL GXBLIN
C
      GO TO 25
c   * ----------group 19  section 2 ---- start of sweep.
  192 CONTINUE
      IF(SURFT) CALL GXSURFT
      IF(MOVBFC) CALL GXBFGR
      IF(CHMKIN) CALL GXCHKI(0,0.0)
      IF(.NOT.DONE192.AND.IDISPA/=0.AND.ISWEEP==1) THEN ! dumping on start of sweep 1
        I1=ITWO(3,2,PNAM(1:2)=='SW')
        IF(PNAM(I1:I1+1)=='IN') THEN ! PNAM ends in IN so dump initial field
          IF(STEADY.OR.(.NOT.STEADY.AND.(I1==3.OR.(I1==2.AND.
     1                                              ISTEP==1)))) THEN
            ISTP0=FSTEP ; ISWP0=FSWEEP ; FSTEP=0 ; FSWEEP=0
            ISTP1=ISTEP ; ISWP1=ISWEEP ; ISTEP=0 ; ISWEEP=0
            TIM0=TIM; TIM=0.0
            CALL DUMPS(PNAM(1:1),GNAM(1:1),0,1,0,1)
            FSTEP = ISTP0 ; FSWEEP=ISWP0 ; DONE192=.TRUE.
            ISTEP = ISTP1 ; ISWEEP=ISWP1 ; TIM=TIM0
          ENDIF
        ENDIF
      ENDIF
      GO TO 25
C   * ----------GROUP 19  SECTION 3 ---- START OF IZ SLAB.
C.... Modify porosities
  193 IF(IPORIA/=0) THEN
        CALL GXPORA(SETPOR)
        IF(BFC) THEN
          CALL BGEOM(1)
          CALL BGEOM(2)
        ENDIF
      ENDIF
C
      IF(LBMARK/=0) L0MARK= L0F(LBMARK)
      GO TO 25
C
C   * ------------------- SECTION 4 ---- Start of iterations over slab.
C
 194  CONTINUE
      IF(SURFT) CALL GXSURFT
      IF(KEBUOY) CALL GXKEGB
C.... Set DOSKIN to true for WALL patches and correct G-function
      IF(TURB) THEN
        IF(NMWALL>0) CALL GXWFUN
      ENDIF
C.... GXBLIN compute PASQUILL-F profiles for buoyancy forces
      IF(PASQBUOY) CALL GXBLIN
      IF(IENUTA==10.OR.IENUTA==11) CALL GXKW_WILCOX
      IF(IENUTA==15.OR.IENUTA==16) CALL GXKW_WILCOX
      IF(IENUTA.GE.17.AND.IENUTA.LE.20) CALL GXKW_MENTER
Call realisable k-e model after upgrading of G-function in E2conn for walls
      IF(IENUTA.EQ.14) CALL GX_REALISABLE_KE
C.... Calculate phase-average z-direction velocities
C     Copy coding if velocities for other directions are required
      IF(.NOT.ONEPHS) THEN
        IF(LBWAVE>0) THEN
          L0WAVE=L0F(LBWAVE)
          L0R1=L0F(9); L0R2=L0F(10)
          L0W1=L0F(7); L0W2=L0F(8)
          RH1=1.0/RHO1A; RH2=1.0/RHO2A
          DO 1940 I=1,NXNY
            F(L0WAVE+I)=(F(L0W1+I)*F(L0R1+I)*RH1 +
     1                   F(L0W2+I)*F(L0R2+I)*RH2) /
     1                  (F(L0R1+I)*RH1 + F(L0R2+I)*RH2)
 1940     CONTINUE
          if(dbgrnd) call prn('wave',lbwave)
        ENDIF
C
        IF(STORE(P2)) THEN
          IF(STORE(VPOR)) THEN
C.... Extra second-phase pressure proportional to volume fraction
C     divided by porosity.
            CALL FN15(P2,R2,VPOR,0.0,PHS2PA)
          ELSE
            IF(PHS2P) THEN
C.... Extra second-phase pressure proportional to:
C     a*(r2-r2crit)**b if r2 gt r2crit.
              L0R2=L0F(R2)
              L0P2=L0F(P2)
              RLX=ABS(DTFALS(P2))
              DO 1941 I=1,NXNY
                IF(F(L0R2+I)>PHS2PA) THEN
                  GP2=PHS2PB*(F(L0R2+I)-PHS2PA)**PHS2PC
                ELSE
                  GP2=0.0
                ENDIF
                F(L0P2+I)=RLX*GP2+(1.0-RLX)*F(L0P2+I)
 1941         CONTINUE
            ELSE
C.... Extra second-phase pressure proportional to volume fraction
              CALL FN2(P2,R2,0.0,PHS2PA)
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      IF(LSG62) CALL GXNOX
      IF(LSG63) CALL GXNOX
C.... Calculate derivatives of phase 1 temperature
      IF(ITEM1>0.OR.TEMP1>0) THEN
        LBSCAL= ITWO(ITEM1,TEMP1,ITEM1>0); L0SCAL=L0F(LBSCAL)
        IF(ISDTDX.AND.NX>1) THEN
          L0DTDX= L0F(LBDTDX)
          CALL FN1(-L0DTDX,0.0)
          CALL FNDFDX(L0SCAL,L0DTDX,IZSTEP,.TRUE.)
        ENDIF
        IF(ISDTDY.AND.NY>1) THEN
          L0DTDY= L0F(LBDTDY)
          CALL FN1(-L0DTDY,0.0)
          CALL FNDFDY(L0SCAL,L0DTDY,IZSTEP,.TRUE.)
        ENDIF
        IF(ISDTDZ.AND.NZ/=1.AND..NOT.PARAB) THEN
          L0DTDZ= L0F(LBDTDZ)
          CALL FN1(-L0DTDZ,0.0)
          CALL FNDFDZ(L0SCAL,L0DTDZ,IZSTEP,.TRUE.)
        ENDIF
      ENDIF
C.... Calculate derivatives of phase 2 temperature
      IF(ITEM2>0.OR.TEMP2>0) THEN
        LBSCAL= ITWO(ITEM2,TEMP2,ITEM2>0); L0SCAL=L0F(LBSCAL)
        IF(ISDTX2.AND.NX>1) THEN
          L0DTDX= L0F(LBDTX2)
          CALL FN1(L0DTDX,0.0)
          CALL FNDFDX(L0SCAL,L0DTDX,IZSTEP,.TRUE.)
        ENDIF
        IF(ISDTY2.AND.NY>1) THEN
          L0DTDY= L0F(LBDTY2)
          CALL FN1(L0DTDY,0.0)
          CALL FNDFDY(L0SCAL,L0DTDY,IZSTEP,.TRUE.)
        ENDIF
        IF(ISDTZ2.AND.NZ/=1.AND..NOT.PARAB) THEN
          L0DTDZ= L0F(LBDTZ2)
          CALL FN1(L0DTDZ,0.0)
          CALL FNDFDZ(L0SCAL,L0DTDZ,IZSTEP,.TRUE.)
        ENDIF
      ENDIF
      GO TO 25
 1911 CONTINUE
C   * ------------------- SECTION 11---- After calculation of convection
C                                   fluxes for scalars, and of volume
C                                   fractions, but before calculation of
C                                   scalars or velocities
      IF(.NOT.ONEPHS) THEN
        IF(NEZ(CPIP))  CALL GXINTP
        IF(NEZ(CLIFT)) CALL GXLIFT
      ENDIF
      GO TO 25
  199 CONTINUE
C   * ------------------- SECTION 9 ---- Start of solution sequence for
C                                                          a variable
      IF(FSQSOR.AND.INDVAR==LBNAME('FSQ')) CALL GXFSQS
      GO TO 25
 1910 CONTINUE
C   * ------------------- SECTION 10---- Finish of solution sequence for
C                                                          a variable
      IF(CHMKIN) CALL GXCHKI(0,0.0)
      GO TO 25
C
C   * ------------------- SECTION 5 ---- Finish of iterations over slab.
  195 IF(POTVEL) THEN   ! Compute velocities from potential differences.
        IF(ISWEEP/=1) THEN
          CALL GXPOTV(EPOR,NPOR,HPOR,NZ,DZG,BFC)
          IF(POTCMP.AND.ISWEEP>2)
     1             CALL GXPOTC(EPOR,NPOR,HPOR,PORIA,PORIB,NZ) ! allow for
        ENDIF                                                 ! compressibility flow.
      ENDIF
      IF(PSIVEL) CALL GXPSIV(BFC)  ! deduce velocities from stream function psi
      IF(FDFSOL) CALL GXFDFS       ! fully-developed flow
      IF(CHMKIN) CALL GXCHKI(0,0.0) ! chemkin
C....                                Turbulence-energy generation rates
      IF(LBGENK/=0) CALL FN21(LBGENK,GEN1,VIST,0.0,1.0) ! phase 1
      IF(LBGNK2/=0) CALL FN21(LBGNK2,GEN2,VIST,0.0,1.0) ! phase 2
      IF(LBVLSQ/=0) CALL FNVLSQ(LBVLSQ,1) ! sum of squares of phase1 velocity
      IF(LBVABS/=0) THEN          ! absolute velocity, i.e. SQRT(VLSQ)
        IF(LBVLSQ/=0) THEN
          CALL FN30A(LBVABS,LBVLSQ)  ! VABS is SQRT(VLSQ)
        ELSE
          CALL FNVLSQ(LBVABS,1)     ! put VLSQ in VABS, then take square root
          CALL FN30A(LBVABS,LBVABS)  ! VABS is SQRT(VLSQ)
        ENDIF
        IF(ISCREY) THEN
          L0VOL =L0F(VOL)    ; L0VIST=L0F(VIST)  ; L0VISL=L0F(VISL)
          L0CREY=L0F(LBCREY) ; L0VABS=L0F(LBVABS)
          DO I=1,NXNY
            IF(.NOT.SLD(I)) F(L0CREY+I)=F(L0VABS+I)*F(L0VOL+I)**0.3333/
     1                        (F(L0VIST+I)+F(L0VISL+I))
          ENDDO
        ENDIF
      ENDIF
      IF(LBVAB2/=0) THEN          ! absolute velocity, phase 2
          CALL FNVLSQ(LBVAB2,2)     ! put VLSQ in VABS, then take square root
          CALL FN30A(LBVAB2,LBVAB2)  ! VABS is SQRT(VLSQ)
      ENDIF
      IF(ISTREY) THEN
        L0VIST=L0F(VIST); L0VISL=L0F(VISL); L0TREY=L0F(LBTREY)
        DO I=1,NXNY
          IF(.NOT.SLD(I)) F(L0TREY+I) = F(L0VIST+I)/F(L0VISL+I)
        ENDDO
      ENDIF
      IF(ESCRS) CALL GXSCRI
C....Mach number calculation for Compressible Flow; it must be
c    done here as parab=t field printout precedes gr 19 sc 6
      IF(PARAB) THEN
        IF(LBMACH/=0.OR.LBPTOT/=0)
     1                            CALL GXMACH(LBMACH,LBPTOT,LBVABS,0.,1)
        IF(.NOT.ONEPHS) THEN
          IF(LBMAC2/=0.OR.LBPTO2/=0)
     1                            CALL GXMACH(LBMAC2,LBPTO2,LBVAB2,0.,2)
          IF(SLIPVL) THEN  ! conditionally compute the slip velocities
            IF(LBSLU/=0) CALL FN10(LBSLU,U1,U2,0.0,1.0,-1.0)
            IF(LBSLV/=0) CALL FN10(LBSLV,V1,V2,0.0,1.0,-1.0)
            IF(LBSLW/=0) CALL FN10(LBSLW,W1,W2,0.0,1.0,-1.0)
          ENDIF
        ENDIF
      ENDIF
      GO TO 25
C
C GROUP 19  SECTION 6  Finish of iz slab
  196 IF(PARAB) THEN
        IF(IDISPA>0) CALL GXPARA
      ELSEIF(ARRIVT) THEN
        IF(ISWEEP==LSWEEP) CALL GXARRT
      ENDIF
C....Mach number calculation for Compressible Flow
      IF(.NOT.PARAB) THEN
        IF(LBMACH/=0.OR.LBPTOT/=0)
     1                            CALL GXMACH(LBMACH,LBPTOT,LBVABS,0.,1)
        IF(.NOT.ONEPHS) THEN
          IF(LBMAC2/=0.OR.LBPTO2/=0)
     1                            CALL GXMACH(LBMAC2,LBPTO2,LBVAB2,0.,2)
C.... For two-phase flow, conditionally compute the slip velocities
          IF(SLIPVL) THEN
            IF(LBSLU/=0) CALL FN10(LBSLU,U1,U2,0.0,1.0,-1.0)
            IF(LBSLV/=0) CALL FN10(LBSLV,V1,V2,0.0,1.0,-1.0)
            IF(LBSLW/=0) CALL FN10(LBSLW,W1,W2,0.0,1.0,-1.0)
          ENDIF
        ENDIF
      ENDIF
      IF(LSG62) CALL GXNOX
      IF(LSG63) CALL GXNOX
      IF(BLIN) CALL GXBLIN   ! Boundary-layer inlet profile facility
      IF(IZ0) CALL DMPDSP
      ENDIF
      GO TO 25
C
C GROUP 19  SECTION 7  FINISH OF SWEEP.
  197 CONTINUE
C
C.... Sequence permitting intervention to terminate execution
C     when a file named STOPJOB exists, for example created by
C     editing in a second window
C     Note that the name stopjob will not suffice on case-sensitive
C     systems such as UNIX
C     Note that this may not work on all computer systems, in which
C     case the coding from here to STOPJOB END below should be
C     deactivated.
      IERR=-1
      ISTOP=5
      IF(MOD(ISWEEP,ISTOP)==0) THEN
        IF(MIMD.AND.NPROC>1) THEN
          IF(MYID==0) THEN
            IF(LUNIT(26)/=0) CALL OPENZZ(26,IERR)
          ENDIF
c... pass first proc IERR value to others
          CALL RTFSII(IERR)
        ELSE
          IF(LUNIT(26)/=0) CALL OPENZZ(26,IERR)
        ENDIF
      ENDIF
      IF(IERR==0) THEN
        ENUFSW=.TRUE.
        LSTEP=ISTEP
        LSWEEP=ISWEEP+1
        LUNIT(26)=0
        CALL WRIT40('file stopjob found; execution  stopped  ')
        CALL WRYT40('file stopjob found; execution  will be  ')
        CALL WRYT40('stopped after completion of print-out   ')
      ENDIF
c.... stopjob end
      IF(CHMKIN) CALL GXCHKI(0,0.0)
C.... Print-out of Y+, stress and Stanton number at WALL-type patches
      IF(TURB) THEN
        IF(NMWALL>0) THEN
          IF(.NOT.NULLPR.AND.(YPLS .OR. WALPRN)) CALL GXYPLS
        ENDIF
      ENDIF
c
C.... "Time-out" sequence activated by setting the PIL variable ISG20
C     equal to the maximum number of seconds which the run may take.
      IF(ISG20>0) THEN
        CALL SECONDS(NUMSEC)
        IF(NUMSEC>=ISG20) THEN
          CALL WRITST
          CALL WRIT40('Execution "timed out" at numsec seconds')
          CALL WRIT1I('  numsec',ISG20)
          CALL WRITST
          ENUFSW=.TRUE. ; LSTEP=ISTEP
          IF(ISWEEP/=LSWEEP) LSWEEP=ISWEEP+1
        ENDIF
      ENDIF
c
c.... print-out elicited by SPEDAT('PRINT', ...  statements in Q1
      IF(ISWEEP==FSWEEP) THEN
        IF(ISTEP==1) THEN
c.... get number of special-print commands
          NUMSPPR=0
          CALL GETSDI('PRINT','NUMBER',NUMSPPR)
        ENDIF
      ENDIF
      IF(TIMPRN.AND.SWPRIN.AND..NOT.NULLPR) THEN
C
        IF(LBVOXY>0) L0VOXY=M0F1(LBVOXY)
        IF(LBVOYZ>0) L0VOYZ=M0F1(LBVOYZ)
        IF(LBVOZX>0) L0VOZX=M0F1(LBVOZX)
C
        IF(LBVOXY>0.OR.LBVOYZ>0.OR.LBVOZX>0) THEN
          DO IZZ = 1,NZ
            DO IX=1,NX
              DO IY=1,NY
                I3D=IY+(IX-1)*NY+(IZZ-1)*NFM
                IF(LBVOXY/=0) F(L0VOXY+I3D)=VORTXY(1,IX,IY,IZZ)
                IF(LBVOYZ/=0) F(L0VOYZ+I3D)=VORTYZ(1,IX,IY,IZZ)
                IF(LBVOZX/=0) F(L0VOZX+I3D)=VORTZX(1,IX,IY,IZZ)
              ENDDO
            ENDDO
          ENDDO
        ENDIF
      ENDIF
      IF(NUMSPPR>0) THEN
        IF(TIMPRN.AND.SWPRIN.AND..NOT.NULLPR) THEN
          call writ40('Special print-out via spedat commands')
c.... how many commands?
          DO I=1,NUMSPPR
c.... get the command
            CALL GETSDC('PRINT','COMMAND'//DIGIT(I),COMMAND)
c.... carry out the commands
            IF(COMMAND(1:3)=='POW') THEN
c.... do this before the istep test so that result is included in phi
              CALL GETSDR('PRINT','POWER',POWER)
              IF(MOD(ISTEP,NTPRIN)==0) THEN
                CALL WRIT40(COMMAND(10:13)//
     1                      'holds "power-lawed" values of '
     1                      //COMMAND(5:8))
                CALL WRIT1R('power   ',POWER)
              ENDIF
              CALL PUT_POWER( COMMAND(5:8) , COMMAND(10:13),POWER )
              CALL WRITBL
            ELSEIF(MOD(ISTEP,NTPRIN)==0) THEN
              CALL WRIT40('Spedat print command is '//command)
              IF(COMMAND(1:3)=='AVE') THEN
                CALL AVERAGE(AVE,COMMAND(5:8))
                CALL WRIT40
     1          ('mass-average value of         '//COMMAND(5:8))
                CALL WRIT1R('average ',AVE)
                CALL WRITBL
              ELSEIF(COMMAND(1:5)=='TOTAL') THEN
                CALL SUMALL(SUM,COMMAND(7:10))
                CALL WRIT40('total whole-domain amount of       '//
     1                      COMMAND(7:10))
                CALL WRIT1R('TOTAL   ',SUM)
                CALL WRITBL
              ELSEIF(COMMAND(1:6)=='MINMAX') THEN
                CALL WRIT40('hilo3d called for '//COMMAND(8:11))
                CALL HILO3D(LBNAME(COMMAND(8:11)))
                CALL WRIT2R('MinValue',RLOW3D,'MaxValue',HI3D)
                CALL WRITBL
              ELSEIF(COMMAND(1:5)=='VALUE') THEN
                CALL GETSDI('PRINT','IXLOC',IXLOC)
                CALL GETSDI('PRINT','IYLOC',IYLOC)
                CALL GETSDI('PRINT','IZLOC',IZLOC)
                CALL GETVAL(COMMAND(7:10),IXLOC,IYLOC,IZLOC)
              ENDIF
            ENDIF
          ENDDO
        ENDIF
      ENDIF
C.... Dump fields after each IDISPA sweeps, eg for PHOTON viewing
      IF(IDISPA/=0) THEN
        IF(PNAM(1:2)=='SW') THEN
          IF(STEADY) THEN
            CALL DUMPS(PNAM(1:1),GNAM(1:1),ISWEEP,IDISPA,IDISPB,IDISPC)
          ELSE
            STEADY=.TRUE.
            CALL DUMPS(PNAM(1:1),GNAM(1:1),ISWEEP,IDISPA,1,LSWEEP)
            STEADY=.FALSE.
            DONE192=.FALSE.
          ENDIF
          IF(CALFOR) THEN
            IDSPB=ITWO(1,IDISPB,IDISPB==0)
            IDSPC=ITWO(LSWEEP,IDISPC,IDISPC==0)
            IDSPA=ITWO(1,IDISPA,IDISPA==0)
            IF(ISWEEP>=IDSPB.AND.ISWEEP<=IDSPC.AND.
     1                                       MOD(ISWEEP,IDSPA)==0) THEN
              CALL FORCES(1) ! Integrate
              CALL FORCES(3) ! Output to table file
            ENDIF
          ENDIF
        ENDIF
      ENDIF
C.... Angled inlet and outlet                            GXIO
      IF(LGXIO) CALL GXIO
      IF(BLIN) CALL GXBLIN   ! Boundary-layer inlet profile facility
      GO TO 25
C
C GROUP 19  SECTION 8  FINISH OF TIME STEP.
  198 CONTINUE
C.... Save fields and cell-corner coordinates to a series of files,
C     for examination via VR-Viewer or PHOTON
C.... PNAM and GNAM are the names of the field and (if BFC)
C     grid files, and IDISPA is the frequency of dumping.
C     They are EQUIVALENCEd in the included Common file GRDLOC to the
C     PIL variables CSG1 and CSG2, and can therefore be set in the
C     satellite, eq by the line:
C                 CSG1=PHI;CSG2=XYZ;IDISPA=1
C     which will create files p1da or p1, p2da or p2, etc, & x1da or x1,
C     etc according to whether phida and xyzda are true (the default)
C     or false in the PREFIX file
C
!!!      IF(LSG87) CALL GXTCALC(0,0.0)
      IF(PNAM(1:2)/='SW'.AND.IDISPA/=0) THEN
        LGTEMP=BFC ; BFC=BFC.OR.NOTBFC
        CALL DUMPS(PNAM(1:1),GNAM(1:1),ISTEP,IDISPA,IDISPB,IDISPC)
        BFC=LGTEMP
        IF(CALFOR) THEN
          IDSPB=ITWO(1,IDISPB,IDISPB==0)
          IDSPC=ITWO(LSTEP,IDISPC,IDISPC==0)
          IDSPA=ITWO(1,IDISPA,IDISPA==0)
          IF(ISTEP>=IDSPB.AND.ISTEP<=IDSPC.AND.MOD(ISTEP,IDSPA)==0) THEN
            CALL FORCES(1) ! Integrate
            CALL FORCES(3) ! Output to table file
          ENDIF
        ENDIF
C... Save monitor plot image for each time step
        GXMON=(TSTSWP/=0)
        CALL GETSDL('GXMONI','TRANSIENT',GXMON)
        IF(GXMON) THEN
          IFST=ITWO(FSTEP,IDISPB,IDISPB==0) ! first step to dump
          ILST=ITWO(LSTEP,IDISPC,IDISPC==0) ! last step to dump
          GXMON=ISTEP>=IFST.AND.ISTEP<=ILST.AND.
     1                                          MOD(ISTEP,IDISPA)==0
        ENDIF
        IF(GXMON) THEN
          IF(PLTCLSX) THEN
            COMMAND='gxmoni'
            IF    (ISTEP>=100000) THEN
              WRITE (COMMAND(7:12),FMT='(BN,I6)') ISTEP
            ELSEIF(ISTEP>=10000) THEN
              WRITE (COMMAND(7:11),FMT='(BN,I5)') ISTEP
            ELSEIF(ISTEP>=1000) THEN
              WRITE (COMMAND(7:10),FMT='(BN,I4)') ISTEP
            ELSEIF(ISTEP>=100) THEN
              WRITE (COMMAND(7:9),FMT='(BN,I3)') ISTEP
            ELSEIF(ISTEP>=10) THEN
              WRITE (COMMAND(7:8),FMT='(BN,I2)') ISTEP
            ELSE
              WRITE (COMMAND(7:7),FMT='(BN,I1)') ISTEP
            ENDIF
            LL=LENGZZ(COMMAND)
            CALL DUMPZZ(COMMAND(1:LL),LL)
          ELSE
!DEC$ IF .NOT.DEFINED(LINUX) .AND. DEFINED(_X64_)
              !DEC$ ATTRIBUTES C, ALIAS : 'AllStepsONC' :: SDASCON
            !cs interface code
!!            IF(ISTEP==IFST.OR.ISTEP==ILST) THEN
              CALL SDASCON() ! toggle dumping on
 
              !------------------------------------------------------
              !this block is necessary for the super fast cases
              !it pauses earth immediately after sending signal to take screenshots
              !otherwise the data in the gui may skip a few sweeps before it starts the screenshot r
              !corresponding unpause call in abstractions.screenShotAll(special==1)
              !DEC$ ATTRIBUTES C, ALIAS : 'earTogglePause' :: ETP
              CALL ETP()
              !------------------------------------------------------
 
!!            ENDIF
!DEC$ ENDIF
          ENDIF
        ELSE
!DEC$ IF .NOT.DEFINED(LINUX) .AND. DEFINED(_X64_)
          IF(.NOT.PLTCLSX) THEN
!!          IF(ISTEP==IFST.OR.ISTEP==ILST) THEN
              !DEC$ ATTRIBUTES C, ALIAS : 'dumpAllStepsOFFC' :: SDASCOFF
            CALL SDASCOFF() ! toggle dumping  off
!!          ENDIF
          ENDIF
!DEC$ ENDIF
        ENDIF
      ENDIF
C
      IF(MOVBFC) THEN
C....                      Print-out versus time by use of GRAPH
        TIME(ISTEP)=TIM
C....                                        set ordinate values
        OR(ISTEP,1)=VARONE(LBNAME('VOLU'),1,1)
        OR(ISTEP,2)=VARONE(P1,1,1)
        IF(ISTEP/=LSTEP) GO TO 25
C....                                                 call GRAPH
C....                                        set ordinate labels
        LABOR(1)='VOLU' ; LABOR(2)='PRES'
        IOSC(1)=1 ; IOSC(2)=2 ; IOSC(3)=3 ; IOSC(4)=4
        NABDIM=100 ; NORDIM=10 ; NOR=2
        ABSIZ=0.6 ; ORSIZ=0.4
C
        CALL GRAPH(TIME,NABDIM,ISTEP,'TIME',OR,NORDIM,NOR,
     1              LABOR,IOSC,ABSIZ,ORSIZ,3,LUPR1)
C....                                     end of GRAPH sequence
      ENDIF
C
      IF(.NOT.STEADY) THEN
        IF(.NOT.PARAB.AND.IDISPA>0.AND.PNAM=='    ') CALL GXPARA
        IF(IDISPD>0) CALL DMPDSP
      ENDIF
      GO TO 25
C***************************************************************
C
C
C--- GROUP 20. Preliminary print-out
C
   20 IF(ECHO) THEN ! echo all 25 data groups to RESULT file
        CALL DATPRN(Y,Y,Y,Y,  Y,Y,Y,Y,  Y,Y,Y,Y,  Y,Y,Y,Y,
     1              Y,Y,Y,Y,  Y,Y,Y,Y,  Y)
      ELSE          ! echo only data group 1 to RESULT
        CALL DATPRN(Y,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N)
      ENDIF
      GO TO 25
C***************************************************************
C
C--- GROUP 21. Special print-out to screen
   21 CONTINUE
      GO TO 25
C***************************************************************
C   * Make changes for these groups 22 and 23 only in group 19.
C***************************************************************
C
C
C--- GROUP 24. Dumps for restarts
C
   24 CONTINUE
C   * Insert CALL DUMP where required in group 19.
C
   25 CONTINUE
      if(test) call fcheck(25,'grex3 ')
C***************************************************************
C
      IF(SPCLGR) THEN
C.......................................... special GROUNDS
C.... S2SR, for surface-to-surface radiation
        IF(S2SR) THEN
          RADCVD = .FALSE.
          CALL GETSPD('CVD','RADCVD',3,RDUM,IDUM,RADCVD,' ',IERR)
C.... Only call GXS2SR GR.9, SEC 7 when INDVAR is temperature
          IF(INDVAR==ITEM1.OR..NOT.(IGR==9.AND.ISC==7)) THEN
            IF(RADCVD) THEN
              CALL GXS2SR_CVD
            ELSE
              CALL WRIT40('Please set:')
              CALL WRIT40('SPEDAT(SET,CVD,RADCVD,L,T)')
              CALL WRIT40('in Q1')
              CALL SET_ERR(587,'Old S2SR not supported',1)
            ENDIF
          ENDIF
        ENDIF
        IF(NAMGRD=='NONE') THEN
C....NONE
        ELSEIF(NAMGRD=='FURN') THEN
C....FURN
          CALL FURNGR
C.... ESTER, the aluminium-smelter module
        ELSEIF(NAMGRD=='ESTR') THEN
          CALL ESTRGR
C....HOTBOX
        ELSEIF(NAMGRD=='HTBX') THEN
          CALL HTBXGR
C....TACT
        ELSEIF(NAMGRD=='TACT') THEN
          CALL TACTGR
C....COMPRESSIBLE
        ELSEIF(NAMGRD=='CONV') THEN
          CALL GXCONV
C....CVD
        ELSEIF(NAMGRD=='CVD ') THEN
          CALL GXCVD
C....MICA special GROUND's
        ELSEIF(NAMGRD=='MICA') THEN
          CALL MICAGR
C... F1 - calculates drag and lift for formula-1 car
        ELSEIF(NAMGRD=='F1') THEN
          CALL F1GRD
C... Water cooled steam condenser WATSTCON
        ELSEIF(NAMGRD=='WATS') THEN
          CALL GXCOND
        ELSEIF(NAMGRD=='WAVE') THEN
          CALL GXWAVE
        ENDIF
      ENDIF
      if(test) call fcheck(26,'grex3 ')
      NAMSUB='grex3'
      if(dbgloc) then
        call banner(2,namsub,120121)
      endif
      END
c