subroutine mygrad_spec(n_gauss, gradient, residual, params, dim_v)
implicit none
integer, intent(in) :: n_gauss, dim_v
real(xp), intent(in), dimension(3*n_gauss) :: params
real(xp), intent(in), dimension(:), allocatable :: residual
real(xp), intent(inout), dimension(3*n_gauss) :: gradient
integer :: i, k
real(xp) :: g
real(xp), dimension(:,:), allocatable :: dF_over_dB
allocate(dF_over_dB(3*n_gauss, dim_v))
g = 0._xp
dF_over_dB = 0._xp
gradient = 0._xp
do i=1, n_gauss
do k=1, dim_v
dF_over_dB(1+(3*(i-1)),k) = dF_over_dB(1+(3*(i-1)),k) +&
exp( ( -(real(k,xp) - params(2+(3*(i-1))))**2._xp) / (2._xp * params(3+(3*(i-1)))**2._xp))
dF_over_dB(2+(3*(i-1)),k) = dF_over_dB(2+(3*(i-1)),k) +&
params(1+(3*(i-1))) * ( real(k,xp) - params(2+(3*(i-1))) ) / (params(3+(3*(i-1)))**2._xp) *&
exp( ( -(real(k,xp) - params(2+(3*(i-1))))**2._xp) / (2._xp * params(3+(3*(i-1)))**2._xp))
dF_over_dB(3+(3*(i-1)),k) = dF_over_dB(3+(3*(i-1)),k) +&
params(1+(3*(i-1))) * ( real(k,xp) - params(2+(3*(i-1))) )**2._xp / (params(3+(3*(i-1)))**3._xp) *&
exp( ( -(real(k,xp) - params(2+(3*(i-1))))**2._xp) / (2._xp * params(3+(3*(i-1)))**2._xp))
enddo
enddo
do i=1, dim_v
do k=1, 3*n_gauss
gradient(k) = gradient(k) + dF_over_dB(k,i) * residual(i)
end do
end do
end subroutine mygrad_spec