c

C#### JCL 21.03.18 use ISTEP==FSTEP not ISTEP==1
C#### JCL 08.01.18 correction for undefined variable
C#### JCL 08.11.17 include satear to clear undefined variables
C.... file name GXSURFTS.F              220318
C**** SUBROUTINE GXSURFT  is called from group 13 of GREX3, and is
C     entered only when the patch name begins with the characters 'STEN'.
C
C.... The coding provided here presumes that CO in the relevant COVAL
C     statement is FIXFLUX, and that the VAL is:
C     GRND3, or (what is equivalent) DENSITY,     (ISC = 15)
c     the type is east, north or high
C
C     The arguments BFC, CARTES and PARAB are logicals in the SATELLITE,
C     signifying, when set to T, the body-fitted coordinate, rectangular
C     cartesian coordinate and parabolic flow respectively.
C
C
      SUBROUTINE GXSURFT
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'satear'
      INCLUDE 'grdear'
      INCLUDE 'grdbfc'
      INCLUDE 'prpcmn'
      COMMON
C####     1      /IDATA/NX,NY,NZ,IFL2(117) /LDATA/LDFL1(7),XCYCLE,LDFL2(76)
C####     1      /RDATA/TINY,RDAT1(84) 
     1      /GENL/ LGN1(51),SOLIDM,SOLIDN,LGN3(7)
     1      /GENI/ NXNY,IGNI1(8),NFM,IGFIL(39),ITEM1,ITEM2,IGFIL2(4),
     1             IPRPS,IFIL1(4)
     1      /FNDRI/IDF(11),IJPW(6),IJPS(6),IDSG(6),IDR(6),JDR(6),KDR(6)
     1      /NAMFN/NAMFUN,NAMSUB
C####      INTEGER DEN1,DEN2
C####      LOGICAL BFC,CARTES,PARAB,FLUID1,QNE,XCYCLE,LDFL1,LDFL2,SLD,
      LOGICAL FLUID1,QNE,SLD,
     1        TESTSOL,LGN1,SOLIDM,SOLIDN,LGN3,GRN,dbSURFT
      COMMON /GEODMN0/ I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
     1                I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
      COMMON/SLDPRP/SOLPRP,PORPRP,VACPRP,SOLCRT,SOLSAV
      COMMON /FACES/L0FACE,L0FACZ
      CHARACTER DIRCTN, NAMFUN*6, NAMSUB*6
      PARAMETER (NLG=100, NIG=200, NRG=200, NCG=100)
C
      COMMON/LGRND/LG(NLG)/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
      LOGICAL LG
      CHARACTER*4 CG

      SAVE L0YNY,L0YNX,L0YNZ,L0XNY,L0XNX,L0XNZ,L0ZNY,L0ZNX,L0ZNZ,
     1     L0WK1,L0WK2,L0WK3,L0WORK,L0NOMY0,L0NOMX0,L0NOMZ0,L0SURT0,
     1     L0KAPA0,ISURT,IKAPA,INOMX,INOMY,INOMZ,ISURN,ITM
      SAVE UUDF
C
      dbSURFT=.false.
!!!      dbSURFT=.true.
      if(dbSURFT) then
        write(14,*) 'gxSURFT entered for '
        call writ2i('indvar  ',indvar,'isc     ',isc)
        call writ1r('SURFTa   ',SURFTa)
      endif
      NAMSUB = 'GXSURFT'
      IF(IGR==1.AND.ISC==3) THEN
C***********************************************************************
C   * -----------GROUP 1  SECTION  3 ---------------------------
        CALL GXMAKE(L0WORK,NXNY*NZ,'WORK')
        IF(NY>1) CALL GXMAKE(L0YNY,NXNY,'YNY')
        IF(NX>1) CALL GXMAKE(L0XNX,NXNY,'XNX')
        IF(NZ>1) CALL GXMAKE(L0ZNZ,NXNY,'ZNZ')
        IF(NY>1.AND.NX>1)THEN
          CALL GXMAKE(L0XNY,NXNY,'XNY')
          CALL GXMAKE(L0YNX,NXNY,'YNX')
        ENDIF
        IF(NY>1.AND.NZ>1)THEN
          CALL GXMAKE(L0YNZ,NXNY,'YNZ')
          CALL GXMAKE(L0ZNY,NXNY,'ZNY')
        ENDIF
        IF(NZ>1.AND.NX>1)THEN
          CALL GXMAKE(L0ZNX,NXNY,'ZNX')
          CALL GXMAKE(L0XNZ,NXNY,'XNZ')
        ENDIF
        CALL GXMAKE(L0WK1,NXNY,'WK1')
        CALL GXMAKE(L0WK2,NXNY,'WK2')
        CALL GXMAKE(L0WK3,NXNY,'WK3')
        ISURT=LBNAME('SURT')
        IF(ISURT<=0) CALL GXMAKE(L0SURT0,NXNY*NZ,'SURT')
        IKAPA=LBNAME('KAPA')
        IF(IKAPA<=0) CALL GXMAKE(L0KAPA0,NXNY*NZ,'KAPA')
        IF(NX>1) THEN
          INOMX=LBNAME('NOMX')
          IF(INOMX<=0) CALL GXMAKE(L0NOMX0,NXNY*NZ,'NOMX')
        ENDIF
        IF(NY>1) THEN
          INOMY=LBNAME('NOMY')
          IF(INOMY<=0) CALL GXMAKE(L0NOMY0,NXNY*NZ,'NOMY')
        ENDIF
        IF(NZ>1) THEN
          INOMZ=LBNAME('NOMZ')
          IF(INOMZ<=0) CALL GXMAKE(L0NOMZ0,NXNY*NZ,'NOMZ') 
        ENDIF
        ISURN=LBNAME('SURN')
        UUDF=0.3; CALL GETSDR('SURFTENS','UUDF',UUDF)
!! ITM - 1: number of smoothing      
        ITM=4; CALL GETSDI('SURFTENS','NSMOOTH',ITM)
      ELSEIF(IGR==13) THEN
C***********************************************************************
C--- GROUP 13. Boundary conditions and special sources
        IF(INDVAR.LT.5) THEN
c.... which velocity?
c                                             U1 or U2
          DIRCTN = 'E'
          IDIREC=1
          IPLPHI=NY
          IPLSLD=NY
          TESTSOL=SOLIDM
        ELSEIF(INDVAR.LT.7) THEN
c                                             V1 or V2
          DIRCTN = 'N'
          IDIREC=2
          IPLPHI=1
          IPLSLD=1
          TESTSOL=SOLIDM
        ELSE
c                                             W1 or W2
          DIRCTN = 'H'
          IDIREC=3
          IPLPHI=NFM
          IPLSLD=NXNY
          TESTSOL=SOLIDM.OR.SOLIDN
        ENDIF
        if(dbSURFT) then
          call writ40('direction is '//dirctn)
        endif
c
        IF(ISC.LT.16) THEN ! VAL < GRND4
          
          if(dbSURFT) then
            call writ40('in gxSURFTso for value')
            call writ2l('bfc     ',bfc,'cartes  ',cartes)
          endif

C---- VAL= surtension  (ie GRND3)
          IF(ISC.EQ.15) THEN
C.... 'EAST,NORTH or HIGH' force = sigma*(2*rho*/(RHO1+RHO2))*KAPPA*(CFnext-CF)
C..
C.... Preliminaries
C....                                
            if(dbSURFT) call writ40('surface tension option')
            IF(FLUID1(INDVAR)) THEN
              IF(ISURT>0) THEN
                L0FSURT=L0F(ISURT)
                IF(NZ>1)L0FSURTH=L0F(HIGH(ISURT))
                IF(IZ==NZ.AND.NZ>1)L0FSURTH=L0FSURT
              ELSE
                L0FSURT=L0SURT0+(IZ-1)*NXNY
                L0FSURTH=L0SURT0+(MIN(IZ+1,NZ)-1)*NXNY
              ENDIF 
              IF(IKAPA>0) THEN
                L0KAPA=L0F(IKAPA)
                IF(NZ>1)L0FKAPAH=L0F(HIGH(IKAPA))
C#### JCL 08.01.18 L0FKAPA is undefined 
C####                IF(IZ==NZ.AND.NZ>1)L0FKAPAH=L0FKAPA             
                IF(IZ==NZ.AND.NZ>1)L0FKAPAH=L0KAPA             
              ELSE
                L0KAPA=L0KAPA0+(IZ-1)*NXNY
                L0KAPAH=L0KAPA0+(MIN(IZ+1,NZ)-1)*NXNY                
              ENDIF
              L0FKAPA=L0F(EASP2)   
            ELSE
!! no provision for phase 2 (IPSA)           
            ENDIF
C.... Put staggered KAPA in EASP2
            IF(.NOT.PARAB) THEN
              NFM0=NFM
C#### JCL 22.03.18 use ISTEP==FSTEP not ISTEP==1
C####              IF(ISTEP==1.AND.ISWEEP==1)THEN
              IF(ISTEP==FSTEP.AND.ISWEEP==1)THEN
                CALL FN0(EASP2,-L0KAPA)                    
              ELSE
                IF(IDIREC==3.AND.IKAPA<=0)NFM=NXNY                    
                CALL FNAV(EASP2,-L0KAPA,DIRCTN)
                NFM=NFM0              
              ENDIF
              IF(IDIREC==3.AND.ISURT<=0)NFM=NXNY                    
              CALL FNAV(-L0WK1,-L0FSURT,DIRCTN)
              NFM=NFM0
            ENDIF
            if(dbSURFT) then
              call writ40('staggered kappa')
              call prn('EASP2',EASP2)
            endif
c     therefore no further staggering is needed below
c
            
            L0FVAL=L0F(VAL)
            IF(IDIREC==1)THEN
              L0FDSC=L0F(DXG2D)
            ELSEIF(IDIREC==2)THEN
              L0FDSC=L0F(DYG2D)
            ENDIF   
            
C.... Multiply surface tension 2.0*surfta/(RHOL+RHOG))* kappa*(cfnext-cf)
cc    SURFTA=SIGMA
            RHOL= F(INDPRTB(IPRPSA,0)+1)
            RHOG= F(INDPRTB(IPRPSB,0)+1)         
            IF(IPRPS.NE.0) THEN
              SURFTD=2.0*SURFTA/(RHOL+RHOG)
              DO IX=IXF,IXL
                DO IY=IYF,IYL
                  I=IY+(IX-1)*NY
                  IF(IDIREC==1)THEN
                    I1=I+NY
                    IF(IX==NX)I1=I
                    DERIV=(F(L0FSURT+I1)-F(L0FSURT+I))/F(L0FDSC+I) 
                  ELSEIF(IDIREC==2)THEN
                    I1=I+1
                    IF(IY==NY)I1=I                
                    DERIV=(F(L0FSURT+I1)-F(L0FSURT+I))/F(L0FDSC+I) 
                  ELSEIF(IDIREC==3)THEN
                    I1=I+NXNY                                
                    IF(ISURT>0)I1=I+NFM
                    DSCZ=DZG
                    IF(IZ==NZ)THEN
                      I1=I     
                      DSCZ=DZGL
                    ENDIF
                    DERIV=(F(L0FSURT+I1)-F(L0FSURT+I))/DSCZ
                  ENDIF   
                  CFILTAVE=F(L0WK1+I)
                  UUPHI= CFILTAVE-0.5
                  DELTAFUNC=HEAV2(UUPHI,UUDF) 
                  F(L0FVAL+I)=SURFTD*F(L0FKAPA+I)*DERIV*DELTAFUNC
                ENDDO
              ENDDO
 
              if(dbSURFT) call prn('vall',val)
            ENDIF
          ENDIF ! ISC 15
        ENDIF !ISC <16
C.... Set VAL to zero in solids
        IF(TESTSOL) THEN
          L0FVAL= L0F(VAL)
          DO IX= IXF,IXL
            IADD= (IX-1)*NY
            DO IY= IYF,IYL
              I= IY+IADD
              IF(SLD(I)) THEN
                F(L0FVAL+I)= 0.0
              ELSE
                IF(SLD(I+IPLSLD)) F(L0FVAL+I)= 0.0
              ENDIF
            ENDDO
          ENDDO 
        ENDIF
      ELSEIF(IGR==19) THEN ! Group 19
C***********************************************************************
        IF(ISC==2) THEN
C***********************************************************************
c   * ----------group 19  section 2 ---- start of sweep.
          
          DO IZZ = 1,NZ
            L0SURN=L0F(ANYZ(ISURN,IZZ))
            IF(ISURT>0) THEN
              L0SURT=L0F(ANYZ(ISURT,IZZ))
            ELSE
              L0SURT=L0SURT0+(IZZ-1)*NXNY
            ENDIF
            DO I=1,NXNY
              F(L0SURT+I)=F(L0SURN+I)
            ENDDO
          ENDDO
          
!! smoothing surn
          IT=1
!! ITM - 1: number of smoothing      
!          ITM=4

          DO WHILE (IT0) THEN
                L0SURT=L0F(ANYZ(ISURT,IZZ))
                L0SURTZP1=L0F(ANYZ(ISURT,MIN(IZZ+1,NZ)))
                L0SURTZM1=L0F(ANYZ(ISURT,MAX(IZZ-1,1)))            
              ELSE
                L0SURT=L0SURT0+(IZZ-1)*NXNY
                L0SURTZP1=L0SURT0+(MIN(IZZ+1,NZ)-1)*NXNY
                L0SURTZM1=L0SURT0+(MAX(IZZ-1,1 )-1)*NXNY
              ENDIF
              L0IPRPS=L0F(ANYZ(IPRPS,IZZ))
              IF(NX>1) IAWB=I2DAWB+(IZZ-1)*NY  ! west domain boundary area
              IF(NY>1) IASB=I2DASB+(IZZ-1)*NX  ! south domain boundary area
              IF(NZ>1) IALB=I2DALB             ! low domain boundary area
            
              DO IX=1,NX
                DO IY=1,NY
                  I3D=IY+(IX-1)*NY
                  I3DH=I3D+(IZZ-1)*NXNY
                  IF(NX>1)THEN
                    IF(IX==NX)THEN
                      FFE=F(L0SURT+I3D)
                      AFE=F(I3DAEX+I3DH)
                      IF(F(L0IPRPS+I3D)>SOLPRP)AFE=0.0                        
                    ELSE                   
                      FFE=F(L0SURT+I3D+NY)
                      AFE=F(I3DAEX+I3DH+NY)   
                      IF(F(L0IPRPS+I3D+NY)>SOLPRP)AFE=0.0      
                    ENDIF
                    IF(IX==1)THEN
                      FFW=F(L0SURT+I3D)
                      AFW=F(IAWB+I3D) 
                      IF(F(L0IPRPS+I3D)>SOLPRP)AFW=0.0
                    ELSE
                      FFW=F(L0SURT+I3D-NY)
                      AFW=F(I3DAEX+I3DH-NY) 
                      IF(F(L0IPRPS+I3D-NY)>SOLPRP)AFW=0.0                              
                    ENDIF
                  ELSE
                    FFW=0.0; FFE=0.0; AFW=0.0; AFE=0.0
                  ENDIF
                  IF(NY>1)THEN
                    IF(IY==1)THEN
                      FFS=F(L0SURT+I3D) 
                      AFS=F(IASB+I3D) 
                      IF(F(L0IPRPS+I3D)>SOLPRP)AFS=0.0
                    ELSE 
                      FFS=F(L0SURT+I3D-1) 
                      AFS=F(I3DANY+I3DH-1)   
                      IF(F(L0IPRPS+I3D-1)>SOLPRP)AFS=0.0                        
                    ENDIF
                    IF(IY==NY)THEN
                      FFN=F(L0SURT+I3D)
                      AFN=F(I3DANY+I3DH)   
                      IF(F(L0IPRPS+I3D)>SOLPRP)AFN=0.0                        
                    ELSE 
                      FFN=F(L0SURT+I3D+1)
                      AFN=F(I3DANY+I3DH+1)  
                      IF(F(L0IPRPS+I3D+1)>SOLPRP)AFN=0.0
                    ENDIF   
                  ELSE
                    FFS=0.0; FFN=0.0; AFS=0.0; AFN=0.0
                  ENDIF
                  IF(NZ>1)THEN
                    IF(IZZ==1)THEN
                      FFB=F(L0SURT+I3D)
                      AFB=F(IALB+I3D)
                      IF(F(L0IPRPS+I3D)>SOLPRP)AFB=0.0
                    ELSE
                      FFB=F(L0SURTZM1+I3D)
                      AFB=F(I3DAHZ+I3DH-NXNY)       
                      IF(F(L0IPRPS+I3D-NFM)>SOLPRP)AFT=0.0                        
                    ENDIF
                    IF(IZZ==NZ)THEN
                      FFT=F(L0SURT+I3D)
                      AFT=F(I3DAHZ+I3DH)   
                      IF(F(L0IPRPS+I3D)>SOLPRP)AFB=0.0                        
                    ELSE
                      FFT=F(L0SURTZP1+I3D)
                      AFT=F(I3DAHZ+I3DH+NXNY) 
                      IF(F(L0IPRPS+I3D+NFM)>SOLPRP)AFT=0.0
                    ENDIF
                  ELSE
                    FFB=0.0; FFT=0.0; AFB=0.0; AFT=0.0
                  ENDIF
                        
                FTEM=FFE*AFE+FFW*AFW+FFS*AFS+FFN*AFN+FFB*AFB+FFT*AFT
                FTEM=FTEM/(AFE+AFW+AFS+AFN+AFB+AFT+tiny)
           
                  F(L0WORK+I3DH)=FTEM
                ENDDO
              ENDDO
            ENDDO
     
            DO IZZ = 1,NZ
              IF(ISURT>0) THEN
                L0SURT=L0F(ANYZ(ISURT,IZZ))
              ELSE
                L0SURT=L0SURT0+(IZZ-1)*NXNY
              ENDIF
              DO I=1,NXNY
                I3DH=I+(IZZ-1)*NXNY       
                F(L0SURT+I)=F(L0WORK+I3DH)
              ENDDO
            ENDDO   
          ENDDO
!! End Of smoothing
!! compute kappa
!!
!!     For surface tension we will need a table SURT  (Filtered SURN)
!!    And three tables (in 3D), 2 in 2D and 1 in 1 D  will contain at beginning the derivatives of the filtered surn
!!    And on output the forces in the three direction to be addd as source force in the U1, V1 and W1 equations
!!   If we suppose that we have the filtered value of surn then the surface tension coding is as follow:
    
! let's compute Normx,Normy,Normz (gradient of the smoothed colorfunction)
C#### JCL 22.03.18 use ISTEP==FSTEP not ISTEP==1
C####          IF(ISWEEP==1.AND.ISTEP==1)THEN
          IF(ISWEEP==1.AND.ISTEP==FSTEP)THEN
          ELSE
            L0FACZ0=L0FACZ
            DO IZZ=1,NZ
              L0FACZ=L0FACE+(IZZ-1)*NXNY
              IF(ISURT>0) THEN
                L0FSURT=L0F(ANYZ(ISURT,IZZ))
                L0FSURTP=L0F(ANYZ(ISURT,MIN(IZZ+1,NZ)))
                L0FSURTM=L0F(ANYZ(ISURT,MAX(IZZ-1,1)))           
              ELSE
                L0FSURT=L0SURT0+(IZZ-1)*NXNY
                L0FSURTP=L0SURT0+(MIN(IZZ+1,NZ)-1)*NXNY
                L0FSURTM=L0SURT0+(MAX(IZZ-1,1)-1)*NXNY       
              ENDIF
              IF(NY>1)THEN
                IF(INOMY>0) THEN
                  L0FNOMY=L0F(ANYZ(INOMY,IZZ))               
                ELSE
                  L0FNOMY=L0NOMY0+(IZZ-1)*NXNY
                ENDIF 
                CALL FNDFDY(L0FSURT,L0FNOMY,IZZ,.FALSE.)              
              ENDIF
              IF(NX>1)THEN
                IF(INOMX>0) THEN
                  L0FNOMX=L0F(ANYZ(INOMX,IZZ))               
                ELSE
                  L0FNOMX=L0NOMX0+(IZZ-1)*NXNY
                ENDIF
                CALL FNDFDX(L0FSURT,L0FNOMX,IZZ,.FALSE.)
              ENDIF
              IF(NZ>1)THEN
                IF(INOMZ>0) THEN
                  L0FNOMZ=L0F(ANYZ(INOMZ,IZZ))    
                  L0FNOMZ1=L0F(ANYZ(INOMZ,1))           
                  L0FNOMZNZ1=L0F(ANYZ(INOMZ,NZ-1))
                ELSE
                  L0FNOMZ=L0NOMZ0+(IZZ-1)*NXNY    
                  L0FNOMZ1=L0NOMZ0
                  L0FNOMZNZ1=L0NOMZ0+(NZ-1-1)*NXNY
                ENDIF
                IF(IZZ>1.AND.IZZ0) THEN
            L0FSURT=L0F(ISURT)              
          ELSE
            L0FSURT=L0SURT0+(IZ-1)*NXNY
          ENDIF

          IF(IKAPA>0) THEN
            L0FKAPA=L0F(IKAPA)              
          ELSE
            L0FKAPA=L0KAPA0+(IZ-1)*NXNY
          ENDIF
          IF(NY>1)THEN
            IF(INOMY>0) THEN
              L0FNOMY=L0F(INOMY) 
              L0FNOMYP=L0F(ANYZ(INOMY,MIN(IZ+1,NZ)))
              L0FNOMYM=L0F(ANYZ(INOMY,MAX(IZ-1,1)))                     
            ELSE
              L0FNOMY=L0NOMY0+(IZ-1)*NXNY
              L0FNOMYP=L0NOMY0+(MIN(IZ+1,NZ)-1)*NXNY
              L0FNOMYM=L0NOMY0+(MAX(IZ-1,1)-1)*NXNY                     
            ENDIF 
          ENDIF
          IF(NX>1)THEN
            IF(INOMX>0) THEN
              L0FNOMX=L0F(INOMX)     
              L0FNOMXP=L0F(ANYZ(INOMX,MIN(IZ+1,NZ)))
              L0FNOMXM=L0F(ANYZ(INOMX,MAX(IZ-1,1)))                         
            ELSE
              L0FNOMX=L0NOMX0+(IZ-1)*NXNY
              L0FNOMXP=L0NOMX0+(MIN(IZ+1,NZ)-1)*NXNY
              L0FNOMXM=L0NOMX0+(MAX(IZ-1,1)-1)*NXNY                     
            ENDIF 
          ENDIF
          IF(NZ>1)THEN
            IF(INOMZ>0) THEN
              L0FNOMZ=L0F(INOMZ)   
              L0FNOMZP=L0F(ANYZ(INOMZ,MIN(IZ+1,NZ)))
              L0FNOMZM=L0F(ANYZ(INOMZ,MAX(IZ-1,1)))                         
            ELSE
              L0FNOMZ=L0NOMZ0+(IZ-1)*NXNY
              L0FNOMZP=L0NOMZ0+(MIN(IZ+1,NZ)-1)*NXNY
              L0FNOMZM=L0NOMZ0+(MAX(IZ-1,1)-1)*NXNY                                   
            ENDIF 
          ENDIF
C#### JCL 22.03.18 use ISTEP==FSTEP not ISTEP==1
C####          IF(ISWEEP==1.AND.ISTEP==1)THEN
          IF(ISWEEP==1.AND.ISTEP==FSTEP)THEN
            DO I=1,NXNY
              F(L0FKAPA+I)=0.0
            ENDDO
          ELSE
c  
! compute the derivatives in X,Y & Z of the Normal in Y

            IF(NY>1)THEN
              CALL FNDFDY(L0FNOMY,L0YNY,IZ,.FALSE.)   
              IF(NX>1)CALL FNDFDX(L0FNOMY,L0YNX,IZ,.FALSE.)
              IF(NZ>1)THEN
                IF(INOMY<=0)THEN
                  IJPS(5)=NXNY; IJPS(6)=-NXNY
                ENDIF               
                CALL FNDFDZ(L0FNOMY,L0YNZ,IZ,.FALSE.) 
                IF(INOMY<=0)THEN
                  IJPS(5)=NFM; IJPS(6)=-NFM
                ENDIF                        
              ENDIF
            ENDIF
            
! compute the derivatives in X,Y & Z of the Normal in X
            IF(NX>1)THEN
              IF(NY>1)CALL FNDFDY(L0FNOMX,L0XNY,IZ,.FALSE.)   
              CALL FNDFDX(L0FNOMX,L0XNX,IZ,.FALSE.)
              IF(NZ>1)THEN                  
                IF(INOMX<=0)THEN
                  IJPS(5)=NXNY; IJPS(6)=-NXNY
                ENDIF               
                CALL FNDFDZ(L0FNOMX,L0XNZ,IZ,.FALSE.)
                IF(INOMX<=0)THEN
                  IJPS(5)=NFM; IJPS(6)=-NFM
                ENDIF
              ENDIF
            ENDIF
             
! compute the derivatives in X,Y & Z of the Normal in Y 
            IF(NZ>1)THEN
              IF(NY>1)CALL FNDFDY(L0FNOMZ,L0ZNY,IZ,.FALSE.)                 
              IF(NX>1)CALL FNDFDX(L0FNOMZ,L0ZNX,IZ,.FALSE.)
              IF(INOMZ<=0)THEN
                IJPS(5)=NXNY; IJPS(6)=-NXNY
              ENDIF               
              CALL FNDFDZ(L0FNOMZ,L0ZNZ,IZ,.FALSE.)  
              IF(INOMZ<=0)THEN
                IJPS(5)=NFM; IJPS(6)=-NFM
              ENDIF                            
            ENDIF
c    compute DIV (N) = GRADX(Nx)+GradY(Ny)+GradZ(NZ) : L0WK1
             !Y = X1 + X2 + X3 
            CALL FN1(-L0WK1,0.0)
            IF(NY>1)CALL FN34(-L0WK1,-L0YNY,1.0)   !!Y = Y + A * X
            IF(NX>1)CALL FN34(-L0WK1,-L0XNX,1.0)
            IF(NZ>1)CALL FN34(-L0WK1,-L0ZNZ,1.0)
             
c  compute VECMAG = (NY**2 + NX**2 +NZ**2) : L0WK2
            CALL FN1(-L0WK2,0.0)
             
            IF(NY>1)CALL FN49(-L0WK2,-L0FNOMY) !Y = Y + X**2 
            IF(NX>1)CALL FN49(-L0WK2,-L0FNOMX) !Y = Y + X**2
            IF(NZ>1)CALL FN49(-L0WK2,-L0FNOMZ) !Y = Y + X**2  
cc            CALL FN33(-L0WK2,AEPS) !Y = Y + A
            
cc  term11 = NORMY**2*DYNORMY    :  L0YNY    
cc  term22 = NORMX**2*DXNORMX    :  L0XNX
cc  term33 = NORMZ**2*DZNORMZ    :  L0ZNZ
            CALL FN1(-L0WK3,0.0)
            IF(NY>1)THEN
              CALL FN37(-L0YNY,-L0FNOMY,2.0) !Y = Y * X**A
              CALL FN34(-L0WK3,-L0YNY,1.0)     !Y = Y + X
            ENDIF
             
            IF(NX>1)THEN
              CALL FN37(-L0XNX,-L0FNOMX,2.0) !Y = Y * X**A
              CALL FN34(-L0WK3,-L0XNX,1.0)     !Y = Y + X
            ENDIF
            IF(NZ>1)THEN
              CALL FN37(-L0ZNZ,-L0FNOMZ,2.0) !Y = Y * X**A
              CALL FN34(-L0WK3,-L0ZNZ,1.0)     !Y = Y + X
            ENDIF

cc  term12=Normx*NormY*(DNORMXY+DNORMYX) 
            IF(NX>1.AND.NY>1)THEN              
cc term12 = DNORMXY+DNORMYX
              CALL FN34(-L0XNY,-L0YNX,1.0) !Y = Y + A * X               
cc term12= L0XNY*L0GRSX*L0GRSY
              CALL FN77(-L0XNY,-L0FNOMX,-L0FNOMY) !Y = Y * X1 * X2
              CALL FN34(-L0WK3,-L0XNY,1.0)     !Y = Y + X
            ENDIF  
cc  term13=NormY*NormZ*(DNORMYZ+DNORMZY) 
            IF(NY>1.AND.NZ>1)THEN              
cc term13 = DNORMZY+DNORMYZ
              CALL FN34(-L0ZNY,-L0YNZ,1.0) !Y = Y + A * X               
cc term13= L0ZNY*L0GRSZ*L0GRSY
              CALL FN77(-L0ZNY,-L0FNOMZ,-L0FNOMY) !Y = Y * X1 * X2
              CALL FN34(-L0WK3,-L0ZNY,1.0)     !Y = Y + X               
            ENDIF 
               
cc  term23=NormX*NormZ*(DNORMXZ+DNORMZX) 
            IF(NX>1.AND.NZ>1)THEN              
cc term23 = DNORMZX+DNORMXZ
              CALL FN34(-L0ZNX,-L0XNZ,1.0) !Y = Y + A * X               
cc term23= L0ZNX**L0GRSZ*L0GRSX
              CALL FN77(-L0ZNX,-L0FNOMZ,-L0FNOMX) !Y = Y * X1 * X2
              CALL FN34(-L0WK3,-L0ZNX,1.0)     !Y = Y + X               
            ENDIF               
               
cc L0WK3= TERM= term11 + term 22 + term33 + TERM12 + TERM13 + TERM23
cc L0WK3 =  L0WK3/ L0WK2 = TERM= TERM/VECMAG
            CALL FN27(-L0WK3,-L0WK2) !Y = Y / X
              
cc sqrt (vecmag) L0WK2= sqrt(L0WK2)
            CALL FN30(-L0WK2) !Y = SQRT(Y)
               
cc  KAPPA = (TERM-DIVN)/SQRT(VECMAG)
            AEPS=1.0E-14

            DO I=1,NXNY
              F(L0FKAPA+I)=0.0        
              IF(F(L0WK2+I)>AEPS.AND.F(L0FSURT+I)>0.04
     1                            .AND.F(L0FSURT+I)<0.96)THEN
                F(L0FKAPA+I)=(F(L0WK3+I)-F(L0WK1+I))/F(L0WK2+I)
              ENDIF
            ENDDO
c            CALL FN67(-L0FKAPA,-L0WK3,-L0WK1,-L0WK2,1.0,-1.0) !Y = (B1*X1 + B2*X2)/X3
            
          ENDIF
        ENDIF
      ENDIF
      NAMSUB= 'gxSURFT'
      if(dbSURFT) write(14,*) 'leaving gxSURFT '
      END      
C***************************************************************
      FUNCTION HEAV2(PHI,E)
      PARAMETER (PI=3.14159265)      
      IF(PHI.GT.E) HEAV2=0.
      IF(PHI.LT.-E) HEAV2=0.
      IF(ABS(PHI).LE.E) HEAV2=0.5*(1.0+cos(PI*PHI/E))
      END