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 solves for sedimentation using an explicit substepping approach. It
       6             : !!  is faster and handles large cfl and irregular grids better than the normal PPM
       7             : !!  solver (versol), but it is more diffusive.
       8             : !!
       9             : !! @author Andy Ackerman, Chuck Bardeen 
      10             : !! version Aug 2010
      11           0 : subroutine versub(carma, cstate, pcmax, cvert, itbnd, ibbnd, ftop, fbot, cvert_tbnd, cvert_bbnd, &
      12           0 :   vertadvu, vertadvd, vertdifu, vertdifd, rc)
      13             : 
      14             :   ! types
      15             :   use carma_precision_mod
      16             :   use carma_enums_mod
      17             :   use carma_constants_mod
      18             :   use carma_types_mod
      19             :   use carmastate_mod
      20             :   use carma_mod
      21             : 
      22             :   implicit none
      23             : 
      24             :   type(carma_type), intent(in)         :: carma           !! the carma object
      25             :   type(carmastate_type), intent(inout) :: cstate          !! the carma state object
      26             :   real(kind=f), intent(in)             :: pcmax(NZ)       !! maximum particle concentration (#/x/y/z)
      27             :   real(kind=f), intent(inout)          :: cvert(NZ)       !! quantity being transported (#/x/y/z)
      28             :   integer, intent(in)                  :: itbnd           !! top boundary condition
      29             :   integer, intent(in)                  :: ibbnd           !! bottom boundary condition
      30             :   real(kind=f), intent(in)             :: ftop            !! flux at top boundary
      31             :   real(kind=f), intent(in)             :: fbot            !! flux at bottom boundary
      32             :   real(kind=f), intent(in)             :: cvert_tbnd      !! quantity at top boundary
      33             :   real(kind=f), intent(in)             :: cvert_bbnd      !! quantity at bottom boundary
      34             :   real(kind=f), intent(in)             :: vertadvu(NZP1)  !! upward vertical transport rate into level k from level k-1 [cm/s]
      35             :   real(kind=f), intent(in)             :: vertadvd(NZP1)  !! downward vertical transport rate into level k from level k-1 [cm/s]
      36             :   real(kind=f), intent(in)             :: vertdifu(NZP1)  !! upward vertical diffusion rate into level k from level k-1 [cm/s]
      37             :   real(kind=f), intent(in)             :: vertdifd(NZP1)  !! downward vertical diffusion rate into level k from level k-1 [cm/s]
      38             :   integer, intent(inout)               :: rc              !! return code, negative indicates failure
      39             : 
      40             :   !  Declare local variables
      41             :   integer                        :: iz
      42             :   integer                        :: istep
      43             :   integer                        :: nstep_sed
      44           0 :   real(kind=f)                   :: fvert(NZ)
      45           0 :   real(kind=f)                   :: up(NZP1)
      46           0 :   real(kind=f)                   :: dn(NZP1)
      47             :   real(kind=f)                   :: cfl_max
      48             :   real(kind=f)                   :: fvert_1
      49             :   real(kind=f)                   :: fvert_nz
      50             : 
      51             :   ! Determine the total upward and downward velocities.
      52           0 :   up(:) = vertadvu(:) + vertdifu(:)
      53           0 :   dn(:) = vertadvd(:) + vertdifd(:)
      54             : 
      55             :   ! Compute the maximum CFL for each bin that has a significant concentration
      56             :   ! of particles.
      57           0 :   cfl_max = 0._f
      58             : 
      59           0 :   do iz = 1, NZ
      60           0 :     if (pcmax(iz) > SMALL_PC) then
      61           0 :       cfl_max = max(cfl_max, max(abs(up(iz)), abs(up(iz+1)), abs(dn(iz)), abs(dn(iz+1))) * dtime / dz(iz))
      62             :     end if
      63             :   end do
      64             : 
      65             :   ! Use the maximum CFL determined above to figure out how much substepping is
      66             :   ! needed to sediment explicitly without violating the CFL anywhere in the column.
      67           0 :   if (cfl_max >= 0._f) then
      68           0 :     nstep_sed = int(1._f + cfl_max)
      69             :   else
      70             :     nstep_sed = 0
      71             :   endif
      72             :   
      73             :   ! If velocities are in both directions, then more steps are needed to make sure
      74             :   ! that no more than half of the concentration can be transported in either direction.
      75           0 :   if (maxval(up(:) * dn(:)) > 0._f) then
      76           0 :     nstep_sed = nstep_sed * 2
      77             :   end if
      78             : 
      79             :   ! Determine the top and bottom boundary fluxes, keeping in mind that
      80             :   ! the velocities and grid coordinates are reversed in sigma or hybrid
      81             :   ! coordinates
      82           0 :   if ((igridv .eq. I_SIG) .or. (igridv .eq. I_HYBRID)) then
      83           0 :     if (itbnd .eq. I_FLUX_SPEC) then
      84           0 :       fvert_nz = -fbot
      85             :     else
      86           0 :       fvert_nz = cvert_bbnd*dn(NZ+1)
      87             :     end if
      88             : 
      89           0 :     if (ibbnd .eq. I_FLUX_SPEC) then
      90           0 :       fvert_1 = -ftop
      91             :     else
      92           0 :       fvert_1 = cvert_tbnd*up(1)
      93             :     end if
      94             :   
      95             :   else
      96           0 :     if (itbnd .eq. I_FLUX_SPEC) then
      97           0 :       fvert_nz = ftop
      98             :     else
      99           0 :       fvert_nz = cvert_tbnd*dn(NZ+1)
     100             :     end if
     101             : 
     102           0 :     if (ibbnd .eq. I_FLUX_SPEC) then
     103           0 :       fvert_1 = fbot
     104             :     else
     105           0 :       fvert_1 = cvert_bbnd*up(1)
     106             :     end if
     107             :   endif
     108             : 
     109             :   ! Sediment the particles using multiple iterations to satisfy the CFL.
     110           0 :   do istep = 1, nstep_sed
     111             : 
     112             :     ! Determine the net particle flux at each gridbox. The first and last levels
     113             :     ! need special treatment to handle to bottom and top boundary conditions.
     114           0 :     fvert(1) = (-cvert(1)*dn(1) + fvert_1 + cvert(2)*dn(2) - cvert(1)*up(2))
     115             :     
     116           0 :     do iz = 2, NZ-1
     117           0 :       fvert(iz) = (-cvert(iz)*dn(iz) + cvert(iz-1)*up(iz) + cvert(iz+1)*dn(iz+1) - cvert(iz)*up(iz+1))
     118             :     end do
     119             : 
     120           0 :     fvert(NZ) = (-cvert(NZ)*dn(NZ) + cvert(NZ-1)*up(NZ) + fvert_nz - cvert(NZ)*up(NZ+1))
     121             :     
     122             :     ! Now update the actual concentrations.
     123           0 :     cvert(:) = cvert(:) + fvert(:) * dtime / nstep_sed / dz(:)
     124             :   enddo
     125             : 
     126           0 :   return
     127           0 : end subroutine versub