Use modules to compile on the SiCortex machine (from Steve Brandt).

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@202 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files changed, 158 insertions, 145 deletions

diff --git a/src/EHFinder_IsoSurface_optimized.F90 b/src/EHFinder_IsoSurface_optimized.F90
index 90ba5d6..18641a9 100755
--- a/src/EHFinder_IsoSurface_optimized.F90
+++ b/src/EHFinder_IsoSurface_optimized.F90
@@ -6,150 +6,8 @@
 #include "cctk_Arguments.h"
 #include "cctk_Functions.h"
 
-subroutine EHFinder_IsoSurface(CCTK_ARGUMENTS)
-
-  use EHFinder_mod
-  use IsoSurface_mod
-
-  implicit none
-
-  DECLARE_CCTK_PARAMETERS
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_FUNCTIONS
-
-  CCTK_REAL, dimension(:), pointer :: x_coord, y_coord, z_coord
-  CCTK_INT, dimension(:), pointer :: t_index
-  CCTK_INT :: l,contor,p,nsize
-  
-! Find index ranges for points excluding ghost zones and symmetry points for
-! the 3D grid functions.
-#include "include/physical_part.h"
-
-  call CCTK_INFO ('Entered IsoSurface')
-  
-  nx = cctk_lsh(1); ny = cctk_lsh(2); nz = cctk_lsh(3)
-
-  nsize=1000
-  allocate(x_coord(nsize),y_coord(nsize),z_coord(nsize),t_index(5*nsize))
-
-  do p = 1, eh_number_level_sets
-    call triangularization(l,contor,nx,ny,nz,x_coord,y_coord,z_coord, &
-                           t_index,f(:,:,:,p),x(:,1,1),y(1,:,1),z(1,1,:))
-
-    print*, 'contor = ', t_index(1:contor)
-    print*,'The surface consists of', contor/3,'triangles and',l,'coordinates'
-
-    call interpolation_area(cctkGH, l, x_coord,y_coord,z_coord, contor, &
-                            t_index, eh_area(1,p) )
-
-    print*,'The area is ', eh_area(1,p)   
-  end do
-
-  
-
-end subroutine EHFinder_IsoSurface
-
-
-subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,area)
-  
-  use EHFinder_mod
-  implicit none
-  DECLARE_CCTK_PARAMETERS
-!  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_FUNCTIONS 
-
-  CCTK_POINTER, intent(IN) :: cctkGH
-  CCTK_INT, intent(IN) :: npoints, l
-  CCTK_REAL, dimension(:), pointer :: x_coord,y_coord,z_coord
-  CCTK_INT, dimension(:), pointer :: t_index
-  CCTK_REAL, intent(OUT) :: area
-  character(len=200) :: area_interp
-  character(len=128) :: warn_message
-  CCTK_INT :: area_interp_len, local_contor
-  character(len=7) :: area_order
-  CCTK_INT :: interp_handle, table_handle, coord_system_handle
-  CCTK_POINTER, dimension(3) :: interp_coords
-  CCTK_POINTER, dimension(6) :: out_array
-  CCTK_INT, dimension(6) :: in_array
-  CCTK_REAL, dimension(:), allocatable :: metric_xx, metric_yy, metric_zz, metric_xy, metric_yz, metric_xz
-  
-  CCTK_INT, dimension(6), parameter :: out_types = (/ CCTK_VARIABLE_REAL, &
-                                                      CCTK_VARIABLE_REAL, &
-                                                      CCTK_VARIABLE_REAL, &
-                                                      CCTK_VARIABLE_REAL, &
-                                                      CCTK_VARIABLE_REAL, &
-                                                      CCTK_VARIABLE_REAL /)
-  ! Convert the interpolator parameter into a fortran string.
-  call CCTK_FortranString ( area_interp_len, area_interpolator, &
-                                             area_interp )
-
-!allocation of metric arrays
-  allocate(metric_xx(npoints), metric_yy(npoints), metric_zz(npoints), metric_xy(npoints), metric_xz(npoints), metric_yz(npoints))
-
-! Get an interpolation handle.
-  call CCTK_InterpHandle ( interp_handle, area_interp(1:area_interp_len) )
-
-  if ( interp_handle .lt. 0 ) then
-    warn_message = 'Cannot get handle for interpolation. '
-    warn_message = trim(warn_message)//'Forgot to activate an implementation '
-    warn_message = trim(warn_message)//'providing interpolation operators?'
-    call CCTK_WARN( 0, trim(warn_message) )
-  end if
-
-! Write the interpolation order parameter into a string
-  write(area_order,'(a6,i1)') 'order=',area_interpolation_order
-
-! Create a table from the interpolation order string.
-  call Util_TableCreateFromString ( table_handle, area_order )
-  if ( table_handle .lt. 0 ) then
-    call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' )
-  end if
-
-! Get a handle for the coordinate system.
-  call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' )
-  if ( coord_system_handle .lt. 0) then
-    warn_message = 'Cannot get handle for cart3d coordinate system. '
-    warn_message = trim(warn_message)//'Forgot to activate an implementation '
-    warn_message = trim(warn_message)//'providing coordinates?'
-    call CCTK_WARN( 0, trim(warn_message) )
-  endif
-
-! Get the pointers to the interpolation points.
-  interp_coords(1) = CCTK_PointerTo(x_coord)
-  interp_coords(2) = CCTK_PointerTo(y_coord)
-  interp_coords(3) = CCTK_PointerTo(z_coord)
-
-! Get the pointers to the interpolation return arrays.
-  out_array(1) = CCTK_PointerTo(metric_xx)
-  out_array(2) = CCTK_PointerTo(metric_xy)
-  out_array(3) = CCTK_PointerTo(metric_xz)
-  out_array(4) = CCTK_PointerTo(metric_yy)
-  out_array(5) = CCTK_PointerTo(metric_yz)
-  out_array(6) = CCTK_PointerTo(metric_zz)
-!  out_array(7) = CCTK_PointerTo(psii)
-! Get the indices for the grid functions to be interpolated.
-  call CCTK_VarIndex ( in_array(1), 'admbase::gxx' )
-  call CCTK_VarIndex ( in_array(2), 'admbase::gxy' )
-  call CCTK_VarIndex ( in_array(3), 'admbase::gxz' )
-  call CCTK_VarIndex ( in_array(4), 'admbase::gyy' )
-  call CCTK_VarIndex ( in_array(5), 'admbase::gyz' )
-  call CCTK_VarIndex ( in_array(6), 'admbase::gzz' )
-
-  call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, &
-                                 table_handle, coord_system_handle, &
-                                 npoints, CCTK_VARIABLE_REAL, &
-                                 interp_coords, 6, in_array(1:6), &
-                                 6, out_types(1:6), out_array(1:6) )
-  
-  call area_computation(x_coord, y_coord, z_coord, t_index, &
-                        metric_xx, metric_xy, metric_xz, &
-                        metric_yy, metric_yz, metric_zz, &
-                        npoints, l, area )
-
-  deallocate(metric_xx, metric_yy, metric_zz, metric_xz, metric_yz, metric_xy)
-
-  end subroutine interpolation_area
-  
+  module areac
+  contains
   subroutine area_computation(x_coord, y_coord, z_coord, t_index, &
                               g11, g12, g13, &
                               g22, g23, g33, &
@@ -203,8 +61,11 @@ subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,a
     area=area+ds
   end do
   end subroutine area_computation
+  end module
 
-  
+
+ module triangular
+ contains
 subroutine triangularization(n_coords,n_triangles,Nx,Ny,Nz,x_coord,y_coord,z_coord,t_index,grid_value,grid_x,grid_y,grid_z)
 !subroutine triangularization(l,contor,nx,ny,nz,x_coord,y_coord,z_coord,t_index,f(:,:,:,1),x(:,1,1),y(1,:,1),z(1,1,:))
 
@@ -822,3 +683,155 @@ integer, dimension (16,16):: b9(16,16)=reshape&
   close(12)
   print*,'t_index = ',t_index(1:n_triangles)
   end subroutine triangularization
+  end module
+
+module interpa
+contains
+subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,area)
+  
+  use areac
+  use EHFinder_mod
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+!  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS 
+
+  CCTK_POINTER, intent(IN) :: cctkGH
+  CCTK_INT, intent(IN) :: npoints, l
+  CCTK_REAL, dimension(:), pointer :: x_coord,y_coord,z_coord
+  CCTK_INT, dimension(:), pointer :: t_index
+  CCTK_REAL, intent(OUT) :: area
+  character(len=200) :: area_interp
+  character(len=128) :: warn_message
+  CCTK_INT :: area_interp_len, local_contor
+  character(len=7) :: area_order
+  CCTK_INT :: interp_handle, table_handle, coord_system_handle
+  CCTK_POINTER, dimension(3) :: interp_coords
+  CCTK_POINTER, dimension(6) :: out_array
+  CCTK_INT, dimension(6) :: in_array
+  CCTK_REAL, dimension(:), allocatable :: metric_xx, metric_yy, metric_zz, metric_xy, metric_yz, metric_xz
+  
+  CCTK_INT, dimension(6), parameter :: out_types = (/ CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL /)
+  ! Convert the interpolator parameter into a fortran string.
+  call CCTK_FortranString ( area_interp_len, area_interpolator, &
+                                             area_interp )
+
+!allocation of metric arrays
+  allocate(metric_xx(npoints), metric_yy(npoints), metric_zz(npoints), metric_xy(npoints), metric_xz(npoints), metric_yz(npoints))
+
+! Get an interpolation handle.
+  call CCTK_InterpHandle ( interp_handle, area_interp(1:area_interp_len) )
+
+  if ( interp_handle .lt. 0 ) then
+    warn_message = 'Cannot get handle for interpolation. '
+    warn_message = trim(warn_message)//'Forgot to activate an implementation '
+    warn_message = trim(warn_message)//'providing interpolation operators?'
+    call CCTK_WARN( 0, trim(warn_message) )
+  end if
+
+! Write the interpolation order parameter into a string
+  write(area_order,'(a6,i1)') 'order=',area_interpolation_order
+
+! Create a table from the interpolation order string.
+  call Util_TableCreateFromString ( table_handle, area_order )
+  if ( table_handle .lt. 0 ) then
+    call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' )
+  end if
+
+! Get a handle for the coordinate system.
+  call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' )
+  if ( coord_system_handle .lt. 0) then
+    warn_message = 'Cannot get handle for cart3d coordinate system. '
+    warn_message = trim(warn_message)//'Forgot to activate an implementation '
+    warn_message = trim(warn_message)//'providing coordinates?'
+    call CCTK_WARN( 0, trim(warn_message) )
+  endif
+
+! Get the pointers to the interpolation points.
+  interp_coords(1) = CCTK_PointerTo(x_coord)
+  interp_coords(2) = CCTK_PointerTo(y_coord)
+  interp_coords(3) = CCTK_PointerTo(z_coord)
+
+! Get the pointers to the interpolation return arrays.
+  out_array(1) = CCTK_PointerTo(metric_xx)
+  out_array(2) = CCTK_PointerTo(metric_xy)
+  out_array(3) = CCTK_PointerTo(metric_xz)
+  out_array(4) = CCTK_PointerTo(metric_yy)
+  out_array(5) = CCTK_PointerTo(metric_yz)
+  out_array(6) = CCTK_PointerTo(metric_zz)
+!  out_array(7) = CCTK_PointerTo(psii)
+! Get the indices for the grid functions to be interpolated.
+  call CCTK_VarIndex ( in_array(1), 'admbase::gxx' )
+  call CCTK_VarIndex ( in_array(2), 'admbase::gxy' )
+  call CCTK_VarIndex ( in_array(3), 'admbase::gxz' )
+  call CCTK_VarIndex ( in_array(4), 'admbase::gyy' )
+  call CCTK_VarIndex ( in_array(5), 'admbase::gyz' )
+  call CCTK_VarIndex ( in_array(6), 'admbase::gzz' )
+
+  call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, &
+                                 table_handle, coord_system_handle, &
+                                 npoints, CCTK_VARIABLE_REAL, &
+                                 interp_coords, 6, in_array(1:6), &
+                                 6, out_types(1:6), out_array(1:6) )
+  
+  call area_computation(x_coord, y_coord, z_coord, t_index, &
+                        metric_xx, metric_xy, metric_xz, &
+                        metric_yy, metric_yz, metric_zz, &
+                        npoints, l, area )
+
+  deallocate(metric_xx, metric_yy, metric_zz, metric_xz, metric_yz, metric_xy)
+
+  end subroutine interpolation_area
+  end module
+
+subroutine EHFinder_IsoSurface(CCTK_ARGUMENTS)
+
+  use interpa
+  use EHFinder_mod
+  use IsoSurface_mod
+  use triangular
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_REAL, dimension(:), pointer :: x_coord, y_coord, z_coord
+  CCTK_INT, dimension(:), pointer :: t_index
+  CCTK_INT :: l,contor,p,nsize
+  
+! Find index ranges for points excluding ghost zones and symmetry points for
+! the 3D grid functions.
+#include "include/physical_part.h"
+
+  call CCTK_INFO ('Entered IsoSurface')
+  
+  nx = cctk_lsh(1); ny = cctk_lsh(2); nz = cctk_lsh(3)
+
+  nsize=1000
+  allocate(x_coord(nsize),y_coord(nsize),z_coord(nsize),t_index(5*nsize))
+
+  do p = 1, eh_number_level_sets
+    call triangularization(l,contor,nx,ny,nz,x_coord,y_coord,z_coord, &
+                           t_index,f(:,:,:,p),x(:,1,1),y(1,:,1),z(1,1,:))
+
+    print*, 'contor = ', t_index(1:contor)
+    print*,'The surface consists of', contor/3,'triangles and',l,'coordinates'
+
+    call interpolation_area(cctkGH, l, x_coord,y_coord,z_coord, contor, &
+                            t_index, eh_area(1,p) )
+
+    print*,'The area is ', eh_area(1,p)   
+  end do
+
+  
+
+end subroutine EHFinder_IsoSurface
+
+