des_reaction_model.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: DES_REACTION_MODEL                                     !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!                                                                      !
!  Author: J.Musser                                   Date: 16-Jun-10  !
!                                                                      !
!  Comments:                                                           !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE des_reaction_model

      USE compar
      USE constant
      USE des_rxns
      USE des_thermo
      USE discretelement
      USE geometry
      USE indices
      USE param, only: dimension_n_s
      USE param1, only: zero
      USE run, only: any_species_eq, species_eq
      USE physprop, only: nmax
      USE run, only: dt
      USE run, only: solve_ros
      USE functions

      IMPLICIT NONE

! Passed variables
!-----------------------------------------------
! None

! Local variables
!-----------------------------------------------
! Loop counter
      INTEGER :: NN
! total rate of consumption/production of species (g/sec)
      DOUBLE PRECISION  :: SUM_DES_Rs(1:max_pip)

      DOUBLE PRECISION :: PIx4o3
      DOUBLE PRECISION :: o3 = 1.0d0/3.0d0

      DOUBLE PRECISION :: lDT, lOoDT
! Logical for Adams-Bashfort integration.
      LOGICAL,SAVE:: FIRST_PASS = .true.

!---------------------------------------------------------------------//

      IF(.NOT.any_species_eq) RETURN

      pix4o3 = pi*4.0d0/3.0d0

      ldt = merge(dt, dtsolid, des_explicitly_coupled)
      loodt = -1.0d0/ldt

! Bound the amount of mass loss.
      FORALL(nn=1:dimension_n_s)
         WHERE(particle_state(:max_pip) == normal_particle)         &
            des_r_s(:max_pip,nn) = max(des_r_s(:max_pip,nn),        &
            des_x_s(:max_pip,nn)*pmass(:max_pip)*loodt)
      END FORALL

! First-order method: Euler
      IF(intg_euler) THEN
         WHERE(particle_state(:max_pip) == normal_particle)            &
            sum_des_rs(:max_pip) = sum(des_r_s(:max_pip,:),dim=2)

         WHERE(particle_state(:max_pip) == normal_particle)            &
            pmass(:max_pip) = pmass(:max_pip) + ldt*sum_des_rs(:max_pip)

         FORALL(nn=1:dimension_n_s)
            WHERE(particle_state(:max_pip) == normal_particle)         &
               des_x_s(:max_pip,nn) = max(des_x_s(:max_pip,nn) + ldt*  &
               (des_r_s(:max_pip,nn) - des_x_s(:max_pip,nn)*           &
               sum_des_rs(:max_pip))/pmass(:max_pip), zero)
         END FORALL

      ELSE
         IF(first_pass) THEN
         ENDIF
      ENDIF

      DO nn=1,max_pip
         IF(is_normal(nn)) THEN
            IF(solve_ros(pijk(nn,5))) THEN
               ro_sol(nn)= pmass(nn)/pvol(nn)
            ELSE
               des_radius(nn) = (pmass(nn)/(pix4o3*ro_sol(nn)))**o3
               pvol(nn) = pmass(nn)/ro_sol(nn)
            ENDIF
         ENDIF
      ENDDO

! Clear the necessary variables.
      des_r_s = zero

! Flag that the first pass is over
      first_pass = .false.

      RETURN
      END SUBROUTINE des_reaction_model