phys/module_cam_molec_diff.F

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