LCOV - code coverage report
Current view: | top level - chemistry/modal_aero - modal_aero_gasaerexch.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 615 | 785 | 78.3 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 4 | 4 | 100.0 %

              Line data    Source code
           1             : ! modal_aero_gasaerexch.F90
           2             :
           3             :
           4             : !----------------------------------------------------------------------
           5             : !----------------------------------------------------------------------
           6             : !BOP
           7             : !
           8             : ! !MODULE: modal_aero_gasaerexch --- does modal aerosol gas-aerosol exchange
           9             : !
          10             : ! !INTERFACE:
          11             :    module modal_aero_gasaerexch
          12             :
          13             : ! !USES:
          14             :   use shr_kind_mod,    only:  r8 => shr_kind_r8
          15             :   use chem_mods,       only:  gas_pcnst
          16             :   use modal_aero_data, only:  nspec_max, nsoa, npoa, soa_multi_species
          17             :   use ref_pres,        only:  top_lev => clim_modal_aero_top_lev
          18             :   use ppgrid,          only:  pcols, pver
          19             :   use modal_aero_data, only:  ntot_amode, numptr_amode, sigmag_amode
          20             :   use modal_aero_data, only: lptr2_soa_g_amode, lptr2_soa_a_amode, lptr2_pom_a_amode
          21             :
          22             :   implicit none
          23             :   private
          24             :   save
          25             :
          26             : ! !PUBLIC MEMBER FUNCTIONS:
          27             :   public modal_aero_gasaerexch_sub, modal_aero_gasaerexch_init
          28             :
          29             : ! !PUBLIC DATA MEMBERS:
          30             :   integer, parameter :: pcnstxx = gas_pcnst
          31             :   integer, protected, public :: maxspec_pcage != nspec_max
          32             :
          33             :   integer, protected, public :: modefrm_pcage
          34             :   integer, protected, public :: nspecfrm_pcage
          35             :   integer :: modetoo_pcage
          36             :
          37             :   integer, protected, allocatable, public :: lspecfrm_pcage(:)
          38             :   integer, protected, allocatable, public :: lspectoo_pcage(:)
          39             :
          40             :   real(r8), parameter, public :: n_so4_monolayers_pcage = 8.0_r8
          41             :
          42             : ! number of so4(+nh4) monolayers needed to "age" a carbon particle
          43             :
          44             :   real(r8), parameter, public :: &
          45             :               dr_so4_monolayers_pcage = n_so4_monolayers_pcage * 4.76e-10_r8
          46             : ! thickness of the so4 monolayers (m)
          47             : ! for so4(+nh4), use bi-sulfate mw and 1.77 g/cm3,
          48             : !    --> 1 mol so4(+nh4)  = 65 cm^3 --> 1 molecule = (4.76e-10 m)^3
          49             : ! aging criterion is approximate so do not try to distinguish
          50             : !    sulfuric acid, bisulfate, ammonium sulfate
          51             :
          52             :   real(r8), protected, allocatable, public :: soa_equivso4_factor(:)
          53             : ! this factor converts an soa volume to a volume of so4(+nh4)
          54             : ! having same hygroscopicity as the soa
          55             :
          56             :   real (r8) :: fac_m2v_nh4, fac_m2v_so4
          57             :   real (r8), allocatable :: fac_m2v_soa(:)
          58             :
          59             :   real (r8), allocatable :: fac_m2v_pcarbon(:)
          60             :
          61             : ! !DESCRIPTION: This module implements ...
          62             : !
          63             : ! !REVISION HISTORY:
          64             : !
          65             : !   RCE 07.04.13:  Adapted from MIRAGE2 code
          66             : !
          67             : !EOP
          68             : !----------------------------------------------------------------------
          69             : !BOC
          70             :
          71             : ! list private module data here
          72             :
          73             : !EOC
          74             : !----------------------------------------------------------------------
          75             :
          76             :
          77             :   contains
          78             :
          79             :
          80             : !----------------------------------------------------------------------
          81             : !----------------------------------------------------------------------
          82             : !BOP
          83             : ! !ROUTINE:  modal_aero_gasaerexch_sub --- ...
          84             : !
          85             : ! !INTERFACE:
          86     1490712 : subroutine modal_aero_gasaerexch_sub(                            &
          87             :                         lchnk,    ncol,     nstep,               &
          88             :                         loffset,  deltat,                        &
          89             :                         t,        pmid,     pdel,                &
          90             :                         qh2o,               troplev,             &
          91     1489176 :                         q,                  qqcw,                &
          92     1489176 :                         dqdt_other,         dqqcwdt_other,       &
          93     1489176 :                         dgncur_a,           dgncur_awet,         &
          94             :                         sulfeq         )
          95             :
          96             : ! !USES:
          97             : use modal_aero_data,   only:  alnsg_amode,lmassptr_amode,cnst_name_cw
          98             : use modal_aero_data,   only:  lptr_so4_a_amode,lptr_nh4_a_amode
          99             : use modal_aero_data,   only:  modeptr_pcarbon,nspec_amode,specmw_amode,specdens_amode
         100             : use modal_aero_rename, only:  modal_aero_rename_sub
         101             : use rad_constituents,  only: rad_cnst_get_info
         102             : use constituents,      only: pcnst, cnst_mw
         103             :
         104             : use cam_history,       only:  outfld, fieldname_len
         105             : use chem_mods,         only:  adv_mass
         106             : use constituents,      only:  pcnst, cnst_name, cnst_get_ind
         107             : use mo_tracname,       only:  solsym
         108             : use physconst,         only:  gravit, mwdry, rair
         109             : use cam_abortutils,    only:  endrun
         110             : use spmd_utils,        only:  iam, masterproc
         111             : use phys_control,      only:  cam_chempkg_is
         112             :
         113             : implicit none
         114             :
         115             : ! !PARAMETERS:
         116             :    integer,  intent(in)    :: lchnk                ! chunk identifier
         117             :    integer,  intent(in)    :: ncol                 ! number of atmospheric column
         118             :    integer,  intent(in)    :: nstep                ! model time-step number
         119             :    integer,  intent(in)    :: loffset              ! offset applied to modal aero "ptrs"
         120             :    integer,  intent(in)    :: troplev(pcols)       ! tropopause vertical index
         121             :    real(r8), intent(in)    :: deltat               ! time step (s)
         122             :
         123             :    real(r8), intent(inout) :: q(ncol,pver,pcnstxx) ! tracer mixing ratio (TMR) array
         124             :                                                    ! *** MUST BE  #/kmol-air for number
         125             :                                                    ! *** MUST BE mol/mol-air for mass
         126             :                                                    ! *** NOTE ncol dimension
         127             :    real(r8), intent(inout) :: qqcw(ncol,pver,pcnstxx)
         128             :                                                    ! like q but for cloud-borner tracers
         129             :    real(r8), intent(in)    :: dqdt_other(ncol,pver,pcnstxx)
         130             :                                                    ! TMR tendency from other continuous
         131             :                                                    ! growth processes (aqchem, soa??)
         132             :                                                    ! *** NOTE ncol dimension
         133             :    real(r8), intent(in)    :: dqqcwdt_other(ncol,pver,pcnstxx)
         134             :                                                    ! like dqdt_other but for cloud-borner tracers
         135             :    real(r8), intent(in)    :: t(pcols,pver)        ! temperature at model levels (K)
         136             :    real(r8), intent(in)    :: pmid(pcols,pver)     ! pressure at model levels (Pa)
         137             :    real(r8), intent(in)    :: pdel(pcols,pver)     ! pressure thickness of levels (Pa)
         138             :    real(r8), intent(in)    :: qh2o(pcols,pver)     ! water vapor mixing ratio (kg/kg)
         139             :    real(r8), intent(in)    :: dgncur_a(pcols,pver,ntot_amode)
         140             :    real(r8), intent(in)    :: dgncur_awet(pcols,pver,ntot_amode)
         141             :    real(r8), pointer       :: sulfeq(:,:,:)
         142             :
         143             :                                  ! dry & wet geo. mean dia. (m) of number distrib.
         144             :
         145             : ! !DESCRIPTION:
         146             : ! computes TMR (tracer mixing ratio) tendencies for gas condensation
         147             : !    onto aerosol particles
         148             : !
         149             : ! this version does condensation of H2SO4, NH3, and MSA, both treated as
         150             : ! completely non-volatile (gas --> aerosol, but no aerosol --> gas)
         151             : !    gas H2SO4 goes to aerosol SO4
         152             : !    gas MSA (if present) goes to aerosol SO4
         153             : !       aerosol MSA is not distinguished from aerosol SO4
         154             : !    gas NH3 (if present) goes to aerosol NH4
         155             : !       if gas NH3 is not present, then ????
         156             : !
         157             : !
         158             : ! !REVISION HISTORY:
         159             : !   RCE 07.04.13:  Adapted from MIRAGE2 code
         160             : !
         161             : !EOP
         162             : !----------------------------------------------------------------------
         163             : !BOC
         164             :
         165             : ! local variables
         166             :    integer, parameter :: jsrflx_gaexch = 1
         167             :    integer, parameter :: jsrflx_rename = 2
         168             :    integer, parameter :: ldiag1=-1, ldiag2=-1, ldiag3=-1, ldiag4=-1
         169             :    integer, parameter :: method_soa = 2
         170             : !     method_soa=0 is no uptake
         171             : !     method_soa=1 is irreversible uptake done like h2so4 uptake
         172             : !     method_soa=2 is reversible uptake using subr modal_aero_soaexch
         173             :
         174             :    integer :: i, iq, itmpa
         175             :    integer :: idiagss
         176     2978352 :    integer :: ido_so4a(ntot_amode), ido_nh4a(ntot_amode)
         177     2978352 :    integer ::  ido_soaa(ntot_amode,nsoa)
         178             :    integer :: j, jac, jsrf, jsoa
         179             :    integer :: k,p
         180             :    integer :: l, l2, lb, lsfrm, lstoo
         181             :    integer :: l_so4g, l_nh4g, l_msag
         182     2978352 :    integer :: l_soag(nsoa)
         183             :    integer :: n, nn, niter, niter_max, ntot_soamode
         184             :
         185     2978352 :    logical :: is_dorename_atik, dorename_atik(ncol,pver)
         186             :
         187             :    character(len=fieldname_len+3) :: fieldname
         188             :    character(len=100) :: msg ! string for endrun calls
         189             :    character(len=32) :: spec_type
         190             :
         191     2978352 :    real (r8) :: avg_uprt_nh4, avg_uprt_so4, avg_uprt_soa(nsoa)
         192             :    real (r8) :: deltatxx
         193     2978352 :    real (r8) :: dqdt_nh4(ntot_amode), dqdt_so4(ntot_amode)
         194     2978352 :    real (r8) :: dqdt_soa(ntot_amode,nsoa)
         195     2978352 :    real (r8) :: dqdt_soag(nsoa)
         196             :    real (r8) :: fac_volsfc_pcarbon
         197     2978352 :    real (r8) :: fgain_nh4(ntot_amode), fgain_so4(ntot_amode)
         198     2978352 :    real (r8) :: fgain_soa(ntot_amode,nsoa)
         199             :    real (r8) :: g0_soa(nsoa)
         200     2978352 :    real(r8)  :: mw_poa_host(npoa)          ! molec wght of poa used in host code
         201     2978352 :    real(r8)  :: mw_soa_host(nsoa)          ! molec wght of poa used in host code
         202             :    real (r8) :: pdel_fac
         203             :    real (r8) :: qmax_nh4, qnew_nh4, qnew_so4
         204     2978352 :    real (r8) :: qold_nh4(ntot_amode), qold_so4(ntot_amode)
         205     2978352 :    real (r8) :: qold_poa(ntot_amode,npoa)
         206     2978352 :    real (r8) :: qold_soa(ntot_amode,nsoa)
         207     2978352 :    real (r8) :: qold_soag(nsoa)
         208             :    real (r8) :: sum_dqdt_msa, sum_dqdt_so4
         209     2978352 :    real (r8) :: sum_dqdt_soa(nsoa)
         210             :    real (r8) :: sum_dqdt_nh4, sum_dqdt_nh4_b
         211             :    real (r8) :: sum_uprt_msa, sum_uprt_nh4, sum_uprt_so4
         212     2978352 :    real (r8) :: sum_uprt_soa(nsoa)
         213             :    real (r8) :: tmp1, tmp2, tmpa
         214             :    real (r8) :: tmp_kxt, tmp_pxt
         215             :    real (r8) :: tmp_so4a_bgn, tmp_so4a_end
         216             :    real (r8) :: tmp_so4g_avg, tmp_so4g_bgn, tmp_so4g_equ
         217     2978352 :    real (r8) :: uptkrate(ntot_amode,pcols,pver)
         218     2978352 :    real (r8) :: uptkratebb(ntot_amode)
         219     2978352 :    real (r8) :: uptkrate_soa(ntot_amode,nsoa)
         220             :                 ! gas-to-aerosol mass transfer rates (1/s)
         221             :    real (r8) :: vol_core, vol_shell
         222             :    real (r8) :: xferfrac_pcage, xferfrac_max
         223             :    real (r8) :: xferrate
         224             :
         225             :    logical  :: do_msag         ! true if msa gas is a species
         226             :    logical  :: do_nh4g         ! true if nh3 gas is a species
         227             :    logical  :: do_soag_any         ! true if soa gas is a species
         228     2978352 :    logical  :: do_soag(nsoa)       ! true if soa gas is a species
         229             :
         230             :    logical  :: dotend(pcnstxx)          ! identifies species directly involved in
         231             :                                         !    gas-aerosol exchange (gas condensation)
         232             :    logical  :: dotendqqcw(pcnstxx)      ! like dotend but for cloud-borner tracers
         233             :    logical  :: dotendrn(pcnstxx), dotendqqcwrn(pcnstxx)
         234             :                                         ! identifies species involved in renaming
         235             :                                         !    after "continuous growth"
         236             :                                         !    (gas-aerosol exchange and aqchem)
         237             :
         238             :    integer, parameter :: nsrflx = 2     ! last dimension of qsrflx
         239     2978352 :    real(r8) :: dqdt(ncol,pver,pcnstxx)  ! TMR "delta q" array - NOTE dims
         240     2978352 :    real(r8) :: dqqcwdt(ncol,pver,pcnstxx) ! like dqdt but for cloud-borner tracers
         241             :    real(r8) :: qsrflx(pcols,pcnstxx,nsrflx)
         242             :                               ! process-specific column tracer tendencies
         243             :                               ! (1=renaming, 2=gas condensation)
         244             :    real(r8) :: qconff(pcols,pver),qevapff(pcols,pver)
         245             :    real(r8) :: qconbb(pcols,pver),qevapbb(pcols,pver)
         246             :    real(r8) :: qconbg(pcols,pver),qevapbg(pcols,pver)
         247             :    real(r8) :: qcon(pcols,pver),qevap(pcols,pver)
         248             :
         249             :    real(r8) :: qqcwsrflx(pcols,pcnstxx,nsrflx)
         250             :
         251             : !  following only needed for diagnostics
         252     2978352 :    real(r8) :: qold(ncol,pver,pcnstxx)  ! NOTE dims
         253             :    real(r8) :: qnew(ncol,pver,pcnstxx)  ! NOTE dims
         254             :    real(r8) :: qdel(ncol,pver,pcnstxx)  ! NOTE dims
         255             :    real(r8) :: dumavec(1000), dumbvec(1000), dumcvec(1000)
         256     2978352 :    real(r8) :: qqcwold(ncol,pver,pcnstxx)
         257     2978352 :    real(r8) :: dqdtsv1(ncol,pver,pcnstxx)
         258     2978352 :    real(r8) :: dqqcwdtsv1(ncol,pver,pcnstxx)
         259             :
         260             :
         261             : !----------------------------------------------------------------------
         262             :
         263             : ! set gas species indices
         264     1489176 :    call cnst_get_ind( 'H2SO4', l_so4g, .false. )
         265     1489176 :    call cnst_get_ind( 'NH3',   l_nh4g, .false. )
         266     1489176 :    if ( .not. cam_chempkg_is('geoschem_mam4') ) then
         267     1489176 :       call cnst_get_ind( 'MSA',   l_msag, .false. )
         268             :    else
         269           0 :       l_msag = 0
         270             :    endif
         271     1489176 :    l_so4g = l_so4g - loffset
         272     1489176 :    l_nh4g = l_nh4g - loffset
         273     1489176 :    l_msag = l_msag - loffset
         274     1489176 :    if ((l_so4g <= 0) .or. (l_so4g > pcnstxx)) then
         275             :       write( *, '(/a/a,2i7)' )   &
         276           0 :          '*** modal_aero_gasaerexch_sub -- cannot find H2SO4 species',   &
         277           0 :          '    l_so4g, loffset =', l_so4g, loffset
         278           0 :       call endrun( 'modal_aero_gasaerexch_sub error' )
         279             :    end if
         280     1489176 :    do_nh4g = .false.
         281     1489176 :    do_msag = .false.
         282     1489176 :    if ((l_nh4g > 0) .and. (l_nh4g <= pcnstxx)) do_nh4g = .true.
         283     1489176 :    if ((l_msag > 0) .and. (l_msag <= pcnstxx)) do_msag = .true.
         284             :
         285     1489176 :    do_soag_any = .false.
         286     2978352 :    do_soag(:) = .false.
         287     2978352 :    do jsoa = 1, nsoa
         288     1489176 :       l_soag(jsoa) = lptr2_soa_g_amode(jsoa) - loffset
         289     1489176 :       if ((method_soa == 1) .or. (method_soa == 2)) then
         290     1489176 :          if ((l_soag(jsoa) > 0) .and. (l_soag(jsoa) <= pcnstxx)) then
         291     1489176 :             do_soag_any = .true.
         292     1489176 :             do_soag(jsoa) = .true.
         293             :          end if
         294             :       else if (method_soa /= 0) then
         295             :          write(*,'(/a,1x,i10)') '*** modal_aero_gasaerexch_sub - bad method_soa =', method_soa
         296             :          call endrun( 'modal_aero_gasaerexch_sub error' )
         297             :       end if
         298             :    end do ! jsoa
         299             :
         300             : ! set tendency flags
         301     1489176 :    dotend(:) = .false.
         302     1489176 :    dotendqqcw(:) = .false.
         303     7445880 :    ido_so4a(:) = 0
         304     7445880 :    ido_nh4a(:) = 0
         305     8935056 :    ido_soaa(:,:) = 0
         306             :
         307     1489176 :    dotend(l_so4g) = .true.
         308     1489176 :    if ( do_nh4g ) dotend(l_nh4g) = .true.
         309     1489176 :    if ( do_msag ) dotend(l_msag) = .true.
         310     2978352 :    do jsoa = 1, nsoa
         311     2978352 :       if ( do_soag(jsoa) ) dotend(l_soag(jsoa)) = .true.
         312             :    end do
         313             :
         314     1489176 :    ntot_soamode = 0
         315     7445880 :    do n = 1, ntot_amode
         316     5956704 :       l = lptr_so4_a_amode(n)-loffset
         317     5956704 :       if ((l > 0) .and. (l <= pcnstxx)) then
         318     4467528 :          dotend(l) = .true.
         319     4467528 :          ido_so4a(n) = 1
         320     4467528 :          if ( do_nh4g ) then
         321           0 :             l = lptr_nh4_a_amode(n)-loffset
         322           0 :             if ((l > 0) .and. (l <= pcnstxx)) then
         323           0 :                dotend(l) = .true.
         324           0 :                ido_nh4a(n) = 1
         325             :             end if
         326             :          end if
         327             :       end if
         328             :
         329    13402584 :       do jsoa = 1, nsoa
         330    11913408 :          if ( do_soag(jsoa) ) then
         331     5956704 :             l = lptr2_soa_a_amode(n,jsoa)-loffset
         332     5956704 :             if ((l > 0) .and. (l <= pcnstxx)) then
         333     2978352 :                dotend(l) = .true.
         334     2978352 :                ido_soaa(n,jsoa) = 1
         335     2978352 :                ntot_soamode = n
         336             :             end if
         337             :          end if
         338             :       end do ! jsoa
         339             :    end do ! n
         340             :
         341             :
         342     1489176 :    if ( do_soag_any ) ntot_soamode = max( ntot_soamode, modefrm_pcage )
         343             :
         344     1489176 :    if (modefrm_pcage > 0) then
         345     1489176 :       ido_so4a(modefrm_pcage) = 2
         346     1489176 :       if (ido_nh4a(modetoo_pcage) == 1) ido_nh4a(modefrm_pcage) = 2
         347     2978352 :       do jsoa = 1, nsoa
         348     2978352 :          if (ido_soaa(modetoo_pcage,jsoa) == 1) ido_soaa(modefrm_pcage,jsoa) = 2
         349             :       end do
         350     5956704 :       do iq = 1, nspecfrm_pcage
         351     4467528 :          lsfrm = lspecfrm_pcage(iq)-loffset
         352     4467528 :          lstoo = lspectoo_pcage(iq)-loffset
         353     5956704 :          if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
         354     4467528 :             dotend(lsfrm) = .true.
         355     4467528 :             if ((lstoo > 0) .and. (lstoo <= pcnst)) then
         356     4467528 :                dotend(lstoo) = .true.
         357             :             end if
         358             :          end if
         359             :       end do
         360             :
         361             :
         362     1489176 :       n = modeptr_pcarbon
         363     1489176 :       fac_volsfc_pcarbon = exp( 2.5_r8*(alnsg_amode(n)**2) )
         364     1489176 :       xferfrac_max = 1.0_r8 - 10.0_r8*epsilon(1.0_r8)   ! 1-eps
         365             :    end if
         366             :
         367             :
         368             : ! zero out tendencies and other
         369 71735685840 :    dqdt(:,:,:) = 0.0_r8
         370 71735685840 :    dqqcwdt(:,:,:) = 0.0_r8
         371     1489176 :    qsrflx(:,:,:) = 0.0_r8
         372     1489176 :    qqcwsrflx(:,:,:) = 0.0_r8
         373             :
         374             : !-------Initialize evap/cond diagnostics (ncols x pver)-----------
         375     1489176 :    qconff(:,:) = 0.0_r8
         376     1489176 :    qevapff(:,:) = 0.0_r8
         377     1489176 :    qconbb(:,:) = 0.0_r8
         378     1489176 :    qevapbb(:,:) = 0.0_r8
         379     1489176 :    qconbg(:,:) = 0.0_r8
         380     1489176 :    qevapbg(:,:) = 0.0_r8
         381     1489176 :    qcon(:,:) = 0.0_r8
         382     1489176 :    qevap(:,:) = 0.0_r8
         383             : !---------------------------------------------------
         384             :
         385             : ! compute gas-to-aerosol mass transfer rates
         386             :    call gas_aer_uptkrates( ncol,       loffset,                &
         387             :                            q,          t,          pmid,       &
         388     1489176 :                            dgncur_awet,            uptkrate    )
         389             :
         390             :
         391             : ! use this for tendency calcs to avoid generating very small negative values
         392     1489176 :    deltatxx = deltat * (1.0_r8 + 1.0e-15_r8)
         393             :
         394             :
         395     1489176 :    jsrf = jsrflx_gaexch
         396   139982544 :    do k=top_lev,pver
         397  2314006344 :       do i=1,ncol
         398             :
         399             : !   fgain_so4(n) = fraction of total h2so4 uptake going to mode n
         400             : !   fgain_nh4(n) = fraction of total  nh3  uptake going to mode n
         401  4348047600 :          sum_uprt_so4 = 0.0_r8
         402  4348047600 :          sum_uprt_nh4 = 0.0_r8
         403  4348047600 :          sum_uprt_soa = 0.0_r8
         404 10870119000 :          do n = 1, ntot_amode
         405  8696095200 :             uptkratebb(n) = uptkrate(n,i,k)
         406  8696095200 :             if (ido_so4a(n) > 0) then
         407  8696095200 :                fgain_so4(n) = uptkratebb(n)
         408  8696095200 :                sum_uprt_so4 = sum_uprt_so4 + fgain_so4(n)
         409  8696095200 :                if (ido_so4a(n) == 1) then
         410  6522071400 :                   qold_so4(n) = q(i,k,lptr_so4_a_amode(n)-loffset)
         411             :                else
         412  2174023800 :                   qold_so4(n) = 0.0_r8
         413             :                end if
         414             :             else
         415           0 :                fgain_so4(n) = 0.0_r8
         416           0 :                qold_so4(n) = 0.0_r8
         417             :             end if
         418             :
         419  8696095200 :             if (ido_nh4a(n) > 0) then
         420             :                !   2.08 factor is for gas diffusivity (nh3/h2so4)
         421             :                !   differences in fuch-sutugin and accom coef ignored
         422           0 :                fgain_nh4(n) = uptkratebb(n)*2.08_r8
         423           0 :                sum_uprt_nh4 = sum_uprt_nh4 + fgain_nh4(n)
         424           0 :                if (ido_nh4a(n) == 1) then
         425           0 :                   qold_nh4(n) = q(i,k,lptr_nh4_a_amode(n)-loffset)
         426             :                else
         427           0 :                   qold_nh4(n) = 0.0_r8
         428             :                end if
         429             :             else
         430  8696095200 :                fgain_nh4(n) = 0.0_r8
         431  8696095200 :                qold_nh4(n) = 0.0_r8
         432             :             end if
         433             :
         434 17392190400 :             do j = 1, npoa
         435  8696095200 :                l = lptr2_pom_a_amode(n,j)-loffset
         436 17392190400 :                if (l > 0) then
         437  4348047600 :                   qold_poa(n,j) = q(i,k,l)
         438             :                else
         439  4348047600 :                   qold_poa(n,j) = 0.0_r8
         440             :                end if
         441             :             end do
         442             :
         443  8696095200 :             itmpa = 0
         444 17392190400 :             do jsoa = 1, nsoa
         445  8696095200 :                if (ido_soaa(n,jsoa) > 0) then
         446             :                   ! 0.81 factor is for gas diffusivity (soa/h2so4)
         447             :                   ! (differences in fuch-sutugin and accom coef ignored)
         448  6522071400 :                   fgain_soa(n,jsoa) = uptkratebb(n)*0.81_r8
         449  6522071400 :                   sum_uprt_soa(jsoa) = sum_uprt_soa(jsoa) + fgain_soa(n,jsoa)
         450  6522071400 :                   if (ido_soaa(n,jsoa) == 1) then
         451  4348047600 :                      l = lptr2_soa_a_amode(n,jsoa)-loffset
         452  4348047600 :                      qold_soa(n,jsoa) = q(i,k,l)
         453  4348047600 :                      itmpa = itmpa + 1
         454             :                   else
         455  2174023800 :                      qold_soa(n,jsoa) = 0.0_r8
         456             :                   end if
         457             :                else
         458  2174023800 :                   fgain_soa(n,jsoa) = 0.0_r8
         459  2174023800 :                   qold_soa(n,jsoa) = 0.0_r8
         460             :                end if
         461 17392190400 :                uptkrate_soa(n,jsoa) = fgain_soa(n,jsoa)
         462             :             end do ! jsoa
         463             :             ! in previous code versions with nsoa=1,
         464             :             !    qold_poa was non-zero (i.e., loaded from q) only when ido_soaa(n)=1
         465             :             ! thus qold_poa=0 for the primary carbon mode which has ido_soaa=2
         466             :             ! this is probably not how it should be
         467 15218166600 :             if (itmpa == 0) qold_poa(n,:) = 0.0_r8
         468             :
         469             :          end do ! n
         470             :
         471  2174023800 :          if (sum_uprt_so4 > 0.0_r8) then
         472 10870119000 :             do n = 1, ntot_amode
         473 10870119000 :                fgain_so4(n) = fgain_so4(n) / sum_uprt_so4
         474             :             end do
         475             :          end if
         476             : !       at this point (sum_uprt_so4 <= 0.0) only when all the fgain_so4 are zero
         477  2174023800 :          if (sum_uprt_nh4 > 0.0_r8) then
         478           0 :             do n = 1, ntot_amode
         479           0 :                fgain_nh4(n) = fgain_nh4(n) / sum_uprt_nh4
         480             :             end do
         481             :          end if
         482             :
         483  4348047600 :          do jsoa = 1, nsoa
         484  4348047600 :             if (sum_uprt_soa(jsoa) > 0.0_r8) then
         485 10870119000 :                do n = 1, ntot_amode
         486 10870119000 :                   fgain_soa(n,jsoa) = fgain_soa(n,jsoa) / sum_uprt_soa(jsoa)
         487             :                end do
         488             :             end if
         489             :          end do
         490             :
         491             : !   uptake amount (fraction of gas uptaken) over deltat
         492  2174023800 :          avg_uprt_so4 = (1.0_r8 - exp(-deltatxx*sum_uprt_so4))/deltatxx
         493  2174023800 :          avg_uprt_nh4 = (1.0_r8 - exp(-deltatxx*sum_uprt_nh4))/deltatxx
         494             :
         495  4348047600 :          do jsoa = 1, nsoa
         496  4348047600 :             avg_uprt_soa(jsoa) = (1.0_r8 - exp(-deltatxx*sum_uprt_soa(jsoa)))/deltatxx
         497             :          end do
         498             :
         499             : !   sum_dqdt_so4 = so4_a tendency from h2so4 gas uptake (mol/mol/s)
         500             : !   sum_dqdt_msa = msa_a tendency from msa   gas uptake (mol/mol/s)
         501             : !   sum_dqdt_nh4 = nh4_a tendency from nh3   gas uptake (mol/mol/s)
         502             : !   sum_dqdt_soa = soa_a tendency from soa   gas uptake (mol/mol/s)
         503  2174023800 :          sum_dqdt_so4 = q(i,k,l_so4g) * avg_uprt_so4
         504  2174023800 :          if ( do_msag ) then
         505           0 :             sum_dqdt_msa = q(i,k,l_msag) * avg_uprt_so4
         506             :          else
         507             :             sum_dqdt_msa = 0.0_r8
         508             :          end if
         509  2174023800 :          if ( do_nh4g ) then
         510           0 :             sum_dqdt_nh4 = q(i,k,l_nh4g) * avg_uprt_nh4
         511             :          else
         512             :             sum_dqdt_nh4 = 0.0_r8
         513             :          end if
         514             :
         515  4348047600 :          do jsoa = 1, nsoa
         516  4348047600 :             if ( do_soag(jsoa) ) then
         517  2174023800 :                sum_dqdt_soa(jsoa) = q(i,k,l_soag(jsoa)) * avg_uprt_soa(jsoa)
         518             :             else
         519           0 :                sum_dqdt_soa(jsoa) = 0.0_r8
         520             :             end if
         521             :          end do
         522             :
         523  2174023800 :          if ( associated(sulfeq) .and. (k <= troplev(i)) ) then
         524             :             !   compute TMR tendencies for so4 interstial aerosol due to reversible gas uptake
         525             :             !   only above the tropopause
         526             :
         527           0 :             tmp_kxt = deltatxx*sum_uprt_so4  ! sum over modes of uptake_rate*deltat
         528           0 :             tmp_pxt = 0.0_r8
         529           0 :             do n = 1, ntot_amode
         530           0 :                if (ido_so4a(n) <= 0) cycle
         531           0 :                tmp_pxt = tmp_pxt + uptkratebb(n)*sulfeq(i,k,n)
         532             :             end do
         533           0 :             tmp_pxt = max( 0.0_r8, tmp_pxt*deltatxx )  ! sum over modes of uptake_rate*sulfeq*deltat
         534           0 :             tmp_so4g_bgn = q(i,k,l_so4g)
         535             :             ! calc avg h2so4(g) over deltat
         536           0 :             if (tmp_kxt >= 1.0e-5_r8) then
         537             :                ! exponential decay towards equilibrium value solution
         538           0 :                tmp_so4g_equ = tmp_pxt/tmp_kxt
         539           0 :                tmp_so4g_avg = tmp_so4g_equ + (tmp_so4g_bgn-tmp_so4g_equ)*(1.0_r8-exp(-tmp_kxt))/tmp_kxt
         540             :             else
         541             :                ! first order approx for tmp_kxt small
         542           0 :                tmp_so4g_avg = tmp_so4g_bgn*(1.0_r8-0.5_r8*tmp_kxt) + 0.5_r8*tmp_pxt
         543             :             end if
         544           0 :             sum_dqdt_so4 = 0.0_r8
         545           0 :             do n = 1, ntot_amode
         546           0 :                if (ido_so4a(n) <= 0) cycle
         547             :                ! calc change to so4(a) in mode n
         548           0 :                if (ido_so4a(n) == 1) then
         549           0 :                   l = lptr_so4_a_amode(n)-loffset
         550           0 :                   tmp_so4a_bgn = q(i,k,l)
         551             :                else
         552             :                   tmp_so4a_bgn = 0.0_r8
         553             :                end if
         554           0 :                tmp_so4a_end = tmp_so4a_bgn + deltatxx*uptkratebb(n)*(tmp_so4g_avg-sulfeq(i,k,n))
         555           0 :                tmp_so4a_end = max( 0.0_r8, tmp_so4a_end )
         556           0 :                dqdt_so4(n) = (tmp_so4a_end - tmp_so4a_bgn)/deltatxx
         557           0 :                sum_dqdt_so4 = sum_dqdt_so4 + dqdt_so4(n)
         558             :             end do
         559             :             ! do not allow msa condensation in stratosphere
         560             :             ! ( Note that the code for msa has never been used.
         561             :             !   The plan was to simulate msa(g), treat it as non-volatile (like h2so4(g)),
         562             :             !   and treat condensed msa as sulfate, so just one additional tracer. )
         563           0 :             if ( do_msag ) sum_dqdt_msa = 0.0_r8
         564             :
         565             :          else
         566             :             !   compute TMR tendencies for so4 interstial aerosol due to simple gas uptake
         567 10870119000 :             do n = 1, ntot_amode
         568 10870119000 :                dqdt_so4(n) = fgain_so4(n)*(sum_dqdt_so4 + sum_dqdt_msa)
         569             :             end do
         570             :          end if
         571             :
         572             :          !   compute TMR tendencies for nh4 interstial aerosol due to simple gas uptake
         573             :          !   but force nh4/so4 molar ratio <= 2
         574  2174023800 :          sum_dqdt_nh4_b = 0.0_r8
         575 10870119000 :          dqdt_nh4(:) = 0._r8
         576  2174023800 :          if ( do_nh4g ) then
         577           0 :             do n = 1, ntot_amode
         578           0 :                dqdt_nh4(n) = fgain_nh4(n)*sum_dqdt_nh4
         579           0 :                qnew_nh4 = qold_nh4(n) + dqdt_nh4(n)*deltat
         580           0 :                qnew_so4 = qold_so4(n) + dqdt_so4(n)*deltat
         581           0 :                qmax_nh4 = 2.0_r8*qnew_so4
         582           0 :                if (qnew_nh4 > qmax_nh4) then
         583           0 :                   dqdt_nh4(n) = (qmax_nh4 - qold_nh4(n))/deltatxx
         584             :                end if
         585           0 :                sum_dqdt_nh4_b = sum_dqdt_nh4_b + dqdt_nh4(n)
         586             :             end do
         587             :          end if
         588             :
         589  2174023800 :          if (( do_soag_any ) .and. (method_soa > 1)) then
         590             : !   compute TMR tendencies for soag and soa interstial aerosol
         591             : !   using soa parameterization
         592  2174023800 :             niter_max = 1000
         593 13044142800 :             dqdt_soa(:,:) = 0.0_r8
         594  4348047600 :             dqdt_soag(:) = 0.0_r8
         595  4348047600 :             do jsoa = 1, nsoa
         596  4348047600 :                qold_soag(jsoa) = q(i,k,l_soag(jsoa))
         597             :             end do
         598             :
         599             :             ! get molecular weight from the host model
         600 10870119000 :             do n = 1, ntot_amode
         601 43480476000 :               do l = 1, nspec_amode(n)
         602 32610357000 :                   call rad_cnst_get_info(0, n, l, spec_type=spec_type )
         603  8696095200 :                   select case( spec_type )
         604             :                    case('s-organic')
         605  8696095200 :                      mw_soa_host(:) = specmw_amode(l,n)
         606             :                    case('p-organic')
         607 36958404600 :                      mw_poa_host(:) = specmw_amode(l,n)
         608             :                    end select
         609             :                end do
         610             :             end do
         611             :
         612             :             call modal_aero_soaexch( deltat, t(i,k), pmid(i,k), &
         613             :                  niter, niter_max, ntot_amode, ntot_soamode, npoa, nsoa, &
         614             :                  mw_poa_host, mw_soa_host, &
         615             :                  qold_soag, qold_soa, qold_poa, uptkrate_soa, &
         616  2174023800 :                  dqdt_soag, dqdt_soa )
         617  4348047600 :             sum_dqdt_soa(:) = -dqdt_soag(:)
         618             :
         619             :          else if ( do_soag_any ) then
         620             : !   compute TMR tendencies for soa interstial aerosol
         621             : !   due to simple gas uptake
         622             :
         623             :             do jsoa = 1, nsoa
         624             :                do n = 1, ntot_amode
         625             :                   dqdt_soa(n,jsoa) = fgain_soa(n,jsoa)*sum_dqdt_soa(jsoa)
         626             :                end do
         627             :             end do
         628             :          else ! method_soa is neither 1 nor 2, no uptake
         629           0 :             dqdt_soa(:,:) = 0.0_r8
         630             :          end if
         631             :
         632  2174023800 :          pdel_fac = pdel(i,k)/gravit
         633 10870119000 :          do n = 1, ntot_amode
         634  8696095200 :             if (ido_so4a(n) == 1) then
         635  6522071400 :                l = lptr_so4_a_amode(n)-loffset
         636  6522071400 :                dqdt(i,k,l) = dqdt_so4(n)
         637  6522071400 :                qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_so4(n)*pdel_fac
         638             :             end if
         639             :
         640  8696095200 :             if ( do_nh4g ) then
         641           0 :                if (ido_nh4a(n) == 1) then
         642           0 :                   l = lptr_nh4_a_amode(n)-loffset
         643           0 :                   dqdt(i,k,l) = dqdt_nh4(n)
         644           0 :                   qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_nh4(n)*pdel_fac
         645             :                end if
         646             :             end if
         647             :
         648 19566214200 :             do jsoa = 1, nsoa
         649 17392190400 :                if ( do_soag(jsoa) ) then
         650  8696095200 :                   if (ido_soaa(n,jsoa) == 1) then
         651  4348047600 :                      l = lptr2_soa_a_amode(n,jsoa)-loffset
         652  4348047600 :                      dqdt(i,k,l) = dqdt_soa(n,jsoa) !calculated by  modal_aero_soaexch for method_soa=2
         653  4348047600 :                      qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_soa(n,jsoa)*pdel_fac
         654             : !------- Add code for condensation/evaporation diagnostics---
         655  4348047600 :                      if (nsoa.eq.15) then !check for current SOA package
         656           0 :                         if(jsoa.ge.1.and.jsoa.le.5) then ! Fossil SOA species
         657           0 :                            if (dqdt_soa(n,jsoa).ge.0.0_r8) then
         658           0 :                               qconff(i,k)=qconff(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
         659           0 :                            elseif(dqdt_soa(n,jsoa).lt.0.0_r8) then
         660           0 :                               qevapff(i,k)=qevapff(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
         661             :                            endif
         662             :
         663           0 :                         elseif(jsoa.ge.6.and.jsoa.le.10) then ! Biomass SOA species
         664           0 :                            if (dqdt_soa(n,jsoa).ge.0.0_r8) then
         665           0 :                               qconbb(i,k)=qconbb(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
         666           0 :                            elseif(dqdt_soa(n,jsoa).lt.0.0_r8) then
         667           0 :                               qevapbb(i,k)=qevapbb(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
         668             :                            endif
         669             :
         670           0 :                         elseif(jsoa.ge.11.and.jsoa.le.15) then ! Biomass SOA species
         671           0 :                            if (dqdt_soa(n,jsoa).ge.0.0_r8) then
         672           0 :                               qconbg(i,k)=qconbg(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
         673           0 :                            elseif(dqdt_soa(n,jsoa).lt.0.0_r8) then
         674           0 :                               qevapbg(i,k)=qevapbg(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
         675             :                            endif
         676             :
         677             :                         endif ! jsoa
         678             :                      endif !nsoa
         679  4348047600 :                      if (nsoa.eq.5) then !check for current SOA package
         680           0 :                            if (dqdt_soa(n,jsoa).ge.0.0_r8) then
         681           0 :                               qcon(i,k)=qcon(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
         682           0 :                            elseif(dqdt_soa(n,jsoa).lt.0.0_r8) then
         683           0 :                               qevap(i,k)=qevap(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
         684             :                            endif
         685             :                      endif !nsoa
         686             : !---------------------------------------------------------------------------------------------------------------------
         687             :                   end if
         688             :                end if
         689             :             end do
         690             :          end do ! n
         691             :
         692             : !   compute TMR tendencies for h2so4, nh3, and msa gas
         693             : !   due to simple gas uptake
         694  2174023800 :          l = l_so4g
         695  2174023800 :          dqdt(i,k,l) = -sum_dqdt_so4
         696  2174023800 :          qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
         697             :
         698  2174023800 :          if ( do_msag ) then
         699           0 :             l = l_msag
         700           0 :             dqdt(i,k,l) = -sum_dqdt_msa
         701           0 :             qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
         702             :          end if
         703             :
         704  2174023800 :          if ( do_nh4g ) then
         705           0 :             l = l_nh4g
         706           0 :             dqdt(i,k,l) = -sum_dqdt_nh4_b
         707           0 :             qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
         708             :          end if
         709             :
         710  4348047600 :          do jsoa = 1, nsoa
         711  4348047600 :             if ( do_soag(jsoa) ) then
         712  2174023800 :                l = l_soag(jsoa)
         713  2174023800 :                dqdt(i,k,l) = -sum_dqdt_soa(jsoa)
         714             : ! dqdt for gas is negative of the sum of dqdt for aerosol soa species in each mode: Manish
         715  2174023800 :                qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
         716             :             end if
         717             :          end do
         718             :
         719             : !   compute TMR tendencies associated with primary carbon aging
         720  2312517168 :          if (modefrm_pcage > 0) then
         721  2174023800 :             n = modeptr_pcarbon
         722  2174023800 :             tmpa = 0.0_r8
         723  4348047600 :             do jsoa = 1, nsoa
         724  4348047600 :                tmpa = tmpa + dqdt_soa(n,jsoa)*fac_m2v_soa(jsoa)*soa_equivso4_factor(jsoa)
         725             :             end do
         726             :             vol_shell = deltat *   &
         727  2174023800 :                  ( dqdt_so4(n)*fac_m2v_so4 + dqdt_nh4(n)*fac_m2v_nh4 + tmpa )
         728  2174023800 :             vol_core = 0.0_r8
         729  6522071400 :             do l = 1, nspec_amode(n)
         730             :                vol_core = vol_core + &
         731  6522071400 :                     q(i,k,lmassptr_amode(l,n)-loffset)*fac_m2v_pcarbon(l)
         732             :             end do
         733             : !   ratio1 = vol_shell/vol_core =
         734             : !      actual hygroscopic-shell-volume/carbon-core-volume after gas uptake
         735             : !   ratio2 = 6.0_r8*dr_so4_monolayers_pcage/(dgncur_a*fac_volsfc_pcarbon)
         736             : !      = (shell-volume corresponding to n_so4_monolayers_pcage)/core-volume
         737             : !      The 6.0/(dgncur_a*fac_volsfc_pcarbon) = (mode-surface-area/mode-volume)
         738             : !   Note that vol_shell includes both so4+nh4 AND soa as "equivalent so4",
         739             : !      The soa_equivso4_factor accounts for the lower hygroscopicity of soa.
         740             : !
         741             : !   Define xferfrac_pcage = min( 1.0, ratio1/ratio2)
         742             : !   But ratio1/ratio2 == tmp1/tmp2, and coding below avoids possible overflow
         743             : !
         744  2174023800 :             tmp1 = vol_shell*dgncur_a(i,k,n)*fac_volsfc_pcarbon
         745  2174023800 :             tmp2 = max( 6.0_r8*dr_so4_monolayers_pcage*vol_core, 0.0_r8 )
         746  2174023800 :             if (tmp1 >= tmp2) then
         747             :                xferfrac_pcage = xferfrac_max
         748             :             else
         749  2173297671 :                xferfrac_pcage = min( tmp1/tmp2, xferfrac_max )
         750             :             end if
         751             :
         752  2174023800 :             if (xferfrac_pcage > 0.0_r8) then
         753  8696095200 :                do iq = 1, nspecfrm_pcage
         754  6522071400 :                   lsfrm = lspecfrm_pcage(iq)-loffset
         755  6522071400 :                   lstoo = lspectoo_pcage(iq)-loffset
         756  6522071400 :                   xferrate = (xferfrac_pcage/deltat)*q(i,k,lsfrm)
         757  6522071400 :                   dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferrate
         758  6522071400 :                   qsrflx(i,lsfrm,jsrf) = qsrflx(i,lsfrm,jsrf) - xferrate*pdel_fac
         759  8696095200 :                   if ((lstoo > 0) .and. (lstoo <= pcnst)) then
         760  6522071400 :                      dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferrate
         761  6522071400 :                      qsrflx(i,lstoo,jsrf) = qsrflx(i,lstoo,jsrf) + xferrate*pdel_fac
         762             :                   end if
         763             :                end do
         764             :
         765  2174023800 :                if (ido_so4a(modetoo_pcage) > 0) then
         766  2174023800 :                   l = lptr_so4_a_amode(modetoo_pcage)-loffset
         767  2174023800 :                   dqdt(i,k,l) = dqdt(i,k,l) + dqdt_so4(modefrm_pcage)
         768  2174023800 :                   qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_so4(modefrm_pcage)*pdel_fac
         769             :                end if
         770             :
         771  2174023800 :                if (ido_nh4a(modetoo_pcage) > 0) then
         772           0 :                   l = lptr_nh4_a_amode(modetoo_pcage)-loffset
         773           0 :                   dqdt(i,k,l) = dqdt(i,k,l) + dqdt_nh4(modefrm_pcage)
         774           0 :                   qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_nh4(modefrm_pcage)*pdel_fac
         775             :                end if
         776             :
         777  4348047600 :                do jsoa = 1, nsoa
         778  4348047600 :                   if (ido_soaa(modetoo_pcage,jsoa) > 0) then
         779  2174023800 :                      l = lptr2_soa_a_amode(modetoo_pcage,jsoa)-loffset
         780  2174023800 :                      dqdt(i,k,l) = dqdt(i,k,l) + dqdt_soa(modefrm_pcage,jsoa)
         781  2174023800 :                      qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_soa(modefrm_pcage,jsoa)*pdel_fac
         782             :                   end if
         783             :                end do
         784             :
         785             :             end if
         786             :
         787             :          end if
         788             :
         789             :       end do   ! "i = 1, ncol"
         790             :    end do     ! "k = top_lev, pver"
         791             :
         792             : ! set "temporary testing arrays"
         793 71735685840 :    qold(:,:,:) = q(:,:,:)
         794 71735685840 :    qqcwold(:,:,:) = qqcw(:,:,:)
         795 71735685840 :    dqdtsv1(:,:,:) = dqdt(:,:,:)
         796 71735685840 :    dqqcwdtsv1(:,:,:) = dqqcwdt(:,:,:)
         797             :
         798             :
         799             : !
         800             : ! do renaming calcs
         801             : !
         802     1489176 :    dotendrn(:) = .false.
         803     1489176 :    dotendqqcwrn(:) = .false.
         804  2314006344 :    dorename_atik(1:ncol,:) = .true.
         805     1489176 :    is_dorename_atik = .true.
         806             :    call modal_aero_rename_sub(                              &
         807             :         'modal_aero_gasaerexch_sub',            &
         808             :         lchnk,             ncol,      nstep,    &
         809             :         loffset,           deltat,              &
         810             :         pdel,              troplev,             &
         811             :         dotendrn,          q,                   &
         812             :         dqdt,              dqdt_other,          &
         813             :         dotendqqcwrn,      qqcw,                &
         814             :         dqqcwdt,           dqqcwdt_other,       &
         815             :         is_dorename_atik,  dorename_atik,       &
         816             :         jsrflx_rename,     nsrflx,              &
         817     1489176 :         qsrflx,            qqcwsrflx            )
         818             :
         819             :
         820             : !  This applies dqdt tendencies for all species
         821             : !  apply the dqdt to update q (and same for qqcw)
         822             : !
         823    47653632 :    do l = 1, pcnstxx
         824    46164456 :       if ( dotend(l) .or. dotendrn(l) ) then
         825  2939633424 :          do k = top_lev, pver
         826 48594133224 :             do i = 1, ncol
         827 48562860528 :                q(i,k,l) = q(i,k,l) + dqdt(i,k,l)*deltat
         828             :             end do
         829             :          end do
         830             :       end if
         831    47653632 :       if ( dotendqqcw(l) .or. dotendqqcwrn(l) ) then
         832  1959755616 :          do k = top_lev, pver
         833 32396088816 :             do i = 1, ncol
         834 32375240352 :                qqcw(i,k,l) = qqcw(i,k,l) + dqqcwdt(i,k,l)*deltat
         835             :             end do
         836             :          end do
         837             :       end if
         838             :    end do
         839             :
         840             : ! diagnostics start -------------------------------------------------------
         841             : !!$   if (ldiag3 > 0) then
         842             : !!$   if (icol_diag > 0) then
         843             : !!$      i = icol_diag
         844             : !!$      write(*,'(a,3i5)') 'gasaerexch ppp nstep,lat,lon', nstep, latndx(i), lonndx(i)
         845             : !!$      write(*,'(2i5,3(2x,a))') 0, 0, 'ppp', 'pdel for all k'
         846             : !!$      write(*,'(1p,7e12.4)') (pdel(i,k), k=top_lev,pver)
         847             : !!$
         848             : !!$      write(*,'(a,3i5)') 'gasaerexch ddd nstep,lat,lon', nstep, latndx(i), lonndx(i)
         849             : !!$      do l = 1, pcnstxx
         850             : !!$         lb = l + loffset
         851             : !!$
         852             : !!$         if ( dotend(l) .or. dotendrn(l) ) then
         853             : !!$            write(*,'(2i5,3(2x,a))') 1, l, 'ddd1', cnst_name(lb),    'qold for all k'
         854             : !!$            write(*,'(1p,7e12.4)') (qold(i,k,l), k=top_lev,pver)
         855             : !!$            write(*,'(2i5,3(2x,a))') 1, l, 'ddd2', cnst_name(lb),    'qnew for all k'
         856             : !!$            write(*,'(1p,7e12.4)') (q(i,k,l), k=top_lev,pver)
         857             : !!$            write(*,'(2i5,3(2x,a))') 1, l, 'ddd3', cnst_name(lb),    'dqdt from conden for all k'
         858             : !!$            write(*,'(1p,7e12.4)') (dqdtsv1(i,k,l), k=top_lev,pver)
         859             : !!$            write(*,'(2i5,3(2x,a))') 1, l, 'ddd4', cnst_name(lb),    'dqdt from rename for all k'
         860             : !!$            write(*,'(1p,7e12.4)') ((dqdt(i,k,l)-dqdtsv1(i,k,l)), k=top_lev,pver)
         861             : !!$            write(*,'(2i5,3(2x,a))') 1, l, 'ddd5', cnst_name(lb),    'dqdt other for all k'
         862             : !!$            write(*,'(1p,7e12.4)') (dqdt_other(i,k,l), k=top_lev,pver)
         863             : !!$         end if
         864             : !!$
         865             : !!$         if ( dotendqqcw(l) .or. dotendqqcwrn(l) ) then
         866             : !!$            write(*,'(2i5,3(2x,a))') 2, l, 'ddd1', cnst_name_cw(lb), 'qold for all k'
         867             : !!$            write(*,'(1p,7e12.4)') (qqcwold(i,k,l), k=top_lev,pver)
         868             : !!$            write(*,'(2i5,3(2x,a))') 2, l, 'ddd2', cnst_name_cw(lb), 'qnew for all k'
         869             : !!$            write(*,'(1p,7e12.4)') (qqcw(i,k,l), k=top_lev,pver)
         870             : !!$            write(*,'(2i5,3(2x,a))') 2, l, 'ddd3', cnst_name_cw(lb), 'dqdt from conden for all k'
         871             : !!$            write(*,'(1p,7e12.4)') (dqqcwdtsv1(i,k,l), k=top_lev,pver)
         872             : !!$            write(*,'(2i5,3(2x,a))') 2, l, 'ddd4', cnst_name_cw(lb), 'dqdt from rename for all k'
         873             : !!$            write(*,'(1p,7e12.4)') ((dqqcwdt(i,k,l)-dqqcwdtsv1(i,k,l)), k=top_lev,pver)
         874             : !!$            write(*,'(2i5,3(2x,a))') 2, l, 'ddd5', cnst_name_cw(lb), 'dqdt other for all k'
         875             : !!$            write(*,'(1p,7e12.4)') (dqqcwdt_other(i,k,l), k=top_lev,pver)
         876             : !!$         end if
         877             : !!$
         878             : !!$      end do
         879             : !!$
         880             : !!$      write(*,'(a,3i5)') 'gasaerexch fff nstep,lat,lon', nstep, latndx(i), lonndx(i)
         881             : !!$      do l = 1, pcnstxx
         882             : !!$         lb = l + loffset
         883             : !!$         if ( dotend(l) .or. dotendrn(l) .or. dotendqqcw(l) .or. dotendqqcwrn(l) ) then
         884             : !!$            write(*,'(i5,2(2x,a,2l3))') l, &
         885             : !!$               cnst_name(lb), dotend(l), dotendrn(l), &
         886             : !!$               cnst_name_cw(lb), dotendqqcw(l), dotendqqcwrn(l)
         887             : !!$         end if
         888             : !!$      end do
         889             : !!$
         890             : !!$   end if
         891             : !!$   end if
         892             : ! diagnostics end ---------------------------------------------------------
         893             :
         894             : !-----Outfld for condensation/evaporation------------------------------
         895     1489176 :    if (nsoa.eq.5) then !check for current SOA package
         896           0 :       call outfld(trim('qcon_gaex'), qcon(:,:), pcols, lchnk )
         897           0 :       call outfld(trim('qevap_gaex'), qevap(:,:), pcols, lchnk )
         898             :    endif
         899             : !-----------------------------------------------------------------------
         900     1489176 :    if (nsoa.eq.15) then !check for current SOA package
         901           0 :       call outfld(trim('qconff_gaex'), qconff(:,:), pcols, lchnk )
         902           0 :       call outfld(trim('qevapff_gaex'), qevapff(:,:), pcols, lchnk )
         903           0 :       call outfld(trim('qconbb_gaex'), qconbb(:,:), pcols, lchnk )
         904           0 :       call outfld(trim('qevapbb_gaex'), qevapbb(:,:), pcols, lchnk )
         905           0 :       call outfld(trim('qconbg_gaex'), qconbg(:,:), pcols, lchnk )
         906           0 :       call outfld(trim('qevapbg_gaex'), qevapbg(:,:), pcols, lchnk )
         907             :    endif
         908             : !-----------------------------------------------------------------------
         909             :
         910             : !   do history file column-tendency fields
         911    47653632 :    do l = 1, pcnstxx
         912    46164456 :       lb = l + loffset
         913   139982544 :       do jsrf = 1, 2
         914   323151192 :          do jac = 1, 2
         915   276986736 :             if (jac == 1) then
         916    92328912 :                if (jsrf == jsrflx_gaexch) then
         917    46164456 :                   if ( .not. dotend(l) ) cycle
         918    19359288 :                   fieldname = trim(cnst_name(lb)) // '_sfgaex1'
         919    46164456 :                else if (jsrf == jsrflx_rename) then
         920    46164456 :                   if ( .not. dotendrn(l) ) cycle
         921    20848464 :                   fieldname = trim(cnst_name(lb)) // '_sfgaex2'
         922             :                else
         923             :                   cycle
         924             :                end if
         925   671375952 :                do i = 1, ncol
         926   671375952 :                   qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf)*(adv_mass(l)/mwdry)
         927             :                end do
         928    40207752 :                call outfld( fieldname, qsrflx(:,l,jsrf), pcols, lchnk )
         929             :             else
         930    92328912 :                if (jsrf == jsrflx_gaexch) then
         931             :                   cycle
         932    46164456 :                else if (jsrf == jsrflx_rename) then
         933    46164456 :                   if ( .not. dotendqqcwrn(l) ) cycle
         934    20848464 :                   fieldname = trim(cnst_name_cw(lb)) // '_sfgaex2'
         935             :                else
         936             :                   cycle
         937             :                end if
         938   348120864 :                do i = 1, ncol
         939   348120864 :                   qqcwsrflx(i,l,jsrf) = qqcwsrflx(i,l,jsrf)*(adv_mass(l)/mwdry)
         940             :                end do
         941    20848464 :                call outfld( fieldname, qqcwsrflx(:,l,jsrf), pcols, lchnk )
         942             :             end if
         943             :          end do ! jac = ...
         944             :       end do ! jsrf = ...
         945             :    end do ! l = ...
         946             :
         947     1489176 :    return
         948     1489176 :    end subroutine modal_aero_gasaerexch_sub
         949             :
         950             :
         951             : !----------------------------------------------------------------------
         952             : !----------------------------------------------------------------------
         953     1489176 : subroutine gas_aer_uptkrates( ncol,       loffset,                &
         954     1489176 :                               q,          t,          pmid,       &
         955     1489176 :                               dgncur_awet,            uptkrate    )
         956             :
         957             : !
         958             : !                         /
         959             : !   computes   uptkrate = | dx  dN/dx  gas_conden_rate(Dp(x))
         960             : !                         /
         961             : !   using Gauss-Hermite quadrature of order nghq=2
         962             : !
         963             : !       Dp = particle diameter (cm)
         964             : !       x = ln(Dp)
         965             : !       dN/dx = log-normal particle number density distribution
         966             : !       gas_conden_rate(Dp) = 2 * pi * gasdiffus * Dp * F(Kn,ac)
         967             : !           F(Kn,ac) = Fuchs-Sutugin correction factor
         968             : !           Kn = Knudsen number
         969             : !           ac = accomodation coefficient
         970             : !
         971             :
         972     1489176 : use physconst, only: mwdry, rair
         973             :
         974             : implicit none
         975             :
         976             :
         977             :    integer,  intent(in) :: ncol                 ! number of atmospheric column
         978             :    integer,  intent(in) :: loffset
         979             :    real(r8), intent(in) :: q(ncol,pver,pcnstxx) ! Tracer array (mol,#/mol-air)
         980             :    real(r8), intent(in) :: t(pcols,pver)        ! Temperature in Kelvin
         981             :    real(r8), intent(in) :: pmid(pcols,pver)     ! Air pressure in Pa
         982             :    real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode)
         983             :
         984             :    real(r8), intent(out) :: uptkrate(ntot_amode,pcols,pver)
         985             :                             ! gas-to-aerosol mass transfer rates (1/s)
         986             :
         987             :
         988             : ! local
         989             :    integer, parameter :: nghq = 2
         990             :    integer :: i, iq, k, l1, l2, la, n
         991             :
         992             :    ! Can use sqrt here once Lahey is gone.
         993             :    real(r8), parameter :: tworootpi = 3.5449077_r8
         994             :    real(r8), parameter :: root2 = 1.4142135_r8
         995             :    real(r8), parameter :: beta = 2.0_r8
         996             :
         997             :    real(r8) :: aircon
         998             :    real(r8) :: const
         999             :    real(r8) :: dp, dum_m2v
        1000             :    real(r8) :: dryvol_a(pcols,pver)
        1001             :    real(r8) :: gasdiffus, gasspeed
        1002             :    real(r8) :: freepathx2, fuchs_sutugin
        1003             :    real(r8) :: knudsen
        1004             :    real(r8) :: lndp, lndpgn, lnsg
        1005             :    real(r8) :: num_a
        1006             :    real(r8) :: rhoair
        1007             :    real(r8) :: sumghq
        1008             :    real(r8), save :: xghq(nghq), wghq(nghq) ! quadrature abscissae and weights
        1009             :
        1010             :    data xghq / 0.70710678_r8, -0.70710678_r8 /
        1011             :    data wghq / 0.88622693_r8,  0.88622693_r8 /
        1012             :
        1013             :
        1014             : ! outermost loop over all modes
        1015     7445880 :    do n = 1, ntot_amode
        1016             :
        1017             : ! 22-aug-2007 rc easter - get number from q array rather
        1018             : !    than computing a "bounded" number conc.
        1019             : !! compute dry volume = sum_over_components{ component_mass / density }
        1020             : !!    (m3-AP/mol-air)
        1021             : !! compute it for all i,k to improve accessing q array
        1022             : !      dryvol_a(1:ncol,:) = 0.0_r8
        1023             : !      do l1 = 1, nspec_amode(n)
        1024             : !         l2 = lspectype_amode(l1,n)
        1025             : !! dum_m2v converts (kmol-AP/kmol-air) to (m3-AP/kmol-air)
        1026             : !! [m3-AP/kmol-AP]= [kg-AP/kmol-AP]  / [kg-AP/m3-AP]
        1027             : !         dum_m2v = specmw_amode(l2) / specdens_amode(l2)
        1028             : !         la = lmassptr_amode(l1,n)
        1029             : !         dryvol_a(1:ncol,:) = dryvol_a(1:ncol,:)    &
        1030             : !                            + max(0.0_r8,q(1:ncol,:,la))*dum_m2v
        1031             : !      end do
        1032             :
        1033             : ! loops k and i
        1034   561419352 :       do k=top_lev,pver
        1035  9256025376 :       do i=1,ncol
        1036             :
        1037  8696095200 :          rhoair = pmid(i,k)/(rair*t(i,k))   ! (kg-air/m3)
        1038             : !        aircon = 1.0e3*rhoair/mwdry        ! (mol-air/m3)
        1039             :
        1040             : !!   "bounded" number conc. (#/m3)
        1041             : !        num_a = dryvol_a(i,k)*v2ncur_a(i,k,n)*aircon
        1042             :
        1043             : !   number conc. (#/m3) -- note q(i,k,numptr) is (#/kmol-air)
        1044             : !   so need aircon in (kmol-air/m3)
        1045  8696095200 :          aircon = rhoair/mwdry              ! (kmol-air/m3)
        1046  8696095200 :          num_a = q(i,k,numptr_amode(n)-loffset)*aircon
        1047             :
        1048             : !   gasdiffus = h2so4 gas diffusivity from mosaic code (m^2/s)
        1049             : !               (pmid must be Pa)
        1050  8696095200 :          gasdiffus = 0.557e-4_r8 * (t(i,k)**1.75_r8) / pmid(i,k)
        1051             : !   gasspeed = h2so4 gas mean molecular speed from mosaic code (m/s)
        1052  8696095200 :          gasspeed  = 1.470e1_r8 * sqrt(t(i,k))
        1053             : !   freepathx2 = 2 * (h2so4 mean free path)  (m)
        1054  8696095200 :          freepathx2 = 6.0_r8*gasdiffus/gasspeed
        1055             :
        1056  8696095200 :          lnsg   = log( sigmag_amode(n) )
        1057  8696095200 :          lndpgn = log( dgncur_awet(i,k,n) )   ! (m)
        1058  8696095200 :          const  = tworootpi * num_a * exp(beta*lndpgn + 0.5_r8*(beta*lnsg)**2)
        1059             :
        1060             : !   sum over gauss-hermite quadrature points
        1061  8696095200 :          sumghq = 0.0_r8
        1062 26088285600 :          do iq = 1, nghq
        1063 17392190400 :             lndp = lndpgn + beta*lnsg**2 + root2*lnsg*xghq(iq)
        1064 17392190400 :             dp = exp(lndp)
        1065             :
        1066             : !   knudsen number
        1067 17392190400 :             knudsen = freepathx2/dp
        1068             : !  Changed by Manish Shrivastava on 7/17/2013 to use accom=1; because we do not know better
        1069             : !   following assumes accomodation coefficient = ac = 1. instead 0.65 ! answer change needs to be tested
        1070             : !   (Adams & Seinfeld, 2002, JGR, and references therein)
        1071             : !           fuchs_sutugin = (0.75*ac*(1. + knudsen)) /
        1072             : !                           (knudsen*(1.0 + knudsen + 0.283*ac) + 0.75*ac)
        1073             :             fuchs_sutugin = (0.4875_r8*(1._r8 + knudsen)) /   &
        1074 17392190400 :                             (knudsen*(1.184_r8 + knudsen) + 0.4875_r8)
        1075 26088285600 :             sumghq = sumghq + wghq(iq)*dp*fuchs_sutugin/(dp**beta)
        1076             :          end do
        1077  9250068672 :          uptkrate(n,i,k) = const * gasdiffus * sumghq
        1078             :
        1079             :       end do   ! "do i = 1, ncol"
        1080             :       end do   ! "do k = 1, pver"
        1081             :
        1082             :    end do   ! "do n = 1, ntot_soamode"
        1083             :
        1084             :
        1085     1489176 :    return
        1086             :    end subroutine gas_aer_uptkrates
        1087             :
        1088             : !----------------------------------------------------------------------
        1089             :
        1090  2174023800 :       subroutine modal_aero_soaexch( dtfull, temp, pres, &
        1091             :           niter, niter_max, ntot_amode, ntot_soamode, ntot_poaspec, ntot_soaspec, &
        1092  2174023800 :           mw_poa_host, mw_soa_host, &
        1093  2174023800 :           g_soa_in, a_soa_in, a_poa_in, xferrate_in, &
        1094  2174023800 :           g_soa_tend, a_soa_tend )
        1095             : !         g_soa_tend, a_soa_tend, g0_soa, idiagss )
        1096             :
        1097             : !-----------------------------------------------------------------------
        1098             : !
        1099             : ! Purpose:
        1100             : !
        1101             : ! calculates condensation/evaporation of "soa gas"
        1102             : ! to/from multiple aerosol modes in 1 grid cell
        1103             : !
        1104             : ! key assumptions
        1105             : ! (1) ambient equilibrium vapor pressure of soa gas
        1106             : !     is given by p0_soa_298 and delh_vap_soa
        1107             : ! (2) equilibrium vapor pressure of soa gas at aerosol
        1108             : !     particle surface is given by raoults law in the form
        1109             : !     g_star = g0_soa*[a_soa/(a_soa + a_opoa)]
        1110             : ! (3) (oxidized poa)/(total poa) is equal to frac_opoa (constant)
        1111             : !
        1112             : !
        1113             : ! Author: R. Easter and R. Zaveri
        1114             : ! Additions to run with multiple BC, SOA and POM's: Shrivastava et al., 2015
        1115             : !-----------------------------------------------------------------------
        1116             :
        1117             :       use mo_constants, only: rgas ! Gas constant (J/K/mol)
        1118             :
        1119             :       implicit none
        1120             :
        1121             :       real(r8), intent(in)  :: dtfull                     ! full integration time step (s)
        1122             :       real(r8), intent(in)  :: temp                       ! air temperature (K)
        1123             :       real(r8), intent(in)  :: pres                       ! air pressure (Pa)
        1124             :       integer,  intent(out) :: niter                      ! number of iterations performed
        1125             :       integer,  intent(in)  :: niter_max                  ! max allowed number of iterations
        1126             :       integer,  intent(in)  :: ntot_amode                 ! number of modes
        1127             :       integer,  intent(in)  :: ntot_soamode               ! number of modes having soa
        1128             :       integer,  intent(in)  :: ntot_poaspec               ! number of poa species
        1129             :       integer,  intent(in)  :: ntot_soaspec               ! number of soa species
        1130             :       real(r8), intent(in)  :: mw_poa_host(ntot_poaspec)  ! molec wght of poa used in host code
        1131             :       real(r8), intent(in)  :: mw_soa_host(ntot_soaspec)  ! molec wght of poa used in host code
        1132             :       real(r8), intent(in)  :: g_soa_in(ntot_soaspec)               ! initial soa gas mixrat (mol/mol at host mw)
        1133             :       real(r8), intent(in)  :: a_soa_in(ntot_amode,ntot_soaspec)    ! initial soa aerosol mixrat (mol/mol at host mw)
        1134             :       real(r8), intent(in)  :: a_poa_in(ntot_amode,ntot_poaspec)    ! initial poa aerosol mixrat (mol/mol at host mw)
        1135             :       real(r8), intent(in)  :: xferrate_in(ntot_amode,ntot_soaspec) ! gas-aerosol mass transfer rate (1/s)
        1136             :       real(r8), intent(out) :: g_soa_tend(ntot_soaspec)             ! soa gas mixrat tendency (mol/mol/s at host mw)
        1137             :       real(r8), intent(out) :: a_soa_tend(ntot_amode,ntot_soaspec)  ! soa aerosol mixrat tendency (mol/mol/s at host mw)
        1138             : !     integer,  intent(in)  :: idiagss
        1139             :
        1140             :       integer :: ll
        1141             :       integer :: m,k
        1142             :
        1143  4348047600 :       logical :: skip_soamode(ntot_amode)   ! true if this mode does not have soa
        1144             :
        1145             :       real(r8), parameter :: a_min1 = 1.0e-20_r8
        1146             :       real(r8), parameter :: g_min1 = 1.0e-20_r8
        1147             :       real(r8), parameter :: alpha = 0.05_r8     ! parameter used in calc of time step
        1148             :       real(r8), parameter :: dtsub_fixed = -1.0_r8  ! fixed sub-step for time integration (s)
        1149             :
        1150  4348047600 :       real(r8) :: a_ooa_sum_tmp(ntot_soamode)          ! total ooa (=soa+opoa) in a mode
        1151  4348047600 :       real(r8) :: a_opoa(ntot_soamode)                 ! oxidized-poa aerosol mixrat (mol/mol at actual mw)
        1152  4348047600 :       real(r8) :: a_soa(ntot_soamode,ntot_soaspec)     ! soa aerosol mixrat (mol/mol at actual mw)
        1153  4348047600 :       real(r8) :: a_soa_tmp(ntot_soamode,ntot_soaspec) ! temporary soa aerosol mixrat (mol/mol)
        1154  4348047600 :       real(r8) :: beta(ntot_soamode,ntot_soaspec)      ! dtcur*xferrate
        1155  4348047600 :       real(r8) :: delh_vap_soa(ntot_soaspec)           ! delh_vap_soa = heat of vaporization for gas soa (J/mol)
        1156  4348047600 :       real(r8) :: del_g_soa_tmp(ntot_soaspec)
        1157             :       real(r8) :: dtcur                                ! current time step (s)
        1158             :       real(r8) :: dtmax                                ! = (dtfull-tcur)
        1159  4348047600 :       real(r8) :: g0_soa(ntot_soaspec)                 ! ambient soa gas equilib mixrat (mol/mol at actual mw)
        1160  4348047600 :       real(r8) :: g_soa(ntot_soaspec)                  ! soa gas mixrat (mol/mol at actual mw)
        1161  4348047600 :       real(r8) :: g_star(ntot_soamode,ntot_soaspec)    ! soa gas mixrat that is in equilib
        1162             :                                                        ! with each aerosol mode (mol/mol)
        1163  4348047600 :       real(r8) :: mw_poa(ntot_poaspec)                 ! actual molec wght of poa
        1164  4348047600 :       real(r8) :: mw_soa(ntot_soaspec)                 ! actual molec wght of soa
        1165  4348047600 :       real(r8) :: opoa_frac(ntot_poaspec)              ! fraction of poa that is opoa
        1166  4348047600 :       real(r8) :: phi(ntot_soamode,ntot_soaspec)       ! "relative driving force"
        1167  4348047600 :       real(r8) :: p0_soa(ntot_soaspec)                 ! soa gas equilib vapor presssure (atm)
        1168  4348047600 :       real(r8) :: p0_soa_298(ntot_soaspec)             ! p0_soa_298 = soa gas equilib vapor presssure (atm) at 298 k
        1169  4348047600 :       real(r8) :: sat(ntot_soamode,ntot_soaspec)       ! sat(m,ll) = g0_soa(ll)/a_ooa_sum_tmp(m) = g_star(m,ll)/a_soa(m,ll)
        1170             :                                                        !    used by the numerical integration scheme -- it is not a saturation rato!
        1171             :       real(r8) :: tcur                                 ! current integration time (from 0 s)
        1172             :       real(r8) :: tmpa, tmpb, tmpf
        1173  4348047600 :       real(r8) :: tot_soa(ntot_soaspec)                ! g_soa + sum( a_soa(:) )
        1174  2174023800 :       real(r8) :: xferrate(ntot_amode,ntot_soaspec)    ! gas-aerosol mass transfer rate (1/s)
        1175             :
        1176             : ! Changed by Manish Shrivastava
        1177  4348047600 :       opoa_frac(:) = 0.0_r8 !POA does not form solution with SOA for all runs; set opoa_frac=0.0_r8  by Manish Shrivastava
        1178  4348047600 :       mw_poa(:) = 250.0_r8
        1179  4348047600 :       mw_soa(:) = 250.0_r8
        1180             :
        1181             :       ! New SOA properties added by Manish Shrivastava on 09/27/2012
        1182  2174023800 :       if (ntot_soaspec ==1) then
        1183  4348047600 :          p0_soa_298(:) = 9.7831E-11_r8
        1184  4348047600 :          delh_vap_soa(:) = 131.0e3_r8
        1185  4348047600 :          opoa_frac(:) = 0.0_r8
        1186           0 :       elseif (ntot_soaspec ==2) then
        1187             :          ! same for anthropogenic and biomass burning species
        1188           0 :          p0_soa_298 (1) = 1.0e-10_r8
        1189           0 :          p0_soa_298 (2) = 1.0e-10_r8
        1190           0 :          delh_vap_soa(:) = 156.0e3_r8
        1191           0 :       elseif(ntot_soaspec ==5) then
        1192             :          ! 5 volatility bins for each of the a combined SOA classes ( including biomass burning, fossil fuel, biogenic)
        1193           0 :          p0_soa_298 (1) = 9.7831E-13_r8 !soaff0 C*=0.01ug/m3
        1194           0 :          p0_soa_298 (2) = 9.7831E-12_r8 !soaff1 C*=0.10ug/m3
        1195           0 :          p0_soa_298 (3) = 9.7831E-11_r8 !soaff2 C*=1.0ug/m3
        1196           0 :          p0_soa_298 (4) = 9.7831E-10_r8 !soaff3 C*=10.0ug/m3
        1197           0 :          p0_soa_298 (5) = 9.7831E-9_r8  !soaff4 C*=100.0ug/m3
        1198             :
        1199           0 :          delh_vap_soa(1) = 153.0e3_r8
        1200           0 :          delh_vap_soa(2) = 142.0e3_r8
        1201           0 :          delh_vap_soa(3) = 131.0e3_r8
        1202           0 :          delh_vap_soa(4) = 120.0e3_r8
        1203           0 :          delh_vap_soa(5) = 109.0e3_r8
        1204           0 :       elseif(ntot_soaspec ==15) then
        1205             :          !
        1206             :          ! 5 volatility bins for each of the 3 SOA classes ( biomass burning, fossil fuel, biogenic)
        1207             :          ! SOA species 1-5 are for anthropogenic while 6-10 are for biomass burning SOA
        1208             :          ! SOA species 11-15 are for biogenic SOA, based on Cappa et al., Reference needs to be updated
        1209             :          ! For MW=250.0
        1210           0 :          p0_soa_298 (1) = 9.7831E-13_r8 !soaff0 C*=0.01ug/m3
        1211           0 :          p0_soa_298 (2) = 9.7831E-12_r8 !soaff1 C*=0.10ug/m3
        1212           0 :          p0_soa_298 (3) = 9.7831E-11_r8 !soaff2 C*=1.0ug/m3
        1213           0 :          p0_soa_298 (4) = 9.7831E-10_r8 !soaff3 C*=10.0ug/m3
        1214           0 :          p0_soa_298 (5) = 9.7831E-9_r8  !soaff4 C*=100.0ug/m3
        1215           0 :          p0_soa_298 (6) = 9.7831E-13_r8 !soabb0 C*=0.01ug/m3
        1216           0 :          p0_soa_298 (7) = 9.7831E-12_r8 !soabb1 C*=0.10ug/m3
        1217           0 :          p0_soa_298 (8) = 9.7831E-11_r8 !soabb2 C*=1.0ug/m3
        1218           0 :          p0_soa_298 (9) = 9.7831E-10_r8 !soabb3 C*=10.0ug/m3
        1219           0 :          p0_soa_298 (10) = 9.7831E-9_r8  !soabb4 C*=100.0ug/m3
        1220           0 :          p0_soa_298 (11) = 9.7831E-13_r8 !soabg0 C*=0.01ug/m3
        1221           0 :          p0_soa_298 (12) = 9.7831E-12_r8 !soabg1 C*=0.1ug/m3
        1222           0 :          p0_soa_298 (13) = 9.7831E-11_r8 !soabg2 C*=1.0ug/m3
        1223           0 :          p0_soa_298 (14) = 9.7831E-10_r8 !soabg3 C*=10.0ug/m3
        1224           0 :          p0_soa_298 (15) = 9.7831E-9_r8  !soabg4 C*=100.0ug/m3
        1225             :
        1226             :          !
        1227             :          ! have to be adjusted to 15 species, following the numbers by Epstein et al., 2012
        1228             :          !
        1229           0 :          delh_vap_soa(1) = 153.0e3_r8
        1230           0 :          delh_vap_soa(2) = 142.0e3_r8
        1231           0 :          delh_vap_soa(3) = 131.0e3_r8
        1232           0 :          delh_vap_soa(4) = 120.0e3_r8
        1233           0 :          delh_vap_soa(5) = 109.0e3_r8
        1234           0 :          delh_vap_soa(6) = 153.0e3_r8
        1235           0 :          delh_vap_soa(7) = 142.0e3_r8
        1236           0 :          delh_vap_soa(8) = 131.0e3_r8
        1237           0 :          delh_vap_soa(9) = 120.0e3_r8
        1238           0 :          delh_vap_soa(10) = 109.0e3_r8
        1239           0 :          delh_vap_soa(11) = 153.0e3_r8
        1240           0 :          delh_vap_soa(12) = 142.0e3_r8
        1241           0 :          delh_vap_soa(13) = 131.0e3_r8
        1242           0 :          delh_vap_soa(14) = 120.0e3_r8
        1243           0 :          delh_vap_soa(15) = 109.0e3_r8
        1244             :       endif
        1245             :
        1246             :       !BSINGH - Initialized g_soa_tend and a_soa_tend to circumvent the undefined behavior (04/16/12)
        1247  4348047600 :       g_soa_tend(:)   = 0.0_r8
        1248 13044142800 :       a_soa_tend(:,:) = 0.0_r8
        1249             :
        1250             :       ! determine which modes have non-zero transfer rates
        1251             :       !    and are involved in the soa gas-aerosol transfer
        1252             :       ! for diameter = 1 nm and number = 1 #/cm3, xferrate ~= 1e-9 s-1
        1253 10870119000 :       do m = 1, ntot_soamode
        1254  8696095200 :          skip_soamode(m) = .true.
        1255 19566214200 :          do ll = 1, ntot_soaspec
        1256  8696095200 :             xferrate(m,ll) = xferrate_in(m,ll)
        1257 17392190400 :             skip_soamode(m) = .false.
        1258             :          end do
        1259             :       end do
        1260             :
        1261             :       ! convert incoming mixing ratios from mol/mol at the "host-code" molec. weight (12.0 in cam5)
        1262             :       !    to mol/mol at the "actual" molec. weight (currently assumed to be 250.0)
        1263             :       ! also
        1264             :       !    force things to be non-negative
        1265             :       !    calc tot_soa(ll)
        1266             :       !    calc a_opoa (always slightly >0)
        1267  4348047600 :       do ll = 1, ntot_soaspec
        1268  2174023800 :          tmpf = mw_soa_host(ll)/mw_soa(ll)
        1269  2174023800 :          g_soa(ll) = max( g_soa_in(ll), 0.0_r8 ) * tmpf
        1270  2174023800 :          tot_soa(ll) = g_soa(ll)
        1271 13044142800 :          do m = 1, ntot_soamode
        1272  8696095200 :             if ( skip_soamode(m) ) cycle
        1273  8696095200 :             a_soa(m,ll) = max( a_soa_in(m,ll), 0.0_r8 ) * tmpf
        1274 10870119000 :             tot_soa(ll) = tot_soa(ll) + a_soa(m,ll)
        1275             :          end do
        1276             :       end do
        1277             :
        1278 10870119000 :       do m = 1, ntot_soamode
        1279  8696095200 :          if ( skip_soamode(m) ) cycle
        1280  8696095200 :          a_opoa(m) = 0.0_r8
        1281 19566214200 :          do ll = 1, ntot_poaspec
        1282 17392190400 :             a_opoa(m) = opoa_frac(ll)*a_poa_in(m,ll)
        1283             :          end do
        1284             :       end do
        1285             :
        1286             :       ! calc ambient equilibrium soa gas
        1287  4348047600 :       do ll = 1, ntot_soaspec
        1288  2174023800 :          p0_soa(ll) = p0_soa_298(ll) * &
        1289  2174023800 :               exp( -(delh_vap_soa(ll)/rgas)*((1.0_r8/temp)-(1.0_r8/298.0_r8)) )
        1290  4348047600 :          g0_soa(ll) = 1.01325e5_r8*p0_soa(ll)/pres
        1291             :       end do
        1292             :
        1293  2174023800 :       niter = 0
        1294  2174023800 :       tcur = 0.0_r8
        1295  2174023800 :       dtcur = 0.0_r8
        1296 13044142800 :       phi(:,:) = 0.0_r8
        1297 13044142800 :       g_star(:,:) = 0.0_r8
        1298             :
        1299             : !     if (idiagss > 0) then
        1300             : !        write(luna,'(a,1p,10e11.3)') 'p0, g0_soa', p0_soa, g0_soa
        1301             : !        write(luna,'(3a)') &
        1302             : !           'niter, tcur,   dtcur,    phi(:),                       ', &
        1303             : !           'g_star(:),                    ', &
        1304             : !           'a_soa(:),                     g_soa'
        1305             : !        write(luna,'(3a)') &
        1306             : !           '                         sat(:),                       ', &
        1307             : !           'sat(:)*a_soa(:)               ', &
        1308             : !           'a_opoa(:)'
        1309             : !        write(luna,'(i3,1p,20e10.2)') niter, tcur, dtcur, &
        1310             : !           phi(:), g_star(:), a_soa(:), g_soa
        1311             : !     end if
        1312             :
        1313             :
        1314             : ! integration loop -- does multiple substeps to reach dtfull
        1315             : time_loop: &
        1316  8438889498 :       do while (tcur < dtfull-1.0e-3_r8 )
        1317             :
        1318  6264865698 :       niter = niter + 1
        1319  6264865698 :       if (niter > niter_max) exit
        1320             :
        1321 31324328490 :       tmpa = 0.0_r8  ! time integration parameter for all soa species
        1322 31324328490 :       do m = 1, ntot_soamode
        1323 25059462792 :          if ( skip_soamode(m) ) cycle
        1324 56383791282 :          a_ooa_sum_tmp(m) = a_opoa(m) + sum( a_soa(m,1:ntot_soaspec) )
        1325             :       end do
        1326 12529731396 :       do ll = 1, ntot_soaspec
        1327             :          tmpb = 0.0_r8  ! time integration parameter for a single soa species
        1328 31324328490 :          do m = 1, ntot_soamode
        1329 25059462792 :             if ( skip_soamode(m) ) cycle
        1330 25059462792 :             sat(m,ll) = g0_soa(ll)/max( a_ooa_sum_tmp(m), a_min1 )
        1331 25059462792 :             g_star(m,ll) = sat(m,ll)*a_soa(m,ll)
        1332 25059462792 :             phi(m,ll) = (g_soa(ll) - g_star(m,ll))/max( g_soa(ll), g_star(m,ll), g_min1 )
        1333 31324328490 :             tmpb = tmpb + xferrate(m,ll)*abs(phi(m,ll))
        1334             :          end do
        1335 12529731396 :          tmpa = max( tmpa, tmpb )
        1336             :       end do
        1337             :
        1338             :       if (dtsub_fixed > 0.0_r8) then
        1339             :          dtcur = dtsub_fixed
        1340             :          tcur = tcur + dtcur
        1341             :       else
        1342  6264865698 :          dtmax = dtfull-tcur
        1343  6264865698 :          if (dtmax*tmpa <= alpha) then
        1344             : ! here alpha/tmpa >= dtmax, so this is final substep
        1345             :             dtcur = dtmax
        1346             :             tcur = dtfull
        1347             :          else
        1348  4090843651 :             dtcur = alpha/tmpa
        1349  4090843651 :             tcur = tcur + dtcur
        1350             :          end if
        1351             :       end if
        1352             :
        1353             : ! step 1 - for modes where soa is condensing, estimate "new" a_soa(m,ll)
        1354             : !    using an explicit calculation with "old" g_soa
        1355             : !    and g_star(m,ll) calculated using "old" a_soa(m,ll)
        1356             : ! do this to get better estimate of "new" a_soa(m,ll) and sat(m,ll)
        1357 31324328490 :       do m = 1, ntot_soamode
        1358 25059462792 :          if ( skip_soamode(m) ) cycle
        1359 50118925584 :          do ll = 1, ntot_soaspec
        1360             :             ! first ll loop calcs a_soa_tmp(m,ll) & a_ooa_sum_tmp
        1361 25059462792 :             a_soa_tmp(m,ll) = a_soa(m,ll)
        1362 25059462792 :             beta(m,ll) = dtcur*xferrate(m,ll)
        1363 25059462792 :             del_g_soa_tmp(ll) = g_soa(ll) - g_star(m,ll)
        1364 50118925584 :             if (del_g_soa_tmp(ll) > 0.0_r8) then
        1365 19816107746 :                a_soa_tmp(m,ll) = a_soa(m,ll) + beta(m,ll)*del_g_soa_tmp(ll)
        1366             :             end if
        1367             :          end do
        1368 50118925584 :          a_ooa_sum_tmp(m) = a_opoa(m) + sum( a_soa_tmp(m,1:ntot_soaspec) )
        1369 56383791282 :          do ll = 1, ntot_soaspec
        1370             :             ! second ll loop calcs sat & g_star
        1371 50118925584 :             if (del_g_soa_tmp(ll) > 0.0_r8) then
        1372 19816107746 :                sat(m,ll) = g0_soa(ll)/max( a_ooa_sum_tmp(m), a_min1 )
        1373 19816107746 :                g_star(m,ll) = sat(m,ll)*a_soa_tmp(m,ll)   ! this just needed for diagnostics
        1374             :             end if
        1375             :          end do
        1376             :       end do
        1377             :
        1378             : ! step 2 - implicit in g_soa and semi-implicit in a_soa,
        1379             : !    with g_star(m,ll) calculated semi-implicitly
        1380 14703755196 :       do ll = 1, ntot_soaspec
        1381             :          tmpa = 0.0_r8
        1382             :          tmpb = 0.0_r8
        1383 31324328490 :          do m = 1, ntot_soamode
        1384 25059462792 :             if ( skip_soamode(m) ) cycle
        1385 25059462792 :             tmpa = tmpa + a_soa(m,ll)/(1.0_r8 + beta(m,ll)*sat(m,ll))
        1386 31324328490 :             tmpb = tmpb + beta(m,ll)/(1.0_r8 + beta(m,ll)*sat(m,ll))
        1387             :          end do
        1388             :
        1389  6264865698 :          g_soa(ll) = (tot_soa(ll) - tmpa)/(1.0_r8 + tmpb)
        1390  6264865698 :          g_soa(ll) = max( 0.0_r8, g_soa(ll) )
        1391 37589194188 :          do m = 1, ntot_soamode
        1392 25059462792 :             if ( skip_soamode(m) ) cycle
        1393 25059462792 :             a_soa(m,ll) = (a_soa(m,ll) + beta(m,ll)*g_soa(ll))/   &
        1394 31324328490 :                        (1.0_r8 + beta(m,ll)*sat(m,ll))
        1395             :          end do
        1396             :       end do
        1397             :
        1398             : !     if (idiagss > 0) then
        1399             : !        write(luna,'(i3,1p,20e10.2)') niter, tcur, dtcur, &
        1400             : !           phi(:), g_star(:), a_soa(:), g_soa
        1401             : !        write(luna,'(23x,1p,20e10.2)') &
        1402             : !           sat(:), sat(:)*a_soa(:), a_opoa(:)
        1403             : !     end if
        1404             :
        1405             : !     if (niter > 9992000) then
        1406             : !        write(luna,'(a)') '*** to many iterations'
        1407             : !        exit
        1408             : !     end if
        1409             :
        1410             :       end do time_loop
        1411             :
        1412             :
        1413             : ! calculate outgoing tendencies (at the host-code molec. weight)
        1414             : ! (a_soa & g_soa are at actual mw, but a_soa_in & g_soa_in are at host-code mw)
        1415  4348047600 :       do ll = 1, ntot_soaspec
        1416  2174023800 :          tmpf = mw_soa(ll)/mw_soa_host(ll)
        1417  2174023800 :          g_soa_tend(ll) = (g_soa(ll)*tmpf - g_soa_in(ll))/dtfull
        1418 13044142800 :          do m = 1, ntot_soamode
        1419  8696095200 :             if ( skip_soamode(m) ) cycle
        1420 10870119000 :             a_soa_tend(m,ll) = (a_soa(m,ll)*tmpf - a_soa_in(m,ll))/dtfull
        1421             :          end do
        1422             :       end do
        1423             :
        1424             :
        1425  2174023800 :       return
        1426             :
        1427             :       end subroutine modal_aero_soaexch
        1428             :
        1429             : !----------------------------------------------------------------------
        1430             :
        1431             : !----------------------------------------------------------------------
        1432             :
        1433        1536 :       subroutine modal_aero_gasaerexch_init
        1434             :
        1435             : !-----------------------------------------------------------------------
        1436             : !
        1437             : ! Purpose:
        1438             : !    set do_adjust and do_aitken flags
        1439             : !    create history fields for column tendencies associated with
        1440             : !       modal_aero_calcsize
        1441             : !
        1442             : ! Author: R. Easter
        1443             : !
        1444             : !-----------------------------------------------------------------------
        1445             :
        1446             : use modal_aero_data
        1447             : use modal_aero_rename
        1448             :
        1449             : use cam_abortutils, only: endrun
        1450             : use cam_history,    only: addfld, add_default, fieldname_len, horiz_only
        1451             : use constituents,   only: pcnst, cnst_get_ind, cnst_name
        1452             : use spmd_utils,     only: masterproc
        1453             : use phys_control,   only: phys_getopts
        1454             :
        1455             : implicit none
        1456             :
        1457             : !-----------------------------------------------------------------------
        1458             : ! arguments
        1459             :
        1460             : !-----------------------------------------------------------------------
        1461             : ! local
        1462             :    integer  :: ipair, iq, iqfrm, iqfrm_aa, iqtoo, iqtoo_aa
        1463             :    integer  :: jac,jsoa,p
        1464             :    integer  :: l, l1, l2, lsfrm, lstoo, lunout
        1465             :    integer  :: l_so4g, l_nh4g, l_msag
        1466             :    integer  :: m, mfrm, mtoo
        1467             :    integer  :: n, nacc, nait
        1468             :    integer  :: nchfrm, nchfrmskip, nchtoo, nchtooskip, nspec
        1469             :
        1470             :    logical  :: do_msag, do_nh4g
        1471        3072 :    logical  :: do_soag_any, do_soag(nsoa)
        1472             :    logical  :: dotend(pcnst), dotendqqcw(pcnst)
        1473             :
        1474             :    real(r8) :: tmp1, tmp2
        1475             :
        1476             :    character(len=fieldname_len)   :: tmpnamea, tmpnameb
        1477             :    character(len=fieldname_len+3) :: fieldname
        1478             :    character(128)                 :: long_name
        1479             :    character(128)                 :: msg
        1480             :    character(8)                   :: unit
        1481             :
        1482             :    logical                        :: history_aerosol      ! Output the MAM aerosol tendencies
        1483             :    logical                        :: history_aerocom    ! Output the aerocom history
        1484             :    !-----------------------------------------------------------------------
        1485             :
        1486        1536 :         call phys_getopts( history_aerosol_out        = history_aerosol   )
        1487             :
        1488        1536 :         maxspec_pcage = nspec_max
        1489        4608 :         allocate(lspecfrm_pcage(maxspec_pcage))
        1490        3072 :         allocate(lspectoo_pcage(maxspec_pcage))
        1491        4608 :         allocate(soa_equivso4_factor(nsoa))
        1492        3072 :         allocate(fac_m2v_soa(nsoa))
        1493        4608 :         allocate(fac_m2v_pcarbon(nspec_max))
        1494        1536 :         lunout = 6
        1495             : !
        1496             : !   define "from mode" and "to mode" for primary carbon aging
        1497             : !
        1498             : !   skip (turn off) aging if either is absent,
        1499             : !      or if accum mode so4 is absent
        1500             : !
        1501        1536 :         modefrm_pcage = -999888777
        1502        1536 :         modetoo_pcage = -999888777
        1503        1536 :         if ((modeptr_pcarbon <= 0) .or. (modeptr_accum <= 0)) goto 15000
        1504        1536 :         l = lptr_so4_a_amode(modeptr_accum)
        1505        1536 :         if ((l < 1) .or. (l > pcnst)) goto 15000
        1506             :
        1507        1536 :         modefrm_pcage = modeptr_pcarbon
        1508        1536 :         modetoo_pcage = modeptr_accum
        1509             :
        1510             : !
        1511             : !   define species involved in each primary carbon aging pairing
        1512             : !       (include aerosol water)
        1513             : !
        1514             : !
        1515        1536 :         mfrm = modefrm_pcage
        1516        1536 :         mtoo = modetoo_pcage
        1517             :
        1518        1536 :         if (mfrm < 10) then
        1519             :             nchfrmskip = 1
        1520           0 :         else if (mfrm < 100) then
        1521             :             nchfrmskip = 2
        1522             :         else
        1523           0 :             nchfrmskip = 3
        1524             :         end if
        1525        1536 :         if (mtoo < 10) then
        1526             :             nchtooskip = 1
        1527           0 :         else if (mtoo < 100) then
        1528             :             nchtooskip = 2
        1529             :         else
        1530           0 :             nchtooskip = 3
        1531             :         end if
        1532        1536 :         nspec = 0
        1533             :
        1534        7680 : aa_iqfrm: do iqfrm = -1, nspec_amode(mfrm)
        1535             :
        1536        6144 :             if (iqfrm == -1) then
        1537        1536 :                 lsfrm = numptr_amode(mfrm)
        1538        1536 :                 lstoo = numptr_amode(mtoo)
        1539        4608 :             else if (iqfrm == 0) then
        1540             : !   bypass transfer of aerosol water due to primary-carbon aging
        1541             :                 cycle aa_iqfrm
        1542             : !               lsfrm = lwaterptr_amode(mfrm)
        1543             : !               lstoo = lwaterptr_amode(mtoo)
        1544             :             else
        1545        3072 :                 lsfrm = lmassptr_amode(iqfrm,mfrm)
        1546        3072 :                 lstoo = 0
        1547             :             end if
        1548        4608 :             if ((lsfrm < 1) .or. (lsfrm > pcnst)) cycle aa_iqfrm
        1549             :
        1550        4608 :             if (lsfrm>0 .and. iqfrm>0 ) then
        1551        3072 :                 nchfrm = len( trim( cnst_name(lsfrm) ) ) - nchfrmskip
        1552             :
        1553             : ! find "too" species having same lspectype_amode as the "frm" species
        1554             : ! AND same cnst_name (except for last 1/2/3 characters which are the mode index)
        1555        9216 :                 do iqtoo = 1, nspec_amode(mtoo)
        1556             : !                   if ( lspectype_amode(iqtoo,mtoo) .eq.   &
        1557             : !                        lspectype_amode(iqfrm,mfrm) ) then
        1558        9216 :                         lstoo = lmassptr_amode(iqtoo,mtoo)
        1559        9216 :                         nchtoo = len( trim( cnst_name(lstoo) ) ) - nchtooskip
        1560        9216 :                         if (cnst_name(lsfrm)(1:nchfrm) == cnst_name(lstoo)(1:nchtoo)) then
        1561             :                             exit
        1562             :                         else
        1563        6144 :                             lstoo = 0
        1564             :                         end if
        1565             : !                   end if
        1566             :                 end do
        1567             :             end if
        1568             :
        1569        4608 :             if ((lstoo < 1) .or. (lstoo > pcnst)) lstoo = 0
        1570        4608 :             nspec = nspec + 1
        1571        4608 :             lspecfrm_pcage(nspec) = lsfrm
        1572        7680 :             lspectoo_pcage(nspec) = lstoo
        1573             :         end do aa_iqfrm
        1574             :
        1575        1536 :         nspecfrm_pcage = nspec
        1576             :
        1577             : !
        1578             : !   output results
        1579             : !
        1580        1536 :         if ( masterproc ) then
        1581             :
        1582           2 :         write(lunout,9310)
        1583             :
        1584           2 :           mfrm = modefrm_pcage
        1585           2 :           mtoo = modetoo_pcage
        1586           2 :           write(lunout,9320) 1, mfrm, mtoo
        1587             :
        1588           8 :           do iq = 1, nspecfrm_pcage
        1589           6 :             lsfrm = lspecfrm_pcage(iq)
        1590           6 :             lstoo = lspectoo_pcage(iq)
        1591           8 :             if (lstoo .gt. 0) then
        1592           6 :                 write(lunout,9330) lsfrm, cnst_name(lsfrm),   &
        1593          12 :                         lstoo, cnst_name(lstoo)
        1594             :             else
        1595           0 :                 write(lunout,9340) lsfrm, cnst_name(lsfrm)
        1596             :             end if
        1597             :           end do
        1598             :
        1599           2 :         write(lunout,*)
        1600             :
        1601             :         end if ! ( masterproc )
        1602             :
        1603             : 9310    format( / 'subr. modal_aero_gasaerexch_init - primary carbon aging pointers' )
        1604             : 9320    format( 'pair', i3, 5x, 'mode', i3, ' ---> mode', i3 )
        1605             : 9330    format( 5x, 'spec', i3, '=', a, ' ---> spec', i3, '=', a )
        1606             : 9340    format( 5x, 'spec', i3, '=', a, ' ---> LOSS' )
        1607             :
        1608             :
        1609             : 15000 continue
        1610             :
        1611             : ! set tendency flags and gas species indices and flags
        1612        1536 :       dotend(:) = .false.
        1613             :
        1614        1536 :       call cnst_get_ind( 'H2SO4', l_so4g, .false. )
        1615        1536 :       if ((l_so4g <= 0) .or. (l_so4g > pcnst)) then
        1616             :          write( *, '(/a/a,2i7)' )   &
        1617           0 :             '*** modal_aero_gasaerexch_init -- cannot find H2SO4 species',   &
        1618           0 :             '    l_so4g=', l_so4g
        1619           0 :          call endrun( 'modal_aero_gasaerexch_init error' )
        1620             :       end if
        1621        1536 :       dotend(l_so4g) = .true.
        1622             :
        1623        1536 :       call cnst_get_ind( 'NH3',   l_nh4g, .false. )
        1624        1536 :       do_nh4g = .false.
        1625        1536 :       if ((l_nh4g > 0) .and. (l_nh4g <= pcnst)) then
        1626           0 :          do_nh4g = .true.
        1627           0 :          dotend(l_nh4g) = .true.
        1628             :       end if
        1629             :
        1630        1536 :       call cnst_get_ind( 'MSA',   l_msag, .false. )
        1631        1536 :       do_msag = .false.
        1632        1536 :       if ((l_msag > 0) .and. (l_msag <= pcnst)) then
        1633           0 :          do_msag = .true.
        1634           0 :          dotend(l_msag) = .true.
        1635             :       end if
        1636             :
        1637        1536 :       do_soag_any = .false.
        1638        3072 :       do_soag(:) = .false.
        1639        3072 :       do jsoa = 1, nsoa
        1640        1536 :          l = lptr2_soa_g_amode(jsoa)
        1641        3072 :          if ((l > 0) .and. (l <= pcnst)) then
        1642        1536 :             do_soag_any = .true.
        1643        1536 :             do_soag(jsoa) = .true.
        1644        1536 :             dotend(l) = .true.
        1645             :          end if
        1646             :       end do
        1647             :
        1648             :
        1649        7680 :       do n = 1, ntot_amode
        1650        6144 :          l = lptr_so4_a_amode(n)
        1651        6144 :          if ((l > 0) .and. (l <= pcnst)) then
        1652        4608 :             dotend(l) = .true.
        1653        4608 :             if ( do_nh4g ) then
        1654           0 :                l = lptr_nh4_a_amode(n)
        1655           0 :                if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
        1656             :             end if
        1657             :          end if
        1658       13824 :          do jsoa = 1, nsoa
        1659       12288 :             if ( do_soag(jsoa) ) then
        1660        6144 :                l = lptr2_soa_a_amode(n,jsoa)
        1661        6144 :                if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
        1662             :             end if
        1663             :          end do
        1664             :       end do
        1665             :
        1666        1536 :       if (modefrm_pcage > 0) then
        1667        6144 :          do iq = 1, nspecfrm_pcage
        1668        4608 :             lsfrm = lspecfrm_pcage(iq)
        1669        4608 :             lstoo = lspectoo_pcage(iq)
        1670        6144 :             if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
        1671        4608 :                dotend(lsfrm) = .true.
        1672        4608 :                if ((lstoo > 0) .and. (lstoo <= pcnst)) then
        1673        4608 :                   dotend(lstoo) = .true.
        1674             :                end if
        1675             :             end if
        1676             :          end do
        1677             :       end if
        1678             :
        1679             : !---------define history fields for new cond/evap diagnostics----------------------------------------
        1680        1536 :       fieldname=trim('qconff_gaex')
        1681        1536 :       long_name = trim('3D fields for Fossil SOA condensation')
        1682        1536 :       unit = 'kg/kg/s'
        1683        3072 :       call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
        1684        1536 :       if ( history_aerosol ) then
        1685           0 :          call add_default( fieldname,  1, ' ' )
        1686             :       endif
        1687        1536 :       if ( masterproc ) write(*,'(3(a,3x))') 'qconff addfld', fieldname, unit
        1688             :
        1689        1536 :       fieldname=trim('qevapff_gaex')
        1690        1536 :       long_name = trim('3D fields for Fossil SOA evaporation')
        1691        3072 :       call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
        1692        1536 :       if ( history_aerosol ) then
        1693           0 :          call add_default( fieldname,  1, ' ' )
        1694             :       endif
        1695        1536 :       if ( masterproc ) write(*,'(3(a,3x))') 'qevapff addfld', fieldname, unit
        1696             :
        1697        1536 :       fieldname=trim('qconbb_gaex')
        1698        1536 :       long_name = trim('3D fields for Biomass SOA condensation')
        1699        3072 :       call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
        1700        1536 :       if ( history_aerosol ) then
        1701           0 :          call add_default( fieldname,  1, ' ' )
        1702             :       endif
        1703        1536 :       if ( masterproc ) write(*,'(3(a,3x))') 'qconbb addfld', fieldname, unit
        1704             :
        1705        1536 :       fieldname=trim('qevapbb_gaex')
        1706        1536 :       long_name = trim('3D fields for Biomass SOA evaporation')
        1707        3072 :       call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
        1708        1536 :       if ( history_aerosol ) then
        1709           0 :          call add_default( fieldname,  1, ' ' )
        1710             :       endif
        1711        1536 :       if ( masterproc ) write(*,'(3(a,3x))') 'qevapbb addfld', fieldname, unit
        1712             :
        1713        1536 :       fieldname=trim('qconbg_gaex')
        1714        1536 :       long_name = trim('3D fields for Biogenic SOA condensation')
        1715        3072 :       call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
        1716        1536 :       if ( history_aerosol ) then
        1717           0 :          call add_default( fieldname,  1, ' ' )
        1718             :       endif
        1719        1536 :       if ( masterproc ) write(*,'(3(a,3x))') 'qconbg addfld', fieldname, unit
        1720             :
        1721        1536 :       fieldname=trim('qevapbg_gaex')
        1722        1536 :       long_name = trim('3D fields for Biogenic SOA evaporation')
        1723        3072 :       call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
        1724        1536 :       if ( history_aerosol ) then
        1725           0 :          call add_default( fieldname,  1, ' ' )
        1726             :       endif
        1727        1536 :       if ( masterproc ) write(*,'(3(a,3x))') 'qevapbg addfld', fieldname, unit
        1728             :
        1729        1536 :       fieldname=trim('qcon_gaex')
        1730        1536 :       long_name = trim('3D fields for SOA condensation')
        1731        3072 :       call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
        1732        1536 :       if ( history_aerosol ) then
        1733           0 :          call add_default( fieldname,  1, ' ' )
        1734             :       endif
        1735        1536 :       if ( masterproc ) write(*,'(3(a,3x))') 'qcon addfld', fieldname, unit
        1736             :
        1737        1536 :       fieldname=trim('qevap_gaex')
        1738        1536 :       long_name = trim('3D fields for Biogenic SOA evaporation')
        1739        3072 :       call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
        1740        1536 :       if ( history_aerosol ) then
        1741           0 :          call add_default( fieldname,  1, ' ' )
        1742             :       endif
        1743        1536 :       if ( masterproc ) write(*,'(3(a,3x))') 'qevap addfld', fieldname, unit
        1744             : !------------------------------------------------------------------------------
        1745             :
        1746             : !  define history fields for basic gas-aer exchange
        1747             : !  and primary carbon aging from that
        1748       64512 :       do l = 1, pcnst
        1749       62976 :          if ( .not. dotend(l) ) cycle
        1750             :
        1751       19968 :          tmpnamea = cnst_name(l)
        1752       19968 :          fieldname = trim(tmpnamea) // '_sfgaex1'
        1753       19968 :          long_name = trim(tmpnamea) // ' gas-aerosol-exchange primary column tendency'
        1754       19968 :          unit = 'kg/m2/s'
        1755       19968 :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
        1756       19968 :          if ( history_aerosol ) then
        1757           0 :             call add_default( fieldname, 1, ' ' )
        1758             :          endif
        1759       21504 :          if ( masterproc ) write(*,'(3(a,3x))') 'gasaerexch addfld', fieldname, unit
        1760             :
        1761             :       end do   ! l = ...
        1762             : !  define history fields for aitken-->accum renaming
        1763        1536 :       dotend(:) = .false.
        1764        1536 :       dotendqqcw(:) = .false.
        1765        6144 :       do ipair = 1, npair_renamexf
        1766       26112 :          do iq = 1, nspecfrm_renamexf(ipair)
        1767       19968 :             lsfrm = lspecfrma_renamexf(iq,ipair)
        1768       19968 :             lstoo = lspectooa_renamexf(iq,ipair)
        1769       19968 :             if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
        1770       19968 :                dotend(lsfrm) = .true.
        1771       19968 :                if ((lstoo > 0) .and. (lstoo <= pcnst)) then
        1772       19968 :                   dotend(lstoo) = .true.
        1773             :                end if
        1774             :             end if
        1775             :
        1776       19968 :             lsfrm = lspecfrmc_renamexf(iq,ipair)
        1777       19968 :             lstoo = lspectooc_renamexf(iq,ipair)
        1778       24576 :             if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
        1779       19968 :                dotendqqcw(lsfrm) = .true.
        1780       19968 :                if ((lstoo > 0) .and. (lstoo <= pcnst)) then
        1781       19968 :                   dotendqqcw(lstoo) = .true.
        1782             :                end if
        1783             :             end if
        1784             :          end do ! iq = ...
        1785             :       end do ! ipair = ...
        1786             :
        1787       64512 :       do l = 1, pcnst
        1788      190464 :       do jac = 1, 2
        1789      125952 :          if (jac == 1) then
        1790       62976 :             if ( .not. dotend(l) ) cycle
        1791       21504 :             tmpnamea = cnst_name(l)
        1792             :          else
        1793       62976 :             if ( .not. dotendqqcw(l) ) cycle
        1794       21504 :             tmpnamea = cnst_name_cw(l)
        1795             :          end if
        1796             :
        1797       43008 :          fieldname = trim(tmpnamea) // '_sfgaex2'
        1798       43008 :          long_name = trim(tmpnamea) // ' gas-aerosol-exchange renaming column tendency'
        1799       43008 :          unit = 'kg/m2/s'
        1800       43008 :          if ((tmpnamea(1:3) == 'num') .or. &
        1801        9216 :              (tmpnamea(1:3) == 'NUM')) unit = '#/m2/s'
        1802       43008 :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
        1803       43008 :          if ( history_aerosol ) then
        1804           0 :             call add_default( fieldname, 1, ' ' )
        1805             :          endif
        1806      105984 :          if ( masterproc ) write(*,'(3(a,3x))') 'gasaerexch addfld', fieldname, unit
        1807             :       end do   ! jac = ...
        1808             :       end do   ! l = ...
        1809             :
        1810             :
        1811             : ! set for used in aging calcs:
        1812             : !    fac_m2v_so4, fac_m2v_nh4, fac_m2v_soa(:)
        1813             : !    soa_equivso4_factor(:)
        1814        3072 :       soa_equivso4_factor = 0.0_r8
        1815        1536 :       if (modefrm_pcage > 0) then
        1816        1536 :          n = modeptr_accum
        1817        1536 :          l = lptr_so4_a_amode(n) ; l2 = -1
        1818        1536 :          if (l <= 0) call endrun( 'modal_aero_gasaerexch_init error a001 finding accum. so4' )
        1819       10752 :          do l1 = 1, nspec_amode(n)
        1820       10752 :             if (lmassptr_amode(l1,n) == l) then
        1821             : !               l2 = lspectype_amode(l1,n)
        1822        1536 :                l2 = l1
        1823             : !               fac_m2v_so4 = specmw_amode(l2) / specdens_amode(l2)
        1824        1536 :                fac_m2v_so4 = specmw_amode(l1,n) / specdens_amode(l1,n)
        1825             : !               tmp2 = spechygro(l2)
        1826        1536 :                tmp2 = spechygro(l1,n)
        1827             :
        1828             :             end if
        1829             :          end do
        1830        1536 :          if (l2 <= 0) call endrun( 'modal_aero_gasaerexch_init error a002 finding accum. so4' )
        1831             :
        1832        1536 :          l = lptr_nh4_a_amode(n) ; l2 = -1
        1833        1536 :          if (l > 0) then
        1834           0 :             do l1 = 1, nspec_amode(n)
        1835           0 :                if (lmassptr_amode(l1,n) == l) then
        1836             : !                  l2 = lspectype_amode(l1,n)
        1837           0 :                   l2 = l1
        1838             : !                  fac_m2v_nh4 = specmw_amode(l2) / specdens_amode(l2)
        1839           0 :                   fac_m2v_nh4 = specmw_amode(l1,n) / specdens_amode(l1,n)
        1840             :
        1841             :                end if
        1842             :             end do
        1843           0 :             if (l2 <= 0) call endrun( 'modal_aero_gasaerexch_init error a002 finding accum. nh4' )
        1844             :          else
        1845        1536 :             fac_m2v_nh4 = fac_m2v_so4
        1846             :          end if
        1847             :
        1848        3072 :          do jsoa = 1, nsoa
        1849        1536 :             l = lptr2_soa_a_amode(n,jsoa) ; l2 = -1
        1850        1536 :             if (l <= 0) then
        1851           0 :                write( msg, '(a,i4)') 'modal_aero_gasaerexch_init error a001 finding accum. jsoa =', jsoa
        1852           0 :                call endrun( msg )
        1853             :             end if
        1854       10752 :             do l1 = 1, nspec_amode(n)
        1855       10752 :                if (lmassptr_amode(l1,n) == l) then
        1856             : !                  l2 = lspectype_amode(l1,n)
        1857        1536 :                   l2 = l1
        1858             : !                  fac_m2v_soa(jsoa) = specmw_amode(l2) / specdens_amode(l2)
        1859        1536 :                   fac_m2v_soa(jsoa) = specmw_amode(l1,n) / specdens_amode(l1,n)
        1860             : !                  soa_equivso4_factor(jsoa) = spechygro(l2)/tmp2
        1861        1536 :                   soa_equivso4_factor(jsoa) = spechygro(l1,n)/tmp2
        1862             :                end if
        1863             :             end do
        1864        3072 :             if (l2 <= 0) then
        1865           0 :                write( msg, '(a,i4)') 'modal_aero_gasaerexch_init error a002 finding accum. jsoa =', jsoa
        1866           0 :                call endrun( msg )
        1867             :             end if
        1868             :          end do
        1869             :
        1870       10752 :          fac_m2v_pcarbon(:) = 0.0_r8
        1871        1536 :          n = modeptr_pcarbon
        1872        4608 :          do l = 1, nspec_amode(n)
        1873             : !            l2 = lspectype_amode(l,n)
        1874             : !      fac_m2v converts (kmol-AP/kmol-air) to (m3-AP/kmol-air)
        1875             : !           [m3-AP/kmol-AP]    = [kg-AP/kmol-AP]  / [kg-AP/m3-AP]
        1876             : !            fac_m2v_pcarbon(l) = specmw_amode(l2) / specdens_amode(l2)
        1877        4608 :             fac_m2v_pcarbon(l) = specmw_amode(l,n) / specdens_amode(l,n)
        1878             :          end do
        1879             :       end if
        1880             :
        1881             :
        1882        1536 :       return
        1883             :
        1884        3072 :       end subroutine modal_aero_gasaerexch_init
        1885             :
        1886             :
        1887             : !----------------------------------------------------------------------
        1888             :
        1889             : end module modal_aero_gasaerexch

Generated by: LCOV version 1.14