LCOV - code coverage report
Current view: | top level - physics/cam - molec_diff.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 104 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 3 | 0.0 %

              Line data    Source code
           1             : module molec_diff
           2             :
           3             :   !------------------------------------------------------------------------------------------------- !
           4             :   ! Module to compute molecular diffusivity for various constituents                                 !
           5             :   !                                                                                                  !
           6             :   ! Public interfaces :                                                                              !
           7             :   !                                                                                                  !
           8             :   !    init_molec_diff           Initializes time independent coefficients                           !
           9             :   !    init_timestep_molec_diff  Time-step initialization for molecular diffusivity                  !
          10             :   !    compute_molec_diff        Computes constituent-independent terms for moleculuar diffusivity   !
          11             :   !    vd_lu_qdecomp             Computes constituent-dependent terms for moleculuar diffusivity and !
          12             :   !                              updates terms in the triadiagonal matrix used for the implicit      !
          13             :   !                              solution of the diffusion equation                                  !
          14             :   !                                                                                                  !
          15             :   !---------------------------Code history---------------------------------------------------------- !
          16             :   ! Modularized     :  J. McCaa, September 2004                                                      !
          17             :   ! Lastly Arranged :  S. Park,  January.  2010                                                      !
          18             :   !                    M. Mills, November  2011
          19             :   !------------------------------------------------------------------------------------------------- !
          20             :
          21             :   use perf_mod
          22             :   use air_composition, only: mbarv
          23             :   use phys_control,    only: waccmx_is             !WACCM-X runtime switch
          24             :
          25             :   implicit none
          26             :   private
          27             :   save
          28             :
          29             :   public init_molec_diff
          30             :   public compute_molec_diff
          31             :   public vd_lu_qdecomp
          32             :
          33             :   ! ---------- !
          34             :   ! Parameters !
          35             :   ! ---------- !
          36             :
          37             :   integer,  parameter   :: r8 = selected_real_kind(12) ! 8 byte real
          38             :
          39             :   real(r8), parameter   :: km_fac = 3.55E-7_r8         ! Molecular viscosity constant [ unit ? ]
          40             :   real(r8), parameter   :: pwr    = 2._r8/3._r8        ! Exponentiation factor [ no unit ]
          41             :   real(r8), parameter   :: d0     = 1.52E20_r8         ! Diffusion factor [ m-1 s-1 ] molec sqrt(kg/kmol/K) [ unit ? ]
          42             :                                                        ! Aerononmy, Part B, Banks and Kockarts (1973), p39
          43             :                                                        ! Note text cites 1.52E18 cm-1 ...
          44             :
          45             :   real(r8)              :: mw_dry                      ! Molecular weight of dry air
          46             :   real(r8)              :: n_avog                      ! Avogadro's number [ molec/kmol ]
          47             :
          48             :   real(r8), allocatable :: mw_fac(:)                   ! sqrt(1/M_q + 1/M_d) in constituent diffusivity [  unit ? ]
          49             :   real(r8), allocatable :: alphath(:)                  ! Thermal diffusion factor, -0.38 for H, 0 for others
          50             :
          51             :   logical :: waccmx_mode = .false.
          52             :
          53             : contains
          54             :
          55             :   !============================================================================ !
          56             :   !                                                                             !
          57             :   !============================================================================ !
          58             :
          59           0 :   subroutine init_molec_diff( kind, ncnst, mw_dry_in, n_avog_in, &
          60             :                               errstring)
          61             :
          62             :     use constituents,     only: cnst_mw, cnst_get_ind
          63             :
          64             :     integer,  intent(in)  :: kind           ! Kind of reals being passed in
          65             :     integer,  intent(in)  :: ncnst          ! Number of constituents
          66             :     real(r8), intent(in)  :: mw_dry_in      ! Molecular weight of dry air
          67             :     real(r8), intent(in)  :: n_avog_in      ! Avogadro's number [ molec/kmol ]
          68             :
          69             :     character(len=*), intent(out) :: errstring
          70             :
          71             :     ! Local
          72             :
          73             :     integer               :: m              ! Constituent index
          74             :     integer               :: indx_H         ! Constituent index for H
          75             :     integer               :: ierr           ! Allocate error check
          76             :
          77           0 :     errstring = ' '
          78             :
          79           0 :     mw_dry     = mw_dry_in
          80           0 :     n_avog     = n_avog_in
          81             :
          82           0 :     if( kind /= r8 ) then
          83           0 :        errstring = 'inconsistent KIND of reals passed to init_molec_diff'
          84           0 :        return
          85             :     end if
          86             :
          87             :     ! Determine whether WACCM-X is on.
          88           0 :     waccmx_mode = waccmx_is('ionosphere') .or. waccmx_is('neutral')
          89             :
          90             :   ! Molecular weight factor in constitutent diffusivity
          91             :   ! ***** FAKE THIS FOR NOW USING MOLECULAR WEIGHT OF DRY AIR FOR ALL TRACERS ****
          92             :
          93           0 :     allocate(mw_fac(ncnst))
          94           0 :     do m = 1, ncnst
          95           0 :        mw_fac(m) = d0 * mw_dry * sqrt(1._r8/mw_dry + 1._r8/cnst_mw(m)) / n_avog
          96             :     end do
          97             :
          98             :     !--------------------------------------------------------------------------------------------
          99             :     ! For WACCM-X, get H data index and initialize thermal diffusion coefficient
         100             :     !--------------------------------------------------------------------------------------------
         101           0 :     if ( waccmx_mode ) then
         102             :
         103           0 :       call cnst_get_ind('H',  indx_H)
         104             :
         105           0 :       allocate(alphath(ncnst), stat=ierr)
         106           0 :       if ( ierr /= 0 ) then
         107           0 :          errstring = 'allocate failed in init_molec_diff'
         108           0 :          return
         109             :       end if
         110           0 :       alphath(:ncnst) = 0._r8
         111           0 :       alphath(indx_H) = -0.38_r8
         112             :
         113             :     endif
         114             :
         115           0 :   end subroutine init_molec_diff
         116             :
         117             :   !============================================================================ !
         118             :   !                                                                             !
         119             :   !============================================================================ !
         120             :
         121           0 :   subroutine compute_molec_diff(lchnk, pcols, pver, ncnst, ncol, &
         122           0 :        kvm, kvt, tint, rhoi, kq_scal, cnst_mw, &
         123           0 :        mw_fac_out, nbot_molec)
         124             :
         125             :     use cam_thermo,       only: kmvis, kmcnd
         126             :     use air_composition,  only: cpairv
         127             :
         128             :     ! --------------------- !
         129             :     ! Input-Output Argument !
         130             :     ! --------------------- !
         131             :
         132             :     integer,  intent(in)    :: pcols
         133             :     integer,  intent(in)    :: pver
         134             :     integer,  intent(in)    :: ncnst
         135             :     integer,  intent(in)    :: ncol                      ! Number of atmospheric columns
         136             :     integer,  intent(in)    :: lchnk                     ! Chunk identifier
         137             :
         138             :     real(r8), intent(inout) :: kvm(pcols,pver+1)         ! Viscosity ( diffusivity for momentum )
         139             :     real(r8), intent(out)   :: kvt(pcols,pver+1)         ! Kinematic molecular conductivity
         140             :     real(r8), intent(in)    :: tint(pcols,pver+1)        ! Interface temperature [ K ]
         141             :     real(r8), intent(in)    :: rhoi(pcols,pver+1)        ! Density ( rho ) at interfaces
         142             :     real(r8), intent(in)    :: cnst_mw(ncnst)            ! Constituent molecular weight
         143             :
         144             :     real(r8), intent(out)   :: kq_scal(pcols,pver+1)     ! kq_fac*sqrt(T)*m_d/rho for molecular diffusivity
         145             :     real(r8), intent(out)   :: mw_fac_out(pcols,pver+1,ncnst) ! composition dependent mw_fac on interface level
         146             :     integer,  intent(in)    :: nbot_molec
         147             :
         148             :     ! --------------- !
         149             :     ! Local variables !
         150             :     ! --------------- !
         151             :
         152             :     integer                 :: m                          ! Constituent index
         153             :     integer                 :: k                          ! Level index
         154             :
         155           0 :     real(r8)                :: mbarvi(pcols,nbot_molec,ncnst) ! mbarv on interface level
         156             :
         157           0 :     real(r8)                :: mkvisc(ncol)               ! Molecular kinematic viscosity c*tint**(2/3)/rho
         158             :
         159             :     ! ----------------------- !
         160             :     ! Main Computation Begins !
         161             :     ! ----------------------- !
         162             :
         163             :     !
         164             :     ! Need variable mw_fac for kvt and constant otherwise.
         165             :     !
         166             :     ! Then compute molecular kinematic viscosity, heat diffusivity and
         167             :     ! factor for constituent diffusivity.
         168             :     ! This is a key part of the code. For WACCM-X, use constituent
         169             :     ! dependent molecular viscosity and conductivity.
         170             :     !
         171             :
         172           0 :     kvt     = 0._r8
         173           0 :     kq_scal = 0._r8
         174           0 :     if (waccmx_mode) then
         175           0 :       do m = 1, ncnst
         176           0 :          mbarvi(:ncol,1,m) = .75_r8*mbarv(:ncol,1,lchnk)+0.5_r8*mbarv(:ncol,2,lchnk) &
         177           0 :             -.25_r8*mbarv(:ncol,3,lchnk)
         178           0 :         do k = 2, nbot_molec
         179           0 :            mbarvi(:ncol,k,m) = 0.5_r8 * (mbarv(:ncol,k-1,lchnk)+mbarv(:ncol,k,lchnk))
         180           0 :            mw_fac_out(:ncol,k,m) = d0 * mbarvi(:ncol,k,m) * &
         181           0 :                 sqrt(1._r8/mbarvi(:ncol,k,m) + 1._r8/cnst_mw(m)) / n_avog
         182             :         enddo
         183           0 :         mw_fac_out(:ncol,1,m) = 1.5_r8*mw_fac_out(:ncol,2,m)-.5_r8*mw_fac_out(:ncol,3,m)
         184           0 :         do k = nbot_molec+1, pver+1
         185           0 :           mw_fac_out(:ncol,k,m) = mw_fac_out(:ncol,nbot_molec,m)
         186             :         enddo
         187             :       end do
         188             :
         189           0 :       do k = 1, nbot_molec
         190           0 :          mkvisc  = kmvis(:ncol,k,lchnk) / rhoi(:ncol,k)
         191           0 :          kvm(:ncol,k) = kvm(:ncol,k) + mkvisc
         192           0 :          mkvisc  = kmcnd(:ncol,k,lchnk) / rhoi(:ncol,k)
         193           0 :          kvt(:ncol,k) = mkvisc
         194           0 :          kq_scal(:ncol,k) = sqrt(tint(:ncol,k)) / rhoi(:ncol,k)
         195             :       end do
         196             :
         197             :     else
         198           0 :       do m = 1, ncnst
         199           0 :          mw_fac_out(:,:,m) = mw_fac(m)
         200             :       end do
         201             :
         202           0 :       do k = 1, nbot_molec
         203           0 :          mkvisc   = km_fac * tint(:ncol,k)**pwr / rhoi(:ncol,k)
         204           0 :          kvm(:ncol,k) = kvm(:ncol,k) + mkvisc
         205           0 :          kvt(:ncol,k) = mkvisc * cpairv(:ncol,k,lchnk)
         206           0 :          kq_scal(:ncol,k) = sqrt(tint(:ncol,k)) / rhoi(:ncol,k)
         207             :       end do
         208             :     endif
         209             :
         210           0 :   end subroutine compute_molec_diff
         211             :
         212             :   !============================================================================ !
         213             :   !                                                                             !
         214             :   !============================================================================ !
         215             :
         216           0 :   function vd_lu_qdecomp( &
         217             :        pcols , pver   , ncol       , fixed_ubc  , mw     , &
         218           0 :        kv    , kq_scal, mw_facm    , dpidz_sq   , p      , &
         219             :        interface_boundary, molec_boundary, &
         220           0 :        tint  , ztodt  , nbot_molec , &
         221           0 :        lchnk , t          , m      , no_molec_decomp)      result(decomp)
         222             :
         223             :     use coords_1d,           only: Coords1D
         224             :     use linear_1d_operators, only: BoundaryType, TriDiagDecomp
         225             :     use vdiff_lu_solver,     only: fin_vol_lu_decomp
         226             :
         227             :     !------------------------------------------------------------------------------ !
         228             :     ! Add the molecular diffusivity to the turbulent diffusivity for a consitutent. !
         229             :     ! Update the superdiagonal (ca(k)), diagonal (cb(k)) and subdiagonal (cc(k))    !
         230             :     ! coefficients of the tridiagonal diffusion matrix, also ze and denominator.    !
         231             :     !------------------------------------------------------------------------------ !
         232             :
         233             :     ! ---------------------- !
         234             :     ! Input-Output Arguments !
         235             :     ! ---------------------- !
         236             :
         237             :     integer,  intent(in)    :: pcols
         238             :     integer,  intent(in)    :: pver
         239             :     integer,  intent(in)    :: ncol                  ! Number of atmospheric columns
         240             :
         241             :     integer,  intent(in)    :: nbot_molec
         242             :
         243             :     logical,  intent(in)    :: fixed_ubc             ! Fixed upper boundary condition flag
         244             :     real(r8), intent(in)    :: kv(pcols,pver+1)      ! Eddy diffusivity
         245             :     real(r8), intent(in)    :: kq_scal(pcols,pver+1) ! Molecular diffusivity ( kq_fac*sqrt(T)*m_d/rho )
         246             :     real(r8), intent(in)    :: mw                    ! Molecular weight for this constituent
         247             :     real(r8), intent(in)    :: mw_facm(pcols,pver+1) ! composition dependent sqrt(1/M_q + 1/M_d) for this constituent
         248             :     real(r8), intent(in)    :: dpidz_sq(ncol,pver+1) ! (g*rho)**2 (square of vertical derivative of pint)
         249             :     type(Coords1D), intent(in) :: p                  ! Pressure coordinates
         250             :     type(BoundaryType), intent(in) :: interface_boundary ! Boundary on grid edge.
         251             :     type(BoundaryType), intent(in) :: molec_boundary ! Boundary at edge of molec_diff region.
         252             :     real(r8), intent(in)    :: tint(pcols,pver+1)    ! Interface temperature [ K ]
         253             :     real(r8), intent(in)    :: ztodt                 ! 2 delta-t [ s ]
         254             :
         255             :     integer,  intent(in)    :: lchnk                ! Chunk number
         256             :     real(r8), intent(in)    :: t(pcols,pver)        ! temperature
         257             :     integer,  intent(in)    :: m                    ! cnst index
         258             :
         259             :     ! Decomposition covering levels without vertical diffusion.
         260             :     type(TriDiagDecomp), intent(in) :: no_molec_decomp
         261             :
         262             :     ! LU decomposition information for solver.
         263             :     type(TriDiagDecomp) :: decomp
         264             :
         265             :     ! --------------- !
         266             :     ! Local Variables !
         267             :     ! --------------- !
         268             :
         269             :     ! Level index.
         270             :     integer :: k
         271             :
         272             :     ! Molecular diffusivity for constituent.
         273           0 :     real(r8)                :: kmq(ncol,nbot_molec+1)
         274             :
         275             :     ! Term for drift due to molecular separation: (m_i/m - 1) / p
         276           0 :     real(r8)                :: mw_term(ncol,nbot_molec+1)
         277             :
         278             :     ! Diffusion coefficient.
         279           0 :     real(r8)                :: diff_coef(ncol,nbot_molec+1)
         280             :     ! Advection velocity.
         281           0 :     real(r8)                :: advect_v(ncol,nbot_molec+1)
         282             :
         283             :     ! 1/mbar * d(mbar)/dp
         284           0 :     real(r8)                :: gradm(ncol,nbot_molec+1)
         285             :
         286             :     ! alphaTh/T * dT/dp, for now alphaTh is non-zero only for H.
         287           0 :     real(r8)                :: gradt(ncol,nbot_molec+1)
         288             :
         289             :     ! mbarv at interface
         290           0 :     real(r8)                :: mbarvi(ncol)
         291             :
         292             :     ! ----------------------- !
         293             :     ! Main Computation Begins !
         294             :     ! ----------------------- !
         295             :
         296             :     ! --------------------------------------------------------------------- !
         297             :     ! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the !
         298             :     ! tridiagonal diffusion matrix. The diagonal elements  (cb=1+ca+cc) are !
         299             :     ! a combination of ca and cc; they are not required by the solver.      !
         300             :     !---------------------------------------------------------------------- !
         301             :
         302           0 :     call t_startf('vd_lu_qdecomp')
         303             :
         304           0 :     kmq  = 0._r8
         305           0 :     mw_term = 0._r8
         306           0 :     gradm = 0._r8
         307           0 :     gradt = 0._r8
         308             :
         309             :     ! Compute difference between scale heights of constituent and dry air
         310             :
         311           0 :     if ( waccmx_mode ) then
         312             :
         313             :        ! Top level first.
         314           0 :        k = 1
         315           0 :        mbarvi = .75_r8*mbarv(:ncol,k,lchnk)+0.5_r8*mbarv(:ncol,k+1,lchnk) &
         316           0 :             -.25_r8*mbarv(:ncol,k+2,lchnk)
         317           0 :        mw_term(:,k) = (mw/mbarvi - 1._r8) / p%ifc(:,k)
         318           0 :        gradm(:,k) = (mbarv(:ncol,k,lchnk)-mbarvi)/ &
         319           0 :             (p%mid(:,k)-p%ifc(:,k))/ &
         320           0 :             (mbarv(:ncol,k,lchnk)+mbarvi)*2._r8
         321             :
         322           0 :        if (alphath(m) /= 0._r8) then
         323           0 :           gradt(:,k) = alphath(m)*(t(:ncol,k)-tint(:ncol,k))/ &
         324           0 :                (p%mid(:ncol,k)-p%ifc(:ncol,k))/ &
         325           0 :                (t(:ncol,k)+tint(:ncol,k))*2._r8
         326             :        end if
         327             :
         328             :        ! Interior of molecular diffusion region.
         329           0 :        do k = 2, nbot_molec
         330           0 :           mbarvi = 0.5_r8 * (mbarv(:ncol,k-1,lchnk)+mbarv(:ncol,k,lchnk))
         331           0 :           mw_term(:,k) = (mw/mbarvi - 1._r8) / p%ifc(:,k)
         332           0 :           gradm(:,k) = (mbarv(:ncol,k,lchnk)-mbarv(:ncol,k-1,lchnk)) * &
         333           0 :                p%rdst(:,k-1)/mbarvi
         334             :        enddo
         335             :
         336           0 :        if (alphath(m) /= 0._r8) then
         337           0 :           do k = 2, nbot_molec
         338           0 :              gradt(:,k) = alphath(m)*(t(:ncol,k)-t(:ncol,k-1)) &
         339           0 :                   *p%rdst(:,k-1)/tint(:ncol,k)
         340             :           end do
         341             :        end if
         342             :
         343             :        ! Leave nbot_molec+1 terms as zero, because molecular diffusion is
         344             :        ! small at the lower boundary.
         345             :
         346             :     else
         347             :
         348           0 :        do k = 1, nbot_molec
         349           0 :           mw_term(:,k) = (mw/mw_dry - 1._r8) / p%ifc(:ncol,k)
         350             :        enddo
         351             :
         352             :     endif
         353             :
         354             :     !-------------------- !
         355             :     ! Molecular diffusion !
         356             :     !-------------------- !
         357             :
         358             :     ! Start with non-molecular portion of diffusion.
         359             :
         360             :     ! Molecular diffusion coefficient.
         361           0 :     do k = 1, nbot_molec
         362           0 :        kmq(:,k)  = kq_scal(:ncol,k) * mw_facm(:ncol,k)
         363             :     end do
         364             :
         365           0 :     diff_coef = kv(:ncol,:nbot_molec+1) + kmq
         366             :
         367             :     ! "Drift" terms.
         368           0 :     advect_v = kmq*mw_term
         369           0 :     if ( waccmx_mode ) then
         370             :        advect_v = advect_v - kmq*gradt - &
         371           0 :             (kv(:ncol,:nbot_molec+1) + kmq)*gradm
         372             :     end if
         373             :
         374             :     ! Convert from z to pressure representation.
         375           0 :     diff_coef = dpidz_sq(:,:nbot_molec+1) * diff_coef
         376           0 :     advect_v = dpidz_sq(:,:nbot_molec+1) * advect_v
         377             :
         378           0 :     if( fixed_ubc ) then
         379             :        decomp = fin_vol_lu_decomp(ztodt, p, &
         380             :             coef_q_diff=diff_coef, coef_q_adv=advect_v, &
         381             :             upper_bndry=interface_boundary, &
         382             :             lower_bndry=molec_boundary, &
         383           0 :             graft_decomp=no_molec_decomp)
         384             :     else
         385             :        decomp = fin_vol_lu_decomp(ztodt, p, &
         386             :             coef_q_diff=diff_coef, coef_q_adv=advect_v, &
         387             :             lower_bndry=molec_boundary, &
         388           0 :             graft_decomp=no_molec_decomp)
         389             :     end if
         390             :
         391           0 :     call t_stopf('vd_lu_qdecomp')
         392             :
         393           0 :   end function vd_lu_qdecomp
         394             :
         395             : end module molec_diff

Generated by: LCOV version 1.14