array out of bounds in ridge_ana.F90 #64
Description
I think there is a bug in rigde_ana.F90. This code:
DO np = 1, 6
DO j=1-nhalo+1,ncube+nhalo-1,bloc
DO i=1-nhalo+1,ncube+nhalo-1,bloc
aa = terr_dev_halo(i:i+bloc,j:j+bloc,np)
ijaa = MAXLOC( aa )
im = ijaa(1)-1 + i
jm = ijaa(2)-1 + j
if ( terr_dev_halo(im,jm,np) >= thsh ) terr_max_halo(im,jm,np)= terr_dev_halo(im,jm,np)
END DO
END DO
can make terr_dev_halo go out of bounds when bloc>1 (array bounds for terr_dev_halo is DIMENSION(1-nhalo:ncube+nhalo, 1-nhalo:ncube+nhalo, 6).
My temporary fix is:
DO np = 1, 6
DO j=1-nhalo+1,ncube+nhalo-1,bloc
DO i=1-nhalo+1,ncube+nhalo-1,bloc
if (i-bloc<1-nhalo.or.j-bloc<1-nhalo.or.i+bloc>ncube+nhalo.or.j+bloc>ncube+nhalo) then
!out of bounds - skip
else
aa = terr_dev_halo(i:i+bloc,j:j+bloc,np)
ijaa = MAXLOC( aa )
im = ijaa(1)-1 + i
jm = ijaa(2)-1 + j
if ( terr_dev_halo(im,jm,np) >= thsh ) terr_max_halo(im,jm,np)= terr_dev_halo(im,jm,np)
end if
END DO
END DO
Discovered problem on Casper with:
./cube_to_target --jmax_segments 100000 --smoothing_over_ocean --grid_descriptor_file /glade/campaign/cesm/cesmdata/inputdata/share/meshes/ne16np4_ESMFmesh_cdf5_c20211018.nc --intermediate_cs_name /glade/campaign/cgd/amp/pel/topo/cubedata/mars-ncube3000.nc --output_grid ne16np4 --smoothing_scale 1.0 --name_email_of_creator Peter Hjort Lauritzen, pel@ucar.edu --output_data_directory output/