function a_Interpolate_array00_wa(wa_data,alon,alat)
!
! 緯度 alon, 経度 alat における関数値を
! その球面調和変換係数 wa_data から補間計算する
!
real(8), intent(IN) :: wa_data(:,:) ! スペクトルデータ
real(8), intent(IN) :: alon ! 補間する位置(経度)
real(8), intent(IN) :: alat ! 補間する位置(緯度)
real(8) :: a_Interpolate_array00_wa(size(wa_data,2)) ! 補間した値
real(8) :: mu
real(8), dimension(size(wa_data,2)) :: y0, y1, y2, AnmPnm
integer :: k,m
mu = sin(alat)
a_Interpolate_array00_wa = 0.0D0
!---- Σa_n^0 P_n^0 の計算
y2 = 0.0D0
y1 = 0.0D0
do k=nm,1,-1
y0 = alpha(k,0,mu) * y1 + beta(k+1,0)*y2 + wa_data(l_nm(k,0),:)
y2 = y1
y1 = y0
enddo
a_Interpolate_array00_wa = ( beta(1,0) * y2 + mu*sqrt(3.0D0) * y1 + wa_data(l_nm(0,0),:) ) * Pmm(0,mu)
!---- cos mλ Σa_n^m P_n^m の計算
do m=1,nm
y2 = 0.0D0
y1 = 0.0D0
do k=nm,m+1,-1
y0 = alpha(k,m,mu) * y1 + beta(k+1,m) * y2 + wa_data(l_nm(k,m),:)
y2 = y1
y1 = y0
enddo
AnmPnm =(wa_data(l_nm(m,m),:) + beta(m+1,m)*y2 + mu*sqrt(2.0D0*m+3)*y1 ) * Pmm(m,mu)
a_Interpolate_array00_wa = a_Interpolate_array00_wa + AnmPnm * 2.0D0 * cos(m*alon)
end do
!---- sin λ Σa_n^m P_n^m の計算n
do m=1,nm
y2 = 0.0D0
y1 = 0.0D0
do k=nm,m+1,-1
y0 = alpha(k,m,mu) * y1 + beta(k+1,m)*y2 + wa_data(l_nm(k,-m),:)
y2 = y1
y1 = y0
enddo
AnmPnm =(wa_data(l_nm(m,-m),:) + beta(m+1,m)*y2 + mu*sqrt(2.0D0*m+3)*y1 ) * Pmm(m,mu)
a_Interpolate_array00_wa = a_Interpolate_array00_wa - AnmPnm * 2.0D0 * sin(m*alon)
end do
end function a_Interpolate_array00_wa