scalsource Subroutine

subroutine scalsource

  use modglobal,  only : pi,nsv,ib,ie,jb,je,kb,ke,ih,jh,kh,ihc,jhc,khc,xf,zf,xh,zh,dx,dy,imax,itot,jmax,jtot,lchem,&
                         xS,yS,zS,xSb,ySb,zSb,xSe,ySe,zSe,SS,sigS,lscasrc,lscasrcl,lscasrcr,libm,dxfi,dzfi,nscasrc,scasrcp,nscasrcl,scasrcl
  use modfields,  only : svp,svpp
  use modmpi,     only : myid,myidx,myidy,mpierr,MY_REAL,comm3d,MPI_SUM

  implicit none
  integer :: i,j,k,n,ns
  real :: dxi, dyi
  real :: ra2 = 0.
  real :: scalsum = 0.
  real :: scalsumt = 0.
  real :: px = 0., py = 0., pz = 0.0, vx = 0., vy = 0., vz = 0., lsx = 0., lsy = 0., lsz = 0., dot_projection = 0.

  dxi = 1./dx
  dyi = 1./dy
 
  !  Input passive scalar point sources
  if (lscasrc .AND. nsv.gt.0) then 

    do n=1,nsv 
      do ns=1,nscasrc
        xS = scasrcp(ns,1,n)
        yS = scasrcp(ns,2,n)
        zS = scasrcp(ns,3,n)
        SS = scasrcp(ns,4,n)
        sigS = scasrcp(ns,5,n)
        do k=kb,ke
          do j=jb,je
            do i=ib,ie
                
              ra2 = ((i+myidx*imax-0.5)*dx-xS)**2 + ((j+myidy*jmax-0.5)*dy-yS)**2 + (zf(k)-zS)**2

              if (ra2 .LE. 9*sigS**2) then
                svp(i,j,k,n) = svp(i,j,k,n) + dxi * dyi * dzfi(k) * SS*exp(-ra2/(2*sigS**2))
              end if

            end do
          end do
        end do
      end do
    end do

  end if !lscasrc

  !  Input passive scalar line sources
  if (lscasrcl .AND. nsv.gt.0) then 
    
    do n=1,nsv
      do ns=1,nscasrcl
        xSb = scasrcl(ns,1,n)
        ySb = scasrcl(ns,2,n)
        zSb = scasrcl(ns,3,n)
        xSe = scasrcl(ns,4,n)
        ySe = scasrcl(ns,5,n)
        zSe = scasrcl(ns,6,n)
        SS = scasrcl(ns,7,n)
        sigS = scasrcl(ns,8,n)
        
        lsx = xSe-xSb
        lsy = ySe-ySb
        lsz = zSe-zSb
        
        do k=kb,ke
          do j=jb,je
            do i=ib,ie
              
              px = (i+myidx*imax-0.5)*dx
              py = (j+myidy*jmax-0.5)*dy
              pz = zf(k)

              vx = px-xSb
              vy = py-ySb
              vz = pz-zSb

              dot_projection = (vx * lsx + vy * lsy + vz * lsz) / (lsx * lsx + lsy * lsy + lsz * lsz)

              if (dot_projection < 0.0) then
                ra2 = (px-xSb)**2 + (py-ySb)**2 + (pz-zSb)**2
              elseif (dot_projection > 1.0) then
                ra2 = (px-xSe)**2 + (py-ySe)**2 + (pz-zSe)**2
              else
                ra2 = (px-(xSb+dot_projection*lsx))**2 + (py-(ySb+dot_projection*lsy))**2 + (pz-(zSb+dot_projection*lsz))**2
              end if

              if (ra2 .LE. 9*sigS**2) then
                svp(i,j,k,n) = svp(i,j,k,n) + dxi * dyi * dzfi(k) * &
                              sqrt(2.0*pi)*SS*sigS*exp(-ra2/(2*sigS**2)) * erf(sqrt((9*sigS**2-ra2)/(2*sigS**2)))
              end if

            end do
          end do
        end do
      end do
    end do

  end if !lscasrcl

end subroutine scalsource