LCOV - code coverage report
Current view: | top level - chemistry/bulk_aero - mo_aerosols.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 26 | 97 | 26.8 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 1 | 2 | 50.0 %

              Line data    Source code
           1             : module mo_aerosols
           2             :   !-----------------------------------------------------------------
           3             :   !
           4             :   ! this module computes the production of ammonium nitrate
           5             :   ! using the formulation by Seinfeld and Pandis (p531, 1998)
           6             :   ! with the simplification of activity coefficients and
           7             :   ! aerosol molality using the parameterizations
           8             :   ! from Metzger et al. (JGR, ACH-16, 107(D16), 2002)
           9             :   !
          10             :   ! written by Jean-Francois Lamarque (April 2004)
          11             :   ! adapted for CAM (May 2004)
          12             :   !
          13             :   !-----------------------------------------------------------------
          14             :
          15             :   use shr_kind_mod, only : r8 => shr_kind_r8
          16             :   use ppgrid,       only : pver
          17             :   use cam_logfile,  only: iulog
          18             :
          19             :   private
          20             :   public :: aerosols_inti,aerosols_formation
          21             :   public :: has_aerosols
          22             :
          23             :   save
          24             :
          25             :   integer, target  :: spc_ndx(5)
          26             :   integer, pointer :: nh3_ndx, nh4no3_ndx, nh4_ndx
          27             :   integer, pointer :: so4_ndx, hno3_ndx
          28             :   integer  :: xhno3_ndx, xnh4no3_ndx
          29             :   integer  :: nu_i(2)
          30             :   real(r8)     :: zeta_inv
          31             :   real(r8)     :: z_i(2)
          32             :   logical  :: has_aerosols = .true.
          33             :
          34             : contains
          35             :
          36        1536 :   subroutine aerosols_inti()
          37             :
          38             :     use mo_chem_utls, only : get_spc_ndx
          39             :     use cam_history,  only : addfld
          40             :     use spmd_utils,   only : masterproc
          41             :
          42             :     implicit none
          43             :
          44             :     !-----------------------------------------------------------------
          45             :     !   ... local variables
          46             :     !-----------------------------------------------------------------
          47             :     integer :: m
          48             :
          49        1536 :     nh3_ndx    => spc_ndx(1)
          50        1536 :     nh4no3_ndx => spc_ndx(2)
          51        1536 :     so4_ndx    => spc_ndx(3)
          52        1536 :     hno3_ndx   => spc_ndx(4)
          53        1536 :     nh4_ndx    => spc_ndx(5)
          54             :
          55             :     !-----------------------------------------------------------------
          56             :     !   ... set species index
          57             :     !-----------------------------------------------------------------
          58        1536 :     nh3_ndx    = get_spc_ndx( 'NH3'    )
          59        1536 :     nh4no3_ndx = get_spc_ndx( 'NH4NO3' )
          60        1536 :     so4_ndx    = get_spc_ndx( 'SO4'    )
          61        1536 :     hno3_ndx   = get_spc_ndx( 'HNO3'   )
          62        1536 :     nh4_ndx    = get_spc_ndx( 'NH4'    )
          63        1536 :     xnh4no3_ndx = get_spc_ndx( 'XNH4NO3' )
          64        1536 :     xhno3_ndx   = get_spc_ndx( 'XHNO3'   )
          65             :
          66        1536 :     has_aerosols = all( spc_ndx(:) > 0 )
          67        1536 :     if( .not. has_aerosols ) then
          68        1536 :        if (masterproc) then
          69           2 :           write(iulog,*) '-----------------------------------------'
          70           2 :           write(iulog,*) 'mozart will NOT do nh4no3'
          71           2 :           write(iulog,*) 'following species are missing'
          72          12 :           do m = 1,size(spc_ndx)
          73          12 :              if( spc_ndx(m) < 1 ) then
          74          10 :                 write(iulog,*) m
          75             :              end if
          76             :           end do
          77           2 :           write(iulog,*) '-----------------------------------------'
          78             :        endif
          79        1536 :        return
          80             :     else
          81           0 :        if (masterproc) then
          82           0 :           write(iulog,*) '-----------------------------------------'
          83           0 :           write(iulog,*) 'mozart will do nh4no3'
          84           0 :           write(iulog,*) '-----------------------------------------'
          85             :        end if
          86             :     end if
          87             :
          88             :     !
          89             :     ! define parameters
          90             :     !
          91             :     ! ammonium nitrate (NH4NO3)
          92             :     !
          93           0 :     zeta_inv = 1._r8/4._r8
          94           0 :     nu_i(1)  = 4
          95           0 :     z_i (1)  = 1._r8
          96             :     !
          97             :     ! ammonium sulfate
          98             :     !
          99           0 :     nu_i(2)  = 4
         100           0 :     z_i (2)  = 0.5_r8
         101             :
         102           0 :     call addfld ('TSO4_VMR', (/ 'lev' /), 'A', 'mol/mol','total sulfate in mo_aerosols')
         103           0 :     call addfld ('THNO3_VMR',(/ 'lev' /), 'A','mol/mol','total nitric acid in mo_aerosols')
         104             :
         105           0 :     return
         106        1536 :   end subroutine aerosols_inti
         107             :
         108           0 :   subroutine aerosols_formation( ncol, lchnk, tfld, rh, qin)
         109             :
         110        1536 :     use ppgrid, only        : pcols, pver
         111             :     use chem_mods,  only    : gas_pcnst, adv_mass
         112             :     use cam_history, only   : outfld
         113             :
         114             :     implicit none
         115             :     !
         116             :     ! input arguments
         117             :     !
         118             :     !
         119             :     ! input arguments
         120             :     !
         121             :     integer,  intent(in)    :: ncol        ! number columns in chunk
         122             :     integer,  intent(in)    :: lchnk       ! chunk index
         123             :     real(r8), intent(in)    :: tfld(:,:)   ! temperature
         124             :     real(r8), intent(in)    :: rh(:,:)     ! relative humidity
         125             :     real(r8), intent(inout) :: qin(:,:,:)  ! xported species ( vmr )
         126             :
         127             :     !
         128             :     ! local variables
         129             :     !
         130             :     integer :: i,j,k,n
         131             :     integer                           :: domain_number   ! concentration domain
         132             :     real(r8)                          :: sulfate_state   ! fraction of sulfate neutralized by ammonia
         133           0 :     real(r8), dimension(ncol,pver)    :: tso4, &         ! total sulfate
         134           0 :                                          thno3,&         ! total nitric acid
         135           0 :                                          txhno3          ! total nitric acid ( XHNO3 )
         136             :     real(r8)                          :: tnh3, &         ! total ammonia
         137             :                                          fnh3, &         ! free ammonia
         138             :                                          rhd,  &         ! relative humidity of deliquescence
         139             :                                          gamma, &        ! activity coefficient
         140             :                                          ssm_nh4no3, &   ! single solute molality for NH4NO3
         141             :                                          ta, &           ! total ammonia
         142             :                                          tn, &           ! total nitrate
         143             :                                          kp, &           ! equilibrium constant
         144             :                                          nh4no3          ! ammonium nitrate produced
         145             :     real(r8) :: log_t
         146             :     real(r8) :: ti
         147             :     real(r8) :: xnh4no3
         148             :
         149           0 :     do k=1,pver
         150           0 :        do i=1,ncol
         151             :
         152             :           !
         153             :           ! compute total concentrations
         154             :           !
         155           0 :           tnh3      = (qin(i,k,  nh3_ndx)+qin(i,k,nh4_ndx))
         156           0 :           tso4(i,k) = qin(i,k,so4_ndx)
         157             :           !
         158             :           ! define concentration domain
         159             :           !
         160           0 :           if ( tnh3 < tso4(i,k) ) then
         161             :              domain_number = 4
         162             :              sulfate_state = 1.0_r8
         163           0 :           elseif ( tnh3 < 2._r8*tso4(i,k) ) then
         164             :              domain_number = 3
         165             :              sulfate_state = 1.5_r8
         166             :           else
         167           0 :              domain_number = 2
         168           0 :              sulfate_state = 2.0_r8
         169             :           endif
         170             :           !
         171             :           ! define free ammonia (ammonia available for ammonium nitrate production)
         172             :           !
         173           0 :           fnh3 = tnh3 - sulfate_state * tso4(i,k)
         174           0 :           fnh3 = max(0._r8,fnh3)
         175             :           !
         176             :           ! convert initial concentrations to ppbv
         177             :           !
         178           0 :           tso4(i,k)  = tso4(i,k) * 1.e9_r8
         179           0 :           tnh3       = tnh3 * 1.e9_r8
         180           0 :           fnh3       = fnh3 * 1.e9_r8
         181           0 :           thno3(i,k) = (qin(i,k,hno3_ndx)+qin(i,k,nh4no3_ndx)) * 1.e9_r8
         182           0 :           if ( xhno3_ndx > 0 .and. xnh4no3_ndx > 0 ) then
         183           0 :             txhno3(i,k) = (qin(i,k,xhno3_ndx)+qin(i,k,xnh4no3_ndx)) * 1.e9_r8
         184             :           endif
         185             :           !
         186             :           ! compute relative humidity of deliquescence (%) for NH4NO3
         187             :           ! (Seinfeld and Pandis, p532)
         188             :           !
         189           0 :           ti    = 1._r8/tfld(i,k)
         190           0 :           rhd   = 0.01_r8 * exp( 1.6954_r8 + 723.7_r8*ti )
         191           0 :           log_t = log( tfld(i,k)/298._r8 )
         192           0 :           if ( rh(i,k) < rhd ) then
         193             :              !
         194             :              ! crystalline ammonium nitrate
         195             :              !
         196             :              ! compute equilibrium constant
         197             :              !
         198           0 :              kp = exp( 84.6_r8 - 24220._r8*ti - 6.1_r8*log_t )
         199             :              !
         200             :           else
         201             :              !
         202             :              ! aqueous phase ammonium nitrate
         203             :              !
         204             :              ! compute activity coefficients (from Menzger et al.)
         205             :              !
         206           0 :              n = domain_number
         207           0 :              gamma = (rh(i,k)**n/(1000._r8/n*(1._r8-rh(i,k))+n))**zeta_inv
         208             :              !
         209             :              ! compute single solute molality for NH4NO3
         210             :              !
         211           0 :              ssm_nh4no3 = (1000._r8 * 0.81_r8 * nu_i(1) * (1._r8/rh(i,k)-1._r8)/80._r8)**z_i(1)
         212             :              !
         213             :              ! compute equilibrium constant
         214             :              !
         215           0 :              kp = (gamma*ssm_nh4no3)**2 * exp( 53.19_r8 - 15850.62_r8*ti + 11.51_r8*log_t )
         216             :
         217             :           endif
         218             :           !
         219             :           ! calculate production of NH4NO3 (in ppbv) using Seinfeld and Pandis (p534, 1998)
         220             :           !
         221           0 :           ta = fnh3
         222           0 :           tn = thno3(i,k)
         223           0 :           nh4no3 = 0.5_r8 * (ta + tn - sqrt(max(0._r8,(ta+tn)**2 - 4._r8*(ta*tn-kp))))
         224           0 :           nh4no3 = max(0._r8,nh4no3)
         225           0 :           if ( xhno3_ndx > 0 .and. xnh4no3_ndx > 0 ) then
         226           0 :              tn = txhno3(i,k)
         227           0 :              xnh4no3 = 0.5_r8 * (ta + tn - sqrt(max(0._r8,(ta+tn)**2 - 4._r8*(ta*tn-kp))))
         228           0 :              xnh4no3 = max(0._r8,xnh4no3)
         229             :           endif
         230             :           !
         231             :           ! reset concentrations according to equilibrium calculation
         232             :           !
         233           0 :           qin(i,k,nh4no3_ndx)  = nh4no3
         234           0 :           if ( xhno3_ndx > 0 ) then
         235           0 :              qin(i,k,xnh4no3_ndx)  = xnh4no3
         236             :           endif
         237           0 :           qin(i,k,nh3_ndx   )  = max(0._r8,(fnh3-nh4no3))
         238           0 :           qin(i,k,nh4_ndx   )  = max(0._r8,(tnh3-(fnh3-nh4no3)))
         239           0 :           qin(i,k,hno3_ndx  )  = max(0._r8,(thno3(i,k)-nh4no3))
         240           0 :           if ( xhno3_ndx > 0 ) then
         241           0 :              qin(i,k,xhno3_ndx  )  = max(0._r8,(txhno3(i,k)-xnh4no3))
         242             :           endif
         243           0 :           qin(i,k,so4_ndx   )  = tso4(i,k)
         244             :           !
         245             :           ! convert from ppbv to vmr
         246             :           !
         247           0 :           qin(i,k,nh4no3_ndx)  = qin(i,k,nh4no3_ndx) * 1.e-9_r8
         248           0 :           qin(i,k,nh3_ndx   )  = qin(i,k,nh3_ndx   ) * 1.e-9_r8
         249           0 :           qin(i,k,nh4_ndx   )  = qin(i,k,nh4_ndx   ) * 1.e-9_r8
         250           0 :           qin(i,k,hno3_ndx  )  = qin(i,k,hno3_ndx  ) * 1.e-9_r8
         251           0 :           qin(i,k,so4_ndx   )  = qin(i,k,so4_ndx   ) * 1.e-9_r8
         252           0 :           if ( xhno3_ndx > 0 ) then
         253           0 :              qin(i,k,xnh4no3_ndx)  = qin(i,k,xnh4no3_ndx) * 1.e-9_r8
         254             :           endif
         255           0 :           if ( xhno3_ndx > 0 ) then
         256           0 :              qin(i,k,xhno3_ndx  )  = qin(i,k,xhno3_ndx  ) * 1.e-9_r8
         257             :           endif
         258             :
         259             :        end do
         260             :     end do
         261             :     !
         262             :     ! outputs
         263             :     !
         264           0 :     call outfld ('TSO4_VMR' ,tso4 (:ncol,:), ncol, lchnk )
         265           0 :     call outfld ('THNO3_VMR',thno3(:ncol,:), ncol, lchnk )
         266             :
         267           0 :     return
         268           0 :   end subroutine aerosols_formation
         269             :
         270             : end module mo_aerosols

Generated by: LCOV version 1.14