Added flag for handling surface versus surface0, fixed typos, changed code to...

Added flag for handling surface versus surface0, fixed typos, changed code to resize surface arrays rather than deallocate/reallocate entire surface array structure

src/remove_surfaces.f90

@@ -16,9 +16,10 @@
 !------------------------------------------------------------------------------|
 !------------------------------------------------------------------------------|
 
-subroutine remove_surfaces (params,surface)
+subroutine remove_surfaces (params,surface,surf0flag,threadinfo)
 
 use definitions
+use threads
 
 !------------------------------------------------------------------------------|
 !(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
@@ -36,13 +37,15 @@ implicit none
 include 'mpif.h'
 
 type (parameters) params
-type (sheet) surface 
+type (sheet) :: surface(params%nsorig)
+integer :: surf0flag
+type (thread) threadinfo
 
 !------------------------------------------------------------------------------|
 !(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
 !------------------------------------------------------------------------------|
-type (sheet) surfacetemp
-integer :: i,j,iproc,nproc,ierr,err,is
+type (sheet),dimension(:),allocatable :: surfacetemp
+integer :: i,j,iproc,nproc,ierr,err,is,isnew
 
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
@@ -53,9 +56,21 @@ call mpi_comm_rank (mpi_comm_world,iproc,ierr)
 allocate (surfacetemp(params%nsorig),stat=err); if (err.ne.0) call stop_run ('Error alloc surfacetemp in remove_surfaces$')
 
 do is=1,params%nsorig
+
+
+  !write (*,*) 'In first loop with is=',is
+  !write (*,*) 'surface(is)%nsurface: ',surface(is)%nsurface
+  !write (*,*) 'surface(is)%nt: ',surface(is)%nt
+  !write (*,*) 'size(surface): ',size(surface)
+
+
   ! Store number of surface particles and triangles in original surfaces
   surfacetemp(is)%nsurface=surface(is)%nsurface
-  surfacetemp(is)%nt=surface(is)%nt
+  if (surf0flag == 0) surfacetemp(is)%nt=surface(is)%nt
+
+
+  !write (*,*) 'Stored surface particles/triangles'
+
 
   ! Allocate memory for temporary surface fields
   allocate (surfacetemp(is)%x(surfacetemp(is)%nsurface),stat=threadinfo%err)
@@ -74,8 +89,14 @@ do is=1,params%nsorig
   if (params%debug.gt.1) call heap (threadinfo,'surfacetemp(is)%r','remove_surfaces',size(surfacetemp(is)%r),'dp',+1)
   allocate (surfacetemp(is)%s(surfacetemp(is)%nsurface),stat=threadinfo%err)
   if (params%debug.gt.1) call heap (threadinfo,'surfacetemp(is)%s','remove_surfaces',size(surfacetemp(is)%s),'dp',+1)
-  allocate (surfacetemp(is)%icon(3,surfacetemp(is)%nt),stat=threadinfo%err)
-  if (params%debug.gt.1) call heap (threadinfo,'surfacetemp(is)%icon','remove_surfaces',size(surfacetemp(is)%icon),'int',+1)
+  if (surf0flag == 0) then
+    allocate (surfacetemp(is)%icon(3,surfacetemp(is)%nt),stat=threadinfo%err)
+    if (params%debug.gt.1) call heap (threadinfo,'surfacetemp(is)%icon','remove_surfaces',size(surfacetemp(is)%icon),'int',+1)
+  endif
+
+
+  !write (*,*) 'After big allocate'
+
 
   ! Copy data from original surface array to surfacetemp
   surfacetemp(is)%x=surface(is)%x
@@ -86,7 +107,7 @@ do is=1,params%nsorig
   surfacetemp(is)%zn=surface(is)%zn
   surfacetemp(is)%r=surface(is)%r
   surfacetemp(is)%s=surface(is)%s
-  surfacetemp(is)%icon(1:3,1:surfacetemp(is)%nt)=surface(is)%icon(1:3,1:surface(is)%nt)
+  if (surf0flag == 0) surfacetemp(is)%icon(1:3,1:surfacetemp(is)%nt)=surface(is)%icon(1:3,1:surface(is)%nt)
 
   ! Remove original surface fields
   if (params%debug.gt.1) call heap (threadinfo,'surface(is)%x','remove_surfaces',size(surface(is)%x),'dp',-1)
@@ -105,22 +126,24 @@ do is=1,params%nsorig
   deallocate (surface(is)%r)
   if (params%debug.gt.1) call heap (threadinfo,'surface(is)%s','remove_surfaces',size(surface(is)%s),'dp',-1)
   deallocate (surface(is)%s)
-  if (params%debug.gt.1) call heap (threadinfo,'surface(is)%icon','remove_surfaces',size(surface(is)%icon),'int',-1)
-  deallocate (surface(is)%icon)
+  if (surf0flag == 0) then
+    if (params%debug.gt.1) call heap (threadinfo,'surface(is)%icon','remove_surfaces',size(surface(is)%icon),'int',-1)
+    deallocate (surface(is)%icon)
+  endif
 enddo
 
-! Deallocate surface structure
-deallocate (surface)
+! Deallocate surface structure - Commented out for now, probably not possible
+!deallocate (surface)
 
-! Allocate new surface structure
-allocate (surface(params%ns),stat=err); if (err.ne.0) call stop_run ('Error alloc surface in remove_surfaces$')
+! Allocate new surface structure - Commented due to lack of need to reallocate
+!allocate (surface(params%ns),stat=err); if (err.ne.0) call stop_run ('Error alloc surface in remove_surfaces$')
 
 do is=1,params%nsorig
   if (params%surfremindex(is) > 0) then
     isnew=params%surfremindex(is)
     ! Store number of surface particles and triangles in new surfaces
     surface(isnew)%nsurface=surfacetemp(is)%nsurface
-    surface(isnew)%nt=surfacetemp(is)%nt
+    if (surf0flag == 0) surface(isnew)%nt=surfacetemp(is)%nt
 
     ! Allocate memory for new surface fields
     allocate (surface(isnew)%x(surface(isnew)%nsurface),stat=threadinfo%err)
@@ -139,8 +162,10 @@ do is=1,params%nsorig
     if (params%debug.gt.1) call heap (threadinfo,'surface(isnew)%r','remove_surfaces',size(surface(isnew)%r),'dp',+1)
     allocate (surface(isnew)%s(surface(isnew)%nsurface),stat=threadinfo%err)
     if (params%debug.gt.1) call heap (threadinfo,'surface(isnew)%s','remove_surfaces',size(surface(isnew)%s),'dp',+1)
-    allocate (surface(isnew)%icon(3,surface(isnew)%nt),stat=threadinfo%err)
-    if (params%debug.gt.1) call heap (threadinfo,'surface(isnew)%icon','remove_surfaces',size(surface(isnew)%icon),'int',+1)
+    if (surf0flag == 0) then
+      allocate (surface(isnew)%icon(3,surface(isnew)%nt),stat=threadinfo%err)
+      if (params%debug.gt.1) call heap (threadinfo,'surface(isnew)%icon','remove_surfaces',size(surface(isnew)%icon),'int',+1)
+    endif
 
     ! Copy data from surfacetemp arrays to new surface arrays
     surface(isnew)%x=surfacetemp(is)%x
@@ -151,7 +176,7 @@ do is=1,params%nsorig
     surface(isnew)%zn=surfacetemp(is)%zn
     surface(isnew)%r=surfacetemp(is)%r
     surface(isnew)%s=surfacetemp(is)%s
-    surface(isnew)%icon(1:3,1:surface(isnew)%nt)=surfacetemp(is)%icon(1:3,1:surfacetemp(is)%nt)
+    if (surf0flag == 0) surface(isnew)%icon(1:3,1:surface(isnew)%nt)=surfacetemp(is)%icon(1:3,1:surfacetemp(is)%nt)
   endif
 
   ! Remove temporary surface fields
@@ -171,13 +196,43 @@ do is=1,params%nsorig
   deallocate (surfacetemp(is)%r)
   if (params%debug.gt.1) call heap (threadinfo,'surfacetemp(is)%s','remove_surfaces',size(surfacetemp(is)%s),'dp',-1)
   deallocate (surfacetemp(is)%s)
-  if (params%debug.gt.1) call heap (threadinfo,'surfacetemp(is)%icon','remove_surfaces',size(surfacetemp(is)%icon),'int',-1)
-  deallocate (surfacetemp(is)%icon)
+  if (surf0flag == 0) then
+    if (params%debug.gt.1) call heap (threadinfo,'surfacetemp(is)%icon','remove_surfaces',size(surfacetemp(is)%icon),'int',-1)
+    deallocate (surfacetemp(is)%icon)
+  endif
 enddo
 
 ! Deallocate temporary surface structure
 deallocate (surfacetemp)
 
+! Reduce size of unused surfaces to 1 particle and 1 triangle
+do is=params%ns+1,params%nsorig
+  ! Store number of surface particles and triangles in new surfaces
+  surface(is)%nsurface=-1
+  if (surf0flag == 0) surface(is)%nt=-1
+enddo
+
+!  ! Allocate memory for new surface fields
+!  allocate (surface(is)%x(surface(is)%nsurface),stat=threadinfo%err)
+!  if (params%debug.gt.1) call heap (threadinfo,'surface(is)%x','remove_surfaces',size(surface(is)%x),'dp',+1)
+!  allocate (surface(is)%y(surface(is)%nsurface),stat=threadinfo%err)
+!  if (params%debug.gt.1) call heap (threadinfo,'surface(is)%y','remove_surfaces',size(surface(is)%y),'dp',+1)
+!  allocate (surface(is)%z(surface(is)%nsurface),stat=threadinfo%err)
+!  if (params%debug.gt.1) call heap (threadinfo,'surface(is)%z','remove_surfaces',size(surface(is)%z),'dp',+1)
+!  allocate (surface(is)%xn(surface(is)%nsurface),stat=threadinfo%err)
+!  if (params%debug.gt.1) call heap (threadinfo,'surface(is)%xn','remove_surfaces',size(surface(is)%xn),'dp',+1)
+!  allocate (surface(is)%yn(surface(is)%nsurface),stat=threadinfo%err)
+!  if (params%debug.gt.1) call heap (threadinfo,'surface(is)%yn','remove_surfaces',size(surface(is)%yn),'dp',+1)
+!  allocate (surface(is)%zn(surface(is)%nsurface),stat=threadinfo%err)
+!  if (params%debug.gt.1) call heap (threadinfo,'surface(is)%zn','remove_surfaces',size(surface(is)%zn),'dp',+1)
+!  allocate (surface(is)%r(surface(is)%nsurface),stat=threadinfo%err)
+!  if (params%debug.gt.1) call heap (threadinfo,'surface(is)%r','remove_surfaces',size(surface(is)%r),'dp',+1)
+!  allocate (surface(is)%s(surface(is)%nsurface),stat=threadinfo%err)
+!  if (params%debug.gt.1) call heap (threadinfo,'surface(is)%s','remove_surfaces',size(surface(is)%s),'dp',+1)
+!  allocate (surface(is)%icon(3,surface(is)%nt),stat=threadinfo%err)
+!  if (params%debug.gt.1) call heap (threadinfo,'surface(is)%icon','remove_surfaces',size(surface(is)%icon),'int',+1)
+!enddo
+
 end subroutine remove_surfaces
 
 !-------------------------------------------------------------------------------