Fix routines to generate full sphere radial grids

Author: hiromatsui

src/Fortran_libraries/SERIAL_src/SPH_SPECTR_src/half_chebyshev_radial_grid.f90

@@ -44,14 +44,15 @@ subroutine half_chebyshev_distance_shell                          &
       kst = 1
       ked = nlayer_CMB
       do k = kst, ked
-        r_grid(k) = r_CMB * sin(half * pi * dble(k)/dble(nlayer_CMB))
+        r_grid(k)                                                       &
+     &     = r_CMB * cos(half * pi * dble(ked - k)/dble(ked))
       end do
 !
-      kst = nlayer_CMB + 1
+      kst = nlayer_CMB
       ked = min(num_layer, nlayer_CMB + nlayer_CMB/2)
-      do k = kst, ked
-        r_grid(k) = r_CMB + r_CMB * (one - sin(half*pi                  &
-     &                       * dble(k-nlayer_CMB)/dble(nlayer_CMB)) )
+      do k = kst+1, ked
+        r_grid(k)                                                       &
+     &     = r_CMB * (two - cos(half*pi * dble(k-kst)/dble(kst)))
       end do
       dr = r_grid(ked) - r_grid(ked-1)
 !
@@ -93,15 +94,13 @@ subroutine count_half_chebyshev_external(nri, r_CMB, r_max,       &
         r = r_CMB + r_CMB * (one - cos(half*pi*dble(k)/dble(nri)) )
         dr =  r_CMB * (-cos(half*pi*dble(k  )/dble(nri))                &
      &                 +cos(half*pi*dble(k-1)/dble(nri)) )
-!        write(*,*) k, r, dr
       end do
 !
 !
       do
         if(r .ge. r_max) exit
         k = k + 1
         r = r + dr
-!        write(*,*) k, r, dr
       end do
 !
       if(k .le. 1) then
@@ -110,7 +109,7 @@ subroutine count_half_chebyshev_external(nri, r_CMB, r_max,       &
         ngrid_ext = k
       end if
 !
-      nlayer_ICB = 0
+      nlayer_ICB = 1
       nlayer_CMB = nri
       ntot_shell = nlayer_CMB + ngrid_ext
 !

src/Fortran_libraries/UTILS_src/MERGE/parallel_assemble_sph.f90

@@ -89,12 +89,11 @@ subroutine set_mode_table_4_assemble(org_sph, new_sph, j_table)
       type(sph_grids), intent(in) :: new_sph
       type(rj_assemble_tbl), intent(inout) :: j_table
 !
-      integer(kind = kint) :: j_org, j_gl, j_new
+      integer(kind = kint) :: j_org, j_new
       integer(kind = 4) :: l_gl, m_gl
 !
 !
       do j_org = 1, org_sph%sph_rj%nidx_rj(2)
-        j_gl =     org_sph%sph_rj%idx_gl_1d_rj_j(j_org,1)
         l_gl = int(org_sph%sph_rj%idx_gl_1d_rj_j(j_org,2))
         m_gl = int(org_sph%sph_rj%idx_gl_1d_rj_j(j_org,3))
         j_new = find_local_sph_address(new_sph%sph_rj, l_gl, m_gl)