added surface area calculation to nek_expansion. Refs #8

examples/integration_example/integration_example.usr

@@ -20,14 +20,16 @@ c-----------------------------------------------------------------------
       return
       end
 c-----------------------------------------------------------------------
+! Sets the force term in the momentum equation. Use this to set body
+! forces such as gravity. Note that ffx, ffy, and ffz will later be
+! multiplied by the density.
       subroutine userf  (ix,iy,iz,eg)
       include 'SIZE'
       include 'TOTAL'
       include 'NEKUSE'
-      integer e,f,eg
       ffx=0.0
       ffy=0.0
-      ffz=0.0     
+      ffz=0.0
       return
       end
 c-----------------------------------------------------------------------
@@ -305,42 +307,45 @@ C=======================================================================
       parameter (nf_max=100)
 
       include 'SIZE'
-      include 'TOTAL'       
+      include 'TOTAL'
       common/expansion_tdata/n_legendre, m_fourier
-      common/expansion_tcoef/coeff_tij(nl_max,nf_max)  
+      common/expansion_tcoef/coeff_tij(nl_max,nf_max)
       integer e,f
 
-      real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt) 
-      real*8 sint, sint1
+      real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt)
+      real*8 sint, sint1, sarea, sarea1
+      real*8 pi
 
+      pi=4.0*atan(1.0)
       ntot=nx1*ny1*nz1*nelt
 
-      zmin=glmin(zm1,ntot)
-      zmax=glmax(zm1,ntot)
-
       call rzero(fmode,ntot)
-         do i0=1,n_legendre
-           do j0=1,m_fourier
-             call nek_mode(fmode,i0,j0)
-             sint1=0.0 
-             do e=1,nelt
-               do f=1,6
-                 sint=0.0
-                 if (cbc(f,e,1).eq.'W  ') then
-                   call col3(cache,fmode,t,ntot)
-                   call surface_int(sint,sarea,cache,e,f)
-                   sint1=sint1+sint
-                 endif
-               enddo
-             enddo
+      do i0=1,n_legendre
+        do j0=1,m_fourier
+          call nek_mode(fmode,i0,j0)
+          sarea1=0.0
+          sint1= 0.0
+          do e=1,nelt
+            do f=1,6
+              sint=0.0
+              sarea=0.0
+              if (cbc(f,e,1).eq.'W  ') then
+                call col3(cache,fmode,t,ntot)
+                call surface_int(sint,sarea,cache,e,f)
+                sint1=sint1+sint
+                sarea1=sarea1+sarea
+              endif
+            enddo
+          enddo
           call  gop(sint1,wtmp,'+  ',1)
-!                                   
-          coeff_tij(i0,j0)=sint1*2.0/(0.5*(zmax-zmin))
-!         Note that R=0.5 here!!!!
-        enddo 
+          call  gop(sarea1,wtmp,'+  ',1)
+!
+          coeff_tij(i0,j0)=sint1*4.0*pi/sarea1
+!
+        enddo
       enddo
 
-c     For Testing 
+c     For Testing
       call nek_testp()
       return
       end