!> \file scalsource.f90

!> Input point, line or planar scalars
!  This file is part of uDALES.
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  Copyright 2006-2021 the uDALES Team.
!

!! ############  The below part is the scalsource script used by Tom Grylls. Now commented by DMajumdar.
!! If wish to use the below, do not forget to also uncomment Namelist reading of variables xS, yS, zS, SS and sigS
!! and broadcasting them in MPI in the modstartup.f90 file.  ############

! subroutine createscals

!   use modglobal,  only : nsv,ib,ie,jb,je,kb,ke,ih,jh,kh,lscasrcr,cexpnr,ifinput,imax,jmax,jtot,&
!                          libm,nblocks,xSa,ySa,zSa,lscasrc,xS,yS,zS
!   use modfields,  only : scar,scarl
!   use modmpi,     only : myid,MY_REAL,comm3d,mpierr
!   use initfac,    only : block
!   implicit none
!   integer :: i,j,k,n,m,p,il,iu,jl,ju,ku

!   ! if (lscasrc) then

!     ! allocate(xSa(1:nsv))
!     ! allocate(ySa(1:nsv))
!     ! allocate(zSa(1:nsv))

!     ! hard code point source locations
!     ! xSa(1:nsv) = (/ 453., 463., 523., 539., 407. /)
!     ! ySa(1:nsv) = (/ 375., 371., 327., 341., 325. /)
!     ! zSa(1:nsv) = (/ 4.   , 4.   , 4. , 4.  , 4.  /)

!     ! all point sources from position defined in namoptions
!     ! xSa = xS; ySa = yS; zSa = zS

!   ! end if

!   if (lscasrcr .AND. nsv.gt.0) then

!     ! read 2-D field of point sources for vehicular emissions network      
!     open (ifinput, file='scals.inp.'//cexpnr)

!     do j=jb,jtot
!       read(ifinput, *) scar(:,j)
!     end do

!     scarl(ib:ie,jb:je) = scar(ib:ie,jb+myid*jmax:je+myid*jmax)

!     ! Set scarl to 0 if set inside an obstacle
!     if (libm) then
!       do n=1,nblocks
!         il = block(n,1)
!         iu = block(n,2)
!         jl = block(n,3)-myid*jmax
!         ju = block(n,4)-myid*jmax
!         ku = block(n,6)
!         if (ju < jb .or. jl > je) then
!           cycle
!         else
!           if (ju > je) ju=je
!           if (jl < jb) jl=jb
!           ! addition due to all blocks at lowest layer
!           if (ku>0) then
!             scarl(il:iu,jl:ju) = 0.
!           end if
!         end if
!       end do
!     end if

!   end if !lscasrcr

! end subroutine

! subroutine scalsource

!   use modglobal,  only : nsv,ib,ie,jb,je,kb,ke,ih,jh,kh,ihc,jhc,khc,xf,zf,xh,zh,dy,jmax,lchem,&
!                          xS,yS,zS,SS,sigS,lscasrc,lscasrcl,lscasrcr,libm,dxfi,dzfi,nblocks,xSa,ySa,zSa
!   use modfields,  only : svp,svpp,scar,scarl
!   use modmpi,     only : myid,mpierr,MY_REAL,comm3d,MPI_SUM
!   use initfac, only : block

!   implicit none
!   integer :: i,j,k,n,il,iu,jl,ju,kl,ku,ncan
!   real :: xL,zL
!   real :: dyi
!   real :: ra2 = 0.
!   real :: scalsum = 0.
!   real :: scalsumt = 0.
!   real :: Pi = 3.1415927

!   dyi = 1./dy

!   ! 2-D network of point sources at lowest level 
!   if (lscasrcr .AND. nsv.gt.0) then
!       if (lchem) then
!         svp(ib:ie,jb:je,kb+1,1) = svp(ib:ie,jb:je,kb+1,1) + 0.522*scarl 
!         svp(ib:ie,jb:je,kb+1,2) = svp(ib:ie,jb:je,kb+1,2) + 0.2*scarl
!         svp(ib:ie,jb:je,kb+1,4) = svp(ib:ie,jb:je,kb+1,4) + scarl
!       else
!         svp(ib:ie,jb:je,kb+1,1) = svp(ib:ie,jb:je,kb+1,1) + scarl
!       end if
!   end if !lscasrcr

!   scalsum = 0.
 
!   !  Input passive scalar point sources
!   if (lscasrc .AND. nsv.gt.0) then 

!     do n=1,nsv 
!     do k=kb,ke
!       do j=jb,je
!         do i=ib,ie
               
!             ra2 = (xf(i)-xS)**2 + ((j+myid*jmax-0.5)*dy-yS)**2 + (zf(k)-zS)**2

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

!           end if
!         end do
!       end do
!     end do
!     end do

!     ! Set svpp to 0 when set inside an obstacle
!     if (libm) then
!       do n=1,nblocks
!         il = block(n,1)
!         iu = block(n,2)
!         kl = block(n,5)
!         ku = block(n,6)
!         jl = block(n,3)-myid*jmax
!         ju = block(n,4)-myid*jmax
!         if (ju < jb .or. jl > je) then
!           cycle
!         else
!           if (ju > je) ju=je
!           if (jl < jb) jl=jb
!           svp(il:iu,jl:ju,kl:ku,:) = 0.
!         end if
!       end do
!     end if

!     ! Normalise scalar field to 1/s
!     !call MPI_ALLREDUCE(scalsum,scalsumt,1,MY_REAL,MPI_SUM,comm3d,mpierr)
!     !svpp(:,:,:,1) = svpp(:,:,:,1) / scalsumt
!     !svp(:,:,:,1) = svp(:,:,:,1) + svpp(:,:,:,1)

!     scalsum = 0.

!   end if !lscasrc

!   ! Input passive scalar line sources
!   if (lscasrcl .AND. nsv.gt.0) then

!   ncan = count(block(:,6)>0) !tg3315 update due to block at lowest level

!   zL = zf(kb+1)

!     if (nblocks>0) then

!       do n = 1,ncan-1 ! not ncan and ncan+1 because we do not want release in first and last canyon for BCxs==2
!         if (n == ncan+1) then   ! Added to run for pollutant in first canyon

!           !xL = xh(block(1,1) - (block(2,1) - block(1,2)+1)/2)
!           xL = xh(block(1,1)) - 0.5*(xh(block(2,1)) - xh(block(1,2)+1))
! !          ra2 = (i - (block(1,1) - (block(2,1) - block(1,2))/2.0))**2 + (k)**2 !tg3315 commented for chem validation

!         else !cycle through all other canyons

!           !xL = xh(block(n,2) + (block(2,1) - block(1,2)+1)/2)
!           xL = xh(block(n,2)+1) + 0.5*(xh(block(2,1)) - xh(block(1,2)+1))
! !          ra2 =(xf(i) - xL)**2 + zf(k)**2

!        end if

!         do i=ib,ie
!           do k=kb,ke

! !           if (ra2 .LE. 12*sigS**2) then

!               !scalsum = scalsum + dxf(i) * jmax * dy * dzf(k) * (SS/2*Pi*sigS**2) * exp(-ra2/(2*sigS**2))

!               !tg3315 use this if we want to normalise th scalar conc. !sums values in building too...
! !              scalsum = scalsum + dy * (je - jb +1) * ( (SS/4.) * &
! !                        (erf((xh(i+1)-xL)/(sqrt(2.)*sigS)) - erf((xh(i)-xL)/(sqrt(2.)*sigS))) * &
! !                        (erf((zh(k+1)-zh(kb+1))/(sqrt(2.)*sigS)) - erf((zh(k)-zh(kb+1))/(sqrt(2.)*sigS))) + &
! !                        (SS/4.) * &
! !                        (erf((xh(i+1)-xL)/(sqrt(2.)*sigS)) - erf((xh(i)-xL)/(sqrt(2.)*sigS))) * &
! !                        (erf((zh(k+1)+zh(kb+1))/(sqrt(2.)*sigS)) - erf((zh(k)+zh(kb+1))/(sqrt(2.)*sigS))) )
! !                        * dxfi(i) * dzfi(k)

!               svp(i,jb:je,k,1) = svp(i,jb:je,k,1) + ( (SS/4.) * & ! SS in g/ms... no normalisation
!                         (erf((xh(i+1)-xL)/(sqrt(2.)*sigS)) - erf((xh(i)-xL)/(sqrt(2.)*sigS))) * &
!                         (erf((zh(k+1)-zL)/(sqrt(2.)*sigS)) - erf((zh(k)-zL)/(sqrt(2.)*sigS))) + &
!                         (SS/4.) * & ! reflection from ground...
!                         (erf((xh(i+1)-xL)/(sqrt(2.)*sigS)) - erf((xh(i)-xL)/(sqrt(2.)*sigS))) * &
!                         (erf((zh(k)-2*(zh(k)-zh(kb+1))-zL)/(sqrt(2.)*sigS)) - erf((zh(k+1)-2*(zh(k+1)-zh(kb+1))-zL)/(sqrt(2.)*sigS))) ) &
!                         * dxfi(i) * dzfi(k)

!               !svp(i,jb:je,k,3) = svp(i,jb:je,k,3) + ( (SS/4.) * & ! SS in g/ms... no normalisation
!               !          (erf((xh(i+1)-xL)/(sqrt(2.)*sigS)) - erf((xh(i)-xL)/(sqrt(2.)*sigS))) * &
!               !          (erf((zh(k+1)-zL)/(sqrt(2.)*sigS)) - erf((zh(k)-zL)/(sqrt(2.)*sigS))) + &
!               !          (SS/4.) * & ! reflection from ground...
!               !          (erf((xh(i+1)-xL)/(sqrt(2.)*sigS)) - erf((xh(i)-xL)/(sqrt(2.)*sigS))) * &
!               !          (erf((zh(k)-2*(zh(k)-zh(kb+1))-zL)/(sqrt(2.)*sigS)) - erf((zh(k+1)-2*(zh(k+1)-zh(kb+1))-zL)/(sqrt(2.)*sigS))) ) &
!               !          * dxfi(i) * dzfi(k)

! !            end if
!           end do
!         end do
!       end do
!     end if !nblocks

!     ! Set svpp to 0 when set inside an obstacle
!     if (libm) then
!       do n=1,nblocks
!         il = block(n,1)
!         iu = block(n,2)
!         kl = kb
!         ku = block(n,6)
!         jl = block(n,3)-myid*jmax
!         ju = block(n,4)-myid*jmax
!         if (ju < jb .or. jl > je) then
!           cycle
!         else
!           if (ju > je) ju=je
!           if (jl < jb) jl=jb
!           svp(il:iu,jl:ju,kl:ku,:) = 0.
!         end if
!       end do
!     end if !libm

!     ! Normalise scalar field to 1/s
! !    call MPI_ALLREDUCE(scalsum,scalsumt,1,MY_REAL,MPI_SUM,comm3d,mpierr)

! !    write(*,*), 'scalsum', scalsum

! !    if (lchem) then
!       !svpp(:,:,:,1) = svpp(:,:,:,1)
! !      svp(:,:,:,1) = svp(:,:,:,1) + svpp(:,:,:,1)
!       !svp(:,:,:,2) = svp(:,:,:,2) + 0.1518 * svpp(:,:,:,1)
! !    else
! !      svpp(:,:,:,1) = svpp(:,:,:,1)/ scalsumt !tg3315 not normalised 07/11/2017
! !      svp(:,:,:,1) = svp(:,:,:,1) + svpp(:,:,:,1)
!       !svp(:,:,:,2) = svp(:,:,:,2) + 0.1518 * svpp(:,:,:,1)
! !    end if

!     svpp = 0.
!     scalsum = 0.

!   end if !lscasrcl

! end subroutine scalsource

!! ############  ****  ############


!! ############  Updated scalar source modelling. DMajumdar  ############

module modscalsource

  implicit none

contains

subroutine createscals

  use modglobal,  only : nsv,cexpnr,ifinput,lscasrc,nscasrc,scasrcp,lscasrcl,nscasrcl,scasrcl
  use modmpi,     only : myid,MY_REAL,comm3d,mpierr
  implicit none
  integer :: n,m
  character :: cnsv
  character(80) :: chmess

  if (lscasrc .AND. nsv.gt.0 .AND. nscasrc.gt.0) then

    allocate(scasrcp(nscasrc,5,nsv))

    ! read global scalar source locations
    if(myid==0) then
      do m=1,nsv
        write (cnsv, '(i1.1)') m
        open (ifinput,file='scalarsourcep.inp.'//cnsv//'.'//cexpnr)
        read (ifinput,'(a80)') chmess
        read (ifinput,'(a80)') chmess
        do n=1,nscasrc
          read (ifinput,*) &
                scasrcp(n,1,m), &
                scasrcp(n,2,m), &
                scasrcp(n,3,m), &
                scasrcp(n,4,m), &
                scasrcp(n,5,m)             
        end do
        close (ifinput)
      end do

      ! write (6,*) 'nsv, source_n, xS, yS, zS, SS, sigS '
      ! do m=1,nsv
      !   do n=1,nscasrc
      !     write (6,*) &
      !           m , &
      !           n , &
      !           scasrcp(n,1,m), &
      !           scasrcp(n,2,m), &
      !           scasrcp(n,3,m), &
      !           scasrcp(n,4,m), &
      !           scasrcp(n,5,m)               
      !   end do
      ! end do

    end if

    call MPI_BCAST(scasrcp ,5*nscasrc*nsv,MY_REAL ,0,comm3d,mpierr)

  end if

  if (lscasrcl .AND. nsv.gt.0 .AND. nscasrcl.gt.0) then

    allocate(scasrcl(nscasrcl,8,nsv))

    ! read global scalar line source locations
    if(myid==0) then
      do m=1,nsv
        write (cnsv, '(i1.1)') m
        open (ifinput,file='scalarsourcel.inp.'//cnsv//'.'//cexpnr)
        read (ifinput,'(a80)') chmess
        read (ifinput,'(a80)') chmess
        do n=1,nscasrcl
          read (ifinput,*) &
                scasrcl(n,1,m), &
                scasrcl(n,2,m), &
                scasrcl(n,3,m), &
                scasrcl(n,4,m), &
                scasrcl(n,5,m), &
                scasrcl(n,6,m), &
                scasrcl(n,7,m), &
                scasrcl(n,8,m)             
        end do
        close (ifinput)
      end do

      ! write (6,*) 'nsv, source_n, xSb, ySb, zSb, xSe, ySe, zSe, SS, sigS '
      ! do m=1,nsv
      !   do n=1,nscasrcl
      !     write (6,*) &
      !           m , &
      !           n , &
      !           scasrcl(n,1,m), &
      !           scasrcl(n,2,m), &
      !           scasrcl(n,3,m), &
      !           scasrcl(n,4,m), &
      !           scasrcl(n,5,m), &
      !           scasrcl(n,6,m), &
      !           scasrcl(n,7,m), &
      !           scasrcl(n,8,m)               
      !   end do
      ! end do

    end if

    call MPI_BCAST(scasrcl ,8*nscasrcl*nsv,MY_REAL ,0,comm3d,mpierr)

  end if

end subroutine


subroutine scalsource

  use modglobal,  only : pi,nsv,ib,ie,jb,je,kb,ke,zf,dx,dy,imax,jmax,&
                         xS,yS,zS,xSb,ySb,zSb,xSe,ySe,zSe,SS,sigS,lscasrc,lscasrcl,dzfi,nscasrc,scasrcp,nscasrcl,scasrcl
  use modfields,  only : svp
  use modmpi,     only : myidx,myidy

  implicit none
  integer :: i,j,k,n,ns
  real :: dxi, dyi
  real :: ra2 = 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

end module modscalsource
