c
C.... File Name ...... GXGENK.FOR........... 230312
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-centerd 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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      COMMON /CMNGN/ IUG0(2),IVG0(2),IWG0(2),IUIZM0,IUIZP0,IVIZM0,
     1               IVIZP0,IWIZM0,IWIZP0,IDUDI0,IDUDJ0,IDUDK0,IDVDI0,
     1               IDVDJ0,IDVDK0,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
     1       /RSTM2/ J0DUDX,J0DUDY,J0DUDZ,J0DVDX,J0DVDY,J0DVDZ,J0DWDX,
     1               J0DWDY,J0DWDZ,JRSTF(45)
     1       /F0/    LB1(109),KZXCY,LB2(194)
      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
      SAVE ITMP10,ITMP20,ITMP30,IU1C,IV1C,IW1C
      SAVE LSTOU,LSTOV,LSTOW,AVEVEL,DVDL3D,LDBAVR
C
      NAMSUB = 'GXGENF'
      if(indvar.gt.0) then
        dbgloc=debgtz.and.dbgphi(indvar)
      else
        dbgloc=.false.
      endif
      if(flag.or.dbgloc)  then
        call banner(1,namsub,230312)
        call writ2i('igr     ',igr,'isc     ',isc)
      endif
      IF(IGR.EQ.1 .AND. ISC.EQ.2) THEN ! Call from Gr.1; Sec.2 to prepare
        CALL SUB2L( AVEVEL,LSG53, AXISYM,.NOT.CARTES )
        IF(BFC .AND. NX.EQ.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
          CALL SUB3(IUG0(1),IUCRT(1),IVG0(1),IVCRT(1),IWG0(1),IWCRT(1))
          CALL SUB3(IUG0(2),IUCRT(2),IVG0(2),IVCRT(2),IWG0(2),IWCRT(2))
C.... Temporary slab-storage for DUDI,DUDJ,...,DWDK:
          CALL SUB3(IDUDI0,L0F(LD4),IDVDI0,L0F(LD7),IDWDI0,L0F(LD10)  )
          CALL SUB3(IDUDJ0,L0F(LD5),IDVDJ0,L0F(LD8),IDWDJ0,L0F(LD11)  )
          CALL SUB3(IDUDK0,L0F(LD6),IDVDK0,L0F(LD9),IDWDK0,L0F(LDSTAG))
        ELSEIF(CCM) THEN
          CALL SUB2(IUG0(1),LBNAME('UC1'), IUG0(2),LBNAME('UC2'))
          CALL SUB2(IVG0(1),LBNAME('VC1'), IVG0(2),LBNAME('VC2'))
          CALL SUB2(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.
          CALL SUB3L( LSTOU,STORE(3), LSTOV,STORE(5), LSTOW,STORE(7) )
          CALL SUB3( IUG0(1),IUIZM0,  IVG0(1),IVIZM0,  IWG0(1),IWIZM0  )
          CALL SUB3( IUIZM0,L0F(LD4), IVIZM0,L0F(LD5), IWIZM0,L0F(LD6) )
C.... Debug 3D-storage of averaged velocities:
          IU1C= LBNAME('U1C')
          IV1C= LBNAME('V1C')
          IW1C= LBNAME('W1C')
          LDBAVR= IU1C.NE.0 .OR. IV1C.NE.0 .OR. IW1C.NE.0
          CALL SUB3(IUG0(2),0, IVG0(2),0, IWG0(2),0)
          IF(.NOT.ONEPHS) THEN
            CALL SUB3(IUG0(2),IUG0(1), IVG0(2),IVG0(1), IWG0(2),IWG0(1))
          ENDIF
        ENDIF
        IF(BFC.AND.KZXCY.NE.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).NE.0) IUG0(1)= L0F(IUG0(1))
          IF(IVG0(1).NE.0) IVG0(1)= L0F(IVG0(1))
          IF(IWG0(1).NE.0) IWG0(1)= L0F(IWG0(1))
          IF(IUG0(2).NE.0) IUG0(2)= L0F(IUG0(2))
          IF(IVG0(2).NE.0) IVG0(2)= L0F(IVG0(2))
          IF(IWG0(2).NE.0) IWG0(2)= L0F(IWG0(2))
          LSTOU= IUG0(1).NE.0
          LSTOV= IVG0(1).NE.0
          LSTOW= IWG0(1).NE.0
        ENDIF
C.... Define auxiliary logicals:
        CALL SUB3L( LDUDX,.FALSE., LDUDY,.FALSE., LDUDZ,.FALSE. )
        CALL SUB3L( LDVDX,.FALSE., LDVDY,.FALSE., LDVDZ,.FALSE. )
        CALL SUB3L( 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):
        CALL SUB3( ITMP10,L0F(LD1), ITMP20,L0F(LD2), ITMP30,L0F(LD3) )
C.... Store velocity derivatives for the dump into RESULT file:
        DVDL3D= IDUDX.NE.0 .OR. IDUDY.NE.0 .OR. IDUDZ.NE.0 .OR.
     1          IDVDX.NE.0 .OR. IDVDY.NE.0 .OR. IDVDZ.NE.0 .OR.
     1          IDWDX.NE.0 .OR. IDWDY.NE.0 .OR. IDWDZ.NE.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.EQ.1.AND.STRGN1 .OR. IPH.EQ.2.AND.STRGN2) .AND.
     1      INDVAR.GT.10 .AND. INDVAR.NE.ILTLS ) RETURN
C.... If IGEN0.ne.0 calculate generation term, otherwise only DUDX,...
        LDOGEN= IGEN0.NE.0
        IF(GCV) THEN
          CALL BFGENB(IGEN0)
          RETURN
        ENDIF
C.... Calculate Generation for MB-FGE technique
        IF(CCM .AND. NUMBLK.GT.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.GT.1 ) THEN
                  CALL SHINX2(IUIZM0,IUG0(IPH),IUG0(IPH),IUIZP0)
                ELSE
                  CALL AVRGVL(4,IUG0(IPH),IPH,IZSTEP,0)
                ENDIF
                IF(IZSTEP.LT.NZ) CALL AVRGVL(4,IUIZP0,IPH,IZSTEP+1,1)
              ELSE
                CALL AVRGVL(4,IUG0(IPH),IPH,IZSTEP,0)
              ENDIF
            ELSE
              IUG0(IPH)= L0F(U1+IPH-1)
              CALL SUB2( IUIZM0,IUG0(IPH)-NFM, IUIZP0,IUG0(IPH)+NFM )
            ENDIF
          ENDIF
          IF(LSTOV) THEN
            IF(NYNE1 .AND. (LDVDX.OR.LDVDY.OR.LDVDZ)) THEN
              IF(LDVDZ) THEN
                IF(IZSTEP.GT.1) THEN
                  CALL SHINX2(IVIZM0,IVG0(IPH),IVG0(IPH),IVIZP0)
                ELSE
                  CALL AVRGVL(2,IVG0(IPH),IPH,IZSTEP,0)
                ENDIF
                IF(IZSTEP.LT.NZ) CALL AVRGVL(2,IVIZP0,IPH,IZSTEP+1,1)
              ELSE
                CALL AVRGVL(2,IVG0(IPH),IPH,IZSTEP,0)
              ENDIF
            ELSE
              IVG0(IPH)= L0F(V1+IPH-1)
              CALL SUB2( IVIZM0,IVG0(IPH)-NFM, IVIZP0,IVG0(IPH)+NFM )
            ENDIF
          ENDIF
          IF(LSTOW) THEN
            IF(NZNE1 .AND. (LDWDX.OR.LDWDY.OR.LDWDZ)) THEN
              IF(LDWDZ) THEN
                IF(IZSTEP.GT.1 ) THEN
                  CALL SHINX2(IWIZM0,IWG0(IPH),IWG0(IPH),IWIZP0)
                ELSE
                  CALL AVRGVL(6,IWG0(IPH),IPH,IZSTEP,0)
                ENDIF
                IF(IZSTEP.LT.NZ) CALL AVRGVL(6,IWIZP0,IPH,IZSTEP+1,1)
              ELSE
                CALL AVRGVL(6,IWG0(IPH),IPH,IZSTEP,0)
              ENDIF
            ELSE
              IWG0(IPH)= L0F(W1+IPH-1)
            ENDIF
          ENDIF
C.... Put averaged velocities in 3D-storage for DUMP into RESULT file:
          IF(LDBAVR) THEN
            IF(IU1C.NE.0) CALL SHINUM(L0F(IU1C),IUG0(IPH),NXNYST)
            IF(IV1C.NE.0) CALL SHINUM(L0F(IV1C),IVG0(IPH),NXNYST)
            IF(IW1C.NE.0) CALL SHINUM(L0F(IW1C),IWG0(IPH),NXNYST)
          ENDIF
        ELSEIF(BFC) THEN
          IF(FDFSOL .AND. LSTOW .AND. NZ.EQ.1) IWG0(IPH)= L0F(W1+IPH-1)
C.... Calculate and store DUDI,DUDJ,...,DWDK for a current slab:
          CALL UVWIJK(IPH,IZSTEP)
        ENDIF
C.... Calculate generation term for a current slab.
        IF(LDOGEN) CALL ZERNUM(IGEN0,NXNYST)
C.... Term 1 = 2*((DU/DX+J*V/R)**2 + DV/DY**2 +DW/DZ**2)
        IF(LDUDX .OR. LDVDY .OR. LDWDZ) THEN
          CALL SUB3R(AM1,0.0, AM2,0.0, AM3,0.0)
          IF(RSTM) CALL SUB3(ITMP10,J0DUDX,ITMP20,J0DVDY,ITMP30,J0DWDZ)
          IF(LDUDX) THEN
            AM1= 1.0
            CALL GFDUDX(IPH,IZSTEP,ITMP10,AVEVEL)
          ENDIF
          IF(LDVDY) THEN
            AM2= 1.0
            CALL GFDVDY(IPH,IZSTEP,ITMP20,AVEVEL)
          ENDIF
          IF(LDWDZ) THEN
            AM3= 1.0
            CALL GFDWDZ(IPH,IZSTEP,ITMP30,AVEVEL)
          ENDIF
          IF(LDOGEN) CALL GFTRM1(IGEN0,AM1,ITMP10,AM2,ITMP20,AM3,ITMP30)
          IF(DVDL3D) THEN
            IF(IDUDX.NE.0) CALL SHINUM(L0F(IDUDX),ITMP10,NXNYST)
            IF(IDVDY.NE.0) CALL SHINUM(L0F(IDVDY),ITMP20,NXNYST)
            IF(IDWDZ.NE.0) CALL SHINUM(L0F(IDWDZ),ITMP30,NXNYST)
            IF(IDWDZ.NE.0) CALL SHINUM(L0F(IDWDZ),ITMP30,NXNYST)
          ENDIF
        ENDIF
C.... Term 2 = ( DU/DY - J*U/R + DV/DX )**2
        IF(LDUDY .OR. LDVDX) THEN
          CALL SUB2R(AM1,0.0, AM2,0.0)
          IF(RSTM) CALL SUB2(ITMP10,J0DUDY, ITMP20,J0DVDX)
          IF(LDUDY) THEN
            AM1= 1.0
            CALL GFDUDY(IPH,IZSTEP,ITMP10)
          ENDIF
          IF(LDVDX) THEN
            AM2= 1.0
            CALL GFDVDX(IPH,IZSTEP,ITMP20)
          ENDIF
          IF(LDOGEN) CALL GFTRM2(IGEN0,AM1,ITMP10,AM2,ITMP20)
          IF(DVDL3D) THEN
            IF(IDUDY.NE.0) CALL SHINUM(L0F(IDUDY),ITMP10,NXNYST)
            IF(IDVDX.NE.0) CALL SHINUM(L0F(IDVDX),ITMP20,NXNYST)
          ENDIF
        ENDIF
C.... Term 3 = ( DU/DZ + DW/DX )**2
        IF(LDUDZ .OR. LDWDX) THEN
          CALL SUB2R(AM1,0.0, AM2,0.0)
          IF(RSTM) CALL SUB2(ITMP10,J0DUDZ, ITMP20,J0DWDX)
          IF(LDUDZ) THEN
            AM1= 1.0
            CALL GFDUDZ(IPH,IZSTEP,ITMP10)
          ENDIF
          IF(LDWDX) THEN
            AM2= 1.0
            CALL GFDWDX(IPH,IZSTEP,ITMP20)
          ENDIF
          IF(LDOGEN) CALL GFTRM2(IGEN0,AM1,ITMP10,AM2,ITMP20)
          IF(DVDL3D) THEN
            IF(IDUDZ.NE.0) CALL SHINUM(L0F(IDUDZ),ITMP10,NXNYST)
            IF(IDWDX.NE.0) CALL SHINUM(L0F(IDWDX),ITMP20,NXNYST)
          ENDIF
        ENDIF
C.... Term 4 = ( DW/DY + DV/DZ )**2
        IF(LDWDY .OR. LDVDZ) THEN
          CALL SUB2R(AM1,0.0, AM2,0.0)
          IF(RSTM) CALL SUB2(ITMP10,J0DWDY, ITMP20,J0DVDZ)
          IF(LDWDY) THEN
            AM1= 1.0
            CALL GFDWDY(IPH,IZSTEP,ITMP10)
          ENDIF
          IF(LDVDZ) THEN
            AM2= 1.0
            CALL GFDVDZ(IPH,IZSTEP,ITMP20)
          ENDIF
          IF(LDOGEN) CALL GFTRM2(IGEN0,AM1,ITMP10,AM2,ITMP20)
          IF(DVDL3D) THEN
            IF(IDWDY.NE.0) CALL SHINUM(L0F(IDWDY),ITMP10,NXNYST)
            IF(IDVDZ.NE.0) CALL SHINUM(L0F(IDVDZ),ITMP20,NXNYST)
          ENDIF
        ENDIF
      ENDIF
      NAMSUB = 'gxgenf'
      if(flag.or.dbgloc)  call banner(2,namsub,230312)
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... GFTRM1 and GFTRM2 calculate the Generation terms:
c
      SUBROUTINE GFTRM1(IGEN0,A1,ID10,A2,ID20,A3,ID30)
      INCLUDE '/phoenics/d_includ/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 '/phoenics/d_includ/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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(6),IDUDI0,IDUDJ0,IDUDK0,
     1              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
            CALL SUB3(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).NE.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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/cmnrot'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(6),IDUDI0,IDUDJ0,IDUDK0,
     1              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
        CALL SUB3(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=.FALSE.
          IF(LROTOR) THEN
            DO IROT=1,NROTOR
              IF(IZZ.GE.IZ1(IROT).AND.IZZ.LE.IZ2(IROT)) THEN
                DOIT=.TRUE. ! in rotor region
                EXIT
              ENDIF
            ENDDO
          ENDIF
          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
                IY= 1+MOD(IJ-1,NY)
                IF(IY.GE.IY1(IROT).AND.IY.LE.IY2(IROT)) THEN
                  F(IDUDY)= F(IDUDY) -(F(IUC0+IJ)+ANGV(IROT)*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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IUIZM0,IUIZP0,IC1(4),IDUDI0,
     1              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
          CALL SUB3(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,IUG0(IPH),IUIZP0,IDUDZ0,IZZ)
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
      SUBROUTINE GFDVDX(IPH,IZZ,IDVDX0)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(9),IDVDI0,IDVDJ0,IDVDK0,
     1              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
        CALL SUB3(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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(9),IDVDI0,IDVDJ0,IDVDK0,
     1              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
          CALL SUB3(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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(2),IVIZM0,IVIZP0,IC2(5),
     1              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
          CALL SUB3(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,IVG0(IPH),IVIZP0,IDVDZ0,IZZ)
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
      SUBROUTINE GFDWDX(IPH,IZZ,IDWDX0)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC(12),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
        CALL SUB3(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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC(12),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
        CALL SUB3(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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(4),IWIZM0,IWIZP0,IC2(6),
     1              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
          CALL SUB3(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,IWG0(IPH),IWIZP0,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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/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
      CALL SUB2( IDIR,IDIP-1, L0FZST,L0FACZ )
      L0FACZ= L0FACZ + ISLAB*NXNYST
      IF(IDIR.EQ.1) THEN
        IUST0= L0F(V1+IPH-1) + ISLAB*NFM
      ELSEIF(IDIR.EQ.3) THEN
        IUST0= L0F(U1+IPH-1) + ISLAB*NFM
      ELSE
        IUST0= L0F(W1+IPH-1) + ISLAB*NFM
      ENDIF
      XCYCZ= IDIR.EQ.3 .AND. KZXCY.NE.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.EQ.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.EQ.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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/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       /FNDRL/NXNE1,NYNE1,NZNE1 /NAMFN/NAMFUN,NAMSUB
     1       /GENL/  LGL1(14),debgtz,LGL2(45)
     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
      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.GT.0) DBGLOC=DBGLOC.AND.DBGPHI(INDVAR)
      IF(FLAG.OR.DBGLOC)  CALL BANNER(1,NAMSUB,151105)
      IF(CCM) THEN
        CALL UNKEGB
        RETURN
      ENDIF
      IF(IGR.EQ.1 .AND. ISC.EQ.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.LE.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):
        CALL SUB2( JDRDX,LBNAME('DRDX'), JDRDY,LBNAME('DRDY') )
        CALL SUB2( JDRDZ,LBNAME('DRDZ'), JGENB,LBNAME('GENB') )
        CALL SUB3( 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.EQ.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.LT.0.0) VARC3E=.TRUE.
        IF(VARC3E) THEN
          JC3EB=LBNAME('C3EB')
          IF(JC3EB.EQ.0) CALL GXMAKE0(L0C3EB0,NXYZMX,'C3EB')
          C3ERLX=0.3
          CALL GETSDR('KEBUOY','C3ERLX',C3ERLX)
          CALL WRIT1R('KE-C3RLX',C3ERLX)
          CALL SUB3L(SOLU,SOLVE(U1),SOLV,SOLVE(V1),SOLW,SOLVE(W1))
        ENDIF
      ELSEIF(IGR.EQ.19 .AND. ISC.EQ.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.NE.0) L0DRDX= L0F(JDRDX)
        IF(JDRDY.NE.0) L0DRDY= L0F(JDRDY)
        IF(JDRDZ.NE.0) L0DRDZ= L0F(JDRDZ)
        IF(JGENB.NE.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 20 IX= 1,NX
           IADX = (IX-1)*NY
           ANGLE= F(L0XG+IX)
           SINA = SIN(ANGLE)
           COSA = COS(ANGLE)
           DO 20 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
  20      CONTINUE
        ENDIF
C.... The volumetric production of K due to buoyancy forces
        L0ENUT= L0F(VIST)
        DO 30 IJ= 1,NXNYST
  30    F(L0GENB+IJ)= -F(L0ENUT+IJ)*F(L0GENB+IJ)*RECPRT
        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.EQ.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.NE.0) THEN
            L0C3EB = L0F(JC3EB)
          ELSE
            L0C3EB=L0C3EB0+(IZ-1)*NXNYST
          ENDIF
          L0XG= L0F(LXXG)
          DO IX= 1,NX
            IADX = (IX-1)*NY
            IF(.NOT.CARTES) THEN
              ANGLE= F(L0XG+IX)
              SINA = SIN(ANGLE)
              COSA = COS(ANGLE)
            ENDIF
            DO IY= 1,NY
              IJ= IY+IADX
C.... Cartesian velocity vector
              CALL VECTOR(VF,0.0,0.0,0.0)
              IF(SOLW) VF(3)=F(L0W1+IJ)
              IF(CARTES .OR. BFC) THEN
                IF(SOLU) VF(1)=F(L0U1+IJ)
                IF(SOLV) VF(2)=F(L0V1+IJ)
              ELSE
                IF(SOLU) VF(1)=  COSA*F(L0U1+IJ) + SINA*F(L0V1+IJ)
                IF(SOLV) VF(2)= -SINA*F(L0U1+IJ) + COSA*F(L0V1+IJ)
              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
          ENDDO
        ENDIF
C
      ELSEIF(IGR.EQ.13 .AND. ISC.EQ.5) THEN
C-----------------------------------------------------------------------
C     Call from Gr.13; Sec.5 to process CO=GRND4.
C-----------------------------------------------------------------------
        CALL SUB2( L0CO,L0F(CO), L0KE,L0F(KE) )
        IF(INDVAR.EQ.KE) THEN
          DO 40 IJ= 1,NXNYST
  40      F(L0CO+IJ)= AMAX1( FIXFLU, -F(L0GENB+IJ)/(F(L0KE+IJ)+TINY) )
        ELSEIF(INDVAR.EQ.EP) THEN
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
          IF(JC3EB.NE.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.EQ.13 .AND. ISC.EQ.16) THEN
C-----------------------------------------------------------------------
C     Call from Gr.13; Sec.16 to process VAL=GRND4.
C-----------------------------------------------------------------------
        L0VAL= L0F(VAL)
        RFF  = 1./FIXFLU
        IF(INDVAR.EQ.KE) THEN
          DO 60 IJ= 1,NXNYST
  60      F(L0VAL+IJ)= AMAX1( 0.0, RFF*F(L0GENB+IJ) )
        ELSEIF(INDVAR.EQ.EP) THEN
          CALL SUB2( L0KE,L0F(KE), L0EP,L0F(EP) )
          IF(LBEPKE.NE.0) L0EPKE=L0F(LBEPKE)
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
          IF(JC3EB.NE.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.NE.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,191005)
      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.LE.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).LT.0.01*AA2
  10    CONTINUE
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c