usr_dqmom.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: USR_DQMOM                                              !
!  Author: rong fan                                   Date: dd-mmm-yy  !
!                                                                      !
!  Purpose: This routine is called from the time loop and is           !
!           user-definable tosolve teh population equation             !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE usr_dqmom

      USE run, only: call_dqmom, time, dt
      USE scalars, only: nscalar
      USE scalars, only: ijk_index, ystart
      USE fldvar, only: rop_s, ro_s, scalar
      use param, only: dim_scalar2

      USE geometry
      USE indices
      USE compar
      USE usr
      USE functions

      IMPLICIT NONE


      INTEGER :: IJK
! phase index
      double precision t1,t2
! beginning time and stop time
      double precision eps
! error
      double precision h1,hmin
      integer nok, nbad,I,K
      double precision max1, min1

      IF(.NOT.call_dqmom) RETURN

!  Insert user-defined code here
      IF(time<= 1e-15) THEN
         t1 =time
         t2= time
      ELSE
         t1= time-dt
         t2= time
      ENDIF


      eps=1.0e-3
      h1=1.0e-4
      hmin=0
      nok=0
      nbad=0
!      n=2*Nscalar



      DO ijk = ijkstart3, ijkend3

         IF(.NOT.fluid_at(ijk)) cycle

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

         DO i=nscalar+1,2*nscalar
            ystart(i)=scalar(ijk,i-nscalar)
         ENDDO

         ijk_index=ijk
         max1=ystart(1)

         DO k=2,nscalar
            max1=max(max1,ystart(k))
         ENDDO

         min1=ystart(1)

         DO k=2,nscalar
            min1=min(min1,ystart(k))
         ENDDO


         IF(max1>1.0e-3) THEN
            call odeint(ystart,dim_scalar2,t1,t2,eps,h1, hmin,nok,nbad)
         ENDIF

         DO i=1,nscalar
            rop_s(ijk,i) = ystart(i)*ro_s(ijk,i)
            scalar(ijk,i) = ystart(i+nscalar)
         ENDDO

      ENDDO

      RETURN
      END SUBROUTINE usr_dqmom