rh_source Subroutine

  subroutine rh_source(qt,lQgrid,x)
    use mod_global_parameters
    

    double precision, intent(in)    :: qt
    double precision, intent(inout) :: lQgrid
    double precision, intent(in)    :: x(1:ndim)

    double precision                :: tshift,tr,stt1,stt2,stt3
    double precision                :: tl1,tl2,tl3,tl4,tl5,tl6
    double precision                :: pulse1,pulse2,pulse3,pulse4,pulse5,pulse6
    integer                         :: iip1,iip2,iip3,i

    tshift = qt - trelax
    
    if(qt .lt. trelax) then
      tr=zero
      lQgrid=zero
      return
    elseif(qt .lt. trelax+tramp) then
      tr=(qt-trelax)/tramp
    else
      tr=one
    endif

    do i = 4, ntimes
      if(qt .lt. trelax + timearray(i)) then
          tl1 = dexp(-(tshift - timearray(i-3))**2/(0.5d0*dtimearr(i-3))**2)
          tl2 = dexp(-(tshift - timearray(i-2))**2/(0.5d0*dtimearr(i-2))**2)
          tl3 = dexp(-(tshift - timearray(i-1))**2/(0.5d0*dtimearr(i-1))**2)
          tl4 = dexp(-(tshift - timearray(i)  )**2/(0.5d0*dtimearr(i)  )**2)
          tl5 = dexp(-(tshift - timearray(i+1))**2/(0.5d0*dtimearr(i+1))**2)
          tl6 = dexp(-(tshift - timearray(i+2))**2/(0.5d0*dtimearr(i+2))**2)
          ! remapped on the unit cube
          stt1=(xprobmax1-xprobmin1)/(dble(vlim1)-1.d0)
          stt2=(xprobmax2-xprobmin2)/(dble(vlim2)-1.d0)
          stt3=(xprobmax3-xprobmin3)/(dble(vlim3)-1.d0)
          iip1=floor((x(1)-xprobmin1)/stt1)+1
          iip2=floor((x(2)-xprobmin2)/stt2)+1
          iip3=floor((x(3)-xprobmin3)/stt3)+1
          iip1=max(1,iip1)
          iip2=max(1,iip2)
          iip3=max(1,iip3)
          pulse1=va1(iip1+(i-4)*vlim1)*va2(iip2+(i-4)*vlim2)*va3(iip3+(i-4)*vlim3)
          pulse2=va1(iip1+(i-3)*vlim1)*va2(iip2+(i-3)*vlim2)*va3(iip3+(i-3)*vlim3)
          pulse3=va1(iip1+(i-2)*vlim1)*va2(iip2+(i-2)*vlim2)*va3(iip3+(i-2)*vlim3)
          pulse4=va1(iip1+(i-1)*vlim1)*va2(iip2+(i-1)*vlim2)*va3(iip3+(i-1)*vlim3)
          pulse5=va1(iip1+(i)*vlim1)*va2(iip2+(i)*vlim2)*va3(iip3+(i)*vlim3)
          pulse6=va1(iip1+(i+1)*vlim1)*va2(iip2+(i+1)*vlim2)*va3(iip3+(i+1)*vlim3)
          lQgrid = tl1 * pulse1 +  tl2 *pulse2  &
                        +tl3 * pulse3 +  tl4 *pulse4  &
                        +tl5 * pulse5 +  tl6 *pulse6
          exit
      endif 
    enddo 

    lQgrid=tr*lQgrid
  
  end subroutine rh_source