module source_term_module
  use common_module
  implicit none
contains
  subroutine source_term(U0,F0,dF)
    !$ use omp_lib
    implicit none
    real(DP),dimension(xlen,ylen,zlen,8),intent(IN) :: U0
    real(DP),dimension(xlen,ylen,zlen,8,3),intent(IN) :: F0
    real(DP),dimension(xlen,ylen,zlen,8),intent(INOUT) :: dF
    real(DP) :: sc(8),se(8)
    integer :: szi,szj,szk,cnt,s(3)
    szi=xr(2)-xr(1)+1; szj=yr(2)-yr(1)+1; szk=zr(2)-zr(1)+1
    !$omp parallel do private(s,sc,se),default(shared)
    do cnt=1,szi*szj*szk
       call loopcounter(cnt,(/szk,szj,szi/),(/zr(1),yr(1),xr(1)/),s)
       call save_iflux(s(3),s(2),s(1),dF(s(3),s(2),s(1),:)&
            *eq0(dble(initial_flux_cancel-1)))
       sc=source_curve(s(3),s(2),s(1),U0,F0)
       se=source_external(s(3),s(2),s(1),U0,F0)
       dF(s(3),s(2),s(1),:)=dF(s(3),s(2),s(1),:)-sc
       call save_iflux(s(3),s(2),s(1),dF(s(3),s(2),s(1),:)&
            *eq0(dble(initial_flux_cancel-2)))
       dF(s(3),s(2),s(1),:)=dF(s(3),s(2),s(1),:)-se
       call save_iflux(s(3),s(2),s(1),dF(s(3),s(2),s(1),:)&
            *eq0(dble(initial_flux_cancel-3)))
       dF(s(3),s(2),s(1),:)=dF(s(3),s(2),s(1),:)-reffer_iflux(s(3),s(2),s(1))
    end do
    !$omp end parallel do
  contains
    function reffer_iflux(i,j,k)
      implicit none
      integer,intent(IN) :: i,j,k
      real(DP) :: reffer_iflux(8)
      integer :: s,sz(4)
      sz=shape(initial_flux)
      reffer_iflux=0._DP
      do s=1,sz(4)
         reffer_iflux(s)=initial_flux(min(i,sz(1)),min(j,sz(2)),min(k,sz(3)),s)
      end do
      reffer_iflux=reffer_iflux*iflux_damp(tax(mod(n-1,taxlen)+1))
    end function reffer_iflux
    subroutine save_iflux(i,j,k,flux)
      implicit none
      integer,intent(IN) :: i,j,k
      real(DP),intent(IN) :: flux(8)
      real(DP) :: update,iflux(8)
      integer :: sz(4),s
      sz=shape(initial_flux)
      update=eq0(maxval(abs(reffer_iflux(i,j,k))))
      iflux=reffer_iflux(i,j,k)
      do s=1,sz(4)
         initial_flux(min(i,sz(1)),min(j,sz(2)),min(k,sz(3)),s)&
              =iflux(s)+update*flux(s)
      end do
    end subroutine save_iflux

    function source_curve(xi,yj,zk,U0,F0)
      implicit none
      real(DP),dimension(xlen,ylen,zlen,8),intent(IN):: U0
      real(DP),dimension(xlen,ylen,zlen,8,3),intent(IN) :: F0
      real(DP),dimension(8) :: source_curve
      real(DP) :: U8(8),hc(3),pc(3),hfR(3),pfR(3),hfL(3),pfL(3),sfR,sfL
      real(DP) :: dh(3,3)
      real(DP) :: FR,FL,TrPi(3),denomi
      integer,intent(IN) :: xi,yj,zk
      integer axis1,axis2,axis3,i3(3)
      source_curve=0._DP
      if (size(hx) .eq. 1) return

      call cell_center(xi,yj,zk,pc,hc)
      do axis1=1,3
         i3=0**abs(axis1-(/1,2,3/))
         axis2=mod(axis1,3)+1 ; axis3=mod(axis2,3)+1
         call cell_face(xi,yj,zk,axis1, 1,pfR,hfR,sfR)
         call cell_face(xi,yj,zk,axis1,-1,pfL,hfL,sfL)
         FR=F0(xi,yj,zk,Vax(axis1),axis1)
         FL=F0(max(xi-i3(1),1),max(yj-i3(2),1),max(zk-i3(3),1),Vax(axis1),axis1)

         TrPi(axis1)=.5_DP*(FR*sfR*hfR(axis1)+FL*sfL*hfL(axis1))&
              /hc(1)/hc(2)/hc(3)
         denomi=pfR(axis1)-pfL(axis1); denomi=ne0(denomi)/(denomi+eq0(denomi))
         do axis2=1,3
            dh(axis2,axis1)=(hfR(axis2)-hfL(axis2))*denomi/hc(axis2)
         end do
      end do
      do axis1=1,3
         source_curve(Vax(axis1))=sum(TrPi*dh(:,axis1))/hc(axis1)
      end do
    end function source_curve
    function source_external(xi,yj,zk,U0,F0)
      implicit none
      real(DP),dimension(xlen,ylen,zlen,8),intent(IN):: U0
      real(DP),dimension(xlen,ylen,zlen,8,3),intent(IN) :: F0
      real(DP),dimension(8) :: source_external
      real(DP) :: U8(8),F8(8,3),exf(3)
      integer,intent(IN) :: xi,yj,zk
      integer axis1,axis2,axis3,kdlt(3)
      source_external=0._DP
      U8=CT_preprocess(U0,xi,yj,zk)
      do axis1=1,3
         kdlt=0**abs(axis1-(/1,2,3/))!*udim(axis1)
         F8(:,axis1)=(F0(xi,yj,zk,:,axis1)&
                     +F0(max(xi-kdlt(1),1),max(yj-kdlt(2),1),&
                         max(zk-kdlt(3),1),:,axis1))*.5_DP
      end do
      !External Force +++++++++++++++++++++++++++++++++++++++++
      do axis1=1,3
         exf(axis1)=eforce(xi,yj,zk,axis1,U0)
         source_external(Vax(axis1))=U8(Sax(1))*exf(axis1)
         source_external(Sax(2))=&
              source_external(Sax(2))+U8(Vax(axis1))*exf(axis1)
      end do
      !--------------------------------------------------------
      !External Heating +++++++++++++++++++++++++++++++++++++++
      source_external(Sax(2))=source_external(Sax(2))+eheat(xi,yj,zk,U0)
      !--------------------------------------------------------
    end function source_external
    subroutine loopcounter(count,sizes,org,subs)
      implicit none
      integer,intent(IN) :: count,sizes(:),org(:)
      integer,intent(OUT) :: subs(:)
      integer :: dim,s,quo(size(subs))
      dim=size(sizes)
      quo=count
      do s=2,dim; quo(s)=(quo(s-1)-1)/sizes(s-1)+1
      end do
      do s=1,dim; subs(s)=mod(quo(s)-1,sizes(s))+org(s)
      end do
    end subroutine loopcounter
  end subroutine source_term
end module source_term_module