c
C.... File Name ...... GXGENK.FOR........... 120416
C    File includes the following subroutines and functions:
C       GXGENF, GXSQUR, GFTRM1, GFTRM2, GFDUDX, GFDUDY, GFDUDZ, GFDVDX,
C       GFDVDY, GFDVDZ, GFDWDX, GFDWDY, GFDWDZ, AVRGVL, GXKEGB, AXIBFC.
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... GXGENF calculates derivatives of flow velocities and Generation
C     for a current slab. It is called twice; from Gr.1 Sec.2 to
C     define auxiliary variables; and from Gr.19 Sec.4 to make the
C     calculations. If IGEN0=0 when GXGENF is called from Gr.19 Sec.4,
C     it calculates and stores only DUDX,...,DWDZ. This call is normally
C     done for Reynolds stress turbulence model.
C
C.... The dummy IPH is the phase index, which indicates that the Generation
C     is calculated for the first-phase (IPH=1) or the second-phase (IPH=2)
C     The IGEN0 is the slab-wise address of the Generation storage.
C
C.... There is one /LSG/-logicals temporary in use:
C        LSG53=T activates the use of cell-centered averaged velocities
C                for both 'in-line' and 'off-line' derivatives.
C.... To check calculations of the generation function the user may
C     activate DEBUG 3D-storage of averaged velocities, or velocity
C     derivatives by storing appropriate variables in Q1 (see details
C     below).
 
c
      SUBROUTINE GXGENF(IPH,IGEN0)
      INCLUDE 'farray'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      COMMON /CMNGN/ IUG0(2),IVG0(2),IWG0(2),IUIZM0(2),IUIZP0(2),
     1               IVIZM0(2),IVIZP0(2),IWIZM0(2),IWIZP0(2),
     1               IDUDI0,IDUDJ0,IDUDK0,IDVDI0,IDVDJ0,IDVDK0,
     1               IDWDI0,IDWDJ0,IDWDK0
     1       /GENI/  IG1(2),NXNYST,IG2(6),NFM,IG3(23),ILTLS,IG4(26)
     1       /GENL/  LGL1(14),debgtz,LGL2(43),STRGN1,STRGN2
     1       /NAMFN/NAMFUN,NAMSUB
     1       /BFICRT/IUCRT(2),IVCRT(2),IWCRT(2),ITMP(6)
     1       /GENFL/ LDUDX,LDUDY,LDUDZ,LDVDX,LDVDY,LDVDZ,LDWDX,LDWDY,
     1               LDWDZ,COLVEL,AXISYM  /FNDRL/NXNE1,NYNE1,NZNE1
     1       /LBDFDL/IDUDX,IDUDY,IDUDZ,IDVDX,IDVDY,IDVDZ,IDWDX,IDWDY,
     1               IDWDZ,IDSDX,IDSDY,IDSDZ,IDU2X,IDU2Y,IDU2Z,IDV2X,
     1               IDV2Y,IDV2Z,IDW2X,IDW2Y,IDW2Z
     1       /RSTM2/ J0DUDX,J0DUDY,J0DUDZ,J0DVDX,J0DVDY,J0DVDZ,J0DWDX,
     1               J0DWDY,J0DWDZ,JRSTF(45)
     1       /F0/    LB1(109),KZXCY,LB2(194)
     1       /LSG/   LSGF1(52),LSG53,LSGF2(47)
      LOGICAL LSGF1,LSG53,LSGF2
      LOGICAL NXNE1,NYNE1,NZNE1,COLVEL,AXISYM,LDUDX,LDUDY,LDUDZ,LDVDX,
     1        LDVDY,LDVDZ,LDWDX,LDWDY,LDWDZ,LSTOU,LSTOV,LSTOW,AVEVEL,
     1        DVDL3D,LDOGEN,LGL1,STRGN1,STRGN2,LDBAVR,dbgloc,LGL2,debgtz
      CHARACTER*6 NAMFUN,NAMSUB
      DIMENSION IU1C(2),IV1C(2),IW1C(2)
      SAVE ITMP10,ITMP20,ITMP30,IU1C,IV1C,IW1C
      SAVE LSTOU,LSTOV,LSTOW,AVEVEL,DVDL3D,LDBAVR
C
      NAMSUB = 'GXGENF'
      if(indvar>0) then
        dbgloc=debgtz.and.dbgphi(indvar)
      else
        dbgloc=.false.
      endif
      if(flag.or.dbgloc)  then
        call banner(1,namsub,250516)
        call writ2i('igr     ',igr,'isc     ',isc)
      endif
      IF(IGR==1 .AND. ISC==2) THEN ! Call from Gr.1; Sec.2 to prepare
        AVEVEL=LSG53; AXISYM=.NOT.CARTES
        IF(BFC .AND. NX==1) CALL AXIBFC(AXISYM) ! axisymmetry
             ! (if Y-Z grid is axisymmetric, additional term V/R appears as
             !  part of dU/dX):
C.... Define velocities to use for Generation calculation:
        IF(BFC) THEN
          IUG0(1)=IUCRT(1); IVG0(1)=IVCRT(1); IWG0(1)=IWCRT(1)
          IUG0(2)=IUCRT(2); IVG0(2)=IVCRT(2); IWG0(2)=IWCRT(2)
C.... Temporary slab-storage for DUDI,DUDJ,...,DWDK:
          IDUDI0=L0F(LD4); IDVDI0=L0F(LD7); IDWDI0=L0F(LD10)
          IDUDJ0=L0F(LD5); IDVDJ0=L0F(LD8); IDWDJ0=L0F(LD11)
          IDUDK0=L0F(LD6); IDVDK0=L0F(LD9); IDWDK0=L0F(LDSTAG)
        ELSEIF(CCM) THEN
          IUG0(1)=LBNAME('UC1'); IUG0(2)=LBNAME('UC2')
          IVG0(1)=LBNAME('VC1'); IVG0(2)=LBNAME('VC2')
          IWG0(1)=LBNAME('WC1'); IWG0(2)=LBNAME('WC2')
        ELSE
C.... For staggered solver on Cartesian geometry there are two options:
C      1. AVEVEL=F (default) Use face-centered velocities U1,...,W1 for
C                 the 'in-line' derivatives (dU/dX,...) and averaged
C                 cell-centered velocities for the 'off-line' (dU/dY..).
C      2. AVEVEL=T Use averaged cell-centered velocities for the both
C                 'in-line' and 'off-line' derivatives.
C     Storage for averaged velocities at the current and high slabs is
C     allocated in EARTH. LD4,LD5 and LD6 are used as temporary storage
C     for low slab.
          LSTOU=STORE(3); LSTOV=STORE(5); LSTOW=STORE(7)
          IUG0(1)=IUIZM0(1);  IVG0(1)=IVIZM0(1); IWG0(1)=IWIZM0(1)
          IUIZM0(1)=L0F(LD4); IVIZM0(1)=L0F(LD5); IWIZM0(1)=L0F(LD6)
C.... Debug 3D-storage of averaged velocities:
          IU1C(1)= LBNAME('U1C'); IV1C(1)= LBNAME('V1C')
          IW1C(1)= LBNAME('W1C')
          IU1C(2)= LBNAME('U2C'); IV1C(2)= LBNAME('V2C')
          IW1C(2)= LBNAME('W2C')
          LDBAVR= IU1C(1)/=0 .OR. IV1C(1)/=0 .OR. IW1C(1)/=0 .OR.
     1            IU1C(2)/=0 .OR. IV1C(2)/=0 .OR. IW1C(2)/=0
          IUG0(2)=0; IVG0(2)=0; IWG0(2)=0
          IF(.NOT.ONEPHS) THEN
          IUG0(2)=IUIZM0(2);  IVG0(2)=IVIZM0(2); IWG0(2)=IWIZM0(2)
          ENDIF
        ENDIF
        IF(BFC.AND.KZXCY/=0) CALL XCYVEL(IGR,IZ,IUSL0,IVSL0,IWSL0)
C.... COLVEL indicates that velocities are cell-centered
        COLVEL= BFC .OR. CCM
C.... Find zero initial addresses for colocated velocities:
        IF(COLVEL) THEN
          IF(IUG0(1)/=0) IUG0(1)= L0F(IUG0(1))
          IF(IVG0(1)/=0) IVG0(1)= L0F(IVG0(1))
          IF(IWG0(1)/=0) IWG0(1)= L0F(IWG0(1))
          IF(IUG0(2)/=0) IUG0(2)= L0F(IUG0(2))
          IF(IVG0(2)/=0) IVG0(2)= L0F(IVG0(2))
          IF(IWG0(2)/=0) IWG0(2)= L0F(IWG0(2))
          LSTOU= IUG0(1)/=0
          LSTOV= IVG0(1)/=0
          LSTOW= IWG0(1)/=0
        ENDIF
C.... Define auxiliary logicals:
        LDUDX=.FALSE.; LDUDY=.FALSE.; LDUDZ=.FALSE.
        LDVDX=.FALSE.; LDVDY=.FALSE.; LDVDZ=.FALSE.
        LDWDX=.FALSE.; LDWDY=.FALSE.; LDWDZ=.FALSE.
        IF(.NOT.PARAB) THEN
          LDUDX= LSTOU .AND. (DUDX .OR. GENK .OR. RSTM) .AND. NXNE1
          LDUDX= LDUDX .OR. (AXISYM .AND. LSTOV)
          LDUDY= LSTOU .AND. (DUDY .OR. GENK .OR. RSTM) .AND. NYNE1
          LDUDZ= LSTOU .AND. (DUDZ .OR. GENK .OR. RSTM) .AND. NZNE1
          LDVDX= LSTOV .AND. (DVDX .OR. GENK .OR. RSTM) .AND. NXNE1
          LDVDY= LSTOV .AND. (DVDY .OR. GENK .OR. RSTM) .AND. NYNE1
          LDVDZ= LSTOV .AND. (DVDZ .OR. GENK .OR. RSTM) .AND. NZNE1
          LDWDX= LSTOW .AND. (DWDX .OR. GENK .OR. RSTM) .AND. NXNE1
          LDWDY= LSTOW .AND. (DWDY .OR. GENK .OR. RSTM) .AND. NYNE1
          LDWDZ= LSTOW .AND. (DWDZ .OR. GENK .OR. RSTM) .AND. NZNE1
        ELSE
          LDWDX= LSTOW .AND. (DWDX .OR. GENK .OR. RSTM) .AND. NXNE1
          LDWDY= LSTOW .AND. (DWDY .OR. GENK .OR. RSTM) .AND. NYNE1
        ENDIF
C.... Define addresses for temporary storage (LD1,LD2,LD3 are in use):
        ITMP10=L0F(LD1); ITMP20=L0F(LD2); ITMP30=L0F(LD3)
C.... Store velocity derivatives for the dump into RESULT file:
        DVDL3D= IDUDX/=0 .OR. IDUDY/=0 .OR. IDUDZ/=0 .OR.
     1          IDVDX/=0 .OR. IDVDY/=0 .OR. IDVDZ/=0 .OR.
     1          IDWDX/=0 .OR. IDWDY/=0 .OR. IDWDZ/=0 .OR.
     1          IDU2X/=0 .OR. IDU2Y/=0 .OR. IDU2Z/=0 .OR.
     1          IDV2X/=0 .OR. IDV2Y/=0 .OR. IDV2Z/=0 .OR.
     1          IDW2X/=0 .OR. IDW2Y/=0 .OR. IDW2Z/=0
C.... For Reynolds stress model DUDX,...,DWDZ are put into the storage
C     provided in UST191 subroutine:
        DVDL3D= DVDL3D .AND. .NOT.RSTM
      ELSE ! Call from Gr.19; Sec.4 to calculate Generation or DUDX,...,etc.
C-----------------------------------------------------------------------
C.... If GEN1 (or GEN2) is in 3D-storage don't calculate it again when
C     Gr.19 Sec.4 is entered for variables solved whole-field:
        IF( (IPH==1.AND.STRGN1 .OR. IPH==2.AND.STRGN2) .AND.
     1      INDVAR>10 .AND. INDVAR/=ILTLS ) RETURN
C.... If IGEN0/=0 calculate generation term, otherwise only DUDX,...
        LDOGEN= IGEN0/=0
        IF(GCV) THEN
          CALL BFGENB(IGEN0)
          RETURN
        ENDIF
C.... Calculate Generation for MB-FGE technique
        IF(CCM .AND. NUMBLK>0) THEN
          CALL UNGENB(IPH-1,IGEN0)
          RETURN
        ENDIF
C.... Recalculate geometric coefficients for moving BFC-grids:
        IF(MOVBFC) CALL IJKXYZ(IZSTEP)
        IF(.NOT.COLVEL) THEN
C.... Average staggered velocities:
          IF(LSTOU) THEN
            IF(NXNE1 .AND. (LDUDX.OR.LDUDY.OR.LDUDZ)) THEN
              IF(LDUDZ) THEN
                IF(IZSTEP>1 ) THEN
                CALL SHINX2(IUIZM0(IPH),IUG0(IPH),IUG0(IPH),IUIZP0(IPH))
                ELSE
                  CALL AVRGVL(4,IUG0(IPH),IPH,IZSTEP,0)
                ENDIF
                IF(IZSTEP1) THEN
                CALL SHINX2(IVIZM0(IPH),IVG0(IPH),IVG0(IPH),IVIZP0(IPH))
                ELSE
                  CALL AVRGVL(2,IVG0(IPH),IPH,IZSTEP,0)
                ENDIF
                IF(IZSTEP1 ) THEN
                CALL SHINX2(IWIZM0(IPH),IWG0(IPH),IWG0(IPH),IWIZP0(IPH))
                ELSE
c w average for current slab
                  CALL AVRGVL(6,IWG0(IPH),IPH,IZSTEP,0)
                ENDIF
c w average for high slab
                IF(IZSTEP
      SUBROUTINE GFTRM1(IGEN0,A1,ID10,A2,ID20,A3,ID30)
      INCLUDE 'farray'
      COMMON      /GENI/IG1(2),NXNYST,IG2(57)
      I1M0= ID10-IGEN0
      I2M0= ID20-IGEN0
      I3M0= ID30-IGEN0
      DO I= IGEN0+1,IGEN0+NXNYST
        F(I)= 2.*(A1*F(I1M0+I)**2 + A2*F(I2M0+I)**2 + A3*F(I3M0+I)**2)
      ENDDO
      END
C-----------------------------------------------------------------------
c
      SUBROUTINE GFTRM2(IGEN0,A1,ID10,A2,ID20)
      INCLUDE 'farray'
      COMMON      /GENI/IG1(2),NXNYST,IG2(57)
      LOGICAL QGE
      I1M0= ID10-IGEN0
      I2M0= ID20-IGEN0
      DO I= IGEN0+1,IGEN0+NXNYST
        IF(QGE(ABS(A1*F(I1M0+I)),1.E10).OR.QGE(ABS(A2*F(I2M0+I)),1.E10))
     1                                                              THEN
          F(I)=1.E10
        ELSE
          F(I)= F(I) + ( A1*F(I1M0+I) + A2*F(I2M0+I) )**2
        ENDIF
      ENDDO
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
      SUBROUTINE GFDUDX(IPH,IZZ,IDUDX0,AVEVEL)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(12),IDUDI0,IDUDJ0,
     `              IDUDK0,IC2(6)  /GENFL/LGF(9),COLVEL,AXISYM
     1       /FNDRI/IDIDX0,IDJDX0,IDKDX0,IDF(8),IJPW(6),IJPS(6),IDSG(6),
     1              IDR(6),JDR(6),KDR(6) /FNDRL/NXNE1,NYNE1,NZNE1
     1       /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
     1       /RDATA/TINY,RDT1(16),RINNER,RDT2(67)
     1       /IDATA/NX,NY,NZ,IDT1(117) /LDATA/LDT1(19),BFC,LDT2(64)
      LOGICAL LGF,COLVEL,AXISYM,SLD,AVEVEL,LDT1,BFC,LDT2,NXNE1,NYNE1,
     1        NZNE1
C
      IADS= (IZZ-1)*NFM
C.... Main term dU/dX:
      IF(NXNE1) THEN
        IF(COLVEL) THEN
          IF(BFC) THEN
            IADW= (IZZ-1)*NXNYST
            IDX=IDIDX0+IADW; JDX=IDJDX0+IADW; KDX=IDKDX0+IADW
            DO 10 IJ= 1,NXNYST
             F(IDUDX0+IJ)= F(IDUDI0+IJ)*F(IDX+IJ)+F(IDUDJ0+IJ)*F(JDX+IJ)
     1                   + F(IDUDK0+IJ)*F(KDX+IJ)
  10        CONTINUE
          ELSE
            CALL FNDFDX(IUG0(IPH)+IADS,IDUDX0,IZZ,.FALSE.)
          ENDIF
        ELSE
C.... Treatment of staggered velocities on non-BFC grids:
          IF(AVEVEL) THEN
            CALL FNDFDX(IUG0(IPH),IDUDX0,IZZ,.FALSE.)
          ELSE
            CALL AVDVEL(3,4,L0F(3+IPH-1),IDUDX0,IZZ)
          ENDIF
        ENDIF
      ELSE
        CALL ZERNUM(IDUDX0,NXNYST)
      ENDIF
C.... Additional term for polar-cylindrical geometry:
      IF(AXISYM .AND. IVG0(IPH)/=0) THEN
        IVC0= ITWO(IVG0(IPH)+IADS,IVG0(IPH),COLVEL)
        IF(BFC) THEN
C.... For BFC this term appears only if NX=1:
          DO 20 IY= 1,NY
            IF(.NOT.SLD(IY)) THEN
              IDUDX= IDUDX0+IY
              XCP  = COORD1(1,IY,IZZ,1)
              YCP  = COORD1(1,IY,IZZ,2)
              RCP  = SQRT(XCP*XCP+YCP*YCP) + RINNER
              F(IDUDX)= F(IDUDX) + F(IVC0+IY)/(RCP + TINY)
            ENDIF
  20      CONTINUE
        ELSE
          IR0= L0F(LXYR)
          DO 30 IJ= 1,NXNYST
            IF(.NOT.SLD(IJ)) THEN
              IDUDX   = IDUDX0+IJ
              F(IDUDX)= F(IDUDX) + F(IVC0+IJ)/(F(IR0+IJ) + TINY)
            ENDIF
  30      CONTINUE
        ENDIF
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
      SUBROUTINE GFDUDY(IPH,IZZ,IDUDY0)
      INCLUDE 'farray'
      INCLUDE 'cmnrot'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(12),IDUDI0,IDUDJ0,
     1              IDUDK0,IC2(6)  /GENFL/LGF(9),COLVEL,AXISYM
     1       /FNDRI/IDF1(3),IDIDY0,IDJDY0,IDKDY0,IDF2(41)
     1       /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(6),NXNYNZ,IG4(43)
     1       /RDATA/TINY,RDT1(16),RINNER,RDT2(67)
     1       /IDATA/NX,NY,NZ,IDT1(117) /LDATA/LDT1(19),BFC,LDT2(64)
      COMMON /LROT/LROTOR
      LOGICAL LGF,COLVEL,AXISYM,SLD,LDT1,BFC,LDT2,DOIT,LROTOR
C
C.... Main term dU/dY:
      IUC0= ITWO(IUG0(IPH)+(IZZ-1)*NFM,IUG0(IPH),COLVEL)
      IF(BFC) THEN
        IADW= (IZZ-1)*NXNYST
        IDY=IDIDY0+IADW; JDY=IDJDY0+IADW; KDY=IDKDY0+IADW
        DO 10 IJ= 1,NXNYST
          F(IDUDY0+IJ)= F(IDUDI0+IJ)*F(IDY+IJ)+F(IDUDJ0+IJ)*F(JDY+IJ)
     1                + F(IDUDK0+IJ)*F(KDY+IJ)
  10    CONTINUE
      ELSE
        CALL FNDFDY(IUC0,IDUDY0,IZZ,.FALSE.)
      ENDIF
C.... Additional term for polar-cylindrical geometry:
      IF(AXISYM) THEN
        IF(BFC) THEN
C.... For BFC this term appears only if NX=1:
          DO 20 IY= 1,NY
            IF(.NOT.SLD(IY)) THEN
              IDUDY= IDUDY0+IY
              XCP  = COORD1(1,IY,IZZ,1)
              YCP  = COORD1(1,IY,IZZ,2)
              RCP  = SQRT(XCP*XCP+YCP*YCP) + RINNER
              F(IDUDY)= F(IDUDY) - F(IUC0+IY)/(RCP + TINY)
            ENDIF
  20      CONTINUE
        ELSE
          IR0= L0F(LXYR)
          DOIT=LROTOR
          DO 30 IJ= 1,NXNYST
            IF(.NOT.SLD(IJ)) THEN
              IDUDY   = IDUDY0+IJ
              IF(DOIT) THEN  ! in rotor region make U1 in static co-ordinates
                IJK=IJ+(IZZ-1)*NXNYST; IROTC=INROT(IJK)
                IF(IROTC>0) THEN
                  IY= 1+MOD(IJ-1,NY)
                  F(IDUDY)= F(IDUDY)-(F(IUC0+IJ)+ANGV(IROTC)*F(IR0+IJ))/
     1                                                (F(IR0+IJ) + TINY)
                ELSE
                  F(IDUDY)= F(IDUDY)-F(IUC0+IJ)/(F(IR0+IJ) + TINY)
                ENDIF
              ELSE
                F(IDUDY)= F(IDUDY) - F(IUC0+IJ)/(F(IR0+IJ) + TINY)
              ENDIF
            ENDIF
  30      CONTINUE
        ENDIF
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
      SUBROUTINE GFDUDZ(IPH,IZZ,IDUDZ0)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IUIZM0(2),IUIZP0(2),IC1(8),
     1              IDUDI0,IDUDJ0,IDUDK0,IC2(6)
     1       /FNDRI/IDF1(6),IDIDZ0,IDJDZ0,IDKDZ0,IDF2(38)
     1       /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
     1       /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
      LOGICAL LGF,COLVEL,AXISYM,LDT1,BFC,LDT2
C
      IF(COLVEL) THEN
        IF(BFC) THEN
          IADW= (IZZ-1)*NXNYST
          IDZ=IDIDZ0+IADW; JDZ=IDJDZ0+IADW; KDZ=IDKDZ0+IADW
          DO 10 IJ= 1,NXNYST
            F(IDUDZ0+IJ)= F(IDUDI0+IJ)*F(IDZ+IJ)+F(IDUDJ0+IJ)*F(JDZ+IJ)
     1                  + F(IDUDK0+IJ)*F(KDZ+IJ)
  10      CONTINUE
        ELSE
          CALL FNDFDZ(IUG0(IPH)+(IZZ-1)*NFM,IDUDZ0,IZZ,.FALSE.)
        ENDIF
      ELSE
C.... Treatment of staggered velocities on non-BFC grids:
        CALL FUVWDZ(IUIZM0(IPH),IUG0(IPH),IUIZP0(IPH),IDUDZ0,IZZ)
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
      SUBROUTINE GFDVDX(IPH,IZZ,IDVDX0)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(15),IDVDI0,IDVDJ0,
     1              IDVDK0,IC2(3) /GENI/IG1(2),NXNYST,IG2(6),NFM,IG3(50)
     1       /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
     1       /FNDRI/IDIDX0,IDJDX0,IDKDX0,IDF(44)
      LOGICAL LGF,COLVEL,AXISYM,LDT1,BFC,LDT2
C
      IF(BFC) THEN
        IADW= (IZZ-1)*NXNYST
        IDX=IDIDX0+IADW; JDX=IDJDX0+IADW; KDX=IDKDX0+IADW
        DO 10 IJ= 1,NXNYST
          F(IDVDX0+IJ)= F(IDVDI0+IJ)*F(IDX+IJ) + F(IDVDJ0+IJ)*F(JDX+IJ)
     1                + F(IDVDK0+IJ)*F(KDX+IJ)
  10    CONTINUE
      ELSE
        IVC0= ITWO(IVG0(IPH)+(IZZ-1)*NFM,IVG0(IPH),COLVEL)
        CALL FNDFDX(IVC0,IDVDX0,IZZ,.FALSE.)
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
      SUBROUTINE GFDVDY(IPH,IZZ,IDVDY0,AVEVEL)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(15),IDVDI0,IDVDJ0,
     1              IDVDK0,IC2(3) /GENI/IG1(2),NXNYST,IG2(6),NFM,IG3(50)
     1       /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
     1       /FNDRI/IDF1(3),IDIDY0,IDJDY0,IDKDY0,IDF2(41)
      LOGICAL LGF,COLVEL,AXISYM,AVEVEL,LDT1,BFC,LDT2
C
      IF(COLVEL) THEN
        IF(BFC) THEN
          IADW= (IZZ-1)*NXNYST
          IDY=IDIDY0+IADW; JDY=IDJDY0+IADW; KDY=IDKDY0+IADW
          DO 10 IJ= 1,NXNYST
            F(IDVDY0+IJ)= F(IDVDI0+IJ)*F(IDY+IJ)+F(IDVDJ0+IJ)*F(JDY+IJ)
     1                  + F(IDVDK0+IJ)*F(KDY+IJ)
  10      CONTINUE
        ELSE
          CALL FNDFDY(IVG0(IPH)+(IZZ-1)*NFM,IDVDY0,IZZ,.FALSE.)
        ENDIF
      ELSE
C.... Treatment of staggered velocities:
        IF(AVEVEL) THEN
          CALL FNDFDY(IVG0(IPH),IDVDY0,IZZ,.FALSE.)
        ELSE
          CALL AVDVEL(1,2,L0F(5+IPH-1),IDVDY0,IZZ)
        ENDIF
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
      SUBROUTINE GFDVDZ(IPH,IZZ,IDVDZ0)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(4),IVIZM0(2),IVIZP0(2),
     1              IC2(7),IDVDI0,IDVDJ0,IDVDK0,IC3(3)
     1       /FNDRI/IDF1(6),IDIDZ0,IDJDZ0,IDKDZ0,IDF2(38)
     1       /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
     1       /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
      LOGICAL LGF,COLVEL,AXISYM,LDT1,BFC,LDT2
C
      IF(COLVEL) THEN
        IF(BFC) THEN
          IADW= (IZZ-1)*NXNYST
          IDZ=IDIDZ0+IADW; JDZ=IDJDZ0+IADW; KDZ=IDKDZ0+IADW
          DO 10 IJ= 1,NXNYST
            F(IDVDZ0+IJ)= F(IDVDI0+IJ)*F(IDZ+IJ)+F(IDVDJ0+IJ)*F(JDZ+IJ)
     1                  + F(IDVDK0+IJ)*F(KDZ+IJ)
  10      CONTINUE
        ELSE
          CALL FNDFDZ(IVG0(IPH)+(IZZ-1)*NFM,IDVDZ0,IZZ,.FALSE.)
        ENDIF
      ELSE
C.... Treatment of staggered velocities on non-BFC grids:
        CALL FUVWDZ(IVIZM0(IPH),IVG0(IPH),IVIZP0(IPH),IDVDZ0,IZZ)
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
      SUBROUTINE GFDWDX(IPH,IZZ,IDWDX0)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC(18),IDWDI0,IDWDJ0,IDWDK0
     1       /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
     1       /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
     1       /FNDRI/IDIDX0,IDJDX0,IDKDX0,IDF(44)
      LOGICAL LGF,COLVEL,AXISYM,LDT1,BFC,LDT2
C
      IF(BFC) THEN
        IADW= (IZZ-1)*NXNYST
        IDX=IDIDX0+IADW; JDX=IDJDX0+IADW; KDX=IDKDX0+IADW
        DO 10 IJ= 1,NXNYST
          F(IDWDX0+IJ)= F(IDWDI0+IJ)*F(IDX+IJ) + F(IDWDJ0+IJ)*F(JDX+IJ)
     1                + F(IDWDK0+IJ)*F(KDX+IJ)
  10    CONTINUE
      ELSE
        IWC0= ITWO(IWG0(IPH)+(IZZ-1)*NFM,IWG0(IPH),COLVEL)
        CALL FNDFDX(IWC0,IDWDX0,IZZ,.FALSE.)
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
      SUBROUTINE GFDWDY(IPH,IZZ,IDWDY0)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC(18),IDWDI0,IDWDJ0,IDWDK0
     1       /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
     1       /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
     1       /FNDRI/IDF1(3),IDIDY0,IDJDY0,IDKDY0,IDF2(41)
      LOGICAL LGF,COLVEL,AXISYM,LDT1,BFC,LDT2
C
      IF(BFC) THEN
        IADW= (IZZ-1)*NXNYST
        IDY=IDIDY0+IADW; JDY=IDJDY0+IADW; KDY=IDKDY0+IADW
        DO 10 IJ= 1,NXNYST
          F(IDWDY0+IJ)= F(IDWDI0+IJ)*F(IDY+IJ)+F(IDWDJ0+IJ)*F(JDY+IJ)
     1                + F(IDWDK0+IJ)*F(KDY+IJ)
  10    CONTINUE
      ELSE
        IWC0= ITWO(IWG0(IPH)+(IZZ-1)*NFM,IWG0(IPH),COLVEL)
        CALL FNDFDY(IWC0,IDWDY0,IZZ,.FALSE.)
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
      SUBROUTINE GFDWDZ(IPH,IZZ,IDWDZ0,AVEVEL)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(8),IWIZM0(2),IWIZP0(2),
     1              IC2(6),IDWDI0,IDWDJ0,IDWDK0
     1       /FNDRI/IDF1(6),IDIDZ0,IDJDZ0,IDKDZ0,IDF2(38)
     1       /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
     1       /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
      LOGICAL AVEVEL,LGF,COLVEL,AXISYM,LDT1,BFC,LDT2
C
      IF(COLVEL) THEN
        IF(BFC) THEN
          IADW= (IZZ-1)*NXNYST
          IDZ=IDIDZ0+IADW; JDZ=IDJDZ0+IADW; KDZ=IDKDZ0+IADW
          DO 10 IJ= 1,NXNYST
            F(IDWDZ0+IJ)= F(IDWDI0+IJ)*F(IDZ+IJ)+F(IDWDJ0+IJ)*F(JDZ+IJ)
     1                  + F(IDWDK0+IJ)*F(KDZ+IJ)
  10      CONTINUE
        ELSE
          CALL FNDFDZ(IWG0(IPH)+(IZZ-1)*NFM,IDWDZ0,IZZ,.FALSE.)
        ENDIF
      ELSE
C.... Treatment of staggered velocities on non-BFC grids:
        IF(AVEVEL) THEN
          CALL FUVWDZ(IWIZM0(IPH),IWG0(IPH),IWIZP0(IPH),IDWDZ0,IZZ)
        ELSE
          CALL AVDVEL(5,6,L0F(7+IPH-1),IDWDZ0,IZZ)
        ENDIF
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... AVRGVL averages staggered velocities U1,...,W2. Result is stored
C     in temporary storage IUAV0. The later is used in calculations of
C     the generation function for staggered solver on non-BFC geometry.
C     The user might print-out averaged velocities by storing U1C,V1C,
C     etc in Q1 (see GXGENF for details).
C.... IDIP is the direction index; which can be either IDIP=4 for U1;
C     or IDIP=2 for V1; or IDIP=6 for W1.
C.... ISLAB is slab index; ISLAB=0 -averaging is for the current slab;
C     ISLAB=1 -averaging is for the HIGH-slab.
C.... AVRGVL is called from GXGENF.
C
c
      SUBROUTINE AVRGVL(IDIP,IUAV0,IPH,IZZ,ISLAB)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      COMMON /IDATA/NX,NY,NZ,IDT1(117)  /F0/LB1(109),KZXCY,LB2(194)
     1       /GENI/ NXNY,NXM1NY,NXNYST,IG1(6),NFM,IG2(50)
     1       /FNDRI/IDF(11),IJPW(6),IJPS(6),IDSG(6),IDR(6),JDR(6),KDR(6)
     1       /FACES/L0FACE,L0FACZ  /FACDIR/IBND(6),IFCP(6),IWLP(6)
      LOGICAL SLD,LSLDR,BLKM,BLKP,XCYCZ,NEZ
C.... Preliminaries
      IDIR=IDIP-1; L0FZST=L0FACZ
      L0FACZ= L0FACZ + ISLAB*NXNYST
      IF(IDIR==1) THEN
        IUST0= L0F(V1+IPH-1) + ISLAB*NFM
      ELSEIF(IDIR==3) THEN
        IUST0= L0F(U1+IPH-1) + ISLAB*NFM
      ELSE
        IUST0= L0F(W1+IPH-1) + ISLAB*NFM
      ENDIF
      XCYCZ= IDIR==3 .AND. KZXCY/=0
      IF(XCYCZ) XCYCZ= NEZ(F(KZXCY+IZZ))
C.... Loop over the current slab
      DO 20 IX= 1,NX
       IADX= (IX-1)*NY
       IF(XCYCZ .AND. IX==1) IJPS(4)= NXM1NY
       DO 10 IY= 1,NY
        IJ= IY+IADX
        IF( SLD(IJ) ) THEN
          F(IUAV0+IJ)= 0.0
        ELSE
          IUSTC= IUST0+IJ
          BLKM = LSLDR(IJ,IDIP)
          BLKP = LSLDR(IJ,IDIR)
          IF(BLKM .AND. BLKP) THEN
            F(IUAV0+IJ)= 0.0
          ELSEIF(BLKM) THEN
            F(IUAV0+IJ)= F(IUSTC)
          ELSEIF(BLKP) THEN
            F(IUAV0+IJ)= F(IUSTC+IJPS(IDIP))
          ELSE
            F(IUAV0+IJ)= 0.5*( F(IUSTC) + F(IUSTC+IJPS(IDIP)) )
          ENDIF
        ENDIF
  10   CONTINUE
       IF(XCYCZ .AND. IX==1) IJPS(4)= -NY
  20  CONTINUE
      L0FACZ= L0FZST
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... GXKEGB is called from group 1, group 19 section 4, and group 13 of
C     GREX3 by setting the coefficient & value to GRND4 in the COVAL
C     statements for KE and EP for the patch name KEBUOY.
C     See the 'K-Epsilon turbulence model' and 'GRAVitational' entries
C     in the PHOENICS Encyclopaedia for further details.
C
C.... BUOYA, BUOYB and BUOYC signify the resolutes of the gravitational
C     acceleration in the cartesian coordinate direction XC, YC and ZC.
C
C.... The coefficient GCEB determines the magnitude of the buoyancy
C     production term in the EP equation. There is no "standard" value
C     for GCEB, although it is recommended that GCEB should be close to
C     zero for stably-stratified flows, and close to unity for unstably-
C     stratified flows. The default in PHOENICS-VR is to compute GCEB
C     from the function: GCEB = tanh (v/U) where v is the velocity
C     component parallel to the gravity vector, and U is the velocity
C     component perpendicular to the gravity vector.
C     This default setting is GCEB=1.0, but this
C     may be overwritten by setting: SPEDAT(SET,KEBUOY,GCEB,R,0.2), say,
C     in the Q1 file.
C
c
      SUBROUTINE GXKEGB
      INCLUDE 'farray'
      INCLUDE 'satear'
      INCLUDE 'grdear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdbfc'
      COMMON /INTDMN/IDMN,NUMDMN,NXMAX,NYMAX,NZMAX,NXYMAX,NXYZMX,NCLTOT
      COMMON /GENI/ IG1(2),NXNYST,IG2(30),ILTLS,IG3(15),ITEM1,ITEM2,
     1              IG4(4),IPRPS,IG5(4)
     1       /F0/    LB1(109),KZXCY,LB2(194)
     1       /FNDRL/NXNE1,NYNE1,NZNE1 /NAMFN/NAMFUN,NAMSUB
     1       /GENL/  LGL1(14),debgtz,LGL2(45)
     1       /FNDRI/IDF(11),IJPW(6),IJPS(6),IDSG(6),IDR(6),JDR(6),KDR(6)
     1       /LRNTM2/LBWDIS,LBFMU,LBFONE,LBFTWO,LBREYN,LBREYT,LBEPKE
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
      COMMON/BFICRT/IUCRT(2),IVCRT(2),IWCRT(2),ITMP(6)
      DIMENSION GV(3),VF(3),VPG(3),VNG(3)
      LOGICAL NXNE1,NYNE1,NZNE1,GRN,BOUSS,SLDEF,LGL1,debgtz,LGL2,
     1        dbgloc,solu,solv,solw,VARC3E,XCYCZ,BLKM,BLKP,LSLDR
      CHARACTER*6 NAMFUN,NAMSUB
      SAVE BOUSS,JSCAL,JDRDX,JDRDY,JDRDZ,JGENB,L0DRDX,L0DRDY,L0DRDZ,
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
     1     L0GENB,L0SCAL,GCEB,RECPRT,SLDEF,L0C3EB0,SOLU,SOLV,SOLW,
     1     JC3EB,VARC3E,C3ERLX
C
      NAMSUB='GXKEGB'
      DBGLOC=DEBGTZ
      IF(INDVAR>0) DBGLOC=DBGLOC.AND.DBGPHI(INDVAR)
      IF(FLAG.OR.DBGLOC)  CALL BANNER(1,NAMSUB,250516)
      IF(CCM) THEN
        CALL UNKEGB
        RETURN
      ENDIF
      IF(IGR==1 .AND. ISC==2) THEN
C-----------------------------------------------------------------------
C     CALL FROM GR.1; SEC.2 TO PREPARE CALCULATIONS
C-----------------------------------------------------------------------
C.... THE BOUSSINESQ APPROXIMATION IS EMPLOYED WHEN RHO1 IS CONSTANT.
        BOUSS= .NOT.GRN(RHO1)
        SLDEF= .FALSE.
        IF(BOUSS) THEN
          IF(SOLVE(ITEM1)) THEN
            JSCAL= ITEM1
            SLDEF= .TRUE.
          ELSEIF(SOLVE(H1)) THEN
            JSCAL= H1
          ELSE
            CALL WRYT40('GXKEGB REQUIRES TEM1 OR H1 TO BE SOLVED ')
            CALL WRYT40('FOR THE CURRENT GRAVITY &/OR RHO1 MODELS')
            CALL WRIT40('GXKEGB REQUIRES TEM1 OR H1 TO BE SOLVED ')
            CALL WRIT40('FOR THE CURRENT GRAVITY &/OR RHO1 MODELS')
            CALL SET_ERR(222,'ERROR. SEE RESULT FILE.',1)
C            CALL EAROUT(1)
          ENDIF
          RECPRT= 1./PRT(JSCAL)
        ELSE
          IF(DEN1<=0) THEN
            CALL WRIT40('GXKEGB REQUIRES STORE(DEN1) IN Q1!  ')
            CALL WRYT40('GXKEGB REQUIRES STORE(DEN1) IN Q1!  ')
            CALL SET_ERR(223,'ERROR. GXKEGB REQUIRES STORE(DEN1) IN Q1.'
     *                  ,1)
C            CALL EAROUT(1)
          ENDIF
          JSCAL = DEN1
          RECPRT= 1./PRT(H1)
        ENDIF
C.... ALLOCATE STORAGE FOR SCALAR DERIVATIVES (LD1,LD2,LD3 ARE IN USE):
        JDRDX=LBNAME('DRDX'); JDRDY=LBNAME('DRDY')
        JDRDZ=LBNAME('DRDZ'); JGENB=LBNAME('GENB')
        L0DRDX=L0F(LD1); L0DRDY=L0F(LD2); L0DRDZ=L0F(LD3)
C.... IF A 3D STORE HAS NOT BE PROVIDED BY STORE(GENB) IN Q1, CREATE
C     A SLAB-WISE STORE AND INITIALISE CONTENTS TO ZERO.
        IF(JGENB==0) CALL GXMAKE0(L0GENB,NXYMAX,'GENB')
C
        GCEB=1.0
        CALL GETSDR('KEBUOY','GCEB',GCEB)
        CALL WRIT1R('K-E GCEB',GCEB)
C!!! MRM 19.10.05 INTRODUCE HENKES C3E FUNCTION
        VARC3E=.FALSE.
        IF(GCEB<0.0) VARC3E=.TRUE.
        IF(VARC3E) THEN
          JC3EB=LBNAME('C3EB')
          IF(JC3EB==0) CALL GXMAKE0(L0C3EB0,NXYZMX,'C3EB')
          C3ERLX=0.3
          CALL GETSDR('KEBUOY','C3ERLX',C3ERLX)
          CALL WRIT1R('KE-C3RLX',C3ERLX)
          SOLU=SOLVE(U1); SOLV=SOLVE(V1); SOLW=SOLVE(W1)
        ENDIF
      ELSEIF(IGR==19 .AND. ISC==4) THEN
C-----------------------------------------------------------------------
C     Call from Gr.19; Sec.4 to calculate Gb term.
C-----------------------------------------------------------------------
C.... If DRDX,...,GENB are stored in Q1, get slab-wise addresses:
        IF(JDRDX/=0) L0DRDX= L0F(JDRDX)
        IF(JDRDY/=0) L0DRDY= L0F(JDRDY)
        IF(JDRDZ/=0) L0DRDZ= L0F(JDRDZ)
        IF(JGENB/=0) L0GENB= L0F(JGENB)
C.... Calculate derivatives of JSCAL
        L0SCAL= L0F(JSCAL)
        CALL ZERNM3(L0DRDX,L0DRDY,L0DRDZ,NXNYST)
        IF(NXNE1) CALL FNDFDX(L0SCAL,L0DRDX,IZSTEP,SLDEF)
        IF(NYNE1) CALL FNDFDY(L0SCAL,L0DRDY,IZSTEP,SLDEF)
        IF(NZNE1.AND..NOT.PARAB) CALL FNDFDZ(L0SCAL,L0DRDZ,IZSTEP,SLDEF)
        IF(CARTES .OR. BFC) THEN
C... BFC or Cartesian geometries:
          DO 10 IJ= 1,NXNYST
            GDRDX= F(L0DRDX+IJ)
            GDRDY= F(L0DRDY+IJ)
            GDRDZ= F(L0DRDZ+IJ)
            F(L0GENB+IJ)= BUOYA*GDRDX + BUOYB*GDRDY + BUOYC*GDRDZ
  10      CONTINUE
        ELSE
C.... Polar-cylindrical geometry:
          L0XG= L0F(LXXG)
          DO IX= 1,NX
            IADX = (IX-1)*NY
            ANGLE= F(L0XG+IX)
            SINA = SIN(ANGLE); COSA = COS(ANGLE)
            DO IY= 1,NY
              IJ= IY+IADX
              GDRDX= F(L0DRDX+IJ)
              GDRDY= F(L0DRDY+IJ)
              GDRDZ= F(L0DRDZ+IJ)
              F(L0GENB+IJ)= BUOYA*( COSA*GDRDX+SINA*GDRDY)
     1                  + BUOYB*(-SINA*GDRDX+COSA*GDRDY) + BUOYC*GDRDZ
            ENDDO
          ENDDO
        ENDIF
C.... The volumetric production of K due to buoyancy forces
        L0ENUT= L0F(VIST)
        DO IJ= 1,NXNYST
          F(L0GENB+IJ)= -F(L0ENUT+IJ)*F(L0GENB+IJ)*RECPRT
        ENDDO
        IF(BOUSS) THEN
C.... Put expansion coefficient into L0GENB:
C.... FN46(Y,X,A,B) ........... Y = Y*(A+B*X):
          CALL FN46(-L0GENB,LDVO1,0.0,-1.0)
C.... Divide by Cp when enthalpy is the solved for variable
C.... FN27(Y,A) ............... Y = Y / X
          IF(JSCAL==H1) CALL FN27(-L0GENB,LCP1)
        ELSE
          CALL FN27(-L0GENB,DEN1)
        ENDIF
C
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
C.... Compute C3E from velocity function
        IF(VARC3E) THEN
C.... Unit gravity vector
          CALL VECTOR(GV,BUOYA,BUOYB,BUOYC)
          CALL NORM(GV)
          IF(BFC) THEN
            IF(SOLU) L0U1=L0F(IUCRT(1))
            IF(SOLV) L0V1=L0F(IVCRT(1))
            IF(SOLW) L0W1=L0F(IWCRT(1))
          ELSE
            IF(SOLU) L0U1=L0F(U1)
            IF(SOLV) L0V1=L0F(V1)
            IF(SOLW) L0W1=L0F(W1)
          ENDIF
          IF(JC3EB/=0) THEN
            L0C3EB = L0F(JC3EB)
          ELSE
            L0C3EB=L0C3EB0+(IZ-1)*NXNYST
          ENDIF
          L0XG= L0F(LXXG)
          XCYCZ= SOLU .AND. KZXCY/=0
          IF(XCYCZ) XCYCZ= NEZ(F(KZXCY+IZ))
          DO IX= 1,NX
            IADX = (IX-1)*NY
            IF(XCYCZ .AND. IX==1) IJPS(4)= NXNYST-NX
            DO IY= 1,NY
              IJ= IY+IADX
C.... Cartesian velocity vector
              CALL VECTOR(VF,0.0,0.0,0.0)
              IF(SOLW) THEN
                IDIR=5; IDIP=6
                BLKM = LSLDR(IJ,IDIP); BLKP = LSLDR(IJ,IDIR)
                IF(BLKM .AND. BLKP) THEN
                  VF(3)=0.0
                ELSEIF(BLKM) THEN
                  VF(3)=F(L0W1+IJ)
                ELSEIF(BLKP) THEN
                  VF(3)=F(L0W1+IJ+IJPS(IDIP))
                ELSE
                  VF(3)=0.5*(F(L0W1+IJ)+F(L0W1+IJ+IJPS(IDIP)))
                ENDIF
              ENDIF
              IF(SOLU) THEN
                IDIR=3; IDIP=4
                BLKM = LSLDR(IJ,IDIP); BLKP = LSLDR(IJ,IDIR)
                IF(BLKM .AND. BLKP) THEN
                  VF(1)=0.0
                ELSEIF(BLKM) THEN
                  VF(1)=F(L0U1+IJ)
                ELSEIF(BLKP) THEN
                  VF(1)=F(L0U1+IJ+IJPS(IDIP))
                ELSE
                  VF(1)=0.5*(F(L0U1+IJ)+F(L0U1+IJ+IJPS(IDIP)))
                ENDIF
              ENDIF
              IF(SOLV) THEN
                IDIR=1; IDIP=2
                BLKM = LSLDR(IJ,IDIP); BLKP = LSLDR(IJ,IDIR)
                IF(BLKM .AND. BLKP) THEN
                  VF(2)=0.0
                ELSEIF(BLKM) THEN
                  VF(2)=F(L0V1+IJ)
                ELSEIF(BLKP) THEN
                  VF(2)=F(L0V1+IJ+IJPS(IDIP))
                ELSE
                  VF(2)=0.5*(F(L0V1+IJ)+F(L0V1+IJ+IJPS(IDIP)))
                ENDIF
              ENDIF
              IF(.NOT.CARTES) THEN
                UVEL=VF(1); VVEL=VF(2)
                ANGLE= F(L0XG+IX)
                SINA = SIN(ANGLE); COSA = COS(ANGLE)
                IF(SOLU) VF(1)=  COSA*UVEL + SINA*VVEL
                IF(SOLV) VF(2)= -SINA*UVEL + COSA*VVEL
              ENDIF
C.... velocity vectors parallel and normal to g
              VFCPG = DOT(VF,GV)
              CALL MULVEC(GV,VFCPG,VPG)
              CALL SUBVEC(VF,VPG,VNG)
              VFCNG = VECMAG(VNG)
              VRAT = ABS(VFCPG/(VFCNG + TINY))
              C3ENEW = TANH(VRAT)
              F(L0C3EB+IJ)=C3ENEW*C3ERLX+(1.-C3ERLX)*F(L0C3EB+IJ)
            ENDDO
            IF(XCYCZ .AND. IX==1) IJPS(4)= -NY
          ENDDO
        ENDIF
C
      ELSEIF(IGR==13 .AND. ISC==5) THEN
C-----------------------------------------------------------------------
C     Call from Gr.13; Sec.5 to process CO=GRND4.
C-----------------------------------------------------------------------
        L0CO=L0F(CO); L0KE=L0F(KE)
        IF(INDVAR==KE) THEN
          DO 40 IJ= 1,NXNYST
  40      F(L0CO+IJ)= AMAX1( FIXFLU, -F(L0GENB+IJ)/(F(L0KE+IJ)+TINY) )
        ELSEIF(INDVAR==EP) THEN
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
          IF(JC3EB/=0) THEN
            L0C3EB = L0F(JC3EB)
          ELSE
            L0C3EB=L0C3EB0+(IZ-1)*NXNYST
          ENDIF
          DO 50 IJ= 1,NXNYST
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
            GC3EB=GCEB
            IF(VARC3E) GC3EB = F(L0C3EB+IJ)
  50      F(L0CO+IJ)=AMAX1(FIXFLU,-GC3EB*F(L0GENB+IJ)/(F(L0KE+IJ)+TINY))
        ENDIF
      ELSEIF(IGR==13 .AND. ISC==16) THEN
C-----------------------------------------------------------------------
C     Call from Gr.13; Sec.16 to process VAL=GRND4.
C-----------------------------------------------------------------------
        L0VAL= L0F(VAL)
        RFF  = 1./FIXFLU
        IF(INDVAR==KE) THEN
          DO 60 IJ= 1,NXNYST
  60      F(L0VAL+IJ)= AMAX1( 0.0, RFF*F(L0GENB+IJ) )
        ELSEIF(INDVAR==EP) THEN
          L0KE=L0F(KE); L0EP=L0F(EP)
          IF(LBEPKE/=0) L0EPKE=L0F(LBEPKE)
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
          IF(JC3EB/=0) THEN
            L0C3EB = L0F(JC3EB)
          ELSE
            L0C3EB=L0C3EB0+(IZ-1)*NXNYST
          ENDIF
          DO 70 IJ= 1,NXNYST
            GEPDK= F(L0EP+IJ)/(F(L0KE+IJ)+TINY)
            IF(LBEPKE/=0) F(L0EPKE+IJ) = GEPDK
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
            GC3EB=GCEB
            IF(VARC3E) GC3EB = F(L0C3EB+IJ)
            F(L0VAL+IJ)= AMAX1( RFF*GC3EB*GEPDK*F(L0GENB+IJ), 0.0 )
  70      CONTINUE
        ENDIF
      ENDIF
      namsub='gxkegb'
      if(flag.or.dbgloc)  call banner(2,namsub,250516)
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... AXIBFC checks the BFC-grid and returns AXISYM= T if the grid is
C     axisymmetric.
c
      SUBROUTINE AXIBFC(AXISYM)
      COMMON /IDATA/ NX,NY,NZ,IDT1(117) /RDATA/TINY,RDT1(84)
     1       /XYZCRN/X1,X2,X3,X4,X5,X6,X7,X8,Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8,
     1               Z1,Z2,Z3,Z4,Z5,Z6,Z7,Z8
      LOGICAL AXISYM
C
      IF(NY<=2) THEN
        AXISYM= .FALSE.
      ELSE
        AXISYM= .TRUE.
        DO 10 IZZ= 1,NZ
          CALL CELCOR(1,2,IZZ)
          A1= ABS(2.*(X4-X1)/(Y1+Y4+TINY))
          A2= ABS(2.*(X2-X3)/(Y2+Y3+TINY))
          A3= ABS(2.*(X5-X8)/(Y5+Y8+TINY))
          A4= ABS(2.*(X6-X7)/(Y6+Y7+TINY))
          AA1= 0.25*(A1+A2+A3+A4)
          CALL CELCOR(1,NY,IZZ)
          A1= ABS(2.*(X4-X1)/(Y1+Y4+TINY))
          A2= ABS(2.*(X2-X3)/(Y2+Y3+TINY))
          A3= ABS(2.*(X5-X8)/(Y5+Y8+TINY))
          A4= ABS(2.*(X6-X7)/(Y6+Y7+TINY))
          AA2= 0.25*(A1+A2+A3+A4)
          AXISYM= AXISYM .AND. ABS(AA1-AA2)<0.01*AA2
  10    CONTINUE
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c