LCOV - code coverage report
Current view: | top level - chemistry/bulk_aero - mo_setsoa.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 82 | 351 | 23.4 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 3 | 7 | 42.9 %

              Line data    Source code
           1             : module mo_setsoa
           2             :
           3             :   use shr_kind_mod,     only : r8 => shr_kind_r8
           4             :   use cam_logfile,      only : iulog
           5             :   use ppgrid,           only : pcols, pver, begchunk, endchunk
           6             :   use cam_abortutils,   only : endrun
           7             :   use mo_constants,     only : avogadro, Rgas
           8             :
           9             :   implicit none
          10             :   private
          11             :   public :: soa_inti, setsoa, has_soa
          12             :   public :: soa_register
          13             :
          14             :   save
          15             :
          16             :   integer, parameter :: NRX = 10                  ! number of SOA forming reactions
          17             :   integer, parameter :: NPR = 2                   ! number of products (always 2!)
          18             :   integer, target  :: spc_ndx(22)
          19             :   integer, pointer :: soam_ndx, soai_ndx, soab_ndx, soat_ndx, soax_ndx
          20             :   integer, pointer :: sogm_ndx, sogi_ndx, sogb_ndx, sogt_ndx, sogx_ndx
          21             :   integer, pointer :: oc1_ndx, oc2_ndx, c10h16_ndx, isop_ndx, o3_ndx, oh_ndx
          22             :   integer, pointer :: no3_ndx, no_ndx, ho2_ndx, tolo2_ndx, beno2_ndx, xylo2_ndx
          23             :
          24             :   integer :: rxn_soa(nrx),react_ndx(NRX,NPR)
          25             :   real(r8), dimension(NRX,NPR) :: alpha         ! mass-based stoichiometric coefficients
          26             :   real(r8), dimension(NRX,NPR) :: k_om          ! equilibrium gas-particule partition
          27             :   real(r8), dimension(NRX)     :: T1, delH      ! Clausium Clayperson parameters (K, J/mol)
          28             :   real(r8), dimension(nrx,npr) :: fracsog_init ! mass fraction of each SOA class from each reaction
          29             :   real(r8), dimension(nrx,npr) :: fracsoa_init  ! mass fraction of each SOG class from each reaction
          30             :
          31             :   integer :: fracsog_ndx = -1
          32             :   integer :: fracsoa_ndx = -1
          33             :
          34             :   integer, pointer :: soa_ndx
          35             :   integer, pointer ::  bigalk_ndx, toluene_ndx
          36             :
          37             :   real(r8), dimension(6)   :: bulk_yield               ! total yield of condensable compound (ug/m3/ppm)
          38             :   real(r8), dimension(6)   :: fraction                 ! fraction of VOC used in reaction
          39             :
          40             :   real(r8), parameter :: OMscale  = 2.1_r8             ! scaling factor for OM:OC [Turpin and Lim, 2001]
          41             :   logical :: has_soa = .false.
          42             :   logical :: has_soa_equil = .false.
          43             :
          44             : contains
          45             :
          46             : !===============================================================================
          47             : !===============================================================================
          48        1536 : subroutine soa_inti(pbuf2d)
          49             :   use physics_buffer, only : physics_buffer_desc
          50             :   type(physics_buffer_desc), pointer :: pbuf2d(:,:)
          51        1536 :   if ( has_soa_equil) then
          52           0 :      call soa_inti_equil(pbuf2d)
          53             :   else
          54        1536 :      call soa_inti_old()
          55             :   endif
          56        1536 : endsubroutine soa_inti
          57             : !===============================================================================
          58             : !===============================================================================
          59           0 : subroutine setsoa( ncol, lchnk, dt, reaction_rates, tfld, xhnm, vmr, pbuf)
          60        1536 :   use physics_buffer, only : physics_buffer_desc
          61             :   use chem_mods,      only : gas_pcnst, rxntot
          62             :   !-----------------------------------------------------------------------
          63             :   !      ... dummy arguments
          64             :   !-----------------------------------------------------------------------
          65             :   integer, intent(in)      :: ncol                   ! number columns in chunkx
          66             :   integer, intent(in)      :: lchnk                  ! chunk index
          67             :   real(r8), intent(in)     :: dt                     ! time step
          68             :   real(r8), intent(in)     :: reaction_rates(:,:,:)  ! reaction rates
          69             :   real(r8), intent(in)     :: tfld(:,:)              ! temperature (K)
          70             :   real(r8), intent(in)     :: xhnm(:,:)              ! total atms density (molec/cm**3)
          71             :   real(r8), intent(inout)  :: vmr(:,:,:)             ! xported species ( vmr )
          72             :   type(physics_buffer_desc), pointer :: pbuf(:)
          73             :
          74           0 :   if ( has_soa_equil) then
          75           0 :      call setsoa_equil(dt,reaction_rates,tfld,vmr,xhnm,ncol,lchnk,pbuf)
          76             :   else
          77           0 :      call setsoa_old(dt,reaction_rates,tfld,vmr,xhnm,ncol,lchnk)
          78             :   endif
          79             :
          80           0 : end subroutine setsoa
          81             : !===============================================================================
          82             : !===============================================================================
          83        1536 : subroutine soa_inti_old
          84             :
          85           0 :   use mo_chem_utls, only : get_spc_ndx, get_rxt_ndx
          86             :   use cam_history,  only : addfld
          87             :   use ppgrid,       only : pver
          88             :   use spmd_utils,   only : masterproc
          89             :
          90             :   implicit none
          91             :
          92             : !-----------------------------------------------------------------------
          93             : !       ... check if this is an aerosol simulation
          94             : !-----------------------------------------------------------------------
          95        1536 :   if( .not. has_soa ) then
          96             :     return
          97             :   end if
          98             :
          99             :   !-----------------------------------------------------------------------
         100             :   !     ... set reaction indicies
         101             :   !-----------------------------------------------------------------------
         102           0 :   rxn_soa(1) = get_rxt_ndx( 'soa1' )
         103           0 :   if ( rxn_soa(1) <= 0 ) then
         104           0 :      rxn_soa(1) = get_rxt_ndx( 'C10H16_O3' )
         105             :   end if
         106           0 :   rxn_soa(2) = get_rxt_ndx( 'soa2' )
         107           0 :   if ( rxn_soa(2) <= 0 ) then
         108           0 :      rxn_soa(2) = get_rxt_ndx( 'C10H16_OH' )
         109             :   end if
         110           0 :   rxn_soa(3) = get_rxt_ndx( 'soa3' )
         111           0 :   if ( rxn_soa(3) <= 0 ) then
         112           0 :      rxn_soa(3) = get_rxt_ndx( 'C10H16_NO3' )
         113             :   end if
         114           0 :   rxn_soa(4) = get_rxt_ndx( 'soa4' )
         115           0 :   if ( rxn_soa(4) <= 0 ) then
         116           0 :      rxn_soa(4) = get_rxt_ndx( 'TOLUENE_OH' )
         117             :   end if
         118           0 :   rxn_soa(5) = get_rxt_ndx( 'soa4' ) ! TOLUENE is a lumped species and there are two sets of pathways
         119           0 :   if ( rxn_soa(5) <= 0 ) then
         120           0 :      rxn_soa(5) = get_rxt_ndx( 'TOLUENE_OH' )
         121             :   end if
         122           0 :   rxn_soa(6) = get_rxt_ndx( 'soa5' )
         123           0 :   if ( rxn_soa(6) <= 0 ) then
         124           0 :      rxn_soa(6) = get_rxt_ndx( 'BIGALK_OH' )
         125             :   end if
         126           0 :   if( all( rxn_soa(:) < 1 ) ) then
         127           0 :      has_soa = .false.
         128           0 :      return
         129             :   else
         130           0 :      if (masterproc) then
         131           0 :         write(iulog,*) '-----------------------------------------'
         132           0 :         write(iulog,*) 'soa_inti_old: mozart will do soa aerosols'
         133           0 :         write(iulog,*) '-----------------------------------------'
         134             :      endif
         135             :   end if
         136             :
         137             : !
         138             : ! define reactants
         139             : !
         140           0 :   react_ndx(1,1) = c10h16_ndx
         141           0 :   react_ndx(1,2) = o3_ndx
         142           0 :   react_ndx(2,1) = c10h16_ndx
         143           0 :   react_ndx(2,2) = oh_ndx
         144           0 :   react_ndx(3,1) = c10h16_ndx
         145           0 :   react_ndx(3,2) = no3_ndx
         146           0 :   react_ndx(4,1) = toluene_ndx
         147           0 :   react_ndx(4,2) = oh_ndx
         148           0 :   react_ndx(5,1) = toluene_ndx
         149           0 :   react_ndx(5,2) = oh_ndx
         150           0 :   react_ndx(6,1) = bigalk_ndx
         151           0 :   react_ndx(6,2) = oh_ndx
         152             :
         153           0 :   if ( masterproc ) then
         154           0 :      write(iulog,*)'soa_inti ',c10h16_ndx, o3_ndx, oh_ndx, no3_ndx, bigalk_ndx, toluene_ndx
         155           0 :      write(iulog,*)'soa_inti ',soa_ndx, oc1_ndx, oc2_ndx
         156           0 :      write(iulog,*)'soa_inti ',react_ndx
         157             :   endif
         158             : !
         159             : ! define partitioning coefficients for each reaction
         160             : ! bulk yields are from Seinfeld and Pandis (1998)
         161             : !
         162             : ! c10h16 + o3 (from Chung and Seinfeld, JGR, 107, 2002)
         163             : !
         164           0 :   alpha(1,1)    = 0.067_r8
         165           0 :   alpha(1,2)    = 0.354_r8
         166           0 :   k_om (1,1)    = 0.184_r8
         167           0 :   k_om (1,2)    = 0.0043_r8
         168           0 :   fraction(1)   = 1._r8
         169           0 :   bulk_yield(1) = 762._r8
         170             : !
         171             : ! c10h16 + oh (from Chung and Seinfeld, JGR, 107, 2002)
         172             : !
         173           0 :   alpha(2,1)    = 0.067_r8
         174           0 :   alpha(2,2)    = 0.354_r8
         175           0 :   k_om (2,1)    = 0.184_r8
         176           0 :   k_om (2,2)    = 0.0043_r8
         177           0 :   fraction(2)   = 1._r8
         178           0 :   bulk_yield(2) = 762._r8
         179             : !
         180             : ! c10h16 + no3 (from Chung and Seinfeld, JGR, 107, 2002)
         181             : !
         182           0 :   alpha(3,1)    = 1.000_r8
         183           0 :   alpha(3,2)    = 0.000_r8
         184           0 :   k_om (3,1)    = 0.0163_r8
         185           0 :   k_om (3,2)    = 0.0000_r8
         186           0 :   fraction(3)   = 1._r8
         187           0 :   bulk_yield(3) = 762._r8
         188             : !
         189             : ! toluene + oh : toluene (from Odum et al., Environ. Sci. Technol., 1892, 1997)
         190             : !
         191           0 :   alpha(4,1)    = 0.038_r8
         192           0 :   alpha(4,2)    = 0.167_r8
         193           0 :   k_om (4,1)    = 0.042_r8
         194           0 :   k_om (4,2)    = 0.0014_r8
         195           0 :   fraction(4)   = 0.7_r8
         196           0 :   bulk_yield(4) = 424._r8
         197             : !
         198             : ! toluene + oh : m-xylene (from Cocker et al., Atmos. Environ., 6079, 2001)
         199             : !
         200           0 :   alpha(5,1)    = 0.120_r8
         201           0 :   alpha(5,2)    = 0.019_r8
         202           0 :   k_om (5,1)    = 0.060_r8
         203           0 :   k_om (5,2)    = 0.010_r8
         204           0 :   fraction(5)   = 0.2_r8
         205           0 :   bulk_yield(5) = 419._r8
         206             : !
         207             : ! bigalk + oh : only for alkanes >= heptane (assume low-yield aromatics as in Lack et al.)
         208             : !               (from Odum et al., Environ. Sci. Technol., 1892, 1997)
         209             : !
         210           0 :   alpha(6,1)    = 0.071_r8
         211           0 :   alpha(6,2)    = 0.138_r8
         212           0 :   k_om (6,1)    = 0.053_r8
         213           0 :   k_om (6,2)    = 0.0019_r8
         214           0 :   fraction(6)   = 0.1_r8
         215           0 :   bulk_yield(6) = 200._r8
         216             : !
         217           0 :   call addfld( 'SOA_PROD', (/ 'lev' /), 'A', 'kg/kg/s', 'production of SOA' )
         218             :
         219           0 :   return
         220        1536 : end subroutine soa_inti_old
         221             : !
         222           0 : subroutine setsoa_old(dt,reaction_rates,tfld,vmr,xhnm,ncol,lchnk)
         223             : !
         224             : ! secondary organic aerosol for mozart v2.5
         225             : !
         226             : ! based on Lack et al., JGR, 109, D03203, 2004
         227             : !
         228             : ! rewritten by Jean-Francois Lamarque for updated chemical
         229             : ! mechanism (March 2004)
         230             : !
         231             : ! adapted to CAM (May 2004)
         232             : !
         233        1536 :   use ppgrid,           only : pcols, pver
         234             :   use chem_mods,        only : adv_mass, gas_pcnst, rxntot
         235             :   use mo_chem_utls,     only : get_spc_ndx, get_rxt_ndx
         236             :   use cam_history,      only : outfld
         237             :   use cam_abortutils,   only : endrun
         238             : !
         239             :   implicit none
         240             : !
         241             : !-----------------------------------------------------------------------
         242             : !      ... dummy arguments
         243             : !-----------------------------------------------------------------------
         244             :   integer, intent(in)      :: ncol                   ! number columns in chunkx
         245             :   integer, intent(in)      :: lchnk                  ! chunk index
         246             :   real(r8), intent(in)     :: dt                     ! time step
         247             :   real(r8), intent(in)     :: reaction_rates(:,:,:)  ! reaction rates
         248             :   real(r8), intent(inout)  :: vmr(:,:,:)             ! xported species ( vmr )
         249             :   real(r8), intent(in)     :: tfld(:,:)              ! temperature (K)
         250             :   real(r8), intent(in)     :: xhnm(:,:)              ! total atms density (mol/cm**3)
         251             :
         252             : !-----------------------------------------------------------------------
         253             : !      ... local variables
         254             : !-----------------------------------------------------------------------
         255             :   integer :: i,k,n
         256             :   real(r8) :: m_0
         257             :   real(r8) :: mw_soa,yield,prod,soa_mass
         258           0 :   real(r8) :: soa_prod(ncol,pver)
         259             : !
         260             : ! find molecular weight of SOA
         261             : !
         262           0 :   mw_soa = adv_mass(soa_ndx)
         263             : !
         264             :   do k=1,pver
         265             :     do i=1,ncol
         266             : !
         267             : ! calculate initial mass of organic aerosols from OC1 and OC2
         268             : ! and convert to ug/m3
         269             : !
         270             :       m_0 = (vmr(i,k,oc1_ndx)+vmr(i,k,oc2_ndx)) * xhnm(i,k) * adv_mass(oc1_ndx)/avogadro * 1.e12_r8
         271             : !
         272             : ! switch based on a minimum value of m_0.  The bulk approach is
         273             : ! used to initiate the process
         274             : !
         275             :       if ( m_0 <= 0.2_r8 ) then
         276             : !
         277             : ! bulk theory
         278             : !
         279             :         soa_mass = 0._r8
         280             :         do n=1,6
         281             : !
         282             :           if ( rxn_soa(n) <= 0 ) cycle
         283             : !
         284             :           yield = bulk_yield(n)
         285             : !
         286             : ! define chemical production from gas-phase chemistry
         287             : !
         288             :           prod  = reaction_rates(i,k,rxn_soa(n)) * fraction(n) &
         289             :                 * vmr(i,k,react_ndx(n,1)) * vmr(i,k,react_ndx(n,2)) * dt
         290             : !
         291             : ! convert from mixing ratio to ppm
         292             : !
         293             :           prod = 1e6_r8 * prod
         294             : !
         295             : ! collect into total SOA mass
         296             : !
         297             :           soa_mass = soa_mass + yield * prod
         298             : !
         299             :         end do
         300             : !
         301             :       else
         302             : !
         303             : ! partitioning theory
         304             : !
         305             :         soa_mass = 0._r8
         306             :         do n=1,6
         307             : !
         308             :           if ( rxn_soa(n) <= 0 ) cycle
         309             : !
         310             : ! define yield from available m_0
         311             : !
         312             :           yield = soa_yield(m_0,alpha(n,1:2),k_om(n,1:2))
         313             : !
         314             : ! define chemical production from gas-phase chemistry
         315             : !
         316             :           prod  = reaction_rates(i,k,rxn_soa(n)) * fraction(n) &
         317             :                 * vmr(i,k,react_ndx(n,1)) * vmr(i,k,react_ndx(n,2)) * dt
         318             : !
         319             : ! convert from mixing ratio to mass (ug/m3)
         320             : !
         321             :           prod = prod * xhnm(i,k) * mw_soa/avogadro * 1.e12_r8
         322             : !
         323             : ! collect into total SOA mass
         324             : !
         325             :           soa_mass = soa_mass + yield * prod
         326             : !
         327             :         end do
         328             : !
         329             :       endif
         330             : !
         331             : ! convert from ug/m3 to mixing ratio and update vmr
         332             : !
         333             :       vmr(i,k,soa_ndx) = vmr(i,k,soa_ndx) + soa_mass * 1.e-12_r8 * avogadro/(mw_soa*xhnm(i,k))
         334             :       if ( vmr(i,k,soa_ndx) > 1.e0_r8 ) then
         335             :         write(iulog,*)i,k,soa_mass,m_0
         336             :         call endrun('soa_yield: vmr(i,k,soa_ndx) > 1.e0_r8')
         337             :       endif
         338             : !
         339             :       soa_prod(i,k) = soa_mass*1.e-12_r8*avogadro/(28.966_r8*xhnm(i,k)*dt)
         340             :     end do
         341             :   end do
         342             : !
         343             :   call outfld('SOA_PROD',soa_prod(:ncol,:),ncol, lchnk)
         344             :   return
         345             : end subroutine setsoa_old
         346             : !
         347             : real(r8) function soa_yield(m_0,xalpha,xk)
         348             : !
         349             :   implicit none
         350             : !
         351             :   real(r8), intent(in) :: m_0
         352             :   real(r8), intent(in), dimension(2) :: xalpha, xk
         353             : !
         354             :   soa_yield = m_0 * ( ((xalpha(1)*xk(1))/(1._r8+xk(1)*m_0)) &
         355             :                   +   ((xalpha(2)*xk(2))/(1._r8+xk(2)*m_0)) )
         356             : !
         357             :   return
         358           0 : end function soa_yield
         359             : !
         360             :
         361             :   !===============================================================================
         362             :   !===============================================================================
         363        1536 :   subroutine soa_register
         364             :     use physics_buffer, only : pbuf_add_field, dtype_r8
         365             :     use mo_chem_utls,   only : get_spc_ndx
         366             :
         367             :     !-----------------------------------------------------------------------
         368             :     !       ... set species indices
         369             :     !-----------------------------------------------------------------------
         370             :
         371        1536 :     oc1_ndx     => spc_ndx(1)
         372        1536 :     oc2_ndx     => spc_ndx(2)
         373        1536 :     soam_ndx    => spc_ndx(3)
         374        1536 :     soai_ndx    => spc_ndx(4)
         375        1536 :     soat_ndx    => spc_ndx(5)
         376        1536 :     soab_ndx    => spc_ndx(6)
         377        1536 :     soax_ndx    => spc_ndx(7)
         378        1536 :     c10h16_ndx  => spc_ndx(8)
         379        1536 :     isop_ndx    => spc_ndx(9)
         380        1536 :     tolo2_ndx   => spc_ndx(10)
         381        1536 :     beno2_ndx   => spc_ndx(11)
         382        1536 :     xylo2_ndx   => spc_ndx(12)
         383        1536 :     ho2_ndx     => spc_ndx(13)
         384        1536 :     no_ndx      => spc_ndx(14)
         385        1536 :     o3_ndx      => spc_ndx(15)
         386        1536 :     oh_ndx      => spc_ndx(16)
         387        1536 :     no3_ndx     => spc_ndx(17)
         388        1536 :     sogm_ndx    => spc_ndx(18)
         389        1536 :     sogi_ndx    => spc_ndx(19)
         390        1536 :     sogt_ndx    => spc_ndx(20)
         391        1536 :     sogb_ndx    => spc_ndx(21)
         392        1536 :     sogx_ndx    => spc_ndx(22)
         393             :
         394        1536 :     oc1_ndx      = get_spc_ndx( 'OC1' )
         395        1536 :     oc2_ndx      = get_spc_ndx( 'OC2' )
         396        1536 :     c10h16_ndx   = get_spc_ndx( 'C10H16')
         397        1536 :     isop_ndx     = get_spc_ndx( 'ISOP' )
         398        1536 :     tolo2_ndx    = get_spc_ndx( 'TOLO2' )
         399        1536 :     beno2_ndx    = get_spc_ndx( 'BENO2' )
         400        1536 :     xylo2_ndx    = get_spc_ndx( 'XYLO2' )
         401        1536 :     ho2_ndx      = get_spc_ndx( 'HO2' )
         402        1536 :     no_ndx       = get_spc_ndx( 'NO' )
         403        1536 :     o3_ndx       = get_spc_ndx( 'OX' )
         404        1536 :     if( o3_ndx < 1 ) then
         405        1536 :        o3_ndx =  get_spc_ndx( 'O3' )
         406             :     end if
         407        1536 :     oh_ndx       = get_spc_ndx( 'OH' )
         408        1536 :     no3_ndx      = get_spc_ndx( 'NO3' )
         409             :
         410        1536 :     soam_ndx     = get_spc_ndx( 'SOAM' )
         411        1536 :     soai_ndx     = get_spc_ndx( 'SOAI' )
         412        1536 :     soat_ndx     = get_spc_ndx( 'SOAT' )
         413        1536 :     soab_ndx     = get_spc_ndx( 'SOAB' )
         414        1536 :     soax_ndx     = get_spc_ndx( 'SOAX' )
         415             :
         416        1536 :     sogm_ndx     = get_spc_ndx( 'SOGM' )
         417        1536 :     sogi_ndx     = get_spc_ndx( 'SOGI' )
         418        1536 :     sogt_ndx     = get_spc_ndx( 'SOGT' )
         419        1536 :     sogb_ndx     = get_spc_ndx( 'SOGB' )
         420        1536 :     sogx_ndx     = get_spc_ndx( 'SOGX' )
         421             :
         422        1536 :     has_soa_equil = all( spc_ndx(1:22) > 0 )
         423        1536 :     has_soa = has_soa_equil
         424             :
         425        1536 :     if ( has_soa_equil ) then
         426             :        ! fracsog and fracsoa are added to pbuffer for persistence across restarts
         427           0 :        call pbuf_add_field( 'FRACSOG' ,'global',dtype_r8,(/pcols,pver,nrx,npr/), fracsog_ndx)
         428           0 :        call pbuf_add_field( 'FRACSOA' ,'global',dtype_r8,(/pcols,pver,nrx,npr/), fracsoa_ndx)
         429             :     else
         430             :        ! reassign these ndx pointers...
         431        1536 :        soa_ndx     => spc_ndx(1)
         432        1536 :        oc1_ndx     => spc_ndx(2)
         433        1536 :        oc2_ndx     => spc_ndx(3)
         434        1536 :        c10h16_ndx  => spc_ndx(4)
         435        1536 :        o3_ndx      => spc_ndx(5)
         436        1536 :        oh_ndx      => spc_ndx(6)
         437        1536 :        no3_ndx     => spc_ndx(7)
         438        1536 :        bigalk_ndx  => spc_ndx(8)
         439        1536 :        toluene_ndx => spc_ndx(9)
         440             :
         441        1536 :        soa_ndx      = get_spc_ndx( 'SOA' )
         442        1536 :        oc1_ndx      = get_spc_ndx( 'OC1' )
         443        1536 :        oc2_ndx      = get_spc_ndx( 'OC2' )
         444        1536 :        c10h16_ndx   = get_spc_ndx( 'C10H16')
         445        1536 :        o3_ndx       = get_spc_ndx( 'OX' )
         446        1536 :        if( o3_ndx < 1 ) then
         447        1536 :           o3_ndx =  get_spc_ndx( 'O3' )
         448             :        end if
         449        1536 :        oh_ndx       = get_spc_ndx( 'OH' )
         450        1536 :        no3_ndx      = get_spc_ndx( 'NO3' )
         451             :
         452        1536 :        bigalk_ndx   = get_spc_ndx( 'BIGALK' )
         453        1536 :        toluene_ndx  = get_spc_ndx( 'TOLUENE' )
         454             :
         455        1536 :        has_soa = all( spc_ndx(1:7) > 0 )
         456             :     endif
         457             :
         458        1536 :   end subroutine soa_register
         459             :
         460             :   !===============================================================================
         461             :   !===============================================================================
         462           0 :   subroutine soa_inti_equil(pbuf2d)
         463             :
         464        1536 :     use mo_chem_utls, only : get_spc_ndx, get_rxt_ndx
         465             :     use cam_history,  only : addfld
         466             :     use ppgrid,       only : pver
         467             :     use spmd_utils,   only : masterproc
         468             :     use cam_control_mod, only: initial_run
         469             :     use physics_buffer, only : physics_buffer_desc, pbuf_set_field, pbuf_get_chunk, pbuf_get_field
         470             :
         471             :     implicit none
         472             :
         473             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         474           0 :     type(physics_buffer_desc), pointer :: pbuf_ptr(:)
         475             :
         476             :     integer :: astat
         477             :     integer :: i,k,j, c
         478             :
         479           0 :     real(r8), pointer :: fracsog(:,:,:,:)  ! mass fraction of each SOA class from each reaction
         480           0 :     real(r8), pointer :: fracsoa(:,:,:,:)  ! mass fraction of each SOG class from each reaction
         481             :
         482             :     !-----------------------------------------------------------------------
         483             :     !       ... check if this is an aerosol simulation
         484             :     !-----------------------------------------------------------------------
         485           0 :     if( .not. has_soa ) then
         486             :        return
         487             :     end if
         488             :
         489             :     !-----------------------------------------------------------------------
         490             :     !       ... set reaction indicies
         491             :     !-----------------------------------------------------------------------
         492           0 :     rxn_soa(1) = get_rxt_ndx( 'C10H16_O3' )
         493           0 :     rxn_soa(2) = get_rxt_ndx( 'C10H16_OH' )
         494           0 :     rxn_soa(3) = get_rxt_ndx( 'C10H16_NO3' )
         495           0 :     rxn_soa(4) = get_rxt_ndx( 'ISOP_OH' )
         496           0 :     rxn_soa(5) = get_rxt_ndx( 'TOLO2_HO2' )
         497           0 :     rxn_soa(6) = get_rxt_ndx( 'TOLO2_NO' )
         498           0 :     if (rxn_soa(6)<0) then
         499           0 :        rxn_soa(6) = get_rxt_ndx( 'ox_p12' )
         500             :     endif
         501           0 :     rxn_soa(7) = get_rxt_ndx( 'BENO2_HO2' )
         502           0 :     rxn_soa(8) = get_rxt_ndx( 'BENO2_NO' )
         503           0 :     rxn_soa(9) = get_rxt_ndx( 'XYLO2_HO2' )
         504           0 :     rxn_soa(10) = get_rxt_ndx( 'XYLO2_NO' )
         505           0 :     if( all( rxn_soa(:) < 1 ) ) then
         506           0 :        has_soa = .false.
         507           0 :        return
         508             :     else
         509           0 :        if (masterproc) then
         510           0 :           write(iulog,*) '-----------------------------------------'
         511           0 :           write(iulog,*) 'soa_inti_equil: mozart will do soa aerosols'
         512           0 :           write(iulog,*) '-----------------------------------------'
         513             :        endif
         514             :     end if
         515             :
         516             :     !
         517             :     ! define reactants
         518             :     !
         519           0 :     react_ndx(1,1) = c10h16_ndx
         520           0 :     react_ndx(1,2) = o3_ndx
         521           0 :     react_ndx(2,1) = c10h16_ndx
         522           0 :     react_ndx(2,2) = oh_ndx
         523           0 :     react_ndx(3,1) = c10h16_ndx
         524           0 :     react_ndx(3,2) = no3_ndx
         525           0 :     react_ndx(4,1) = isop_ndx
         526           0 :     react_ndx(4,2) = oh_ndx
         527           0 :     react_ndx(5,1) = tolo2_ndx
         528           0 :     react_ndx(5,2) = ho2_ndx
         529           0 :     react_ndx(6,1) = tolo2_ndx
         530           0 :     react_ndx(6,2) = no_ndx
         531           0 :     react_ndx(7,1) = beno2_ndx
         532           0 :     react_ndx(7,2) = ho2_ndx
         533           0 :     react_ndx(8,1) = beno2_ndx
         534           0 :     react_ndx(8,2) = no_ndx
         535           0 :     react_ndx(9,1) = xylo2_ndx
         536           0 :     react_ndx(9,2) = ho2_ndx
         537           0 :     react_ndx(10,1) = xylo2_ndx
         538           0 :     react_ndx(10,2) = no_ndx
         539             :
         540           0 :     if ( masterproc ) then
         541           0 :        print *,'soa_inti ',c10h16_ndx, isop_ndx, tolo2_ndx, beno2_ndx, xylo2_ndx,o3_ndx, oh_ndx, no3_ndx
         542           0 :        print *,'soa_inti ',soam_ndx, soai_ndx, soab_ndx, soat_ndx, soax_ndx, oc1_ndx, oc2_ndx
         543           0 :        print *,'soa_inti ',sogm_ndx, sogi_ndx, sogb_ndx, sogt_ndx, sogx_ndx
         544           0 :        print *,'soa_inti ',react_ndx
         545             :     endif
         546             :     !
         547             :     ! define partitioning coefficients for each reaction
         548             :     ! bulk yields are from Seinfeld and Pandis (1998)
         549             :     !
         550             :     ! c10h16 + o3 (from Chung and Seinfeld, JGR, 107, 2002)
         551             :     !
         552           0 :     alpha(1,1)    = 0.067_r8
         553           0 :     alpha(1,2)    = 0.354_r8
         554           0 :     k_om(1,1)     = 0.184_r8
         555           0 :     k_om(1,2)     = 0.0043_r8
         556           0 :     T1(1)         = 310._r8
         557           0 :     delH(1)       = 42.e3_r8
         558             :
         559             :     ! c10h16 + oh (from Chung and Seinfeld, JGR, 107, 2002)
         560             :     !
         561           0 :     alpha(2,1)    = 0.067_r8
         562           0 :     alpha(2,2)    = 0.354_r8
         563           0 :     k_om(2,1)     = 0.184_r8
         564           0 :     k_om(2,2)     = 0.0043_r8
         565           0 :     T1(2)         = 310._r8
         566           0 :     delH(2)       = 42.e3_r8
         567             :     !
         568             :     ! c10h16 + no3 (from Chung and Seinfeld, JGR, 107, 2002)
         569             :     !
         570           0 :     alpha(3,1)    = 1.000_r8
         571           0 :     alpha(3,2)    = 0.000_r8
         572           0 :     k_om(3,1)     = 0.0163_r8
         573           0 :     k_om(3,2)     = 0.0000_r8
         574           0 :     T1(3)         = 310._r8
         575           0 :     delH(3)       = 42.e3_r8
         576             :     !
         577             :     !! isop + oh (from Henze and Seinfeld, GRL, 2006): low NOx
         578             :     !!
         579             :     !  alpha(4,1)    = 0.232_r8
         580             :     !  alpha(4,2)    = 0.0288_r8
         581             :     !  k_om(4,1)     = 0.00862_r8
         582             :     !  k_om(4,2)     = 1.62_r8
         583             :     !  T1(4)         = 295._r8
         584             :     !  delH(4)       = 42.e3_r8
         585             :     !!
         586             :     ! isop + oh (from Henze and Seinfeld, GRL, 2006): high NOx
         587             :     !
         588           0 :     alpha(4,1)    = 0.264_r8
         589           0 :     alpha(4,2)    = 0.0173_r8
         590           0 :     k_om(4,1)     = 0.00115_r8
         591           0 :     k_om(4,2)     = 1.52_r8
         592           0 :     T1(4)         = 295._r8
         593           0 :     delH(4)       = 42.e3_r8
         594             :     !
         595             :     ! toluene + oh (pers comm Seinfeld and Henze): low NOx (TOLO2 + HO2)
         596             :     !
         597           0 :     alpha(5,1)    = 0.2349_r8
         598           0 :     alpha(5,2)    = 0.0_r8
         599           0 :     k_om(5,1)     = 1000.0_r8
         600           0 :     k_om(5,2)     = 0.0_r8
         601           0 :     T1(5)         = 295._r8
         602           0 :     delH(5)       = 42.e3_r8
         603             :     !
         604             :     ! toluene + oh (pers comm Seinfeld and Henze): high NOx (TOLO2 + NO)
         605             :     !
         606           0 :     alpha(6,1)    = 0.0378_r8
         607           0 :     alpha(6,2)    = 0.0737_r8
         608           0 :     k_om(6,1)     = 0.4300_r8
         609           0 :     k_om(6,2)     = 0.0470_r8
         610           0 :     T1(6)         = 295._r8
         611           0 :     delH(6)       = 42.e3_r8
         612             :     !
         613             :     ! benzene + oh (pers comm Seinfeld and Henze): low NOx (BENO2 + HO2)
         614             :     !
         615           0 :     alpha(7,1)    = 0.2272_r8
         616           0 :     alpha(7,2)    = 0.0_r8
         617           0 :     k_om (7,1)    = 1000.0_r8
         618           0 :     k_om (7,2)    = 0.0_r8
         619           0 :     T1(7)         = 295._r8
         620           0 :     delH(7)       = 42.e3_r8
         621             :     !
         622             :     ! benzene + oh (pers comm Seinfeld and Henze): high NOx (BENO2 + NO)
         623             :     !
         624           0 :     alpha(8,1)    = 0.0442_r8
         625           0 :     alpha(8,2)    = 0.5454_r8
         626           0 :     k_om(8,1)     = 3.3150_r8
         627           0 :     k_om(8,2)     = 0.0090_r8
         628           0 :     T1(8)         = 295._r8
         629           0 :     delH(8)       = 42.e3_r8
         630             :     !
         631             :     ! xylene + oh (pers comm Seinfeld and Henze): low NOx (XYLO2 + HO2)
         632             :     !
         633           0 :     alpha(9,1)    = 0.2052_r8
         634           0 :     alpha(9,2)    = 0.0_r8
         635           0 :     k_om(9,1)     = 1000.0_r8
         636           0 :     k_om(9,2)     = 0.0_r8
         637           0 :     T1(9)         = 295._r8
         638           0 :     delH(9)       = 42.e3_r8
         639             :     !
         640             :     ! xylene + oh (pers comm Seinfeld and Henze): high NOx (XYLO2 + NO)
         641             :     !
         642           0 :     alpha(10,1)    = 0.0212_r8
         643           0 :     alpha(10,2)    = 0.0615_r8
         644           0 :     k_om(10,1)     = 0.7610_r8
         645           0 :     k_om(10,2)     = 0.0290_r8
         646           0 :     T1(10)         = 295._r8
         647           0 :     delH(10)       = 42.e3_r8
         648             :     !
         649           0 :     call addfld( 'SOAM_PROD', (/ 'lev' /), 'A', 'molec/molec/s', 'production of SOAM' )
         650           0 :     call addfld( 'SOAI_PROD', (/ 'lev' /), 'A', 'molec/molec/s', 'production of SOAI' )
         651           0 :     call addfld( 'SOAT_PROD', (/ 'lev' /), 'A', 'molec/molec/s', 'production of SOAT' )
         652           0 :     call addfld( 'SOAB_PROD', (/ 'lev' /), 'A', 'molec/molec/s', 'production of SOAB' )
         653           0 :     call addfld( 'SOAX_PROD', (/ 'lev' /), 'A', 'molec/molec/s', 'production of SOAX' )
         654             :
         655           0 :     call addfld( 'SOAM_dens', (/ 'lev' /), 'A', 'ug/m3', 'density of SOAM' )
         656           0 :     call addfld( 'SOAI_dens', (/ 'lev' /), 'A', 'ug/m3', 'density of SOAI' )
         657           0 :     call addfld( 'SOAT_dens', (/ 'lev' /), 'A', 'ug/m3', 'density of SOAT' )
         658           0 :     call addfld( 'SOAB_dens', (/ 'lev' /), 'A', 'ug/m3', 'density of SOAB' )
         659           0 :     call addfld( 'SOAX_dens', (/ 'lev' /), 'A', 'ug/m3', 'density of SOAX' )
         660             :
         661             :     !
         662             :     !initialize fracsoa for first timestep and store values for future
         663           0 :     fracsoa_init(1:3,:)=0.2_r8
         664           0 :     fracsoa_init(3,2)=0._r8
         665           0 :     fracsoa_init(4,:)=0.5_r8
         666           0 :     fracsoa_init(5:10,:)=0.33_r8
         667           0 :     fracsoa_init(5,2)=0._r8
         668           0 :     fracsoa_init(7,2)=0._r8
         669           0 :     fracsoa_init(9,2)=0._r8
         670             :
         671             :     !initialize fracsog for first timestep and store values for future
         672           0 :     fracsog_init(1:3,:)=0.2_r8
         673           0 :     fracsog_init(3,2)=0._r8
         674           0 :     fracsog_init(4,:)=0.5_r8
         675           0 :     fracsog_init(5:10,:)=0.33_r8
         676           0 :     fracsog_init(5,2)=0._r8
         677           0 :     fracsog_init(7,2)=0._r8
         678           0 :     fracsog_init(9,2)=0._r8
         679             :
         680           0 :     if (initial_run) then
         681           0 :        do c=begchunk, endchunk
         682           0 :           pbuf_ptr=>pbuf_get_chunk(pbuf2d, c)
         683           0 :           call pbuf_get_field(pbuf_ptr, fracsoa_ndx, fracsoa )
         684           0 :           call pbuf_get_field(pbuf_ptr, fracsog_ndx, fracsog )
         685           0 :           do i = 1,pcols
         686           0 :              do k = 1,pver
         687           0 :                 fracsoa( i,k, :,: ) = fracsoa_init(:,:)
         688           0 :                 fracsog( i,k, :,: ) = fracsog_init(:,:)
         689             :              enddo
         690             :           enddo
         691             :        enddo
         692             :     endif
         693             :
         694             :     return
         695           0 :   end subroutine soa_inti_equil
         696             :   !===============================================================================
         697             :   ! clh (08/01/06): added T dependence of gas-aerosol phase partitioning
         698             :   !                 added isoprene as SOA precursor [Henze and Seinfeld, 2006]
         699             :   !                 added anthropogenic SOA according to [Henze et al., 2007]
         700             :   !                 split SOA into classes (SOAM=monoterpenes,SOAI=isoprene,SOAB=benzene,SOAT=toluene
         701             :   !                                         SOAX=xylene)
         702             :   !                 fixed error in yield calculation
         703             :   !                 added scaling for OC to OM in pre-existing aerosol mass
         704             :   !                 modified yield calculation to allow for re-evaporation
         705             :   ! Note: do not need to subtract formed aerosol product from reactants, because gas-phase
         706             :   !       products do not account for entire mass
         707             :   !===============================================================================
         708           0 :   subroutine setsoa_equil(dt,reaction_rates,tfld,vmr,xhnm,ncol,lchnk,pbuf)
         709             :     !
         710             :     ! updated SOA mechanism
         711             :     ! based on Chung and Seinfeld, JGR, 2002
         712             :     !
         713             :     ! implemented in CAM by Colette Heald (summer 2007)
         714             :     !
         715           0 :     use ppgrid,       only : pcols, pver
         716             :     use chem_mods,    only : adv_mass, gas_pcnst, rxntot
         717             :     use cam_history,  only : outfld
         718             :     use physics_buffer, only : physics_buffer_desc, pbuf_get_field
         719             :     !
         720             :     implicit none
         721             :     !
         722             :     !-----------------------------------------------------------------------
         723             :     !      ... dummy arguments
         724             :     !-----------------------------------------------------------------------
         725             :     integer, intent(in)      :: ncol                   ! number columns in chunkx
         726             :     integer, intent(in)      :: lchnk                  ! chunk index
         727             :     real(r8), intent(in)     :: dt                     ! time step
         728             :     real(r8), intent(in)     :: reaction_rates(:,:,:)  ! reaction rates
         729             :     real(r8), intent(inout)  :: vmr(:,:,:)             ! xported species ( vmr )
         730             :     real(r8), intent(in)     :: tfld(:,:)              ! temperature (K)
         731             :     real(r8), intent(in)     :: xhnm(:,:)              ! total atms density (molec/cm**3)
         732             :     type(physics_buffer_desc), pointer :: pbuf(:)
         733             :
         734             :     !-----------------------------------------------------------------------
         735             :     !      ... local variables
         736             :     !-----------------------------------------------------------------------
         737             :     integer :: i,k,n,p,iter
         738             :     real(r8) :: T_fac
         739             :     real(r8) :: delHC, sumorg, poa, mnew, numer, maxM, minM, tol, mw_soa,m_air
         740             :     real(r8) :: k_om_T(NRX,NPR)                   ! equilibrium gas-particule partition (T dependent)
         741             :     real(r8) :: prod(NRX,NPR)                     ! oxidized produced (alpha_i*delHC)
         742             :     real(r8) :: sog0(NRX,NPR)                     ! pre-existing SOG
         743             :     real(r8) :: soa0(NRX,NPR)                     ! pre-exisiting SOA
         744             :     real(r8) :: orggas(NRX,NPR)                   ! intermediate (SOG0+prod of ox)
         745             :     real(r8) :: sog(NRX,NPR)                      ! final SOG
         746             :     real(r8) :: soa(NRX,NPR)                      ! final SOA
         747             :     real(r8), dimension(pcols,pver) :: soam_mass,soai_mass,soat_mass,soab_mass,soax_mass
         748             :     real(r8), dimension(pcols,pver) :: sogm_mass,sogi_mass,sogt_mass,sogb_mass,sogx_mass
         749           0 :     real(r8) :: soam_prod(ncol,pver), soai_prod(ncol,pver),soat_prod(ncol,pver),soab_prod(ncol,pver),soax_prod(ncol,pver)
         750             :
         751           0 :     real(r8), pointer, dimension(:,:,:,:) :: fracsog  ! mass fraction of each SOA class from each reaction
         752           0 :     real(r8), pointer, dimension(:,:,:,:) :: fracsoa  ! mass fraction of each SOG class from each reaction
         753             :
         754           0 :     soam_mass(:,:)=0._r8
         755           0 :     soai_mass(:,:)=0._r8
         756           0 :     soat_mass(:,:)=0._r8
         757           0 :     soab_mass(:,:)=0._r8
         758           0 :     soax_mass(:,:)=0._r8
         759           0 :     sogm_mass(:,:)=0._r8
         760           0 :     sogi_mass(:,:)=0._r8
         761           0 :     sogt_mass(:,:)=0._r8
         762           0 :     sogb_mass(:,:)=0._r8
         763           0 :     sogx_mass(:,:)=0._r8
         764             :
         765           0 :     call pbuf_get_field(pbuf, fracsoa_ndx, fracsoa )
         766           0 :     call pbuf_get_field(pbuf, fracsog_ndx, fracsog )
         767             :
         768           0 :     do i=1,ncol
         769           0 :        do k=1,pver
         770             :           !
         771             :           ! INIALIZATION AND LUMPING
         772             :           !
         773             :           !       calculate mass concentration of air (ug/m3)
         774           0 :           m_air = xhnm(i,k)*28.966_r8/avogadro*1.e12_r8
         775             :           !
         776             :           !       calculate initial mass of POA from OC1 and OC2 (in ug/m3)
         777           0 :           poa = (vmr(i,k,oc1_ndx)+vmr(i,k,oc2_ndx)) * OMscale * xhnm(i,k) * adv_mass(oc1_ndx)/avogadro * 1.e12_r8
         778             :           !
         779             :           !       specify pre-existing SOG/SOA for each class (in ug/m3)
         780             :
         781             :
         782             :           sog0(1,1)=vmr(i,k,sogm_ndx) * xhnm(i,k) * adv_mass(sogm_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 1,1)
         783             :           sog0(1,2)=vmr(i,k,sogm_ndx) * xhnm(i,k) * adv_mass(sogm_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 1,2)
         784             :           sog0(2,1)=vmr(i,k,sogm_ndx) * xhnm(i,k) * adv_mass(sogm_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 2,1)
         785             :           sog0(2,2)=vmr(i,k,sogm_ndx) * xhnm(i,k) * adv_mass(sogm_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 2,2)
         786             :           sog0(3,1)=vmr(i,k,sogm_ndx) * xhnm(i,k) * adv_mass(sogm_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 3,1)
         787             :           sog0(3,2)=0._r8
         788             :           sog0(4,1)=vmr(i,k,sogi_ndx) * xhnm(i,k) * adv_mass(sogi_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 4,1)
         789             :           sog0(4,2)=vmr(i,k,sogi_ndx) * xhnm(i,k) * adv_mass(sogi_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 4,2)
         790             :           sog0(5,1)=vmr(i,k,sogt_ndx) * xhnm(i,k) * adv_mass(sogt_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 5,1)
         791             :           sog0(5,2)=0._r8
         792             :           sog0(6,1)=vmr(i,k,sogt_ndx) * xhnm(i,k) * adv_mass(sogt_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 6,1)
         793             :           sog0(6,2)=vmr(i,k,sogt_ndx) * xhnm(i,k) * adv_mass(sogt_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 6,2)
         794             :           sog0(7,1)=vmr(i,k,sogb_ndx) * xhnm(i,k) * adv_mass(sogb_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 7,1)
         795             :           sog0(7,2)=0._r8
         796             :           sog0(8,1)=vmr(i,k,sogb_ndx) * xhnm(i,k) * adv_mass(sogb_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 8,1)
         797             :           sog0(8,2)=vmr(i,k,sogb_ndx) * xhnm(i,k) * adv_mass(sogb_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 8,2)
         798             :           sog0(9,1)=vmr(i,k,sogx_ndx) * xhnm(i,k) * adv_mass(sogx_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 9,1)
         799             :           sog0(9,2)=0._r8
         800             :           sog0(10,1)=vmr(i,k,sogx_ndx) * xhnm(i,k) * adv_mass(sogx_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 10,1)
         801             :           sog0(10,2)=vmr(i,k,sogx_ndx) * xhnm(i,k) * adv_mass(sogx_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 10,2)
         802             :           !
         803             :           soa0(1,1)=vmr(i,k,soam_ndx) * xhnm(i,k) * adv_mass(soam_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 1,1)
         804             :           soa0(1,2)=vmr(i,k,soam_ndx) * xhnm(i,k) * adv_mass(soam_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 1,2)
         805             :           soa0(2,1)=vmr(i,k,soam_ndx) * xhnm(i,k) * adv_mass(soam_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 2,1)
         806             :           soa0(2,2)=vmr(i,k,soam_ndx) * xhnm(i,k) * adv_mass(soam_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 2,2)
         807             :           soa0(3,1)=vmr(i,k,soam_ndx) * xhnm(i,k) * adv_mass(soam_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 3,1)
         808             :           soa0(3,2)=0._r8
         809             :           soa0(4,1)=vmr(i,k,soai_ndx) * xhnm(i,k) * adv_mass(soai_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 4,1)
         810             :           soa0(4,2)=vmr(i,k,soai_ndx) * xhnm(i,k) * adv_mass(soai_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 4,2)
         811             :           soa0(5,1)=vmr(i,k,soat_ndx) * xhnm(i,k) * adv_mass(soat_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 5,1)
         812             :           soa0(5,2)=0._r8
         813             :           soa0(6,1)=vmr(i,k,soat_ndx) * xhnm(i,k) * adv_mass(soat_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 6,1)
         814             :           soa0(6,2)=vmr(i,k,soat_ndx) * xhnm(i,k) * adv_mass(soat_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 6,2)
         815             :           soa0(7,1)=vmr(i,k,soab_ndx) * xhnm(i,k) * adv_mass(soab_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 7,1)
         816             :           soa0(7,2)=0._r8
         817             :           soa0(8,1)=vmr(i,k,soab_ndx) * xhnm(i,k) * adv_mass(soab_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 8,1)
         818             :           soa0(8,2)=vmr(i,k,soab_ndx) * xhnm(i,k) * adv_mass(soab_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 8,2)
         819             :           soa0(9,1)=vmr(i,k,soax_ndx) * xhnm(i,k) * adv_mass(soax_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 9,1)
         820             :           soa0(9,2)=0._r8
         821             :           soa0(10,1)=vmr(i,k,soax_ndx) * xhnm(i,k) * adv_mass(soax_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 10,1)
         822             :           soa0(10,2)=vmr(i,k,soax_ndx) * xhnm(i,k) * adv_mass(soax_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 10,2)
         823             :           !
         824             :           !
         825             :           !--------------------------------------------------------------
         826             :           ! CHEMISTRY
         827             :           !
         828             :           sumorg=0._r8
         829             :           !
         830             :           do n=1,NRX
         831             :              !
         832             :              !         temperature dependence of paritioning via Clausiaus Clayperon [C&S, 2002]: clh (02/24/06)
         833             :              T_fac = tfld(i,k)/T1(n) * exp(delH(n)/Rgas*(1/tfld(i,k) - 1/T1(n)))
         834             :              k_om_T(n,:)=k_om(n,:)*T_fac
         835             :              !
         836             :              !         find molecular weight of SOA
         837             :              if ( (n >=1) .and. (n < 4) ) then
         838             :                 mw_soa = adv_mass(soam_ndx)
         839             :              else if (n == 4) then
         840             :                 mw_soa = adv_mass(soai_ndx)
         841             :              else if ( (n > 4) .and. (n < 7) ) then
         842             :                 mw_soa = adv_mass(soat_ndx)
         843             :              else if ( (n >= 7) .and. (n < 9) ) then
         844             :                 mw_soa = adv_mass(soab_ndx)
         845             :              else if ( (n >= 9) .and. (n < 11) ) then
         846             :                 mw_soa = adv_mass(soax_ndx)
         847             :              end if
         848             :              !
         849             :              !         define chemical production from gas-phase chemistry (and convert to ug/m3)
         850             :              delHC  = reaction_rates(i,k,rxn_soa(n)) &
         851             :                   * vmr(i,k,react_ndx(n,1)) * vmr(i,k,react_ndx(n,2))* dt &
         852             :                   * xhnm(i,k) * mw_soa/avogadro * 1.e12_r8
         853             :              !
         854             :              !         specify the new total of gas-phase products (before re-partitioning)
         855             :              !         and total up all SOA in gas and aerosol phase
         856             :              do p=1,NPR
         857             :                 orggas(n,p) = sog0(n,p)+alpha(n,p)*delHC
         858             :                 sumorg=sumorg+orggas(n,p)+soa0(n,p)
         859             :                 prod(n,p)=alpha(n,p)*delHC
         860             :              enddo
         861             :              !
         862             :           enddo
         863             :           !
         864             :           !--------------------------------------------------------------
         865             :           !       check to see if no organics no partitioning!
         866             :           !             (set here to previous timestep concentration+oxprod to preserve mass)
         867             :           if (sumorg < 1.e-10_r8) then
         868             :              do n=1,NRX
         869             :                 do p=1,NPR
         870             :                    sog(n,p)=sog0(n,p)
         871             :                    soa(n,p)=soa0(n,p)+prod(n,p)
         872             :                 enddo
         873             :              enddo
         874             :           else
         875             :              ! PARTITION (Need to iteratively solve for MNEW, set tolerances to 0.1 ng/m3 or 1% of MNEW)
         876             :              !
         877             :              numer=0._r8
         878             :              maxM=0._r8
         879             :              !           If POA is essentially zero, equations simplify
         880             :              if (POA < 1.e-10_r8) then
         881             :                 do n=1,NRX
         882             :                    do p=1,NPR
         883             :                       numer = numer + k_om_T(n,p)* (orggas(n,p)+soa0(n,p))
         884             :                    enddo
         885             :                 enddo
         886             :                 !              if numerator is less than 1 then MNEW must be zero
         887             :                 if (numer <= 1._r8) then
         888             :                    mnew=0._r8
         889             :                    iter=0
         890             :                 else
         891             :                    minM = 1.e-40_r8
         892             :                    maxM = poa + sumorg
         893             :                    tol = 1.e-4_r8
         894             :                    mnew = zeroin(minM,maxM,tol,poa,soa0,orggas,k_om_T,sumorg,iter)
         895             :                 end if
         896             :                 !
         897             :                 !           if additional organic mass is less than 1% of POA, or difference between the
         898             :                 !           two is less than 0.1 ng/m3 (the tolerance) then POA is partitioning mass
         899             :              else if ( (sumorg < 0.01_r8*poa) .or. (abs(sumorg-poa) < 1.e-4_r8) ) then
         900             :                 mnew = poa
         901             :                 iter = 0
         902             :                 !           otherwise solve for MNEW iteratively on interval [poa, poa+sumorg]
         903             :              else
         904             :                 maxM = poa+sumorg
         905             :                 minM = poa
         906             :                 tol = 1.e-4_r8
         907             :                 mnew = zeroin(minM,maxM,tol,poa,soa0,orggas,k_om_T,sumorg,iter)
         908             :              end if
         909             :              !
         910             :              !
         911             :              !           Now equilibrium partitioning with new MNEW
         912             :              !           If no MNEW then all the SOA evaporates to gas-phase
         913             :              if (mnew > 0._r8) then
         914             :                 do n=1,NRX
         915             :                    do p=1,NPR
         916             :                       soa(n,p)=k_om_T(n,p)*mnew*(orggas(n,p)+soa0(n,p))/(1._r8+k_om_T(n,p)*mnew)
         917             :                       if (k_om_T(n,p) /= 0._r8) then
         918             :                          sog(n,p)=soa(n,p)/(k_om_T(n,p)*mnew)
         919             :                       else
         920             :                          sog(n,p)=0._r8
         921             :                       end if
         922             :                    enddo
         923             :                 enddo
         924             :
         925             :              else
         926             :                 do n=1,NRX
         927             :                    do p=1,NPR
         928             :                       sog(n,p)=orggas(n,p)+soa0(n,p)
         929             :                       soa(n,p)=1.e-20_r8
         930             :                    enddo
         931             :                 enddo
         932             :              end if
         933             :              !
         934             :           end if
         935             :           !
         936             :           !--------------------------------------------------------------
         937             :           ! LUMP INTO ARRAYS
         938             :           do n=1,NRX
         939             :              do p=1,NPR
         940             :                 if ( (n >=1) .and. (n < 4) ) then
         941             :                    soam_mass(i,k) = soam_mass(i,k) + soa(n,p)
         942             :                    sogm_mass(i,k) = sogm_mass(i,k) + sog(n,p)
         943             :                 else if (n == 4) then
         944             :                    soai_mass(i,k) = soai_mass(i,k) + soa(n,p)
         945             :                    sogi_mass(i,k) = sogi_mass(i,k) + sog(n,p)
         946             :                 else if ( (n > 4) .and. (n < 7) ) then
         947             :                    soat_mass(i,k) = soat_mass(i,k) + soa(n,p)
         948             :                    sogt_mass(i,k) = sogt_mass(i,k) + sog(n,p)
         949             :                 else if ( (n >= 7) .and. (n < 9) ) then
         950             :                    soab_mass(i,k) = soab_mass(i,k) + soa(n,p)
         951             :                    sogb_mass(i,k) = sogb_mass(i,k) + sog(n,p)
         952             :                 else if ( (n >= 9) .and. (n < 11) ) then
         953             :                    soax_mass(i,k) = soax_mass(i,k) + soa(n,p)
         954             :                    sogx_mass(i,k) = sogx_mass(i,k) + sog(n,p)
         955             :                 end if
         956             :              enddo
         957             :           enddo
         958             :           !
         959             :           !       Save mass fraction of each SOA rxn to SOA class
         960             :           !       (but if sumorg essentially zero revert to init fracs OR
         961             :           !        if mnew = 0 (all evap) then fracsoa revert to old, calculate fracsog only)
         962             :           if (sumorg < 1.e-10_r8) then
         963             :              do n=1,NRX
         964             :                 do p=1,NPR
         965             :                    fracsoa(i,k,n,p)=fracsoa_init(n,p)
         966             :                    fracsog(i,k,n,p)=fracsog_init(n,p)
         967             :                 enddo
         968             :              enddo
         969             :           else
         970             :              if (mnew < 1.e-20_r8) then
         971             :                 do n=1,NRX
         972             :                    do p=1,NPR
         973             :                       fracsoa(i,k,n,p)=fracsoa_init(n,p)
         974             :                       if ( (n >=1) .and. (n < 4) ) then
         975             :                          if (sogm_mass(i,k) == 0._r8) then
         976             :                             fracsog(i,k,n,p)=fracsog_init(n,p)
         977             :                          else
         978             :                             fracsog(i,k,n,p)=sog(n,p)/sogm_mass(i,k)
         979             :                          end if
         980             :                       else if (n == 4) then
         981             :                          if (sogi_mass(i,k) == 0._r8) then
         982             :                             fracsog(i,k,n,p)=fracsog_init(n,p)
         983             :                          else
         984             :                             fracsog(i,k,n,p)=sog(n,p)/sogi_mass(i,k)
         985             :                          end if
         986             :                       else if ( (n > 4) .and. (n < 7) ) then
         987             :                          if (sogt_mass(i,k) == 0._r8) then
         988             :                             fracsog(i,k,n,p)=fracsog_init(n,p)
         989             :                          else
         990             :                             fracsog(i,k,n,p)=sog(n,p)/sogt_mass(i,k)
         991             :                          end if
         992             :                       else if ( (n >= 7) .and. (n < 9) ) then
         993             :                          if (sogb_mass(i,k) == 0._r8) then
         994             :                             fracsog(i,k,n,p)=fracsog_init(n,p)
         995             :                          else
         996             :                             fracsog(i,k,n,p)=sog(n,p)/sogb_mass(i,k)
         997             :                          end if
         998             :                       else if ( (n >= 9) .and. (n < 11) ) then
         999             :                          if (sogx_mass(i,k) == 0._r8) then
        1000             :                             fracsog(i,k,n,p)=fracsog_init(n,p)
        1001             :                          else
        1002             :                             fracsog(i,k,n,p)=sog(n,p)/sogx_mass(i,k)
        1003             :                          end if
        1004             :                       end if
        1005             :                       if ( (p==2) .and. (n==3 .or. n==5 .or. n==7 .or. n==9) ) then
        1006             :                          fracsog(i,k,n,p)=0._r8
        1007             :                       end if
        1008             :                    enddo
        1009             :                 enddo
        1010             :              else
        1011             :                 do n=1,NRX
        1012             :                    do p=1,NPR
        1013             :                       if ( (n >=1) .and. (n < 4) ) then
        1014             :                          if (soam_mass(i,k) == 0._r8) then
        1015             :                             fracsoa(i,k,n,p)=fracsoa_init(n,p)
        1016             :                          else
        1017             :                             fracsoa(i,k,n,p)=soa(n,p)/soam_mass(i,k)
        1018             :                          end if
        1019             :                          if (sogm_mass(i,k) == 0._r8) then
        1020             :                             fracsog(i,k,n,p)=fracsog_init(n,p)
        1021             :                          else
        1022             :                             fracsog(i,k,n,p)=sog(n,p)/sogm_mass(i,k)
        1023             :                          end if
        1024             :                       else if (n == 4) then
        1025             :                          if (soai_mass(i,k) == 0._r8) then
        1026             :                             fracsoa(i,k,n,p)=fracsoa_init(n,p)
        1027             :                          else
        1028             :                             fracsoa(i,k,n,p)=soa(n,p)/soai_mass(i,k)
        1029             :                          end if
        1030             :                          if (sogi_mass(i,k) == 0._r8) then
        1031             :                             fracsog(i,k,n,p)=fracsog_init(n,p)
        1032             :                          else
        1033             :                             fracsog(i,k,n,p)=sog(n,p)/sogi_mass(i,k)
        1034             :                          end if
        1035             :                       else if ( (n > 4) .and. (n < 7) ) then
        1036             :                          if (soat_mass(i,k) == 0._r8) then
        1037             :                             fracsoa(i,k,n,p)=fracsoa_init(n,p)
        1038             :                          else
        1039             :                             fracsoa(i,k,n,p)=soa(n,p)/soat_mass(i,k)
        1040             :                          end if
        1041             :                          if (sogt_mass(i,k) == 0._r8) then
        1042             :                             fracsog(i,k,n,p)=fracsog_init(n,p)
        1043             :                          else
        1044             :                             fracsog(i,k,n,p)=sog(n,p)/sogt_mass(i,k)
        1045             :                          end if
        1046             :                       else if ( (n >= 7) .and. (n < 9) ) then
        1047             :                          if (soab_mass(i,k) == 0._r8) then
        1048             :                             fracsoa(i,k,n,p)=fracsoa_init(n,p)
        1049             :                          else
        1050             :                             fracsoa(i,k,n,p)=soa(n,p)/soab_mass(i,k)
        1051             :                          end if
        1052             :                          if (sogb_mass(i,k) == 0._r8) then
        1053             :                             fracsog(i,k,n,p)=fracsog_init(n,p)
        1054             :                          else
        1055             :                             fracsog(i,k,n,p)=sog(n,p)/sogb_mass(i,k)
        1056             :                          end if
        1057             :                       else if ( (n >= 9) .and. (n < 11) ) then
        1058             :                          if (soax_mass(i,k) == 0._r8) then
        1059             :                             fracsoa(i,k,n,p)=fracsoa_init(n,p)
        1060             :                          else
        1061             :                             fracsoa(i,k,n,p)=soa(n,p)/soax_mass(i,k)
        1062             :                          end if
        1063             :                          if (sogx_mass(i,k) == 0._r8) then
        1064             :                             fracsog(i,k,n,p)=fracsog_init(n,p)
        1065             :                          else
        1066             :                             fracsog(i,k,n,p)=sog(n,p)/sogx_mass(i,k)
        1067             :                          end if
        1068             :                       end if
        1069             :                       if ( (p==2) .and. (n==3 .or. n==5 .or. n==7 .or. n==9) ) then
        1070             :                          fracsoa(i,k,n,p)=0._r8
        1071             :                          fracsog(i,k,n,p)=0._r8
        1072             :                       end if
        1073             :                    enddo
        1074             :                 enddo
        1075             :              end if
        1076             :           end if
        1077             :           !
        1078             :           !--------------------------------------------------------------
        1079             :           !
        1080             :           ! calculate NET production in kg/kg/s (subtract initial mass)
        1081             :           !
        1082             :           soam_prod(i,k) = ( soam_mass(i,k)*1.e-12_r8*avogadro/(adv_mass(soam_ndx)*xhnm(i,k)) - &
        1083             :                            vmr(i,k,soam_ndx) )/dt
        1084             :           soai_prod(i,k) = ( soai_mass(i,k)*1.e-12_r8*avogadro/(adv_mass(soai_ndx)*xhnm(i,k)) - &
        1085             :                            vmr(i,k,soai_ndx) )/dt
        1086             :           soat_prod(i,k) = ( soat_mass(i,k)*1.e-12_r8*avogadro/(adv_mass(soat_ndx)*xhnm(i,k)) - &
        1087             :                            vmr(i,k,soat_ndx) )/dt
        1088             :           soab_prod(i,k) = ( soab_mass(i,k)*1.e-12_r8*avogadro/(adv_mass(soab_ndx)*xhnm(i,k)) - &
        1089             :                            vmr(i,k,soab_ndx) )/dt
        1090             :           soax_prod(i,k) = ( soax_mass(i,k)*1.e-12_r8*avogadro/(adv_mass(soax_ndx)*xhnm(i,k)) - &
        1091             :                            vmr(i,k,soax_ndx) )/dt
        1092             :           !
        1093             :           ! convert from ug/m3 to mixing ratio and update vmr
        1094             :           !
        1095             :           vmr(i,k,soam_ndx) = soam_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(soam_ndx)*xhnm(i,k))
        1096             :           vmr(i,k,soai_ndx) = soai_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(soai_ndx)*xhnm(i,k))
        1097             :           vmr(i,k,soat_ndx) = soat_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(soat_ndx)*xhnm(i,k))
        1098             :           vmr(i,k,soab_ndx) = soab_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(soab_ndx)*xhnm(i,k))
        1099             :           vmr(i,k,soax_ndx) = soax_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(soax_ndx)*xhnm(i,k))
        1100             :           vmr(i,k,sogm_ndx) = sogm_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(sogm_ndx)*xhnm(i,k))
        1101             :           vmr(i,k,sogi_ndx) = sogi_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(sogi_ndx)*xhnm(i,k))
        1102             :           vmr(i,k,sogt_ndx) = sogt_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(sogt_ndx)*xhnm(i,k))
        1103             :           vmr(i,k,sogb_ndx) = sogb_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(sogb_ndx)*xhnm(i,k))
        1104             :           vmr(i,k,sogx_ndx) = sogx_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(sogx_ndx)*xhnm(i,k))
        1105             :        enddo
        1106             :     enddo
        1107             :     !
        1108           0 :     call outfld('SOAM_PROD',soam_prod(:ncol,:),ncol, lchnk)
        1109           0 :     call outfld('SOAI_PROD',soai_prod(:ncol,:),ncol, lchnk)
        1110           0 :     call outfld('SOAT_PROD',soat_prod(:ncol,:),ncol, lchnk)
        1111           0 :     call outfld('SOAB_PROD',soab_prod(:ncol,:),ncol, lchnk)
        1112           0 :     call outfld('SOAX_PROD',soax_prod(:ncol,:),ncol, lchnk)
        1113             :
        1114           0 :     call outfld('SOAM_dens',soam_mass(:ncol,:), ncol, lchnk)
        1115           0 :     call outfld('SOAI_dens',soai_mass(:ncol,:), ncol, lchnk)
        1116           0 :     call outfld('SOAT_dens',soat_mass(:ncol,:), ncol, lchnk)
        1117           0 :     call outfld('SOAB_dens',soab_mass(:ncol,:), ncol, lchnk)
        1118           0 :     call outfld('SOAX_dens',soax_mass(:ncol,:), ncol, lchnk)
        1119             :
        1120             :     !
        1121             :     !
        1122           0 :     return
        1123           0 :   end subroutine setsoa_equil
        1124             :
        1125             :   !===============================================================================
        1126             :   !===============================================================================
        1127             :   real(r8) function zeroin(x1,x2,tol,poa,aer,gas,k,totorg,iter)
        1128             :     ! function to iteratively solve function using bilinear method
        1129             :     !
        1130             :     implicit none
        1131             :     !
        1132             :     integer             :: iter
        1133             :     real(r8),intent(in) :: x1, x2                      ! min/max of interval
        1134             :     real(r8),intent(in) :: tol                         ! tolerance (interval of uncertainty)
        1135             :     real(r8),intent(in) :: poa,totorg
        1136             :     real(r8),intent(in) :: aer(NRX,NPR), gas(NRX,NPR)  ! aerosol and gas phase concentrations
        1137             :     real(r8),intent(in) :: k(NRX,NPR)                  ! partitioning coeff
        1138             :     ! local vars
        1139             :     real(r8)            :: xa,xb,xm,fa,fb,fm
        1140             :     !
        1141             :     xa=x1
        1142             :     xb=x2
        1143             :     xm=0._r8
        1144             :     fa=soa_function(xa,poa,aer,gas,k)
        1145             :     fb=soa_function(xb,poa,aer,gas,k)
        1146             :     !
        1147             :     ! check that functions have opposite signs
        1148             :     if (fa >= 0._r8) then
        1149             :        if (fb >=0._r8) then
        1150             :           write(iulog,*) 'ABORT IN ZEROIN: SAME SIGN ON FUNCTION',poa,totorg,x1,x2,fa,fb,aer,gas,k
        1151             :           write(iulog,*) 'ABORT IN ZEROIN: ERROR1: fa, fb ',fa, fb
        1152             :           write(iulog,*) 'ABORT IN ZEROIN: ERROR1: maxval(aer),minval(aer),maxval(gas),minval(gas) ',&
        1153             :                                                    maxval(aer),minval(aer),maxval(gas),minval(gas)
        1154             :           call endrun('ABORT IN ZEROIN: ERROR1')
        1155             :        end if
        1156             :     else
        1157             :        if (fb <=0._r8) then
        1158             :           write(iulog,*) 'ABORT IN ZEROIN: SAME SIGN ON FUNCTION',poa,totorg,x1,x2,fa,fb,aer,gas,k
        1159             :           write(iulog,*) 'ABORT IN ZEROIN: ERROR2: fa, fb ',fa, fb
        1160             :           write(iulog,*) 'ABORT IN ZEROIN: ERROR2: maxval(aer),minval(aer),maxval(gas),minval(gas) ',&
        1161             :                                                    maxval(aer),minval(aer),maxval(gas),minval(gas)
        1162             :           call endrun('ABORT IN ZEROIN: ERROR2')
        1163             :        end if
        1164             :     end if
        1165             :     !
        1166             :     iter=0
        1167             :     do while ((abs(xa-xb) > 2._r8*tol) .and. (abs(xa-xb) > 0.01_r8*xm) )
        1168             :        xm=(xa+xb)/2
        1169             :        fm=soa_function(xm,poa,aer,gas,k)
        1170             :        if (fa >=0._r8) then
        1171             :           if (fm >=0._r8) then
        1172             :              xa=xm
        1173             :              fa=fm
        1174             :           else
        1175             :              xb=xm
        1176             :              fb=fm
        1177             :           end if
        1178             :        else
        1179             :           if (fm < 0._r8) then
        1180             :              xa=xm
        1181             :              fa=fm
        1182             :           else
        1183             :              xb=xm
        1184             :              fb=fm
        1185             :           end if
        1186             :        end if
        1187             :        iter=iter+1
        1188             :     enddo
        1189             :     !
        1190             :     zeroin = (xa+xb)/2
        1191             :     !
        1192             :     return
        1193           0 :   end function zeroin
        1194             :   !===============================================================================
        1195             :   !===============================================================================
        1196             :   real(r8) function soa_function(m0,poa,aer,gas,k)
        1197             :     ! function which calculates SOAeqn (trying to minimize to zero)
        1198             :     !
        1199             :     implicit none
        1200             :     !
        1201             :     real(r8),intent(in) :: m0,poa
        1202             :     real(r8),intent(in) :: aer(NRX,NPR), gas(NRX,NPR)  ! aerosol and gas phase concentrations
        1203             :     real(r8),intent(in) :: k(NRX,NPR)                  ! partitioning coeff
        1204             :     ! local vars
        1205             :     integer             :: n,p
        1206             :     real(r8)            :: value
        1207             :     !
        1208             :     value=0._r8
        1209             :     !
        1210             :     do n=1,NRX
        1211             :        do p=1,NPR
        1212             :           value = value + k(n,p)*(gas(n,p)+aer(n,p))/(1._r8 + k(n,p)*m0)
        1213             :        enddo
        1214             :     enddo
        1215             :     soa_function = value + (poa/m0) - 1._r8
        1216             :     !
        1217             :     !  write(iulog,*) 'clh soa_function out: ',m0,soa_function
        1218             :     return
        1219             :   end function soa_function
        1220             :   !===============================================================================
        1221             :
        1222             : end module mo_setsoa

Generated by: LCOV version 1.14