diff options
Diffstat (limited to 'src/EHFinder_Integrate2.F90')
| -rw-r--r-- | src/EHFinder_Integrate2.F90 | 69 |
1 files changed, 53 insertions, 16 deletions
diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90 index cc77255..81745ea 100644 --- a/src/EHFinder_Integrate2.F90 +++ b/src/EHFinder_Integrate2.F90 @@ -35,6 +35,11 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) CCTK_VARIABLE_REAL, & CCTK_VARIABLE_REAL /) +! If finding of surface failed don't try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + call CCTK_INFO ( 'Finding surface element' ) call CCTK_GroupbboxGN ( ierr, cctkGH, 4, bbox, "ehfinder::surface_arrays" ) if ( ierr .lt. 0 ) then @@ -65,12 +70,12 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get local size for surface arrays" ) end if - + dtheta = ctheta(2,1) - ctheta(1,1) dphi = cphi(1,2) - cphi(1,1) dthetainv = one / dtheta dphiinv = one / dphi - + call CCTK_InterpHandle ( interp_handle, "Lagrange polynomial interpolation" ) if ( interp_handle .lt. 0 ) then call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" ) @@ -80,7 +85,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) if ( table_handle .lt. 0 ) then call CCTK_WARN( 0, "Cannot create parameter table for interpolator" ) end if - + call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" ) if ( coord_system_handle .lt. 0) then call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system. Forgot to activate an implementation providing coordinates ??" ) @@ -207,6 +212,11 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS) character(len=30) :: info_message CCTK_INT, dimension(1) :: lbnd, ubnd, lsh +! If finding of surface failed don't try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + call CCTK_INFO ( 'Integrating area' ) sn = surface_counter - integrate_counter @@ -264,6 +274,11 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS) CCTK_INT, dimension(1) :: lbnd, ubnd, lsh CCTK_REAL :: centroid_x, centroid_y, centroid_z +! If finding of surface failed don't try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + call CCTK_INFO ( 'Integrating centroid' ) sn = surface_counter - integrate_counter @@ -367,6 +382,11 @@ subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS) CCTK_INT :: ithl, ithr, jphl, jphr, itha, jpha CCTK_REAL :: eq_circ_loc, eq_circ, pol_circ_loc, pol_circ +! If finding of surface failed don't try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + sn = surface_counter - integrate_counter ! Find out the lower bounds of the distributed integration grid arrays. @@ -581,13 +601,18 @@ subroutine EHFinder_CopyArea(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - call CCTK_INFO('Copying area data') - ! print*,Xeh_area0 ! print*,Xeh_area20 - eh_area = eh_area2 - +! If finding of surface failed set area to zero. + print*,'Debug2 : ',find_surface_status + if ( find_surface_status .lt. 0 ) then + eh_area = zero + return + else + call CCTK_INFO('Copying area data') + eh_area = eh_area2 + endif end subroutine EHFinder_CopyArea @@ -601,14 +626,21 @@ subroutine EHFinder_CopyCentroid(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - call CCTK_INFO('Copying centroid data') - ! print*,Xeh_centroid0 ! print*,Xeh_centroid20 - eh_centroid_x = eh_centroid2_x - eh_centroid_y = eh_centroid2_y - eh_centroid_z = eh_centroid2_z +! If finding of surface failed set centroids to zero. + if ( find_surface_status .lt. 0 ) then + eh_centroid_x = zero + eh_centroid_y = zero + eh_centroid_z = zero + return + else + call CCTK_INFO('Copying centroid data') + eh_centroid_x = eh_centroid2_x + eh_centroid_y = eh_centroid2_y + eh_centroid_z = eh_centroid2_z + end if end subroutine EHFinder_CopyCentroid @@ -623,9 +655,14 @@ subroutine EHFinder_CopyCircumference(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - call CCTK_INFO('Copying circumference data') - - eh_circ_eq = eh_circ_eq2 - eh_circ_pol = eh_circ_pol2 + if ( find_surface_status .lt. 0 ) then + eh_circ_eq = zero + eh_circ_pol = zero + return + else + call CCTK_INFO('Copying circumference data') + eh_circ_eq = eh_circ_eq2 + eh_circ_pol = eh_circ_pol2 + end if end subroutine EHFinder_CopyCircumference |