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 evaluates particle vertical diffusion coefficients,
       6             : !! dkz(k,i,j) [cm^2 s^-1].
       7             : !!
       8             : !! Method: Uses equation 8.73 from Seinfeld and Pandis [1998] along
       9             : !! with the slip correction factor (bpm) calculated in the fall
      10             : !! velocity setup.
      11             : !!
      12             : !! This routine requires that vertical profiles of temperature <t>,
      13             : !! air density <rhoa>, viscosity <rmu>, and slip correction <bpm> are
      14             : !! defined (i.e., initatm.f and setupvf.f must be called before this).
      15             : !!
      16             : !! NOTE: Eddy diffusion is carried out by the parent model, so the only
      17             : !! diffusion that CARMA does is Brownian diffusion.
      18             : !!
      19             : !! @author Chuck Bardeen
      20             : !! @version Aug-2010
      21     1050624 : subroutine setupbdif(carma, cstate, rc)
      22             : 
      23             :   !     types
      24             :   use carma_precision_mod
      25             :   use carma_enums_mod
      26             :   use carma_constants_mod
      27             :   use carma_types_mod
      28             :   use carmastate_mod
      29             :   use carma_mod
      30             : 
      31             :   implicit none
      32             : 
      33             :   type(carma_type), intent(in)         :: carma   !! the carma object
      34             :   type(carmastate_type), intent(inout) :: cstate  !! the carma state object
      35             :   integer, intent(inout)               :: rc      !! return code, negative indicates failure
      36             : 
      37             :   ! Local declarations  
      38             :   integer                 :: igroup, ibin, iz, k1, k2, nzm1
      39             : 
      40             :   !  Define formats
      41             :   2 format(/,'Brownian diffusion coefficient (prior to interpolation)')
      42             :   3 format(/,'Particle group ',i3,/,' bin   lev  p [dyne/cm2]  T [K]       r [cm]    wet r [cm]    dkz [cm2/s]',/)
      43             :   4 format(i3,4x,i3,5(1pe11.3,4x)) 
      44             : 
      45             :   
      46             :   !  Loop over all groups.
      47     3151872 :   do igroup = 1, NGROUP
      48             :       
      49             :     ! Loop over particle size bins.
      50    45176832 :     do ibin = 1,NBIN
      51             : 
      52             :       ! Loop over all atltitudes.
      53  2985873408 :       do iz = 1, NZ
      54             : 
      55             :         ! Vertical brownian diffusion coefficient
      56  2983772160 :         dkz(iz,ibin,igroup) = (BK*t(iz)*bpm(iz,ibin,igroup)) / (6._f*PI*rmu(iz)*r_wet(iz,ibin,igroup) * rprat(ibin,igroup))
      57             : 
      58             :       enddo
      59             :     enddo
      60             :   enddo
      61             :     
      62             :   ! Print out diffusivities.
      63             : #ifdef CARMA_DEBUG
      64             :   if (do_print_init) then
      65             :    
      66             :     write(LUNOPRT,2)
      67             : 
      68             :     do igroup = 1, NGROUP
      69             :         
      70             :       write(LUNOPRT,3) igroup
      71             :       
      72             :       do ibin = 1,NBIN
      73             :     
      74             :         do iz = NZ, 1, -1
      75             :           write(LUNOPRT,4) ibin,iz,p(iz),t(iz),r(ibin,igroup),r_wet(iz,ibin,igroup),dkz(iz,ibin,igroup)
      76             :         end do
      77             :       enddo
      78             :     enddo
      79             :   
      80             :     write(LUNOPRT,*) ""
      81             :   end if
      82             : #endif
      83             : 
      84             :   !  Interpolate <dkz> from layer mid-pts to layer boundaries.
      85             :   !  <dkz(k)> is the diffusion coefficient at the lower edge of the layer
      86     1050624 :   nzm1 = max(1, NZ-1)
      87             : 
      88             :   ! Set upper boundary before averaging
      89    45176832 :   dkz(NZP1,:,:) = dkz(NZ,:,:)
      90             : 
      91     1050624 :   if (NZ .gt. 1) then
      92    45176832 :     dkz(NZ,:,:) = sqrt(dkz(nzm1,:,:) * dkz(NZ,:,:))
      93             : 
      94     1050624 :     if (NZ .gt. 2) then
      95    72493056 :       do iz = NZ-1, 2, -1
      96  3073075200 :         dkz(iz,:,:) = sqrt(dkz(iz-1,:,:) * dkz(iz,:,:))
      97             :       enddo
      98             :     endif
      99             :   endif
     100             :   
     101             :   ! Scale cartesian diffusivities to the appropriate vertical coordinate system.
     102             :   ! Non--cartesion coordinates are assumed to be positive downward, but
     103             :   ! vertical velocities in this model are always assumed to be positive upward. 
     104     1050624 :   if( igridv .ne. I_CART )then
     105     3151872 :     do igroup=1,NGROUP
     106    45176832 :       do ibin=1,NBIN
     107  3027898368 :         dkz(:,ibin,igroup) = dkz(:,ibin,igroup) / (zmetl(:)**2)
     108             :       enddo
     109             :     enddo
     110             :   endif
     111             : 
     112             :   ! Return to caller with fall velocities evaluated.
     113     1050624 :   return
     114     1050624 : end