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 vertrical transport rates.
       6             : !! Currently treats diffusion only.
       7             : !! Not necessarily generalized for irregular grid.
       8             : !!
       9             : !! @author Eric Jensen
      10             : !! @version Dec-1996
      11   145973900 : subroutine vertdif(carma, cstate, igroup, ibin, itbnd, ibbnd, vertdifu, vertdifd, rc)
      12             : 
      13             :   ! types
      14             :   use carma_precision_mod
      15             :   use carma_enums_mod
      16             :   use carma_constants_mod
      17             :   use carma_types_mod
      18             :   use carmastate_mod
      19             :   use carma_mod
      20             : 
      21             :   implicit none
      22             : 
      23             :   type(carma_type), intent(in)         :: carma           !! the carma object
      24             :   type(carmastate_type), intent(inout) :: cstate          !! the carma state object
      25             :   integer, intent(in)                  :: igroup          !! particle group index
      26             :   integer, intent(in)                  :: ibin            !! particle bin index
      27             :   integer, intent(in)                  :: itbnd           !! top boundary condition
      28             :   integer, intent(in)                  :: ibbnd           !! bottom boundary condition
      29             :   real(kind=f), intent(out)            :: vertdifu(NZP1)  !! upward vertical diffusion rate into level k from level k-1 [cm/s]
      30             :   real(kind=f), intent(out)            :: vertdifd(NZP1)  !! downward vertical diffusion rate into level k from level k-1 [cm/s]
      31             :   integer, intent(inout)               :: rc              !! return code, negative indicates failure
      32             :   
      33             :   ! Local Variables
      34             :   integer      :: k
      35             :   integer      :: nzm1
      36             :   integer      :: itwo
      37             :   real(kind=f) :: dz_avg
      38             :   real(kind=f) :: rhofact
      39             :   real(kind=f) :: xex
      40             :   real(kind=f) :: ttheta
      41             : 
      42             : 
      43             :   ! Set some constants
      44   145973900 :   nzm1 = max( 1, NZ-1 )
      45   145973900 :   itwo = min( 2, NZ   )
      46             : 
      47             :   !  Loop over vertical levels.
      48 10218173000 :   do k = 2, NZ
      49             : 
      50 10072199100 :     dz_avg = dz(k)                            ! layer thickness
      51             : 
      52             :     !  Check the vertical coordinate
      53             : 
      54 10218173000 :     if( igridv .eq. I_CART ) then
      55           0 :       rhofact = log(  rhoa(k)/rhoa(k-1) &
      56           0 :                     * zmet(k-1)/zmet(k) )
      57             :       xex = rhoa(k-1)/rhoa(k) * &
      58           0 :             zmet(k)/zmet(k-1)
      59           0 :       vertdifu(k) = ( rhofact * dkz(k, ibin, igroup) / dz_avg ) / ( 1._f - xex )
      60             : 
      61           0 :       vertdifd(k) = vertdifu(k) * xex
      62             : 
      63             : 
      64             :     !  ...else you're in sigma or hybrid coordinates...
      65 10072199100 :     elseif(( igridv .eq. I_SIG ) .or. ( igridv .eq. I_HYBRID )) then
      66 10072199100 :       vertdifu(k) = dkz(k, ibin, igroup) / dz_avg
      67 10072199100 :       vertdifd(k) = dkz(k, ibin, igroup) / dz_avg
      68             : 
      69             :     !  ...else write an error (maybe redundant)...
      70             :     else
      71           0 :       if (do_print) write(LUNOPRT,*) 'vertdif::ERROR - Invalid vertical grid type (', igridv, ').'
      72           0 :       rc = -1
      73           0 :       return
      74             :     endif
      75             :   enddo
      76             : 
      77             :   ! Fluxes at boundaries specified by user
      78   145973900 :   if( ibbnd .eq. I_FLUX_SPEC ) then
      79           0 :     vertdifu(1) = 0._f
      80           0 :     vertdifd(1) = 0._f
      81             :   endif
      82             : 
      83   145973900 :   if( itbnd .eq. I_FLUX_SPEC ) then
      84           0 :     vertdifu(NZ+1) = 0._f
      85           0 :     vertdifd(NZ+1) = 0._f
      86             :   endif
      87             : 
      88             :   ! Diffusion across boundaries using fixed boundary concentration:
      89   145973900 :   if( ibbnd .eq. I_FIXED_CONC ) then
      90   145973900 :     dz_avg = dz(1)                            ! layer thickness
      91   145973900 :     rhofact = log( rhoa(itwo)/rhoa(1) )
      92   145973900 :     ttheta = rhofact
      93   145973900 :     if( ttheta .ge. 0._f ) then
      94   102515140 :       ttheta = min(ttheta,POWMAX)
      95             :     else
      96    43458760 :       ttheta = max(ttheta,-POWMAX)
      97             :     endif
      98             : 
      99   145973900 :     xex = exp(-ttheta)
     100   145973900 :     if( abs(ONE - xex) .lt. ALMOST_ZERO ) xex = ALMOST_ONE
     101             : 
     102   145973900 :     vertdifu(1) = ( rhofact * dkz(1, ibin, igroup) / dz_avg ) / ( 1._f - xex )
     103   145973900 :     vertdifd(1) = vertdifu(1) * xex
     104             :   endif
     105             : 
     106   145973900 :   if( itbnd .eq. I_FIXED_CONC ) then
     107   145973900 :     dz_avg = dz(NZ)                            ! layer thickness
     108   145973900 :     rhofact = log( rhoa(NZ)/rhoa(nzm1) )
     109   145973900 :     ttheta = rhofact
     110   145973900 :     if( ttheta .ge. 0._f ) then
     111    66839440 :       ttheta = min(ttheta,POWMAX)
     112             :     else
     113    79134460 :       ttheta = max(ttheta,-POWMAX)
     114             :     endif
     115             : 
     116   145973900 :     xex = exp(-ttheta)
     117   145973900 :     if( abs(ONE - xex) .lt. ALMOST_ZERO ) xex = ALMOST_ONE
     118             : 
     119   145973900 :     vertdifu(NZ+1) = ( rhofact * dkz(NZ+1, ibin, igroup) / dz_avg ) / ( 1._f - xex )
     120   145973900 :     vertdifd(NZ+1) = vertdifu(NZ+1) * xex
     121             :   endif
     122             : 
     123             :   ! Return to caller with vertical diffusion rates.
     124             :   return
     125   145973900 : end