Cleaning up and adding comments to make it more readable.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/ADMConstraints/trunk@43 b7a48df3-cbbf-4440-997f-b4b717c9f7fc

1 files changed, 59 insertions, 47 deletions

diff --git a/src/ADMConstraints.F b/src/ADMConstraints.F
index 1a999dd..7c61d21 100644
--- a/src/ADMConstraints.F
+++ b/src/ADMConstraints.F
@@ -47,7 +47,8 @@ c@@*/
       DECLARE_CCTK_PARAMETERS
       DECLARE_CCTK_FUNCTIONS
 
-      integer :: i,j,k,i1,i2,j1,j2,k1,k2
+      integer :: i,j,k
+      integer :: nx,ny,nz
 
 c     Stencil width used for calculating constraints (for outer boundary condition)
       integer, dimension(3),parameter :: sw = 1
@@ -71,33 +72,31 @@ c     Macros from Standard Einstein
 #include "CactusEinstein/Einstein/src/macro/DETG_declare.h"
 #include "CactusEinstein/Einstein/src/macro/UPPERMET_declare.h"
 
+c     Grid parameters.
+
       dx = cctk_delta_space(1)
       dy = cctk_delta_space(2)
       dz = cctk_delta_space(3)
 
-      i1 = 2
-      i2 = cctk_lsh(1)-1
-      j1 = 2
-      j2 = cctk_lsh(2)-1
-      k1 = 2
-      k2 = cctk_lsh(3)-1
-  
-      do i=i1,i2
-        do j=j1,j2
-          do k=k1,k2
-            ham(i,j,k)  = 0.0D0
-            momx(i,j,k) = 0.0D0
-            momy(i,j,k) = 0.0D0
-            momz(i,j,k) = 0.0D0
-          end do
-        end do
-      end do
+      nx = cctk_lsh(1)
+      ny = cctk_lsh(2)
+      nz = cctk_lsh(3)
+
+c     Fill with zeros.
+
+      ham = 0.0D0
+
+      momx = 0.0D0
+      momy = 0.0D0
+      momz = 0.0D0
+
+c     Calculate constraints.
 
-      Pi =  ACOS(-1D0)
+      pi = acos(-1.0D0)
 
-      do k = k1,k2
-         do j = j1,j2
-            do i = i1,i2
+      do k=2,nz-1
+         do j=2,ny-1
+            do i=2,nx-1
 
                ialp  = 1.0D0/alp(i,j,k)
                ialp2 = ialp**2
@@ -105,20 +104,25 @@ c     Macros from Standard Einstein
 c              Calculate the stress energy tensor at this point
 c              ------------------------------------------------
 
+c              Fill with zeroes.
+
                Ttt = 0.0D0; Ttx = 0.0D0; Tty = 0.0D0; Ttz = 0.0D0
-               Txx = 0.0D0; Txy = 0.0D0; Txz = 0.0D0; Tyy = 0.0D0
-               Tyz = 0.0D0; Tzz = 0.0D0
-               
-c              This may be needed for CalcTmunu
-c              --------------------------------
+
+               Txx = 0.0D0; Txy = 0.0D0; Txz = 0.0D0
+               Tyy = 0.0D0; Tyz = 0.0D0; Tzz = 0.0D0
+
+c              Inverse metric may be needed for CalcTmunu.
 
 #include "CactusEinstein/Einstein/src/macro/DETG_guts.h"
 #include "CactusEinstein/Einstein/src/macro/UPPERMET_guts.h"
 
                det = DETG_DETCG
+
                uxx = UPPERMET_UXX; uxy = UPPERMET_UXY; uxz = UPPERMET_UXZ
                uyy = UPPERMET_UYY; uyz = UPPERMET_UYZ; uzz = UPPERMET_UZZ
 
+c              Call macro for stress energy tensor.
+
 #include "CalcTmunu.inc"
 
 c              Calculate the hamiltonian constraint
@@ -154,13 +158,15 @@ c                  = (T_00 - 2 beta^i T_{i0} + beta^i beta^j T_{ij})/alpha^2
 
                ham(i,j,k) = HAMADM_HAMADM - 16.0D0*pi*m_rho
 
+c              From thorn MAHC (this overwrites all we just did).
+
 #ifdef THORN_MAHC
-               if ((some_hydro .eq. MAHCHYDRO).or.
-     .              (some_hydro .eq. MAHCHYDRODUST)) then
-                  ham(i,j,k) = ADM_ham(i,j,k) - 16.0d0 * pi *
-     .                 ( (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k)) *
+               if ((some_hydro.eq.MAHCHYDRO).or.
+     .             (some_hydro.eq.MAHCHYDRODUST)) then
+                  ham(i,j,k) = ADM_ham(i,j,k) - 16.0d0*pi*
+     .                 ((rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k))*
      .                 w_lorentz(i,j,k)**2 - press(i,j,k))
-               endif
+               end if
 #endif
 
 c              Calculate the Momentum constraints
@@ -206,37 +212,43 @@ c                  = - (T_{i0} - beta^j T_{ij})/alpha
                momy(i,j,k) = MOMYADM_MOMYADM - 8.0D0*pi*m_sy
                momz(i,j,k) = MOMZADM_MOMZADM - 8.0D0*pi*m_sz
 
+c              From thorn MAHC (this overwrites all we just did).
+
 #ifdef THORN_MAHC
-               if ((some_hydro .eq. MAHCHYDRO).or.
-     .              (some_hydro .eq. MAHCHYDRODUST)) then
-                  momx(i,j,k) = ADM_momx(i,j,k) -
-     .                 8.0d0 * pi * sx(i,j,k) / sqrt(det)
-                  momy(i,j,k) = ADM_momy(i,j,k) -
-     .                 8.0d0 * pi * sy(i,j,k) / sqrt(det)
-                  momz(i,j,k) = ADM_momz(i,j,k) -
-     .                 8.0d0 * pi * sz(i,j,k) / sqrt(det)
-               endif
+               if ((some_hydro.eq.MAHCHYDRO).or.
+     .             (some_hydro.eq.MAHCHYDRODUST)) then
+                  momx(i,j,k) = ADM_momx(i,j,k)
+     .                 - 8.0d0*pi*sx(i,j,k)/sqrt(det)
+                  momy(i,j,k) = ADM_momy(i,j,k)
+     .                 - 8.0d0*pi*sy(i,j,k)/sqrt(det)
+                  momz(i,j,k) = ADM_momz(i,j,k)
+     .                 - 8.0d0*pi*sz(i,j,k)/sqrt(det)
+               end if
 #endif
-               
-            enddo
-         enddo
-      enddo
+
+            end do
+         end do 
+      end do
 
 #include "CactusEinstein/Einstein/src/macro/DETG_undefine.h"
+#include "CactusEinstein/Einstein/src/macro/UPPERMET_undefine.h"
 #include "CactusEinstein/Einstein/src/macro/HAMADM_undefine.h"
 #include "CactusEinstein/Einstein/src/macro/MOMXADM_undefine.h"
 #include "CactusEinstein/Einstein/src/macro/MOMYADM_undefine.h"
 #include "CactusEinstein/Einstein/src/macro/MOMZADM_undefine.h"
-#include "CactusEinstein/Einstein/src/macro/UPPERMET_undefine.h"
 
-c     Synchronize and apply flat boundary conditions
+c     Apply symmetry boundary conditions.
 
       call CartSymGN(ierr,cctkGH,"admconstraints::admconstraints")
 
+c     Apply flat boundary conditions at outer boundaries.
+
       if (CCTK_Equals(bound,"flat") == 1) then
          call BndFlatGN(ierr,cctkGH,sw,"admconstraints::admconstraints")
       end if
 
+c     Synchronize.
+
       if (constraint_communication.eq.1) then
          call CCTK_SyncGroup(cctkGH,"admconstraints::admconstraints")
       end if