wrfv2_fire/phys/module_cam_molec_diff.F

   1 #define WRF_PORT
   2
   3 module molec_diff
   4
   5   !------------------------------------------------------------------------------------------------- !
   6   ! Module to compute molecular diffusivity for various constituents                                 !
   7   !                                                                                                  !
   8   ! Public interfaces :                                                                              !
   9   !                                                                                                  !
  10   !    init_molec_diff           Initializes time independent coefficients                           !
  11   !    init_timestep_molec_diff  Time-step initialization for molecular diffusivity                  !
  12   !    compute_molec_diff        Computes constituent-independent terms for moleculuar diffusivity   !
  13   !    vd_lu_qdecomp             Computes constituent-dependent terms for moleculuar diffusivity and !
  14   !                              updates terms in the triadiagonal matrix used for the implicit      !
  15   !                              solution of the diffusion equation                                  !
  16   !                                                                                                  !
  17   !---------------------------Code history---------------------------------------------------------- !
  18   ! Modularized     :  J. McCaa, September 2004                                                      !
  19   ! Lastly Arranged :  S. Park,  January.  2010                                                      !
  20   !------------------------------------------------------------------------------------------------- !
  21
  22 #ifndef WRF_PORT
  23     use perf_mod
  24     use cam_logfile,  only : iulog
  25 #else
  26   use module_cam_support,   only: iulog, t_stopf, t_startf
  27 #endif
  28
  29   implicit none
  30   private
  31   save
  32
  33   public init_molec_diff
  34 #ifndef WRF_PORT
  35   public init_timestep_molec_diff
  36 #endif
  37   public compute_molec_diff
  38   public vd_lu_qdecomp
  39
  40   ! ---------- !
  41   ! Parameters !
  42   ! ---------- !
  43
  44   integer,  parameter   :: r8 = selected_real_kind(12) ! 8 byte real
  45
  46   real(r8), parameter   :: km_fac = 3.55E-7_r8         ! Molecular viscosity constant [ unit ? ]
  47   real(r8), parameter   :: pr_num = 1._r8              ! Prandtl number [ no unit ]
  48   real(r8), parameter   :: pwr    = 2._r8/3._r8        ! Exponentiation factor [ unit ? ]
  49   real(r8), parameter   :: d0     = 1.52E20_r8         ! Diffusion factor [ m-1 s-1 ] molec sqrt(kg/kmol/K) [ unit ? ]
  50                                                        ! Aerononmy, Part B, Banks and Kockarts (1973), p39
  51                                                        ! Note text cites 1.52E18 cm-1 ...
  52
  53   real(r8)              :: rair                        ! Gas constant for dry air
  54   real(r8)              :: mw_dry                      ! Molecular weight of dry air
  55   real(r8)              :: n_avog                      ! Avogadro's number [ molec/kmol ]
  56   real(r8)              :: gravit
  57   real(r8)              :: cpair
  58   real(r8)              :: kbtz                        ! Boltzman constant
  59
  60   integer               :: ntop_molec                  ! Top    interface level to which molecular vertical diffusion is applied ( = 1 )
  61   integer               :: nbot_molec                  ! Bottom interface level to which molecular vertical diffusion is applied ( = pver )
  62   real(r8), allocatable :: mw_fac(:)                   ! sqrt(1/M_q + 1/M_d) in constituent diffusivity [  unit ? ]
  63
  64   contains
  65
  66   !============================================================================ !
  67   !                                                                             !
  68   !============================================================================ !
  69
  70   subroutine init_molec_diff( kind, ncnst, rair_in, ntop_molec_in, nbot_molec_in, &
  71                               mw_dry_in, n_avog_in, gravit_in, cpair_in, kbtz_in )
  72
  73     use constituents,     only : cnst_mw
  74     use upper_bc,         only : ubc_init
  75
  76     integer,  intent(in)  :: kind           ! Kind of reals being passed in
  77     integer,  intent(in)  :: ncnst          ! Number of constituents
  78     integer,  intent(in)  :: ntop_molec_in  ! Top interface level to which molecular vertical diffusion is applied ( = 1 )
  79     integer,  intent(in)  :: nbot_molec_in  ! Bottom interface level to which molecular vertical diffusion is applied.
  80     real(r8), intent(in)  :: rair_in
  81     real(r8), intent(in)  :: mw_dry_in      ! Molecular weight of dry air
  82     real(r8), intent(in)  :: n_avog_in      ! Avogadro's number [ molec/kmol ]
  83     real(r8), intent(in)  :: gravit_in
  84     real(r8), intent(in)  :: cpair_in
  85     real(r8), intent(in)  :: kbtz_in        ! Boltzman constant
  86     integer               :: m              ! Constituent index
  87
  88     if( kind .ne. r8 ) then
  89         write(iulog,*) 'KIND of reals passed to init_molec_diff -- exiting.'
  90         stop 'init_molec_diff'
  91     endif
  92
  93     rair       = rair_in
  94     mw_dry     = mw_dry_in
  95     n_avog     = n_avog_in
  96     gravit     = gravit_in
  97     cpair      = cpair_in
  98     kbtz       = kbtz_in
  99     ntop_molec = ntop_molec_in
 100     nbot_molec = nbot_molec_in
 101
 102   ! Initialize upper boundary condition variables
 103
 104     call ubc_init()
 105
 106   ! Molecular weight factor in constitutent diffusivity
 107   ! ***** FAKE THIS FOR NOW USING MOLECULAR WEIGHT OF DRY AIR FOR ALL TRACERS ****
 108
 109     allocate(mw_fac(ncnst))
 110     do m = 1, ncnst
 111        mw_fac(m) = d0 * mw_dry * sqrt(1._r8/mw_dry + 1._r8/cnst_mw(m)) / n_avog
 112     end do
 113
 114   end subroutine init_molec_diff
 115
 116   !============================================================================ !
 117   !                                                                             !
 118   !============================================================================ !
 119 #ifndef WRF_PORT
 120   subroutine init_timestep_molec_diff(state)
 121     !--------------------------- !
 122     ! Timestep dependent setting !
 123     !--------------------------- !
 124     use upper_bc,     only : ubc_timestep_init
 125     use physics_types,only: physics_state
 126     use ppgrid,       only: begchunk, endchunk
 127
 128     type(physics_state), intent(in) :: state(begchunk:endchunk)
 129
 130     call ubc_timestep_init( state)
 131
 132   end subroutine init_timestep_molec_diff
 133 #endif
 134   !============================================================================ !
 135   !                                                                             !
 136   !============================================================================ !
 137
 138   integer function compute_molec_diff( lchnk       ,                                                                      &
 139                                        pcols       , pver       , ncnst      , ncol     , t      , pmid   , pint        , &
 140                                        zi          , ztodt      , kvh        , kvm      , tint   , rhoi   , tmpi2       , &
 141                                        kq_scal     , ubc_t      , ubc_mmr    , dse_top  , cc_top , cd_top , cnst_mw_out , &
 142                                        cnst_fixed_ubc_out , mw_fac_out , ntop_molec_out , nbot_molec_out )
 143
 144     use upper_bc,        only : ubc_get_vals
 145     use constituents,    only : cnst_mw, cnst_fixed_ubc
 146
 147     ! --------------------- !
 148     ! Input-Output Argument !
 149     ! --------------------- !
 150
 151     integer,  intent(in)    :: pcols
 152     integer,  intent(in)    :: pver
 153     integer,  intent(in)    :: ncnst
 154     integer,  intent(in)    :: ncol                      ! Number of atmospheric columns
 155     integer,  intent(in)    :: lchnk                     ! Chunk identifier
 156     real(r8), intent(in)    :: t(pcols,pver)             ! Temperature input
 157     real(r8), intent(in)    :: pmid(pcols,pver)          ! Midpoint pressures
 158     real(r8), intent(in)    :: pint(pcols,pver+1)        ! Interface pressures
 159     real(r8), intent(in)    :: zi(pcols,pver+1)          ! Interface heights
 160     real(r8), intent(in)    :: ztodt                     ! 2 delta-t
 161
 162     real(r8), intent(inout) :: kvh(pcols,pver+1)         ! Diffusivity for heat
 163     real(r8), intent(inout) :: kvm(pcols,pver+1)         ! Viscosity ( diffusivity for momentum )
 164     real(r8), intent(inout) :: tint(pcols,pver+1)        ! Interface temperature
 165     real(r8), intent(inout) :: rhoi(pcols,pver+1)        ! Density ( rho ) at interfaces
 166     real(r8), intent(inout) :: tmpi2(pcols,pver+1)       ! dt*(g*rho)**2/dp at interfaces
 167
 168     real(r8), intent(out)   :: kq_scal(pcols,pver+1)     ! kq_fac*sqrt(T)*m_d/rho for molecular diffusivity
 169     real(r8), intent(out)   :: ubc_mmr(pcols,ncnst)      ! Upper boundary mixing ratios [ kg/kg ]
 170     real(r8), intent(out)   :: cnst_mw_out(ncnst)
 171     logical,  intent(out)   :: cnst_fixed_ubc_out(ncnst)
 172     real(r8), intent(out)   :: mw_fac_out(ncnst)
 173     real(r8), intent(out)   :: dse_top(pcols)            ! dse on top boundary
 174     real(r8), intent(out)   :: cc_top(pcols)             ! Lower diagonal at top interface
 175     real(r8), intent(out)   :: cd_top(pcols)             ! cc_top * dse ubc value
 176     integer,  intent(out)   :: ntop_molec_out
 177     integer,  intent(out)   :: nbot_molec_out
 178
 179     ! --------------- !
 180     ! Local variables !
 181     ! --------------- !
 182
 183     integer                 :: m                          ! Constituent index
 184     integer                 :: i                          ! Column index
 185     integer                 :: k                          ! Level index
 186     real(r8)                :: mkvisc                     ! Molecular kinematic viscosity c*tint**(2/3)/rho
 187     real(r8)                :: ubc_t(pcols)               ! Upper boundary temperature (K)
 188
 189     ! ----------------------- !
 190     ! Main Computation Begins !
 191     ! ----------------------- !
 192
 193   ! Get upper boundary values
 194
 195     call ubc_get_vals( lchnk, ncol, ntop_molec, pint, zi, ubc_t, ubc_mmr )
 196
 197   ! Below are already computed, just need to be copied for output
 198
 199     cnst_mw_out(:ncnst)        = cnst_mw(:ncnst)
 200     cnst_fixed_ubc_out(:ncnst) = cnst_fixed_ubc(:ncnst)
 201     mw_fac_out(:ncnst)         = mw_fac(:ncnst)
 202     ntop_molec_out             = ntop_molec
 203     nbot_molec_out             = nbot_molec
 204
 205   ! Density and related factors for moecular diffusion and ubc.
 206   ! Always have a fixed upper boundary T if molecular diffusion is active. Why ?
 207
 208     tint (:ncol,ntop_molec) = ubc_t(:ncol)
 209     rhoi (:ncol,ntop_molec) = pint(:ncol,ntop_molec) / ( rair * tint(:ncol,ntop_molec) )
 210     tmpi2(:ncol,ntop_molec) = ztodt * ( gravit * rhoi(:ncol,ntop_molec))**2 &
 211                                     / ( pmid(:ncol,ntop_molec) - pint(:ncol,ntop_molec) )
 212
 213   ! Compute molecular kinematic viscosity, heat diffusivity and factor for constituent diffusivity
 214   ! This is a key part of the code.
 215
 216     kq_scal(:ncol,1:ntop_molec-1) = 0._r8
 217     do k = ntop_molec, nbot_molec
 218        do i = 1, ncol
 219           mkvisc       = km_fac * tint(i,k)**pwr / rhoi(i,k)
 220           kvm(i,k)     = kvm(i,k) + mkvisc
 221           kvh(i,k)     = kvh(i,k) + mkvisc * pr_num
 222           kq_scal(i,k) = sqrt(tint(i,k)) / rhoi(i,k)
 223        end do
 224     end do
 225     kq_scal(:ncol,nbot_molec+1:pver+1) = 0._r8
 226
 227   ! Top boundary condition for dry static energy
 228
 229     dse_top(:ncol) = cpair * tint(:ncol,ntop_molec) + gravit * zi(:ncol,ntop_molec)
 230
 231   ! Top value of cc for dry static energy
 232
 233     do i = 1, ncol
 234        cc_top(i) = ztodt * gravit**2 * rhoi(i,ntop_molec) * km_fac * ubc_t(i)**pwr / &
 235                    ( ( pint(i,2) - pint(i,1) ) * ( pmid(i,1) - pint(i,1) ) )
 236     enddo
 237     cd_top(:ncol) = cc_top(:ncol) * dse_top(:ncol)
 238
 239     compute_molec_diff = 1
 240     return
 241   end function compute_molec_diff
 242
 243   !============================================================================ !
 244   !                                                                             !
 245   !============================================================================ !
 246
 247   integer function vd_lu_qdecomp( pcols , pver   , ncol       , fixed_ubc  , mw     , ubc_mmr , &
 248                                   kv    , kq_scal, mw_facm    , tmpi       , rpdel  ,           &
 249                                   ca    , cc     , dnom       , ze         , rhoi   ,           &
 250                                   tint  , ztodt  , ntop_molec , nbot_molec , cd_top )
 251
 252     !------------------------------------------------------------------------------ !
 253     ! Add the molecular diffusivity to the turbulent diffusivity for a consitutent. !
 254     ! Update the superdiagonal (ca(k)), diagonal (cb(k)) and subdiagonal (cc(k))    !
 255     ! coefficients of the tridiagonal diffusion matrix, also ze and denominator.    !
 256     !------------------------------------------------------------------------------ !
 257
 258     ! ---------------------- !
 259     ! Input-Output Arguments !
 260     ! ---------------------- !
 261
 262     integer,  intent(in)    :: pcols
 263     integer,  intent(in)    :: pver
 264     integer,  intent(in)    :: ncol                  ! Number of atmospheric columns
 265
 266     integer,  intent(in)    :: ntop_molec
 267     integer,  intent(in)    :: nbot_molec
 268
 269     logical,  intent(in)    :: fixed_ubc             ! Fixed upper boundary condition flag
 270     real(r8), intent(in)    :: kv(pcols,pver+1)      ! Eddy diffusivity
 271     real(r8), intent(in)    :: kq_scal(pcols,pver+1) ! Molecular diffusivity ( kq_fac*sqrt(T)*m_d/rho )
 272     real(r8), intent(in)    :: mw                    ! Molecular weight for this constituent
 273     real(r8), intent(in)    :: ubc_mmr(pcols)        ! Upper boundary mixing ratios [ kg/kg ]
 274     real(r8), intent(in)    :: mw_facm               ! sqrt(1/M_q + 1/M_d) for this constituent
 275     real(r8), intent(in)    :: tmpi(pcols,pver+1)    ! dt*(g/R)**2/dp*pi(k+1)/(.5*(tm(k+1)+tm(k))**2
 276     real(r8), intent(in)    :: rpdel(pcols,pver)     ! 1./pdel ( thickness bet interfaces )
 277     real(r8), intent(in)    :: rhoi(pcols,pver+1)    ! Density at interfaces [ kg/m3 ]
 278     real(r8), intent(in)    :: tint(pcols,pver+1)    ! Interface temperature [ K ]
 279     real(r8), intent(in)    :: ztodt                 ! 2 delta-t [ s ]
 280
 281     real(r8), intent(inout) :: ca(pcols,pver)        ! -Upper diagonal
 282     real(r8), intent(inout) :: cc(pcols,pver)        ! -Lower diagonal
 283     real(r8), intent(inout) :: dnom(pcols,pver)      ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1)) , 1./(b(k) - c(k)*e(k-1))
 284     real(r8), intent(inout) :: ze(pcols,pver)        ! Term in tri-diag. matrix system
 285
 286     real(r8), intent(out)   :: cd_top(pcols)         ! Term for updating top level with ubc
 287
 288     ! --------------- !
 289     ! Local Variables !
 290     ! --------------- !
 291
 292     integer                 :: i                     ! Longitude index
 293     integer                 :: k, kp1                ! Vertical indicies
 294
 295     real(r8)                :: rghd(pcols,pver+1)    ! (1/H_i - 1/H) * (rho*g)^(-1)
 296     real(r8)                :: kmq(ncol)             ! Molecular diffusivity for constituent
 297     real(r8)                :: wrk0(ncol)            ! Work variable
 298     real(r8)                :: wrk1(ncol)            ! Work variable
 299
 300     real(r8)                :: cb(pcols,pver)        ! - Diagonal
 301     real(r8)                :: kvq(pcols,pver+1)     ! Output vertical diffusion coefficient
 302
 303     ! ----------------------- !
 304     ! Main Computation Begins !
 305     ! ----------------------- !
 306
 307     ! --------------------------------------------------------------------- !
 308     ! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the !
 309     ! tridiagonal diffusion matrix. The diagonal elements  (cb=1+ca+cc) are !
 310     ! a combination of ca and cc; they are not required by the solver.      !
 311     !---------------------------------------------------------------------- !
 312
 313     call t_startf('vd_lu_qdecomp')
 314
 315     kvq(:,:)  = 0._r8
 316     cd_top(:) = 0._r8
 317
 318   ! Compute difference between scale heights of constituent and dry air
 319
 320     do k = ntop_molec, nbot_molec
 321        do i = 1, ncol
 322           rghd(i,k) = gravit / ( kbtz * n_avog * tint(i,k) ) * ( mw - mw_dry )
 323           rghd(i,k) = ztodt * gravit * rhoi(i,k) * rghd(i,k)
 324        enddo
 325     enddo
 326
 327     !-------------------- !
 328     ! Molecular diffusion !
 329     !-------------------- !
 330
 331     do k = nbot_molec - 1, ntop_molec, -1
 332        kp1 = k + 1
 333        kmq(:ncol)  = kq_scal(:ncol,kp1) * mw_facm
 334        wrk0(:ncol) = ( kv(:ncol,kp1) + kmq(:ncol) ) * tmpi(:ncol,kp1)
 335        wrk1(:ncol) = kmq(:ncol) * 0.5_r8 * rghd(:ncol,kp1)
 336      ! Add species separation term
 337        ca(:ncol,k  )  = ( wrk0(:ncol) - wrk1(:ncol) ) * rpdel(:ncol,k)
 338        cc(:ncol,kp1)  = ( wrk0(:ncol) + wrk1(:ncol) ) * rpdel(:ncol,kp1)
 339        kvq(:ncol,kp1) = kmq(:ncol)
 340     end do
 341
 342     if( fixed_ubc ) then
 343         cc(:ncol,ntop_molec) = kq_scal(:ncol,ntop_molec) * mw_facm                 &
 344                              * ( tmpi(:ncol,ntop_molec) + rghd(:ncol,ntop_molec) ) &
 345                              * rpdel(:ncol,ntop_molec)
 346     end if
 347
 348   ! Calculate diagonal elements
 349
 350     do k = nbot_molec - 1, ntop_molec + 1, -1
 351        kp1 = k + 1
 352        cb(:ncol,k) = 1._r8 + ca(:ncol,k) + cc(:ncol,k)                   &
 353                    + rpdel(:ncol,k) * ( kvq(:ncol,kp1) * rghd(:ncol,kp1) &
 354                    - kvq(:ncol,k) * rghd(:ncol,k) )
 355        kvq(:ncol,kp1) = kv(:ncol,kp1) + kvq(:ncol,kp1)
 356     end do
 357
 358     k   = ntop_molec
 359     kp1 = k + 1
 360     if( fixed_ubc ) then
 361         cb(:ncol,k) = 1._r8 + ca(:ncol,k)                                 &
 362                     + rpdel(:ncol,k) * kvq(:ncol,kp1) * rghd(:ncol,kp1)   &
 363                     + kq_scal(:ncol,ntop_molec) * mw_facm                 &
 364                     * ( tmpi(:ncol,ntop_molec) - rghd(:ncol,ntop_molec) ) &
 365                     * rpdel(:ncol,ntop_molec)
 366     else
 367         cb(:ncol,k) = 1._r8 + ca(:ncol,k) &
 368                     + rpdel(:ncol,k) * kvq(:ncol,kp1) * rghd(:ncol,kp1)
 369     end if
 370
 371     k   = nbot_molec
 372     cb(:ncol,k) = 1._r8 + cc(:ncol,k) + ca(:ncol,k) &
 373                 - rpdel(:ncol,k) * kvq(:ncol,k)*rghd(:ncol,k)
 374     do k = 1, nbot_molec + 1, -1
 375        cb(:ncol,k) = 1._r8 + ca(:ncol,k) + cc(:ncol,k)
 376     end do
 377
 378   ! Compute term for updating top level mixing ratio for ubc
 379
 380     if( fixed_ubc ) then
 381         cd_top(:ncol) = cc(:ncol,ntop_molec) * ubc_mmr(:ncol)
 382     end if
 383
 384     !-------------------------------------------------------- !
 385     ! Calculate e(k).                                         !
 386     ! This term is required in solution of tridiagonal matrix !
 387     ! defined by implicit diffusion equation.                 !
 388     !-------------------------------------------------------- !
 389
 390     do k = nbot_molec, ntop_molec + 1, -1
 391        dnom(:ncol,k) = 1._r8 / ( cb(:ncol,k) - ca(:ncol,k) * ze(:ncol,k+1) )
 392        ze(:ncol,k)   = cc(:ncol,k) * dnom(:ncol,k)
 393     end do
 394     k = ntop_molec
 395     dnom(:ncol,k) = 1._r8 / ( cb(:ncol,k) - ca(:ncol,k) * ze(:ncol,k+1) )
 396
 397     vd_lu_qdecomp = 1
 398     call t_stopf('vd_lu_qdecomp')
 399     return
 400
 401   end function vd_lu_qdecomp
 402
 403   end module molec_diff