Line data    Source code

       1             : ! Include shortname defintions, so that the F77 code does not have to be modified to
       2             : ! reference the CARMA structure.
       3             : #include "carma_globaer.h"
       4             : 
       5             : !!  This routine calculates new particle concentrations.
       6             : !!
       7             : !!   The basic form from which the solution is derived is
       8             : !!   ( new_value - old_value ) / dtime = source_term - loss_rate*new_value
       9             : !!
      10             : !!  Modified  Sep-1997  (McKie)
      11             : !!    New particle concentrations due to coagulation processes
      12             : !!    were moved to the csolve routine.  Csolve is called to
      13             : !!    update particle concentrations due to coagulation.
      14             : !!    This new psolve now updates particle concentrations due
      15             : !!    to the faster calcs of the non-coag microphysical processes.
      16             : !!
      17             : !! @author Eric Jensen, Bill McKie
      18             : !! @version Oct-1995, Sep-1997
      19 46919733960 : subroutine psolve(carma, cstate, iz, ibin, ielem, rc)
      20             : 
      21             :   ! types
      22             :   use carma_precision_mod
      23             :   use carma_enums_mod
      24             :   use carma_constants_mod
      25             :   use carma_types_mod
      26             :   use carmastate_mod
      27             :   use carma_mod
      28             : 
      29             :   implicit none
      30             : 
      31             :   type(carma_type), intent(in)         :: carma   !! the carma object
      32             :   type(carmastate_type), intent(inout) :: cstate  !! the carma state object
      33             :   integer, intent(in)                  :: iz      !! z index
      34             :   integer, intent(in)                  :: ibin    !! bin index
      35             :   integer, intent(in)                  :: ielem   !! element index
      36             :   integer, intent(inout)               :: rc      !! return code, negative indicates failure
      37             : 
      38             :   ! Local declarations
      39             :   integer                        :: igroup    ! group index
      40             :   integer                        :: iepart
      41             :   integer                        :: igto
      42             :   integer                        :: iz_no_sed
      43             :   real(kind=f)                   :: ppd       ! particle prodocution rate
      44             :   real(kind=f)                   :: pc_nonuc  ! particles - no nucleation
      45             :   real(kind=f)                   :: pls       ! particle loss rate
      46             :   real(kind=f)                   :: sed_rate
      47             :   real(kind=f)                   :: rnuclgtot
      48             :   real(kind=f)                   :: dsed
      49             : 
      50             : 
      51             :   ! Define current group & particle number concentration element indices
      52 46919733960 :   igroup = igelem(ielem)
      53 46919733960 :   iepart = ienconc(igroup)
      54             : 
      55 46919733960 :   if(do_grow) then
      56             :   
      57             :     ! Compute total production rate
      58 46919733960 :     ppd = rnucpe(ibin,ielem) + rhompe(ibin,ielem) + growpe(ibin,ielem) + evappe(ibin,ielem)
      59             :     
      60             :     ! Sum up nucleation loss rates
      61 >14075*10^7 :     rnuclgtot = sum(rnuclg(ibin,igroup,:))
      62             : 
      63             :     ! Compute total loss rate
      64 46919733960 :     pls = rnuclgtot + growlg(ibin,igroup) + evaplg(ibin,igroup) 
      65             : 
      66             :     ! Figure out the new particle concentration without nucleation.
      67 46919733960 :     pc_nonuc = (pc(iz,ibin,ielem) + dtime * (ppd - rnucpe(ibin,ielem) - rhompe(ibin,ielem))) / (ONE + (pls - rnuclgtot) * dtime)
      68             : 
      69             :     ! Update net particle number concentration during current timestep
      70             :     ! due to production and loss rates.
      71 46919733960 :     pc(iz,ibin,ielem) = (pc(iz,ibin,ielem) + dtime * ppd) / (ONE + pls * dtime)
      72             :     
      73             :     ! Now determine the number of particles produced by nucleation as the difference
      74             :     ! between the actual particle count and that done without nucleation rates.
      75             :     !
      76             :     ! NOTE: This is for statistics and is done as a total for the step, not per substep.
      77 46919733960 :     pc_nucl(iz,ibin,ielem) = pc_nucl(iz,ibin,ielem) + (pc(iz,ibin,ielem) - pc_nonuc)
      78             :   end if 
      79             :   
      80             :   ! Prevent particle concentrations from dropping below SMALL_PC
      81 46919733960 :   call smallconc(carma, cstate, iz, ibin, ielem, rc)
      82             : 
      83             :   !  Return to caller with new particle number concentrations.
      84 46919733960 :   return
      85 46919733960 : end