source_population_eq.f

Go to the documentation of this file.

      SUBROUTINE source_population_eq(x,y,dydx)

      use param1, only: zero, small_number
      USE constant
      USE fldvar
      USE physprop
      USE rdf
      USE scalars

      IMPLICIT NONE

      double precision x,y(*),dydx(*)
      double precision K_v
      double precision m1,m2, dav,theta,c11
      double precision epstotal
      integer L,n,np,i,j,k,IJK

      k_v=pi/6
      epstotal=0.0
      theta=0.0


      DO i=1,nscalar
      a(i)= 0.0
      omega(i)=y(i)/(k_v*(y(nscalar+i)**3))
      ENDDO

      n=2*nscalar
      np=2*nscalar

      ijk=ijk_index

! initinalize variables
      DO i=1,n
       DO j=1,n
         matrix_a(i,j)=zero
         matrix_b(i,j)=zero
         matrix_c(i,j)=zero
         inv_a(i,j)=zero
       ENDDO
      ENDDO

      DO i=1,n
       dydx(i)=zero
      ENDDO

      DO i=0,2*nscalar-1
       s_bar(i)=zero
      ENDDO

! calculate BETA_A(I,J) and A(I)

      DO i=1, nscalar
        epstotal=rop_s(ijk,i)/ro_s(ijk,i)+epstotal
      ENDDO

       If(epstotal>=small_number) THEN
           Do j=1,nscalar
            theta= theta+(rop_s(ijk,j)/ro_s(ijk,j))*theta_m(ijk,j)
           ENDDO
            theta= theta/epstotal
        else
           theta=0.0
        endif


       DO i=1,nscalar
       DO j=1,nscalar

        m1=pi*(scalar(ijk,i)**3)*ro_s(ijk,i)/6.0
        m2=pi*(scalar(ijk,j)**3)*ro_s(ijk,j)/6.0
        dav=(scalar(ijk,i)+scalar(ijk,j))/2.0
        c11=((theta*(m1+m2)**2)/(4*pi*m1*m2))**(0.5)*(4.0/dav)
        beta_a(i,j)=aggregation_eff*pi*(dav**3)*g_0(ijk,i,j)*c11
        a(i)=a(i)+pi*(dav**3)*g_0(ijk,i,j)*c11*omega(j)*breakage_eff

      ENDDO
      ENDDO


! calculate S_bar
!
       DO l=0,2*nscalar-1
        DO i=1,nscalar
         DO j=1,nscalar
           If(y(i)<=1.0e-6) THEN
              IF(y(j)<=1.0e-6) THEN
                s_bar(l)=s_bar(l)+0.0
              ELSE
                s_bar(l)=s_bar(l)+&
                beta_a(i,j)*(1.0/2.0)*omega(i)*omega(j)*&
                ((y(nscalar+j)**3)**(l/3.0))
              ENDIF
           ELSE
              IF(y(j)<=1.0e-6)THEN
                s_bar(l)=s_bar(l)+&
                beta_a(i,j)*(1.0/2.0)*omega(i)*omega(j)*&
                ((y(nscalar+i)**3)**(l/3.0))-&
                beta_a(i,j)*(y(nscalar+i)**l)*omega(i)*omega(j)
              ELSE
                s_bar(l)=s_bar(l)+&
                beta_a(i,j)*(1.0/2.0)*omega(i)*omega(j)*&
                (((y(nscalar+i)**3)+(y(nscalar+j)**3))**(l/3.0))-&
                beta_a(i,j)*(y(nscalar+i)**l)*omega(i)*omega(j)
              ENDIF
            ENDIF
          ENDDO
        ENDDO

          DO i=1,nscalar
              IF(y(i)<1.0e-6) THEN
                s_bar(l)=s_bar(l)+0.0
              ELSE
                s_bar(l)=s_bar(l)+&
                a(i)*(2**((3.0-l)/3.0))*(y(nscalar+i)**l)*omega(i)-&
                a(i)*(y(nscalar+i)**l)*omega(i)
              ENDIF
          ENDDO
       ENDDO

! calcalute inv(A)


       IF(nscalar==2) THEN
         IF(abs(y(3)-y(4))<=1.0e-4) THEN
          y(3)=y(3)*0.999
          y(4)=y(4)*1.001
         ENDIF
       ENDIF

       IF(nscalar==3) THEN
         IF(abs(y(4)-y(5))<=1.0e-4) THEN
          y(4)=y(4)*0.999
          y(5)=y(4)*1.001
         ENDIF
         IF(abs(y(5)-y(6))<=1.0e-4) THEN
          y(5)=y(5)*0.999
          y(6)=y(5)*1.001
         ENDIF
         IF(abs(y(4)-y(6))<=1.0e-4) THEN
          y(4)=y(4)*0.999
          y(6)=y(4)*1.001
         ENDIF
       ENDIF

       DO j=1,nscalar
            matrix_a(1,j)=1.0
            matrix_a(2,j)=0.0
       ENDDO

       DO j=nscalar+1,n
            matrix_a(1,j)=0.0
            matrix_a(2,j)=1.0
       ENDDO

       DO i=3,n
         DO j=1,nscalar
          matrix_a(i,j)=(1.0-i+1.0)*(y(j+nscalar)**((i-1)*1.0))
         ENDDO

         DO j=nscalar+1,n
          matrix_a(i,j)=(i-1.0)*(y(j)**((i-2)*1.0))
         ENDDO
       ENDDO




       DO i=1,n
         inv_a(i,i)=1.0
       ENDDO

       call gaussj(matrix_a,n,np,inv_a,n,n)

       DO i=1,n
         IF(i<=nscalar) THEN
           IF(y(i)>1.0e-3) matrix_c(i,i)=-2*(y(i+nscalar)**3)
         ELSE
           IF(y(i-nscalar)>1.0e-3) matrix_c(i,i)= 4*(y(i)**3)
         ENDIF
       ENDDO

       DO i=1,nscalar
         IF(y(i)>1.0e-3) THEN
           matrix_c(i,i+nscalar)=3*(y(i+nscalar)**2)
           matrix_c(i+nscalar,i)=-3*(y(i+nscalar)**4)
         ENDIF
       ENDDO


       DO i = 1, n
          DO j = 1, n
            DO k = 1, n
            matrix_b(i,j) = matrix_b(i,j) + &
                      matrix_c(i,k)*inv_a(k,j)
            ENDDO
           ENDDO
       ENDDO

       DO i=1,nscalar
         DO j=1,n
          dydx(i)=dydx(i)+k_v*matrix_b(i,j)*s_bar(j-1)
         ENDDO
       ENDDO

       DO i=nscalar+1,n
           IF(y(i-nscalar)>1.0e-6) THEN
           dydx(i)=-dydx(i-nscalar)*y(i)/y(i-nscalar)
          ELSE
           dydx(i)=dydx(i)+0.0
          ENDIF
         DO j=1,n
           IF(y(i-nscalar)>1.0e-6) THEN
          dydx(i)=dydx(i)+(k_v/y(i-nscalar))*(matrix_b(i,j)*s_bar(j-1))
           ELSE
            dydx(i)=dydx(i)+0.0
           ENDIF
         ENDDO
       ENDDO


      return
      END