Build spring clean

aclocal.m4

@@ -428,10 +428,10 @@ AC_MSG_NOTICE([Using PETSC_DIR=$PETSC_DIR])
 PETSC_LINK_LIBS=`make -s -f petsc_makefile_old getlinklibs 2> /dev/null || make -s -f petsc_makefile getlinklibs`
 LIBS="$PETSC_LINK_LIBS $LIBS"
 
-# need to add -Iinclude/ to what we get from petsc, so we can use our own petsc_legacy.h wrapper
+# need to add -I$PWD/include/ to what we get from petsc, so we can use our own petsc_legacy.h wrapper
 PETSC_INCLUDE_FLAGS=`make -s -f petsc_makefile_old getincludedirs 2> /dev/null || make -s -f petsc_makefile getincludedirs`
-CPPFLAGS="$CPPFLAGS $PETSC_INCLUDE_FLAGS -Iinclude/"
-FCFLAGS="$FCFLAGS $PETSC_INCLUDE_FLAGS -Iinclude/"
+CPPFLAGS="$CPPFLAGS $PETSC_INCLUDE_FLAGS -I$PWD/include/"
+FCFLAGS="$FCFLAGS $PETSC_INCLUDE_FLAGS -I$PWD/include/"
 
 # first check we have the right petsc version
 AC_COMPUTE_INT(PETSC_VERSION_MAJOR, "PETSC_VERSION_MAJOR", [#include "petscversion.h"], 

assemble/Free_Surface.F90

@@ -195,7 +195,7 @@ subroutine insert_original_distance_to_bottom(state)
     type(scalar_field) :: original_bottomdist, original_bottomdist_remap
 
     if (.not. has_scalar_field(state, "OriginalDistanceToBottom")) then
-       ewrite(2, *), "Inserting OriginalDistanceToBottom field into state."   
+       ewrite(2, *) "Inserting OriginalDistanceToBottom field into state."   
        bottomdist => extract_scalar_field(state, "DistanceToBottom")
        call allocate(original_bottomdist, bottomdist%mesh, "OriginalDistanceToBottom")
        call zero(original_bottomdist)
@@ -204,7 +204,7 @@ subroutine insert_original_distance_to_bottom(state)
        call deallocate(original_bottomdist)
 
        ! We also cache  the OriginalDistanceToBottom on the pressure mesh
-       ewrite(2, *), "Inserting OriginalDistanceToBottomPressureMesh field into state."   
+       ewrite(2, *) "Inserting OriginalDistanceToBottomPressureMesh field into state."   
        p_mesh => extract_pressure_mesh(state)
        call allocate(original_bottomdist_remap, p_mesh, "OriginalDistanceToBottomPressureMesh")
        call remap_field(original_bottomdist, original_bottomdist_remap)

assemble/Geostrophic_Pressure.F90

@@ -1986,7 +1986,7 @@ subroutine derive_interpolated_p_dirichlet_multiple(base_ps, base_positions, ps,
 
     call clear_boundary_conditions(ps)
 
-    ewrite(2, *), "Adding strong Dirichlet bc for field " // trim(ps(1)%name)
+    ewrite(2, *) "Adding strong Dirichlet bc for field " // trim(ps(1)%name)
     do i = 1, size(surface_eles)
       surface_eles(i) = i
     end do
@@ -1999,7 +1999,7 @@ subroutine derive_interpolated_p_dirichlet_multiple(base_ps, base_positions, ps,
     call insert_surface_field(ps(1), 1, b_ps(1))
     call deallocate(b_ps(1))
     do i = 2, size(ps)
-      ewrite(2, *), "Adding strong Dirichlet bc for field " // trim(ps(i)%name) 
+      ewrite(2, *) "Adding strong Dirichlet bc for field " // trim(ps(i)%name) 
       call add_boundary_condition_surface_elements(ps(i), "InterpolatedBoundary", "dirichlet", surface_eles) 
       ewrite_minmax(b_ps(i))
       call insert_surface_field(ps(i), 1, b_ps(i))
@@ -2165,7 +2165,7 @@ subroutine decompose_p_mean_multiple(matrices, base_ps, positions, ps, solver_pa
       end do
       do i = 1, size(bc_ps, 1)
         call clear_boundary_conditions(bc_ps(i))
-        ewrite(2, *), "Adding strong Dirichlet bc for field " // trim(bc_ps(i)%name)
+        ewrite(2, *) "Adding strong Dirichlet bc for field " // trim(bc_ps(i)%name)
         call add_boundary_condition_surface_elements(bc_ps(i), "ConstantBoundary", "dirichlet", surface_eles) 
         call get_boundary_condition(bc_ps(i), 1, surface_mesh = bc_mesh, &
           & surface_element_list = bc_surface_element_list) 
@@ -2373,7 +2373,7 @@ subroutine decompose_p_optimal_multiple(matrices, base_ps, positions, ps, solver
       end do
       do i = 1, size(bc_ps, 1)
         call clear_boundary_conditions(bc_ps(i))
-        ewrite(2, *), "Adding strong Dirichlet bc for field " // trim(bc_ps(i)%name)
+        ewrite(2, *) "Adding strong Dirichlet bc for field " // trim(bc_ps(i)%name)
         call add_boundary_condition_surface_elements(bc_ps(i), "ConstantBoundary", "dirichlet", surface_eles) 
         call get_boundary_condition(bc_ps(i), 1, surface_mesh = bc_mesh, &
           & surface_element_list = bc_surface_element_list) 

assemble/tests/test_main.cpp

@@ -59,7 +59,11 @@ int main(int argc, char **argv)
 #ifdef HAVE_MPI
   MPI_Init(&argc, &argv);
   // Undo some MPI init shenanigans
-  chdir(getenv("PWD"));
+  int ierr = chdir(getenv("PWD"));
+  if (ierr == -1) {
+        std::cerr << "Unable to switch to directory " << getenv("PWD");
+        abort();
+  }
 #endif
 #ifdef HAVE_PETSC  
   PetscInitialize(&argc, &argv, NULL, PETSC_NULL);

climatology/create_climatology_atlas.cpp

@@ -461,7 +461,8 @@ int main(){
       // cout<<"Reading level "<<k<<endl;
       for(size_t j=0; j<(size_t)jdim0; j++)
   for(size_t i=0; i<(size_t)idim0; i++){
-    fscanf(fp, "%f", &(dat[m][k][j*idim0+i]));
+    int ierr = fscanf(fp, "%f", &(dat[m][k][j*idim0+i]));
+    if (ierr<1) cerr<<"File read error";
   }
     }
     fclose(fp);
@@ -490,7 +491,8 @@ int main(){
       // cout<<"Reading level "<<k<<endl;
       for(size_t j=0; j<(size_t)jdim0; j++)
   for(size_t i=0; i<(size_t)idim0; i++){
-    fscanf(fp, "%f", &(dat[m][k][j*idim0+i]));
+    int ierr = fscanf(fp, "%f", &(dat[m][k][j*idim0+i]));
+    if (ierr<1) cerr<<"File read error";
   }
     }
     fclose(fp);
@@ -528,7 +530,8 @@ int main(){
       for(size_t j=0; j<(size_t)jdim0; j++)
   for(size_t i=0; i<(size_t)idim0; i++){
     float v;
-    fscanf(fp, "%f", &v);
+    int ierr = fscanf(fp, "%f", &v);
+    if (ierr<1) cerr<<"File read error";
     if(k>=24){
       dat[m][k-24][j*idim0+i] = v;
     }
@@ -560,7 +563,8 @@ int main(){
       for(size_t j=0; j<(size_t)jdim0; j++)
   for(size_t i=0; i<(size_t)idim0; i++){
     float v;
-    fscanf(fp, "%f", &v);
+    int ierr = fscanf(fp, "%f", &v);
+    if (ierr<1) cerr<<"File read error";
     if(k>=24){
       dat[m][k-24][j*idim0+i] = v;
     }

configure

@@ -12641,10 +12641,10 @@ $as_echo "$as_me: Using PETSC_DIR=$PETSC_DIR" >&6;}
 PETSC_LINK_LIBS=`make -s -f petsc_makefile_old getlinklibs 2> /dev/null || make -s -f petsc_makefile getlinklibs`
 LIBS="$PETSC_LINK_LIBS $LIBS"
 
-# need to add -Iinclude/ to what we get from petsc, so we can use our own petsc_legacy.h wrapper
+# need to add -I$PWD/include/ to what we get from petsc, so we can use our own petsc_legacy.h wrapper
 PETSC_INCLUDE_FLAGS=`make -s -f petsc_makefile_old getincludedirs 2> /dev/null || make -s -f petsc_makefile getincludedirs`
-CPPFLAGS="$CPPFLAGS $PETSC_INCLUDE_FLAGS -Iinclude/"
-FCFLAGS="$FCFLAGS $PETSC_INCLUDE_FLAGS -Iinclude/"
+CPPFLAGS="$CPPFLAGS $PETSC_INCLUDE_FLAGS -I$PWD/include/"
+FCFLAGS="$FCFLAGS $PETSC_INCLUDE_FLAGS -I$PWD/include/"
 
 # first check we have the right petsc version
 if ac_fn_c_compute_int "$LINENO" ""PETSC_VERSION_MAJOR"" "PETSC_VERSION_MAJOR"        "#include \"petscversion.h\""; then :

diagnostics/Tidal_Diagnostics.F90

@@ -103,7 +103,7 @@ subroutine calculate_free_surface_history(state, s_field)
     real :: spin_up_time, current_time, timestep
     real, dimension(:), allocatable :: saved_snapshots_times
 
-    ewrite(3,*),'in free_surface_history_diagnostics'
+    ewrite(3,*) 'in free_surface_history_diagnostics'
 
     fs_field => extract_scalar_field(state,"FreeSurface",stat)
     call halo_update(fs_field)
@@ -158,13 +158,13 @@ subroutine calculate_free_surface_history(state, s_field)
    ! Note this is after the options check above as we add options to the tree if
    ! they aren't there and they are needed in tidal_harmonics diagnostics
    if(current_time+timestep<spin_up_time) then
-       ewrite(4,*), "Still spinning up."
+       ewrite(4,*) "Still spinning up."
        return
    endif
 
    ! check if we want to save the current timestep at all
    if (mod(timestep_counter,stride)/=0) then
-       ewrite(4,*), "Ignoring this timestep"
+       ewrite(4,*) "Ignoring this timestep"
        return
    end if
 
@@ -181,7 +181,7 @@ subroutine calculate_free_surface_history(state, s_field)
   saved_snapshots_times(new_snapshot_index)=current_time+timestep
   call set_option(trim(base_path) // "saved_snapshots_times", saved_snapshots_times, stat)
   assert(any(stat == (/SPUD_NO_ERROR, SPUD_NEW_KEY_WARNING/)))
-  ewrite(4,*), 'Filling history level: ', min(timestep_counter/stride+1,levels), '/', levels
+  ewrite(4,*) 'Filling history level: ', min(timestep_counter/stride+1,levels), '/', levels
 
   ! lets copy a snapshot of freesurface to s_field(new_snapshot_index)
   hist_fs_field => extract_scalar_field(state,'harmonic'//int2str(new_snapshot_index))
@@ -220,7 +220,7 @@ subroutine calculate_tidal_harmonics(state, s_field)
 
    ! Only if Harmonics weren't already calculated in this timestep
    if (last_update/=timestep) then 
-      ewrite(3,*), "In tidal_harmonics"
+      ewrite(3,*) "In tidal_harmonics"
       allocate(harmonic_fields(get_number_of_harmonic_fields(state)))
       allocate(sigma(nLevels_))
 
@@ -229,9 +229,9 @@ subroutine calculate_tidal_harmonics(state, s_field)
       if (.not. ignoretimestep) then
          ! Initialize the harmonic fields and frequencies if not already done.
          call getHarmonicFields(state, harmonic_fields, nohfs, sigma, M)
-         ewrite(4,*), 'Frequencies to analyse:'
+         ewrite(4,*) 'Frequencies to analyse:'
          do i=1,M
-            ewrite(4,*), sigma(i)
+            ewrite(4,*) sigma(i)
          end do
          ! Calculate harmonics and update (all!) harmonic fields
          call update_harmonic_fields(state, saved_snapshots_times, size(saved_snapshots_times), current_snapshot_index, sigma, M, harmonic_fields, nohfs)
@@ -268,12 +268,12 @@ subroutine getFreeSurfaceHistoryData(state, ignoretimestep, saved_snapshots_time
    call get_option(trim(free_surface_history_path) // "levels", levels)
 
    if( mod(timestep_counter,stride)/=0) then
-        ewrite(4,*),'Do nothing in this timestep.'
+        ewrite(4,*) 'Do nothing in this timestep.'
         ignoretimestep=.true.
         return
    end if
    if(timestep_counter/stride+1 .lt. levels) then
-        ewrite(4,*), 'Do nothing until levels are filled up.'
+        ewrite(4,*) 'Do nothing until levels are filled up.'
         ignoretimestep=.true.
         return
    end if
@@ -350,8 +350,8 @@ subroutine getHarmonicFields(state, harmonic_fields, nohfs, sigma, M)
            end if
        end if
     end do s_field_loop
-    ewrite(4,*), 'Found ',  nohfs, ' constituents to analyse.' 
-    ewrite(4,*), 'Found ', M, ' frequencies to analyse.'
+    ewrite(4,*) 'Found ',  nohfs, ' constituents to analyse.' 
+    ewrite(4,*) 'Found ', M, ' frequencies to analyse.'
 
    if ( M .le. 0 ) then
       FLExit("Internal error in calculate_tidal_harmonics(). No harmonic constituents were found in option tree.")

error_measures/tests/test_main.cpp

@@ -53,7 +53,11 @@ int main(int argc, char **argv) {
 #ifdef HAVE_MPI
   MPI_Init(&argc, &argv);
   // Undo some MPI init shenanigans
-  chdir(getenv("PWD"));
+  int ierr = chdir(getenv("PWD"));
+  if (ierr == -1) {
+        std::cerr << "Unable to switch to directory " << getenv("PWD");
+        abort();
+  }
 #endif
 #ifdef HAVE_PETSC  
   PetscInitialize(&argc, &argv, NULL, PETSC_NULL);

femtools/Coordinates.F90

@@ -1004,7 +1004,7 @@ subroutine higher_order_sphere_projection(positions, s_positions)
     type(scalar_field):: radius, s_radius
     real, dimension(positions%dim):: xyz
 
-    ewrite(1,*), 'In higher_order_sphere_projection'
+    ewrite(1,*) 'In higher_order_sphere_projection'
 
     call allocate(s_radius, s_positions%mesh, "HigherOrderRadius")
     radius=magnitude(positions)

femtools/Diagnostic_variables.F90

@@ -592,13 +592,13 @@ subroutine initialise_diagnostics(filename, state)
         if(stat_mesh(mesh)) then
           column = column + 1
           buffer = field_tag(name = mesh%name, column = column, statistic = "nodes")
-          write(default_stat%diag_unit, "(a)"), trim(buffer)
+          write(default_stat%diag_unit, "(a)") trim(buffer)
           column = column + 1
           buffer = field_tag(name = mesh%name, column = column, statistic = "elements")
-          write(default_stat%diag_unit, "(a)"), trim(buffer)
+          write(default_stat%diag_unit, "(a)") trim(buffer)
           column = column + 1
           buffer = field_tag(name = mesh%name, column = column, statistic = "surface_elements")
-          write(default_stat%diag_unit, "(a)"), trim(buffer)
+          write(default_stat%diag_unit, "(a)") trim(buffer)
         end if
       end do
 

femtools/Exodusii_C_Interface.c

@@ -8,6 +8,8 @@
 #include "exodusII.h"
 #endif
 
+void FLExit(const char*);
+
 /* Open ExodusII File for reading */
 int c_read_ex_open(const char *path, int mode, int *comp_ws, int *io_ws, float *version)
 {

femtools/Read_GMSH.F90

@@ -387,7 +387,7 @@ subroutine read_nodes_coords( fd, filename, gmshFormat, nodes )
 
     ! Skip newline character when in binary mode
     if( gmshFormat == binaryFormat ) then
-       read(fd), newlineChar
+       read(fd) newlineChar
        call ascii_formatting(fd, filename, "read")
     end if
 
@@ -477,7 +477,7 @@ subroutine read_node_column_IDs( fd, filename, gmshFormat, nodes )
 
     ! Skip newline character when in binary mode
     if( gmshFormat == binaryFormat ) then
-      read(fd), newlineChar
+      read(fd) newlineChar
       call ascii_formatting(fd, filename, "read")
     end if
 

femtools/embed_python.c

@@ -33,6 +33,7 @@ USA
 #include "Python.h"
 #endif
 #ifdef HAVE_NUMPY
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
 #include "numpy/arrayobject.h"
 #endif
 
@@ -526,18 +527,18 @@ void set_tensor_field_from_python(char *function, int *function_len, int *dim,
     }
 
     pArray = (PyArrayObject *)
-      PyArray_ContiguousFromObject(pResult, PyArray_DOUBLE, 2, 2);
+      PyArray_ContiguousFromObject(pResult, NPY_DOUBLE, 2, 2);
 
     if (PyErr_Occurred()){
       PyErr_Print();
       *stat=1;
       return;
     }
 
-    if (pArray->dimensions[0] != result_dim[0] || pArray->dimensions[1] != result_dim[1])
+    if (PyArray_DIMS(pArray)[0] != result_dim[0] || PyArray_DIMS(pArray)[1] != result_dim[1])
     {
       fprintf(stderr, "Error: dimensions of array returned from python ([%d, %d]) do not match allocated dimensions of the tensor_field ([%d, %d])).\n", 
-             (int) pArray->dimensions[0], (int) pArray->dimensions[1], result_dim[0], result_dim[1]);
+             (int) PyArray_DIMS(pArray)[0], (int) PyArray_DIMS(pArray)[1], result_dim[0], result_dim[1]);
       *stat=1;
       return;
     }
@@ -547,7 +548,7 @@ void set_tensor_field_from_python(char *function, int *function_len, int *dim,
 
         // Note the transpose for fortran.
         double tmp;
-        tmp = *(double*)(pArray->data + ii * pArray->strides[0] + jj * pArray->strides[1]);
+        tmp = *(double*)(PyArray_DATA(pArray) + ii * PyArray_STRIDES(pArray)[0] + jj * PyArray_STRIDES(pArray)[1]);
         result[i*(result_dim[0] * result_dim[1]) + jj * result_dim[0] + ii] = tmp;
       }
     }

femtools/python_statec.c

@@ -514,7 +514,7 @@ void python_add_array_double_1d(double *arr, int *size, char *name){
   //  the actual data as a byte array(char*)
 
   // Set the array
-  PyObject *a = PyArray_SimpleNewFromData(1, (npy_intp[]){*size}, PyArray_DOUBLE, (char*)arr);
+  PyObject *a = PyArray_SimpleNewFromData(1, (npy_intp[]){*size}, NPY_DOUBLE, (char*)arr);
   PyDict_SetItemString(pDict,name,a);
   Py_DECREF(a);
 #endif
@@ -528,7 +528,7 @@ void python_add_array_double_2d(double *arr, int *sizex, int *sizey, char *name)
 
   // Set the array
   npy_intp dims[] = {*sizey,*sizex};
-  PyObject *a = PyArray_SimpleNewFromData(2, dims, PyArray_DOUBLE, (char*)arr);
+  PyObject *a = PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, (char*)arr);
   PyDict_SetItemString(pDict,name,a);
   char c[200];
   snprintf(c, 200, "%s = numpy.transpose(%s,(1,0))",name,name);
@@ -545,7 +545,7 @@ void python_add_array_double_3d(double *arr, int *sizex, int *sizey, int *sizez,
 
   // Set the array
   npy_intp dims[] = {*sizez,*sizey,*sizex};
-  PyObject *a = PyArray_SimpleNewFromData(3, dims, PyArray_DOUBLE, (char*)arr);
+  PyObject *a = PyArray_SimpleNewFromData(3, dims, NPY_DOUBLE, (char*)arr);
   PyDict_SetItemString(pDict,name,a);
   char c[200];
   snprintf(c, 200, "%s = numpy.transpose(%s,(2,1,0))",name,name);
@@ -561,7 +561,7 @@ void python_add_array_integer_1d(int *arr, int *size, char *name){
   PyObject *pDict = PyModule_GetDict(pMain);
 
   // Set the array
-  PyObject *a = PyArray_SimpleNewFromData(1, (npy_intp[]){*size}, PyArray_INT, (char*)arr);
+  PyObject *a = PyArray_SimpleNewFromData(1, (npy_intp[]){*size}, NPY_INT, (char*)arr);
   PyDict_SetItemString(pDict,name,a);
   Py_DECREF(a);
 #endif
@@ -575,7 +575,7 @@ void python_add_array_integer_2d(int *arr, int *sizex, int *sizey, char *name){
 
   // Set the array
   npy_intp dims[] = {*sizey,*sizex};
-  PyObject *a = PyArray_SimpleNewFromData(2, dims, PyArray_INT, (char*)arr);
+  PyObject *a = PyArray_SimpleNewFromData(2, dims, NPY_INT, (char*)arr);
   PyDict_SetItemString(pDict,name,a);
   char c[200];
   snprintf(c, 200, "%s = numpy.transpose(%s,(1,0))",name,name);
@@ -592,7 +592,7 @@ void python_add_array_integer_3d(int *arr, int *sizex, int *sizey, int *sizez, c
 
   // Set the array
   npy_intp dims[] = {*sizez,*sizey,*sizex};
-  PyObject *a = PyArray_SimpleNewFromData(3, dims, PyArray_INT, (char*)arr);
+  PyObject *a = PyArray_SimpleNewFromData(3, dims, NPY_INT, (char*)arr);
   PyDict_SetItemString(pDict,name,a);
   char c[200];
   snprintf(c, 200, "%s = numpy.transpose(%s,(2,1,0))",name,name);

femtools/tests/test_main.cpp

@@ -59,7 +59,11 @@ int main(int argc, char **argv)
 #ifdef HAVE_MPI
   MPI_Init(&argc, &argv);
   // Undo some MPI init shenanigans
-  chdir(getenv("PWD"));
+  int ierr = chdir(getenv("PWD"));
+  if (ierr == -1) {
+        std::cerr << "Unable to switch to directory " << getenv("PWD");
+        abort();
+  }
 #endif
 #ifdef HAVE_PETSC  
   PetscInitialize(&argc, &argv, NULL, PETSC_NULL);