Moving key routines to a separate source folder - Part 1

src/nek/NEKMOOSE

+      parameter (nl_max=100)
+      parameter (nf_max=100)
+      common/expansion_tdata/n_legendre, m_fourier
+      common/expansion_tcoef/coeff_tij(nl_max,nf_max)
+      common/expansion_fcoef/coeff_fij(nl_max,nf_max)
+      common/expansion_recon/flux_recon(lx1,ly1,lz1,lelt)
+      common/test_passing/ flux_moose, temp_nek
+

src/nek/nek_moose.f

+C=======================================================================
+      subroutine nek_init_step()
+      include 'SIZE'
+      include 'TSTEP'
+      include 'INPUT'
+      include 'CTIMER'
+      real*4 papi_mflops
+      integer*8 papi_flops
+      integer icall, kstep, i, pstep
+      common /cht_coupler/ pstep
+      save icall
+ 
+      if (icall.eq.0) then
+      call nekgsync()
+      if (instep.eq.0) then
+        if(nid.eq.0) write(6,'(/,A,/,A,/)')
+     &     ' nsteps=0 -> skip time loop',
+     &     ' running solver in post processing mode'
+      else
+        if(nio.eq.0) write(6,'(/,A,/)') 'Starting time loop ...'
+      endif
+      isyc  = 0
+      itime = 0
+      if(ifsync) isyc=1
+      itime = 1
+      call nek_comm_settings(isyc,itime)
+      call nek_comm_startstat()
+      istep  = 0
+      icall  = 1
+      endif
+
+      istep=istep+1
+
+      if (lastep .eq. 1) then
+        pstep=2
+      else
+        call nek_advance
+        pstep=2
+      endif
+
+      return
+      end
+C=======================================================================
+      subroutine nek_step()
+
+      include 'SIZE'
+      include 'TSTEP'
+      include 'INPUT'
+      include 'CTIMER'
+      common /cht_coupler/ pstep
+      integer pstep
+
+      pstep=pstep+1
+      call heat(pstep)
+
+      return
+      end
+C=======================================================================
+      subroutine nek_finalize_step()
+      include 'SIZE'
+      include 'TSTEP'
+      include 'INPUT'
+      include 'CTIMER'
+      common /cht_coupler/ pstep
+      integer pstep 
+
+      real*4 papi_mflops
+      integer*8 papi_flops
+      integer icall, kstep, knstep, i
+
+      if (param(103).gt.0)   call q_filter      (param(103))
+      call setup_convect (pstep)  ! Save convective velocity _after_ filter
+
+      call userchk
+      call prepost (.false.,'his')
+      call in_situ_check()
+
+      if (mod(istep,nsteps).eq.0) lastep=1
+      call nek_comm_settings(isyc,0)
+      call comment
+
+      return
+      end
+C=======================================================================
+      subroutine nek_expansion()
+
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKMOOSE' 
+      integer e,f
+
+      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
+
+      call rzero(fmode,ntot)
+      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)
+          call  gop(sarea1,wtmp,'+  ',1)
+!
+          coeff_tij(i0,j0)=sint1*4.0*pi/sarea1
+!
+        enddo
+      enddo
+
+c     For Testing
+      call nek_testp()
+      return
+      end
+C=======================================================================
+      subroutine nek_diag()
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKMOOSE' 
+      common/diags_coeff/diag_c(nl_max,nf_max)
+      integer e,f
+
+      real*8 fmode(lx1,ly1,lz1,lelt)
+      real*8 cache(lx1,ly1,lz1,lelt)
+      real*8 sint, sint1
+
+      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,fmode,ntot)
+                   call surface_int(sint,sarea,cache,e,f)
+                   sint1=sint1+sint
+                 endif
+               enddo
+             enddo
+          call  gop(sint1,wtmp,'+  ',1)
+          diag_c(i0,j0)=sint1*2.0/(0.5*(zmax-zmin))
+          if (nid.eq.0) write(6,*)i0,j0,diag_c(i0,j0)
+        enddo
+      enddo
+
+      return
+      end
+C=======================================================================
+      subroutine nek_testp()
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKMOOSE'      
+      integer e,f
+      real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt)
+      real*8 fun(lx1,ly1,lz1,lelt) 
+      ntot=nx1*ny1*nz1*nelt
+      call rzero(fmode,ntot)
+      call rzero(fun,ntot)   
+         do i0=1,n_legendre
+           do j0=1,m_fourier
+             call nek_mode(fmode,i0,j0)
+             do i=1,ntot
+             fun(i,1,1,1)=fun(i,1,1,1)+fmode(i,1,1,1)*coeff_tij(i0,j0) 
+             enddo  
+          enddo
+         enddo
+      call rzero(cache,ntot) 
+      call sub3(cache,fun,t,ntot)
+      do i=1,ntot
+        c1=cache(i,1,1,1)**2
+        cache(i,1,1,1)=c1 
+      enddo
+ 
+      sint1=0.0
+      sarea1=0.0
+      do e=1,lelt
+        do f=1,6
+          call surface_int(sint,sarea,cache,e,f)
+          if (cbc(f,e,1).eq.'W  ') then
+           sint1=sint1+sint
+           sarea1=sarea1+sarea
+          endif
+         enddo
+      enddo
+      call  gop(sint1,wtmp,'+  ',1)
+      call  gop(sarea1,wtmp,'+  ',1)
+
+      er_avg=sqrt(sint1/sarea1)
+      if (nid.eq.0) write(6,*)"Error: ",er_avg
+      return
+      end
+C=======================================================================
+      subroutine nek_mode(fmode,im,jm)
+      include 'SIZE'
+      include 'TOTAL'
+      integer e,f
+      real*8 fmode(lx1,ly1,lz1,lelt) 
+      ntot=nx1*ny1*nz1*nelt
+      zmin=glmin(zm1,ntot)
+      zmax=glmax(zm1,ntot)
+             do e=1,nelt
+               do f=1,6
+                 if (cbc(f,e,1).eq.'W  ') then
+                   call dsset(nx1,ny1,nz1)
+                   iface  = eface1(f)
+                   js1    = skpdat(1,iface)
+                   jf1    = skpdat(2,iface)
+                   jskip1 = skpdat(3,iface)
+                   js2    = skpdat(4,iface)
+                   jf2    = skpdat(5,iface)
+                   jskip2 = skpdat(6,iface)
+                   do j2=js2,jf2,jskip2
+                     do j1=js1,jf1,jskip1
+                       x=xm1(j1,j2,1,e)
+                       y=ym1(j1,j2,1,e)
+                       z=zm1(j1,j2,1,e)
+                       z_leg=2*((z-zmin)/zmax)-1
+                       theta=atan2(y,x)
+                       fmode(j1,j2,1,e)=
+     &                 pl_leg(z_leg,im-1)*fl_four(theta,jm-1)
+                     enddo
+                   enddo
+                 endif
+              enddo
+            enddo
+      return
+      end 
+c-----------------------------------------------------------------------
+      subroutine flux_reconstruction()
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKMOOSE'  
+      integer e,f
+      real*8 coeff_base, flux_base, flux_moose 
+      real*8 fmode(lx1,ly1,lz1,lelt)
+      ntot=nx1*ny1*nz1*nelt
+
+c     --------------------------------------
+c     The flux from MOOSE must have proper sign
+c     --------------------------------------
+      coeff_base=-1.0 
+      flux_base=0.25 ! from energy conservation in the problem
+c     --------------------------------------
+
+      call rzero(flux_recon,ntot)
+         do i0=1,n_legendre
+           do j0=1,m_fourier
+             call rzero(fmode,ntot)
+             call nek_mode(fmode,i0,j0)
+             do i=1,ntot
+             flux_recon(i,1,1,1)= flux_recon(i,1,1,1)
+     &                  + coeff_base*fmode(i,1,1,1)*coeff_fij(i0,j0)
+             enddo
+          enddo
+         enddo
+
+c---- Below is for testing
+c             do i=1,ntot
+c             flux_recon(i,1,1,1)= 0.0
+c             enddo
+c--------------------------
+
+c---- Renormalization
+      sint1=0.0
+      sarea1=0.0
+
+      do e=1,lelt
+        do f=1,6
+          call surface_int(sint,sarea,flux_recon,e,f)
+          if (cbc(f,e,1).eq.'W  ') then
+           sint1=sint1+sint
+           sarea1=sarea1+sarea
+          endif
+         enddo
+      enddo
+
+      call  gop(sint1,wtmp,'+  ',1)
+      call  gop(sarea1,wtmp,'+  ',1)
+
+      flux_moose=sint1/sarea1
+      do i=1,ntot
+         flux_recon(i,1,1,1)= flux_recon(i,1,1,1)*flux_base/flux_moose
+      enddo
+c-----------------------
+      return
+      end
+c-----------------------------------------------------------------------
+! calculates Legendre polynomials using a recurrence relationship. If
+! n > the maximum Legendre order, the function returns 0.0.
+      function pl_leg(x,n)
+      parameter (nl_max=100) ! should never be bigger than 100
+      real*8 pl,pl_leg
+      real*8 x
+      real*8 pln(0:n)
+      integer n, k
+
+      pln(0) = 1.0
+      pln(1) = x
+
+      if (n.le.1) then
+        pl = pln(n)
+      else if (n.le.nl_max) then
+        do k=1,n-1
+          pln(k+1) = ((2.0*k+1.0)*x*pln(k)-dble(k)*pln(k-1))/(dble(k+1))
+        end do
+        pl = pln(n)
+      else
+        pl = 0.0
+      end if
+
+      pl_leg = pl*sqrt(dble(2*n+1)/2.0)
+
+      return
+      end
+c-----------------------------------------------------------------------
+! calculates Fourier polynomials Fn(x)
+      function fl_four(x,n)
+      real*8 fl_four,pi
+      real*8 x, A
+      integer n, k
+
+      pi=4.0*atan(1.0)
+      A=1.0/sqrt(pi)
+
+      if (n.eq.0) A=1.0/sqrt(2*pi)
+      fl_four=A*cos(n*x)
+
+      return
+      end
+c-----------------------------------------------------------------------