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 fall velocities for particles. Since there are
       6             : !!  several different approaches, this routine dispatches the call to the
       7             : !!  proper subordinate routine based upon the setup routine defined in the
       8             : !!  particle group.
       9             : !!
      10             : !!
      11             : !! @author Andy Ackerman
      12             : !! @version Mar-2010
      13     1050624 : subroutine setupvf(carma, cstate, rc)
      14             : 
      15             :   !     types
      16             :   use carma_precision_mod
      17             :   use carma_enums_mod
      18             :   use carma_constants_mod
      19             :   use carma_types_mod
      20             :   use carmastate_mod
      21             :   use carma_mod
      22             : 
      23             :   implicit none
      24             : 
      25             :   type(carma_type), intent(in)         :: carma   !! the carma object
      26             :   type(carmastate_type), intent(inout) :: cstate  !! the carma state object
      27             :   integer, intent(inout)               :: rc      !! return code, negative indicates failure
      28             : 
      29             :   ! Local declarations  
      30             :   integer                 :: igroup, i, j, k, k1, k2, ibin, iz, nzm1
      31             : 
      32             :   !  Define formats
      33             :   2 format(/,'Fall velocities and Reynolds number (prior to interpolation)')
      34             :   3 format(/,'Particle group ',i3,' using algorithm ',i3,/, &
      35             :          ' bin   lev  p [dyne/cm2]  T [K]       r [cm]    wet r [cm]       bpm', &
      36             :       '   vf [cm/s]        re'/)
      37             :   4 format(i3,4x,i3,7(1pe11.3,4x)) 
      38             :   
      39             :   !  Loop over all groups.
      40     3151872 :   do igroup = 1, NGROUP
      41             :       
      42             :     ! There are different implementations of the fall velocity calculation. Some of
      43             :     ! these routines may be more appropriate for certain types of partciles.
      44     3151872 :     select case(ifallrtn(igroup))
      45             : 
      46             :       case (I_FALLRTN_STD)
      47     2101248 :         call setupvf_std(carma, cstate, igroup, rc)
      48             :       
      49             :       case(I_FALLRTN_STD_SHAPE)
      50           0 :         call setupvf_std_shape(carma, cstate, igroup, rc)
      51             :       
      52             :       case(I_FALLRTN_HEYMSFIELD2010)
      53           0 :         call setupvf_heymsfield2010(carma, cstate, igroup, rc)
      54             :       
      55             :       case default
      56           0 :         if (do_print) write(LUNOPRT,*) "setupvf:: ERROR - Unknown fall velocity routine  (", ifallrtn(igroup), &
      57           0 :           ") for group (", igroup, ")."
      58           0 :         rc = -1
      59     2101248 :         return
      60             :     end select
      61             :   enddo
      62             :     
      63             :   ! Constant value if <ifall> = 0
      64     1050624 :   if (ifall .eq. 0) then
      65           0 :     vf(:,:,:) = vf_const
      66             :   end if
      67             : 
      68             :   ! Print out fall velocities and reynolds' numbers.
      69             : #ifdef CARMA_DEBUG
      70             :   if (do_print_init) then
      71             :    
      72             :     write(LUNOPRT,2)
      73             : 
      74             :     do j = 1, NGROUP
      75             :         
      76             :       write(LUNOPRT,3) j, ifallrtn(j)
      77             :       
      78             :       do i = 1,NBIN
      79             :     
      80             :         do k = NZ, 1, -1
      81             :           write(LUNOPRT,4) i,k,p(k),t(k),r(i,j),r_wet(k,i,j),bpm(k,i,j),vf(k,i,j),re(k,i,j)
      82             :         end do
      83             :       enddo
      84             :     enddo
      85             :   
      86             :     write(LUNOPRT,*) ""
      87             :   end if
      88             : #endif
      89             : 
      90             :   ! Interpolate <vf> from layer mid-pts to layer boundaries.
      91             :   ! <vf(k)> is the fall velocity at the lower edge of the layer
      92     1050624 :   nzm1 = max(1, NZ-1)
      93             : 
      94             :   ! Set upper boundary before averaging
      95    45176832 :   vf(NZP1,:,:) = vf(NZ,:,:)
      96             : 
      97     1050624 :   if (NZ .gt. 1) then
      98    45176832 :     vf(NZ,:,:) = sqrt(vf(nzm1,:,:) * vf(NZ,:,:))
      99             : 
     100     1050624 :     if (NZ .gt. 2) then
     101    32569344 :       do iz = NZ-1, 2, -1
     102  1356355584 :         vf(iz,:,:) = sqrt(vf(iz-1,:,:) * vf(iz,:,:))
     103             :       enddo
     104             :     endif
     105             :   endif
     106             :   
     107             :   ! Scale cartesian fallspeeds to the appropriate vertical coordinate system.
     108             :   ! Non--cartesion coordinates are assumed to be positive downward, but
     109             :   ! vertical velocities in this model are always assumed to be positive upward. 
     110     1050624 :   if( igridv .ne. I_CART )then
     111     3151872 :     do igroup=1,NGROUP
     112    45176832 :       do ibin=1,NBIN
     113  1430949888 :         vf(:,ibin,igroup) = -vf(:,ibin,igroup) / zmetl(:)
     114             :       enddo
     115             :     enddo
     116             :   endif
     117             : 
     118             :   ! Return to caller with fall velocities evaluated.
     119             :   return
     120     1050624 : end