LCOV - code coverage report
Current view: | top level - chemistry/utils - modal_aero_calcsize.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 175 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 5 | 0.0 %

              Line data    Source code
           1             : module modal_aero_calcsize
           2             :
           3             : !   RCE 07.04.13:  Adapted from MIRAGE2 code
           4             :
           5             : use shr_kind_mod,     only: r8 => shr_kind_r8
           6             : use spmd_utils,       only: masterproc
           7             : use physconst,        only: pi, rhoh2o, gravit
           8             :
           9             : use ppgrid,           only: pcols, pver
          10             : use physics_types,    only: physics_state, physics_ptend
          11             : use physics_buffer,   only: physics_buffer_desc, pbuf_get_index, pbuf_old_tim_idx, pbuf_get_field
          12             :
          13             : use phys_control,     only: phys_getopts
          14             : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, rad_cnst_get_aer_props, &
          15             :                             rad_cnst_get_mode_props, rad_cnst_get_mode_num
          16             :
          17             : use cam_logfile,      only: iulog
          18             : use cam_abortutils,   only: endrun
          19             : use cam_history,      only: addfld, add_default, fieldname_len, horiz_only, outfld
          20             : use constituents,     only: pcnst, cnst_name
          21             : use modal_aero_wateruptake, only: modal_strat_sulfate
          22             :
          23             : use ref_pres,         only: top_lev => clim_modal_aero_top_lev
          24             :
          25             : #ifdef MODAL_AERO
          26             :
          27             : ! these are the variables needed for the diagnostic calculation of dry radius
          28             : use modal_aero_data, only: ntot_amode, nspec_amode, nspec_max, &
          29             :                            numptr_amode, &
          30             :                            alnsg_amode, &
          31             :                            voltonumbhi_amode, voltonumblo_amode, &
          32             :                            dgnum_amode, dgnumhi_amode, dgnumlo_amode
          33             :
          34             :
          35             : ! these variables are needed for the prognostic calculations to exchange mass
          36             : ! between modes
          37             : use modal_aero_data,  only: numptrcw_amode, mprognum_amode, qqcw_get_field, lmassptrcw_amode, &
          38             :            lmassptr_amode, modeptr_accum, modeptr_aitken, &
          39             :            specmw_amode, specdens_amode, voltonumb_amode, &
          40             :            cnst_name_cw
          41             :
          42             : use modal_aero_rename, only: lspectooa_renamexf, lspecfrma_renamexf, lspectooc_renamexf, lspecfrmc_renamexf, &
          43             :            modetoo_renamexf, nspecfrm_renamexf, npair_renamexf, modefrm_renamexf
          44             :
          45             :
          46             : #endif
          47             :
          48             :
          49             : implicit none
          50             : private
          51             : save
          52             :
          53             : public modal_aero_calcsize_init, modal_aero_calcsize_sub, modal_aero_calcsize_diag
          54             : public :: modal_aero_calcsize_reg
          55             :
          56             : logical :: do_adjust_default
          57             : logical :: do_aitacc_transfer_default
          58             :
          59             : integer :: dgnum_idx     = -1
          60             : integer :: hygro_idx     = -1
          61             : integer :: dryvol_idx    = -1
          62             : integer :: dryrad_idx    = -1
          63             : integer :: drymass_idx   = -1
          64             : integer :: so4dryvol_idx = -1
          65             : integer :: naer_idx      = -1
          66             : integer :: sulfeq_idx    = -1
          67             :
          68             :
          69             : !===============================================================================
          70             : contains
          71             : !===============================================================================
          72             :
          73           0 : subroutine modal_aero_calcsize_reg()
          74             :   use physics_buffer,   only: pbuf_add_field, dtype_r8
          75             :   use rad_constituents, only: rad_cnst_get_info
          76             :
          77             :   integer :: nmodes
          78             :
          79           0 :   call rad_cnst_get_info(0, nmodes=nmodes)
          80             :
          81           0 :   call pbuf_add_field('DGNUM', 'global',  dtype_r8, (/pcols, pver, nmodes/), dgnum_idx)
          82             :
          83           0 :   call pbuf_add_field('HYGRO',     'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), hygro_idx)
          84           0 :   call pbuf_add_field('DRYVOL',    'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), dryvol_idx)
          85           0 :   call pbuf_add_field('DRYRAD',    'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), dryrad_idx)
          86           0 :   call pbuf_add_field('DRYMASS',   'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), drymass_idx)
          87           0 :   call pbuf_add_field('SO4DRYVOL', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), so4dryvol_idx)
          88           0 :   call pbuf_add_field('NAER',      'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), naer_idx)
          89             :
          90           0 : end subroutine modal_aero_calcsize_reg
          91             :
          92             : !===============================================================================
          93             : !===============================================================================
          94             :
          95           0 : subroutine modal_aero_calcsize_init(pbuf2d)
          96           0 :    use time_manager,  only: is_first_step
          97             :    use physics_buffer,only: pbuf_set_field
          98             :
          99             :    !-----------------------------------------------------------------------
         100             :    !
         101             :    ! Purpose:
         102             :    !    set do_adjust_default and do_aitacc_transfer_default flags
         103             :    !    create history fields for column tendencies associated with
         104             :    !       modal_aero_calcsize
         105             :    !
         106             :    ! Author: R. Easter
         107             :    !
         108             :    !-----------------------------------------------------------------------
         109             :
         110             :    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         111             :
         112             :    ! local
         113             :    integer  :: ipair, iq
         114             :    integer  :: jac
         115             :    integer  :: lsfrm, lstoo
         116             :    integer  :: n, nacc, nait
         117             :    logical  :: history_aerosol
         118             :
         119             :    character(len=fieldname_len)   :: tmpnamea, tmpnameb
         120             :    character(len=fieldname_len+3) :: fieldname
         121             :    character(128)                 :: long_name
         122             :    character(8)                   :: unit
         123             :    !-----------------------------------------------------------------------
         124             :
         125           0 :    call phys_getopts(history_aerosol_out=history_aerosol)
         126             :
         127             :    ! init entities required for both prescribed and prognostic modes
         128             :
         129           0 :    if (is_first_step()) then
         130             :       ! initialize fields in physics buffer
         131           0 :       call pbuf_set_field(pbuf2d, dgnum_idx, 0.0_r8)
         132             :    endif
         133             :
         134             : #ifndef MODAL_AERO
         135           0 :    do_adjust_default          = .false.
         136           0 :    do_aitacc_transfer_default = .false.
         137             : #else
         138             :    !  do_adjust_default allows adjustment to be turned on/off
         139             :    do_adjust_default = .true.
         140             :
         141             :    !  do_aitacc_transfer_default allows aitken <--> accum mode transfer to be turned on/off
         142             :    !  *** it can only be true when aitken & accum modes are both present
         143             :    !      and have prognosed number and diagnosed surface/sigmag
         144             :    nait = modeptr_aitken
         145             :    nacc = modeptr_accum
         146             :    do_aitacc_transfer_default = .false.
         147             :    if ((modeptr_aitken > 0) .and.   &
         148             :       (modeptr_accum  > 0) .and.   &
         149             :       (modeptr_aitken /= modeptr_accum)) then
         150             :       do_aitacc_transfer_default = .true.
         151             :       if (mprognum_amode(nait) <= 0) do_aitacc_transfer_default = .false.
         152             :       if (mprognum_amode(nacc) <= 0) do_aitacc_transfer_default = .false.
         153             :    end if
         154             :
         155             :    if ( .not. do_adjust_default ) return
         156             :
         157             :    !  define history fields for number-adjust source-sink for all modes
         158             :    do n = 1, ntot_amode
         159             :       if (mprognum_amode(n) <= 0) cycle
         160             :
         161             :       do jac = 1, 2
         162             :          if (jac == 1) then
         163             :             tmpnamea = cnst_name(numptr_amode(n))
         164             :          else
         165             :             tmpnamea = cnst_name_cw(numptrcw_amode(n))
         166             :          end if
         167             :          unit = '#/m2/s'
         168             :          fieldname = trim(tmpnamea) // '_sfcsiz1'
         169             :          long_name = trim(tmpnamea) // ' calcsize number-adjust column source'
         170             :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         171             :          if (history_aerosol) then
         172             :             call add_default(fieldname, 1, ' ')
         173             :          end if
         174             :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         175             :
         176             :          fieldname = trim(tmpnamea) // '_sfcsiz2'
         177             :          long_name = trim(tmpnamea) // ' calcsize number-adjust column sink'
         178             :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         179             :          if (history_aerosol) then
         180             :             call add_default(fieldname, 1, ' ')
         181             :          end if
         182             :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         183             :       end do   ! jac = ...
         184             :    end do   ! n = ...
         185             :
         186             :    if ( .not. do_aitacc_transfer_default ) return
         187             :
         188             :    ! check that renaming ipair=1 is aitken-->accum
         189             :    ipair = 1
         190             :    if ((modefrm_renamexf(ipair) .ne. nait) .or.   &
         191             :       (modetoo_renamexf(ipair) .ne. nacc)) then
         192             :       write( 6, '(//2a//)' )   &
         193             :          '*** modal_aero_calcaersize_init error -- ',   &
         194             :          'modefrm/too_renamexf(1) are wrong'
         195             :       call endrun( 'modal_aero_calcaersize_init error' )
         196             :    end if
         197             :
         198             :    ! define history fields for aitken-accum transfer
         199             :    do iq = 1, nspecfrm_renamexf(ipair)
         200             :
         201             :       ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
         202             :       do jac = 1, 2
         203             :
         204             :          ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species
         205             :          ! the lspectooa_renamexf (and lspectooc_renamexf) are accum  species
         206             :          if (jac .eq. 1) then
         207             :             lsfrm = lspecfrma_renamexf(iq,ipair)
         208             :             lstoo = lspectooa_renamexf(iq,ipair)
         209             :          else
         210             :             lsfrm = lspecfrmc_renamexf(iq,ipair)
         211             :             lstoo = lspectooc_renamexf(iq,ipair)
         212             :          end if
         213             :          if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
         214             :
         215             :          if (jac .eq. 1) then
         216             :             tmpnamea = cnst_name(lsfrm)
         217             :             tmpnameb = cnst_name(lstoo)
         218             :          else
         219             :             tmpnamea = cnst_name_cw(lsfrm)
         220             :             tmpnameb = cnst_name_cw(lstoo)
         221             :          end if
         222             :
         223             :          unit = 'kg/m2/s'
         224             :          if ((tmpnamea(1:3) == 'num') .or. &
         225             :             (tmpnamea(1:3) == 'NUM')) unit = '#/m2/s'
         226             :          fieldname = trim(tmpnamea) // '_sfcsiz3'
         227             :          long_name = trim(tmpnamea) // ' calcsize aitken-to-accum adjust column tendency'
         228             :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         229             :          if (history_aerosol) then
         230             :             call add_default(fieldname, 1, ' ')
         231             :          end if
         232             :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         233             :
         234             :          fieldname = trim(tmpnameb) // '_sfcsiz3'
         235             :          long_name = trim(tmpnameb) // ' calcsize aitken-to-accum adjust column tendency'
         236             :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         237             :          if (history_aerosol) then
         238             :             call add_default(fieldname, 1, ' ')
         239             :          end if
         240             :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         241             :
         242             :          fieldname = trim(tmpnamea) // '_sfcsiz4'
         243             :          long_name = trim(tmpnamea) // ' calcsize accum-to-aitken adjust column tendency'
         244             :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         245             :          if (history_aerosol) then
         246             :             call add_default(fieldname, 1, ' ')
         247             :          end if
         248             :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         249             :
         250             :          fieldname = trim(tmpnameb) // '_sfcsiz4'
         251             :          long_name = trim(tmpnameb) // ' calcsize accum-to-aitken adjust column tendency'
         252             :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         253             :          if (history_aerosol) then
         254             :             call add_default(fieldname, 1, ' ')
         255             :          end if
         256             :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         257             :
         258             :       end do   ! jac = ...
         259             :    end do   ! iq = ...
         260             :
         261             : #endif
         262             :
         263           0 : end subroutine modal_aero_calcsize_init
         264             :
         265             : !===============================================================================
         266             :
         267           0 : subroutine modal_aero_calcsize_sub(state, ptend, deltat, pbuf, do_adjust_in, &
         268             :    do_aitacc_transfer_in)
         269             :
         270             :    !-----------------------------------------------------------------------
         271             :    !
         272             :    ! Calculates aerosol size distribution parameters
         273             :    !    mprognum_amode >  0
         274             :    !       calculate Dgnum from mass, number, and fixed sigmag
         275             :    !    mprognum_amode <= 0
         276             :    !       calculate number from mass, fixed Dgnum, and fixed sigmag
         277             :    !
         278             :    ! Also (optionally) adjusts prognostic number to
         279             :    !    be within bounds determined by mass, Dgnum bounds, and sigma bounds
         280             :    !
         281             :    ! Author: R. Easter
         282             :    !
         283             :    !-----------------------------------------------------------------------
         284             :
         285             :    ! arguments
         286             :    type(physics_state), target, intent(in)    :: state       ! Physics state variables
         287             :
         288             :    type(physics_ptend), target, intent(inout)   :: ptend       ! indivdual parameterization tendencies
         289             :
         290             :    real(r8),                    intent(in)    :: deltat      ! model time-step size (s)
         291             :    type(physics_buffer_desc),   pointer       :: pbuf(:)     ! physics buffer
         292             :
         293             :    logical, optional :: do_adjust_in
         294             :    logical, optional :: do_aitacc_transfer_in
         295             :
         296             : #ifdef MODAL_AERO
         297             :
         298             :    ! local
         299             :
         300             :    logical :: do_adjust
         301             :    logical :: do_aitacc_transfer
         302             :
         303             :    integer  :: lchnk                ! chunk identifier
         304             :    integer  :: ncol                 ! number of columns
         305             :
         306             :    real(r8), pointer :: t(:,:)      ! Temperature in Kelvin
         307             :    real(r8), pointer :: pmid(:,:)   ! pressure at model levels (Pa)
         308             :    real(r8), pointer :: pdel(:,:)   ! pressure thickness of levels
         309             :    real(r8), pointer :: q(:,:,:)    ! Tracer MR array
         310             :
         311             :    logical,  pointer :: dotend(:)   ! flag for doing tendency
         312             :    real(r8), pointer :: dqdt(:,:,:) ! TMR tendency array
         313             :
         314             :    real(r8), pointer :: dgncur_a(:,:,:)
         315             :
         316             :    integer  :: i, icol_diag, iduma, ipair, iq
         317             :    integer  :: ixfer_acc2ait, ixfer_ait2acc
         318             :    integer  :: ixfer_acc2ait_sv(pcols,pver), ixfer_ait2acc_sv(pcols,pver)
         319             :    integer  :: j, jac, jsrflx, k
         320             :    integer  :: l, l1, la, lc, lna, lnc, lsfrm, lstoo
         321             :    integer  :: n, nacc, nait
         322             :
         323             :    integer, save  :: idiagaa = 1
         324             :
         325             :    logical  :: dotendqqcw(pcnst)
         326             :    logical  :: noxf_acc2ait(nspec_max)
         327             :
         328             :    character(len=fieldname_len)   :: tmpnamea, tmpnameb
         329             :    character(len=fieldname_len+3) :: fieldname
         330             :
         331             :    real(r8), parameter :: third = 1.0_r8/3.0_r8
         332             :    real(r8), pointer :: fldcw(:,:)
         333             :    real(r8) :: delnum_a2, delnum_c2            !  work variables
         334             :    real(r8) :: delnum_a3, delnum_c3, delnum_t3 !  work variables
         335             :    real(r8) :: deltatinv                     ! 1/deltat
         336             :    real(r8) :: dgncur_c(pcols,pver,ntot_amode)
         337             :    real(r8) :: dgnyy, dgnxx                  ! dgnumlo/hi of current mode
         338             :    real(r8) :: dqqcwdt(pcols,pver,pcnst)     ! cloudborne TMR tendency array
         339             :    real(r8) :: drv_a, drv_c, drv_t           ! dry volume (cm3/mol_air)
         340             :    real(r8) :: drv_t0
         341             :    real(r8) :: drv_a_noxf, drv_c_noxf, drv_t_noxf
         342             :    real(r8) :: drv_a_acc, drv_c_acc
         343             :    real(r8) :: drv_a_accsv(pcols,pver), drv_c_accsv(pcols,pver)
         344             :    real(r8) :: drv_a_aitsv(pcols,pver), drv_c_aitsv(pcols,pver)
         345             :    real(r8) :: drv_a_sv(pcols,pver,ntot_amode), drv_c_sv(pcols,pver,ntot_amode)
         346             :    real(r8) :: dryvol_a(pcols,pver)          ! interstital aerosol dry
         347             :    ! volume (cm^3/mol_air)
         348             :    real(r8) :: dryvol_c(pcols,pver)          ! activated aerosol dry volume
         349             :    real(r8) :: duma, dumb, dumc, dumd        ! work variables
         350             :    real(r8) :: dumfac, dummwdens             ! work variables
         351             :    real(r8) :: frelaxadj                     ! relaxation factor applied
         352             :    ! to size bounds
         353             :    real(r8) :: fracadj                       ! deltat/tadj
         354             :    real(r8) :: num_a0, num_c0, num_t0        ! initial number (#/mol_air)
         355             :    real(r8) :: num_a1, num_c1                ! working number (#/mol_air)
         356             :    real(r8) :: num_a2, num_c2, num_t2        ! working number (#/mol_air)
         357             :    real(r8) :: num_a, num_c, num_t           ! final number (#/mol_air)
         358             :    real(r8) :: num_t_noxf
         359             :    real(r8) :: numbnd                        ! bounded number
         360             :    real(r8) :: num_a_acc, num_c_acc
         361             :    real(r8) :: num_a_accsv(pcols,pver), num_c_accsv(pcols,pver)
         362             :    real(r8) :: num_a_aitsv(pcols,pver), num_c_aitsv(pcols,pver)
         363             :    real(r8) :: num_a_sv(pcols,pver,ntot_amode), num_c_sv(pcols,pver,ntot_amode)
         364             :    real(r8) :: pdel_fac                      !
         365             :    real(r8) :: tadj                          ! adjustment time scale
         366             :    real(r8) :: tadjinv                       ! 1/tadj
         367             :    real(r8) :: v2ncur_a(pcols,pver,ntot_amode)
         368             :    real(r8) :: v2ncur_c(pcols,pver,ntot_amode)
         369             :    real(r8) :: v2nyy, v2nxx, v2nzz           ! voltonumblo/hi of current mode
         370             :    real(r8) :: v2nyyrl, v2nxxrl              ! relaxed voltonumblo/hi
         371             :    real(r8) :: xfercoef
         372             :    real(r8) :: xfercoef_num_acc2ait, xfercoef_vol_acc2ait
         373             :    real(r8) :: xfercoef_num_ait2acc, xfercoef_vol_ait2acc
         374             :    real(r8) :: xferfrac_num_acc2ait, xferfrac_vol_acc2ait
         375             :    real(r8) :: xferfrac_num_ait2acc, xferfrac_vol_ait2acc
         376             :    real(r8) :: xfertend, xfertend_num(2,2)
         377             :
         378             :    integer, parameter :: nsrflx = 4    ! last dimension of qsrflx
         379             :    real(r8) :: qsrflx(pcols,pcnst,nsrflx,2)
         380             :    ! process-specific column tracer tendencies
         381             :    ! 3rd index --
         382             :    !    1="standard" number adjust gain;
         383             :    !    2="standard" number adjust loss;
         384             :    !    3=aitken-->accum renaming; 4=accum-->aitken)
         385             :    ! 4th index --
         386             :    !    1="a" species; 2="c" species
         387             :    !-----------------------------------------------------------------------
         388             :
         389             :    if (present(do_adjust_in)) then
         390             :       do_adjust = do_adjust_in
         391             :    else
         392             :       do_adjust = do_adjust_default
         393             :    end if
         394             :
         395             :    if (present(do_aitacc_transfer_in)) then
         396             :       do_aitacc_transfer = do_aitacc_transfer_in
         397             :    else
         398             :       do_aitacc_transfer = do_aitacc_transfer_default
         399             :    end if
         400             :
         401             :    lchnk = state%lchnk
         402             :    ncol  = state%ncol
         403             :
         404             :    t    => state%t
         405             :    pmid => state%pmid
         406             :    pdel => state%pdel
         407             :    q    => state%q
         408             :
         409             :    dotend => ptend%lq
         410             :    dqdt   => ptend%q
         411             :
         412             :    call pbuf_get_field(pbuf, dgnum_idx, dgncur_a)
         413             :
         414             :    dotendqqcw(:) = .false.
         415             :    dqqcwdt(:,:,:) = 0.0_r8
         416             :    qsrflx(:,:,:,:) = 0.0_r8
         417             :
         418             :    nait = modeptr_aitken
         419             :    nacc = modeptr_accum
         420             :
         421             :    deltatinv = 1.0_r8/(deltat*(1.0_r8 + 1.0e-15_r8))
         422             :    ! tadj = adjustment time scale for number, surface when they are prognosed
         423             :    !           currently set to deltat
         424             :    tadj = deltat
         425             :    tadj = 86400
         426             :    tadj = max( tadj, deltat )
         427             :    tadjinv = 1.0_r8/(tadj*(1.0_r8 + 1.0e-15_r8))
         428             :    fracadj = deltat*tadjinv
         429             :    fracadj = max( 0.0_r8, min( 1.0_r8, fracadj ) )
         430             :
         431             :
         432             :    !
         433             :    !
         434             :    ! the "do 40000" loop does the original (pre jan-2006)
         435             :    !   number adjustment, one mode at a time
         436             :    ! this artificially adjusts number when mean particle size is too large
         437             :    !   or too small
         438             :    !
         439             :    !
         440             :    do n = 1, ntot_amode
         441             :
         442             :
         443             :       ! initialize all parameters to the default values for the mode
         444             :       do k=top_lev,pver
         445             :          do i=1,ncol
         446             :             !    sgcur_a(i,k,n) = sigmag_amode(n)
         447             :             !    sgcur_c(i,k,n) = sigmag_amode(n)
         448             :             dgncur_a(i,k,n) = dgnum_amode(n)
         449             :             dgncur_c(i,k,n) = dgnum_amode(n)
         450             :             v2ncur_a(i,k,n) = voltonumb_amode(n)
         451             :             v2ncur_c(i,k,n) = voltonumb_amode(n)
         452             :             dryvol_a(i,k) = 0.0_r8
         453             :             dryvol_c(i,k) = 0.0_r8
         454             :          end do
         455             :       end do
         456             :
         457             :       ! compute dry volume mixrats =
         458             :       !      sum_over_components{ component_mass mixrat / density }
         459             :       do l1 = 1, nspec_amode(n)
         460             :          ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
         461             :          dummwdens = 1.0_r8 / specdens_amode(l1,n)
         462             :          la = lmassptr_amode(l1,n)
         463             :          do k=top_lev,pver
         464             :             do i=1,ncol
         465             :                dryvol_a(i,k) = dryvol_a(i,k)    &
         466             :                   + max(0.0_r8,q(i,k,la))*dummwdens
         467             :             end do
         468             :          end do
         469             :
         470             :          fldcw => qqcw_get_field(pbuf,lmassptrcw_amode(l1,n),lchnk)
         471             :          do k=top_lev,pver
         472             :             do i=1,ncol
         473             :                dryvol_c(i,k) = dryvol_c(i,k)    &
         474             :                   + max(0.0_r8,fldcw(i,k))*dummwdens
         475             :             end do
         476             :          end do
         477             :       end do
         478             :
         479             :       ! set "short-hand" number pointers
         480             :       lna = numptr_amode(n)
         481             :       lnc = numptrcw_amode(n)
         482             :       fldcw => qqcw_get_field(pbuf,numptrcw_amode(n),lchnk,.true.)
         483             :
         484             :
         485             :       ! go to section for appropriate number/surface diagnosed/prognosed options
         486             :       if (mprognum_amode(n) <= 0) then
         487             :
         488             :          ! option 1 -- number diagnosed (fixed dgnum and sigmag)
         489             :          !    compute number tendencies that will bring numbers to their
         490             :          !    current diagnosed values
         491             :          !
         492             :          if (lna > 0) then
         493             :             dotend(lna) = .true.
         494             :             do k=top_lev,pver
         495             :                do i=1,ncol
         496             :                   dqdt(i,k,lna) = (dryvol_a(i,k)*voltonumb_amode(n)   &
         497             :                      - q(i,k,lna)) * deltatinv
         498             :                end do
         499             :             end do
         500             :          end if
         501             :          if (lnc > 0) then
         502             :             dotendqqcw(lnc) = .true.
         503             :             do k=top_lev,pver
         504             :                do i=1,ncol
         505             :                   dqqcwdt(i,k,lnc) = (dryvol_c(i,k)*voltonumb_amode(n)   &
         506             :                      - fldcw(i,k)) * deltatinv
         507             :                end do
         508             :             end do
         509             :          end if
         510             :       else
         511             :
         512             :
         513             :          !
         514             :          ! option 2 -- number prognosed (variable dgnum, fixed sigmag)
         515             :          !       Compute number tendencies to adjust numbers if they are outside
         516             :          !    the limits determined by current volume and dgnumlo/hi
         517             :          !       The interstitial and activated aerosol fractions can, at times,
         518             :          !    be the lower or upper tail of the "total" distribution.  Thus they
         519             :          !    can be expected to have a greater range of size parameters than
         520             :          !    what is specified for the total distribution (via dgnumlo/hi)
         521             :          !       When both the interstitial and activated dry volumes are positive,
         522             :          !    the adjustment strategy is to (1) adjust the interstitial and activated
         523             :          !    numbers towards relaxed bounds, then (2) adjust the total/combined
         524             :          !    number towards the primary bounds.
         525             :          !
         526             :          ! note
         527             :          !    v2nyy = voltonumblo_amode is proportional to dgnumlo**(-3),
         528             :          !            and produces the maximum allowed number for a given volume
         529             :          !    v2nxx = voltonumbhi_amode is proportional to dgnumhi**(-3),
         530             :          !            and produces the minimum allowed number for a given volume
         531             :          !    v2nxxrl and v2nyyrl are their "relaxed" equivalents.
         532             :          !            Setting frelaxadj=27=3**3 means that
         533             :          !            dgnumlo_relaxed = dgnumlo/3 and dgnumhi_relaxed = dgnumhi*3
         534             :          !
         535             :          ! if do_aitacc_transfer is .true., then
         536             :          !     for n=nacc, multiply v2nyy by 1.0e6 to effectively turn off the
         537             :          !         adjustment when number is too big (size is too small)
         538             :          !     for n=nait, divide   v2nxx by 1.0e6 to effectively turn off the
         539             :          !         adjustment when number is too small (size is too big)
         540             :          !OLD  however, do not change the v2nyyrl/v2nxxrl so that
         541             :          !OLD      the interstitial<-->activated adjustment is not changed
         542             :          !NEW  also change the v2nyyrl/v2nxxrl so that
         543             :          !NEW      the interstitial<-->activated adjustment is turned off
         544             :          !
         545             :       end if
         546             :       frelaxadj = 27.0_r8
         547             :       dumfac = exp(4.5_r8*alnsg_amode(n)**2)*pi/6.0_r8
         548             :       v2nxx = voltonumbhi_amode(n)
         549             :       v2nyy = voltonumblo_amode(n)
         550             :       v2nxxrl = v2nxx/frelaxadj
         551             :       v2nyyrl = v2nyy*frelaxadj
         552             :       dgnxx = dgnumhi_amode(n)
         553             :       dgnyy = dgnumlo_amode(n)
         554             :       if ( do_aitacc_transfer ) then
         555             :          if (n == nait) v2nxx = v2nxx/1.0e6_r8
         556             :          if (n == nacc) v2nyy = v2nyy*1.0e6_r8
         557             :          v2nxxrl = v2nxx/frelaxadj   ! NEW
         558             :          v2nyyrl = v2nyy*frelaxadj   ! NEW
         559             :       end if
         560             :
         561             :       if (do_adjust) then
         562             :          dotend(lna) = .true.
         563             :          dotendqqcw(lnc) = .true.
         564             :       end if
         565             :
         566             :       do  k = top_lev, pver
         567             :          do  i = 1, ncol
         568             :
         569             :             drv_a = dryvol_a(i,k)
         570             :             num_a0 = q(i,k,lna)
         571             :             num_a = max( 0.0_r8, num_a0 )
         572             :             drv_c = dryvol_c(i,k)
         573             :             num_c0 = fldcw(i,k)
         574             :             num_c = max( 0.0_r8, num_c0 )
         575             :
         576             :             if ( do_adjust) then
         577             :
         578             :                !
         579             :                ! do number adjustment for interstitial and activated particles
         580             :                !    adjustments that (1) make numbers non-negative or (2) make numbers
         581             :                !       zero when volume is zero are applied over time-scale deltat
         582             :                !    adjustments that bring numbers to within specified bounds are
         583             :                !       applied over time-scale tadj
         584             :                !
         585             :                if ((drv_a <= 0.0_r8) .and. (drv_c <= 0.0_r8)) then
         586             :                   ! both interstitial and activated volumes are zero
         587             :                   ! adjust both numbers to zero
         588             :                   num_a = 0.0_r8
         589             :                   dqdt(i,k,lna) = -num_a0*deltatinv
         590             :                   num_c = 0.0_r8
         591             :                   dqqcwdt(i,k,lnc) = -num_c0*deltatinv
         592             :                else if (drv_c <= 0.0_r8) then
         593             :                   ! activated volume is zero, so interstitial number/volume == total/combined
         594             :                   ! apply step 1 and 3, but skip the relaxed adjustment (step 2, see below)
         595             :                   num_c = 0.0_r8
         596             :                   dqqcwdt(i,k,lnc) = -num_c0*deltatinv
         597             :                   num_a1 = num_a
         598             :                   numbnd = max( drv_a*v2nxx, min( drv_a*v2nyy, num_a1 ) )
         599             :                   num_a  = num_a1 + (numbnd - num_a1)*fracadj
         600             :                   dqdt(i,k,lna) = (num_a - num_a0)*deltatinv
         601             :
         602             :                else if (drv_a <= 0.0_r8) then
         603             :                   ! interstitial volume is zero, treat similar to above
         604             :                   num_a = 0.0_r8
         605             :                   dqdt(i,k,lna) = -num_a0*deltatinv
         606             :                   num_c1 = num_c
         607             :                   numbnd = max( drv_c*v2nxx, min( drv_c*v2nyy, num_c1 ) )
         608             :                   num_c  = num_c1 + (numbnd - num_c1)*fracadj
         609             :                   dqqcwdt(i,k,lnc) = (num_c - num_c0)*deltatinv
         610             :                else
         611             :                   ! both volumes are positive
         612             :                   ! apply 3 adjustment steps
         613             :                   ! step1:  num_a,c0 --> num_a,c1 forces non-negative values
         614             :                   num_a1 = num_a
         615             :                   num_c1 = num_c
         616             :                   ! step2:  num_a,c1 --> num_a,c2 applies relaxed bounds to the interstitial
         617             :                   !    and activated number (individually)
         618             :                   !    if only only a or c changes, adjust the other in the opposite direction
         619             :                   !    as much as possible to conserve a+c
         620             :                   numbnd = max( drv_a*v2nxxrl, min( drv_a*v2nyyrl, num_a1 ) )
         621             :                   delnum_a2 = (numbnd - num_a1)*fracadj
         622             :                   num_a2 = num_a1 + delnum_a2
         623             :                   numbnd = max( drv_c*v2nxxrl, min( drv_c*v2nyyrl, num_c1 ) )
         624             :                   delnum_c2 = (numbnd - num_c1)*fracadj
         625             :                   num_c2 = num_c1 + delnum_c2
         626             :                   if ((delnum_a2 == 0.0_r8) .and. (delnum_c2 /= 0.0_r8)) then
         627             :                      num_a2 = max( drv_a*v2nxxrl, min( drv_a*v2nyyrl,   &
         628             :                         num_a1-delnum_c2 ) )
         629             :                   else if ((delnum_a2 /= 0.0_r8) .and. (delnum_c2 == 0.0_r8)) then
         630             :                      num_c2 = max( drv_c*v2nxxrl, min( drv_c*v2nyyrl,   &
         631             :                         num_c1-delnum_a2 ) )
         632             :                   end if
         633             :                   ! step3:  num_a,c2 --> num_a,c3 applies stricter bounds to the
         634             :                   !    combined/total number
         635             :                   drv_t = drv_a + drv_c
         636             :                   num_t2 = num_a2 + num_c2
         637             :                   delnum_a3 = 0.0_r8
         638             :                   delnum_c3 = 0.0_r8
         639             :                   if (num_t2 < drv_t*v2nxx) then
         640             :                      delnum_t3 = (drv_t*v2nxx - num_t2)*fracadj
         641             :                      ! if you are here then (num_a2 < drv_a*v2nxx) and/or
         642             :                      !                      (num_c2 < drv_c*v2nxx) must be true
         643             :                      if ((num_a2 < drv_a*v2nxx) .and. (num_c2 < drv_c*v2nxx)) then
         644             :                         delnum_a3 = delnum_t3*(num_a2/num_t2)
         645             :                         delnum_c3 = delnum_t3*(num_c2/num_t2)
         646             :                      else if (num_c2 < drv_c*v2nxx) then
         647             :                         delnum_c3 = delnum_t3
         648             :                      else if (num_a2 < drv_a*v2nxx) then
         649             :                         delnum_a3 = delnum_t3
         650             :                      end if
         651             :                   else if (num_t2 > drv_t*v2nyy) then
         652             :                      delnum_t3 = (drv_t*v2nyy - num_t2)*fracadj
         653             :                      ! if you are here then (num_a2 > drv_a*v2nyy) and/or
         654             :                      !                      (num_c2 > drv_c*v2nyy) must be true
         655             :                      if ((num_a2 > drv_a*v2nyy) .and. (num_c2 > drv_c*v2nyy)) then
         656             :                         delnum_a3 = delnum_t3*(num_a2/num_t2)
         657             :                         delnum_c3 = delnum_t3*(num_c2/num_t2)
         658             :                      else if (num_c2 > drv_c*v2nyy) then
         659             :                         delnum_c3 = delnum_t3
         660             :                      else if (num_a2 > drv_a*v2nyy) then
         661             :                         delnum_a3 = delnum_t3
         662             :                      end if
         663             :                   end if
         664             :                   num_a = num_a2 + delnum_a3
         665             :                   dqdt(i,k,lna) = (num_a - num_a0)*deltatinv
         666             :                   num_c = num_c2 + delnum_c3
         667             :                   dqqcwdt(i,k,lnc) = (num_c - num_c0)*deltatinv
         668             :                end if
         669             :
         670             :             end if ! do_adjust
         671             :
         672             :             !
         673             :             ! now compute current dgn and v2n
         674             :             !
         675             :             if (drv_a > 0.0_r8) then
         676             :                if (num_a <= drv_a*v2nxx) then
         677             :                   dgncur_a(i,k,n) = dgnxx
         678             :                   v2ncur_a(i,k,n) = v2nxx
         679             :                else if (num_a >= drv_a*v2nyy) then
         680             :                   dgncur_a(i,k,n) = dgnyy
         681             :                   v2ncur_a(i,k,n) = v2nyy
         682             :                else
         683             :                   dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third
         684             :                   v2ncur_a(i,k,n) = num_a/drv_a
         685             :                end if
         686             :             end if
         687             :             pdel_fac = pdel(i,k)/gravit   ! = rho*dz
         688             :             jac = 1
         689             :             qsrflx(i,lna,1,jac) = qsrflx(i,lna,1,jac) + max(0.0_r8,dqdt(i,k,lna))*pdel_fac
         690             :             qsrflx(i,lna,2,jac) = qsrflx(i,lna,2,jac) + min(0.0_r8,dqdt(i,k,lna))*pdel_fac
         691             :
         692             :             if (drv_c > 0.0_r8) then
         693             :                if (num_c <= drv_c*v2nxx) then
         694             :                   dgncur_c(i,k,n) = dgnumhi_amode(n)
         695             :                   v2ncur_c(i,k,n) = v2nxx
         696             :                else if (num_c >= drv_c*v2nyy) then
         697             :                   dgncur_c(i,k,n) = dgnumlo_amode(n)
         698             :                   v2ncur_c(i,k,n) = v2nyy
         699             :                else
         700             :                   dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third
         701             :                   v2ncur_c(i,k,n) = num_c/drv_c
         702             :                end if
         703             :             end if
         704             :             jac = 2
         705             :             qsrflx(i,lnc,1,jac) = qsrflx(i,lnc,1,jac) + max(0.0_r8,dqqcwdt(i,k,lnc))*pdel_fac
         706             :             qsrflx(i,lnc,2,jac) = qsrflx(i,lnc,2,jac) + min(0.0_r8,dqqcwdt(i,k,lnc))*pdel_fac
         707             :
         708             :
         709             :             ! save number and dryvol for aitken <--> accum renaming
         710             :             if ( do_aitacc_transfer ) then
         711             :                if (n == nait) then
         712             :                   drv_a_aitsv(i,k) = drv_a
         713             :                   num_a_aitsv(i,k) = num_a
         714             :                   drv_c_aitsv(i,k) = drv_c
         715             :                   num_c_aitsv(i,k) = num_c
         716             :                else if (n == nacc) then
         717             :                   drv_a_accsv(i,k) = drv_a
         718             :                   num_a_accsv(i,k) = num_a
         719             :                   drv_c_accsv(i,k) = drv_c
         720             :                   num_c_accsv(i,k) = num_c
         721             :                end if
         722             :             end if
         723             :             drv_a_sv(i,k,n) = drv_a
         724             :             num_a_sv(i,k,n) = num_a
         725             :             drv_c_sv(i,k,n) = drv_c
         726             :             num_c_sv(i,k,n) = num_c
         727             :
         728             :          end do
         729             :       end do
         730             :
         731             :
         732             :       !
         733             :       ! option 3 -- number and surface prognosed (variable dgnum and sigmag)
         734             :       !             this is not implemented
         735             :       !
         736             :    end do  ! do n = 1, ntot_amode
         737             :
         738             :
         739             :    !
         740             :    !
         741             :    ! the following section (from here to label 49000)
         742             :    !    does aitken <--> accum mode transfer
         743             :    !
         744             :    ! when the aitken mode mean size is too big, the largest
         745             :    !    aitken particles are transferred into the accum mode
         746             :    !    to reduce the aitken mode mean size
         747             :    ! when the accum mode mean size is too small, the smallest
         748             :    !    accum particles are transferred into the aitken mode
         749             :    !    to increase the accum mode mean size
         750             :    !
         751             :    !
         752             :    ixfer_ait2acc_sv(:,:) = 0
         753             :    ixfer_acc2ait_sv(:,:) = 0
         754             :    if ( do_aitacc_transfer ) then
         755             :
         756             :       ! old - on time first step, npair_renamexf will be <= 0,
         757             :       !       in which case need to do modal_aero_rename_init
         758             :       ! new - init is now done through chem_init and things below it
         759             :       if (npair_renamexf .le. 0) then
         760             :          npair_renamexf = 0
         761             :          !        call modal_aero_rename_init
         762             :          if (npair_renamexf .le. 0) then
         763             :             write( 6, '(//a//)' )   &
         764             :                '*** modal_aero_calcaersize_sub error -- npair_renamexf <= 0'
         765             :             call endrun( 'modal_aero_calcaersize_sub error' )
         766             :          end if
         767             :       end if
         768             :
         769             :       ! check that renaming ipair=1 is aitken-->accum
         770             :       ipair = 1
         771             :       if ((modefrm_renamexf(ipair) .ne. nait) .or.   &
         772             :          (modetoo_renamexf(ipair) .ne. nacc)) then
         773             :          write( 6, '(//2a//)' )   &
         774             :             '*** modal_aero_calcaersize_sub error -- ',   &
         775             :             'modefrm/too_renamexf(1) are wrong'
         776             :          call endrun( 'modal_aero_calcaersize_sub error' )
         777             :       end if
         778             :
         779             :       ! set dotend() for species that will be transferred
         780             :       do iq = 1, nspecfrm_renamexf(ipair)
         781             :          lsfrm = lspecfrma_renamexf(iq,ipair)
         782             :          lstoo = lspectooa_renamexf(iq,ipair)
         783             :          if ((lsfrm > 0) .and. (lstoo > 0)) then
         784             :             dotend(lsfrm) = .true.
         785             :             dotend(lstoo) = .true.
         786             :          end if
         787             :          lsfrm = lspecfrmc_renamexf(iq,ipair)
         788             :          lstoo = lspectooc_renamexf(iq,ipair)
         789             :          if ((lsfrm > 0) .and. (lstoo > 0)) then
         790             :             dotendqqcw(lsfrm) = .true.
         791             :             dotendqqcw(lstoo) = .true.
         792             :          end if
         793             :       end do
         794             :
         795             :       ! identify accum species cannot be transferred to aitken mode
         796             :       noxf_acc2ait(:) = .true.
         797             :       do l1 = 1, nspec_amode(nacc)
         798             :          la = lmassptr_amode(l1,nacc)
         799             :          do iq = 1, nspecfrm_renamexf(ipair)
         800             :             if (lspectooa_renamexf(iq,ipair) == la) then
         801             :                noxf_acc2ait(l1) = .false.
         802             :             end if
         803             :          end do
         804             :       end do
         805             :
         806             :       ! v2nzz is voltonumb at the "geometrically-defined" mid-point
         807             :       ! between the aitken and accum modes
         808             :       v2nzz = sqrt(voltonumb_amode(nait)*voltonumb_amode(nacc))
         809             :
         810             :       ! loop over columns and levels
         811             :       do  k = top_lev, pver
         812             :          do  i = 1, ncol
         813             :
         814             :             pdel_fac = pdel(i,k)/gravit   ! = rho*dz
         815             :             xfertend_num(:,:) = 0.0_r8
         816             :
         817             :             ! compute aitken --> accum transfer rates
         818             :             ixfer_ait2acc = 0
         819             :             xfercoef_num_ait2acc = 0.0_r8
         820             :             xfercoef_vol_ait2acc = 0.0_r8
         821             :
         822             :             drv_t = drv_a_aitsv(i,k) + drv_c_aitsv(i,k)
         823             :             num_t = num_a_aitsv(i,k) + num_c_aitsv(i,k)
         824             :             if (drv_t > 0.0_r8) then
         825             :                if (num_t < drv_t*v2nzz) then
         826             :                   ixfer_ait2acc = 1
         827             :                   if (num_t < drv_t*voltonumb_amode(nacc)) then
         828             :                      xferfrac_num_ait2acc = 1.0_r8
         829             :                      xferfrac_vol_ait2acc = 1.0_r8
         830             :                   else
         831             :                      xferfrac_vol_ait2acc = ((num_t/drv_t) - v2nzz)/   &
         832             :                         (voltonumb_amode(nacc) - v2nzz)
         833             :                      xferfrac_num_ait2acc = xferfrac_vol_ait2acc*   &
         834             :                         (drv_t*voltonumb_amode(nacc)/num_t)
         835             :                      if ((xferfrac_num_ait2acc <= 0.0_r8) .or.   &
         836             :                         (xferfrac_vol_ait2acc <= 0.0_r8)) then
         837             :                         xferfrac_num_ait2acc = 0.0_r8
         838             :                         xferfrac_vol_ait2acc = 0.0_r8
         839             :                      else if ((xferfrac_num_ait2acc >= 1.0_r8) .or.   &
         840             :                         (xferfrac_vol_ait2acc >= 1.0_r8)) then
         841             :                         xferfrac_num_ait2acc = 1.0_r8
         842             :                         xferfrac_vol_ait2acc = 1.0_r8
         843             :                      end if
         844             :                   end if
         845             :                   xfercoef_num_ait2acc = xferfrac_num_ait2acc*tadjinv
         846             :                   xfercoef_vol_ait2acc = xferfrac_vol_ait2acc*tadjinv
         847             :                   xfertend_num(1,1) = num_a_aitsv(i,k)*xfercoef_num_ait2acc
         848             :                   xfertend_num(1,2) = num_c_aitsv(i,k)*xfercoef_num_ait2acc
         849             :                end if
         850             :             end if
         851             :
         852             :             ! compute accum --> aitken transfer rates
         853             :             ! accum may have some species (seasalt, dust, poa, lll) that are
         854             :             !    not in aitken mode
         855             :             ! so first divide the accum drv & num into not-transferred (noxf) species
         856             :             !    and transferred species, and use the transferred-species
         857             :             !    portion in what follows
         858             :             ixfer_acc2ait = 0
         859             :             xfercoef_num_acc2ait = 0.0_r8
         860             :             xfercoef_vol_acc2ait = 0.0_r8
         861             :
         862             :             drv_t = drv_a_accsv(i,k) + drv_c_accsv(i,k)
         863             :             num_t = num_a_accsv(i,k) + num_c_accsv(i,k)
         864             :             drv_a_noxf = 0.0_r8
         865             :             drv_c_noxf = 0.0_r8
         866             :             if (drv_t > 0.0_r8) then
         867             :                if (num_t > drv_t*v2nzz) then
         868             :                   do l1 = 1, nspec_amode(nacc)
         869             :
         870             :                      if ( noxf_acc2ait(l1) ) then
         871             :                         ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
         872             :                         dummwdens = 1.0_r8 / specdens_amode(l1,nacc)
         873             :                         la = lmassptr_amode(l1,nacc)
         874             :                         drv_a_noxf = drv_a_noxf    &
         875             :                            + max(0.0_r8,q(i,k,la))*dummwdens
         876             :                         lc = lmassptrcw_amode(l1,nacc)
         877             :
         878             :                         fldcw => qqcw_get_field(pbuf,lmassptrcw_amode(l1,nacc),lchnk)
         879             :                         drv_c_noxf = drv_c_noxf    &
         880             :                            + max(0.0_r8,fldcw(i,k))*dummwdens
         881             :                      end if
         882             :                   end do
         883             :                   drv_t_noxf = drv_a_noxf + drv_c_noxf
         884             :                   num_t_noxf = drv_t_noxf*voltonumblo_amode(nacc)
         885             :                   num_t0 = num_t
         886             :                   drv_t0 = drv_t
         887             :                   num_t = max( 0.0_r8, num_t - num_t_noxf )
         888             :                   drv_t = max( 0.0_r8, drv_t - drv_t_noxf )
         889             :                end if
         890             :             end if
         891             :
         892             :             if (drv_t > 0.0_r8) then
         893             :                if (num_t > drv_t*v2nzz) then
         894             :                   ixfer_acc2ait = 1
         895             :                   if (num_t > drv_t*voltonumb_amode(nait)) then
         896             :                      xferfrac_num_acc2ait = 1.0_r8
         897             :                      xferfrac_vol_acc2ait = 1.0_r8
         898             :                   else
         899             :                      xferfrac_vol_acc2ait = ((num_t/drv_t) - v2nzz)/   &
         900             :                         (voltonumb_amode(nait) - v2nzz)
         901             :                      xferfrac_num_acc2ait = xferfrac_vol_acc2ait*   &
         902             :                         (drv_t*voltonumb_amode(nait)/num_t)
         903             :                      if ((xferfrac_num_acc2ait <= 0.0_r8) .or.   &
         904             :                         (xferfrac_vol_acc2ait <= 0.0_r8)) then
         905             :                         xferfrac_num_acc2ait = 0.0_r8
         906             :                         xferfrac_vol_acc2ait = 0.0_r8
         907             :                      else if ((xferfrac_num_acc2ait >= 1.0_r8) .or.   &
         908             :                         (xferfrac_vol_acc2ait >= 1.0_r8)) then
         909             :                         xferfrac_num_acc2ait = 1.0_r8
         910             :                         xferfrac_vol_acc2ait = 1.0_r8
         911             :                      end if
         912             :                   end if
         913             :                   duma = 1.0e-37_r8
         914             :                   xferfrac_num_acc2ait = xferfrac_num_acc2ait*   &
         915             :                      num_t/max( duma, num_t0 )
         916             :                   xfercoef_num_acc2ait = xferfrac_num_acc2ait*tadjinv
         917             :                   xfercoef_vol_acc2ait = xferfrac_vol_acc2ait*tadjinv
         918             :                   xfertend_num(2,1) = num_a_accsv(i,k)*xfercoef_num_acc2ait
         919             :                   xfertend_num(2,2) = num_c_accsv(i,k)*xfercoef_num_acc2ait
         920             :                end if
         921             :             end if
         922             :
         923             :             ! jump to end-of-loop if no transfer is needed at current i,k
         924             :             if (ixfer_ait2acc+ixfer_acc2ait > 0) then
         925             :                ixfer_ait2acc_sv(i,k) = ixfer_ait2acc
         926             :                ixfer_acc2ait_sv(i,k) = ixfer_acc2ait
         927             :
         928             :                !
         929             :                ! compute new dgncur & v2ncur for aitken & accum modes
         930             :                !
         931             :                ! currently inactive
         932             :                do n = nait, nacc, (nacc-nait)
         933             :                   if (n .eq. nait) then
         934             :                      duma = (xfertend_num(1,1) - xfertend_num(2,1))*deltat
         935             :                      num_a     = max( 0.0_r8, num_a_aitsv(i,k) - duma )
         936             :                      num_a_acc = max( 0.0_r8, num_a_accsv(i,k) + duma )
         937             :                      duma = (drv_a_aitsv(i,k)*xfercoef_vol_ait2acc -   &
         938             :                         (drv_a_accsv(i,k)-drv_a_noxf)*xfercoef_vol_acc2ait)*deltat
         939             :                      drv_a     = max( 0.0_r8, drv_a_aitsv(i,k) - duma )
         940             :                      drv_a_acc = max( 0.0_r8, drv_a_accsv(i,k) + duma )
         941             :                      duma = (xfertend_num(1,2) - xfertend_num(2,2))*deltat
         942             :                      num_c     = max( 0.0_r8, num_c_aitsv(i,k) - duma )
         943             :                      num_c_acc = max( 0.0_r8, num_c_accsv(i,k) + duma )
         944             :                      duma = (drv_c_aitsv(i,k)*xfercoef_vol_ait2acc -   &
         945             :                         (drv_c_accsv(i,k)-drv_c_noxf)*xfercoef_vol_acc2ait)*deltat
         946             :                      drv_c     = max( 0.0_r8, drv_c_aitsv(i,k) - duma )
         947             :                      drv_c_acc = max( 0.0_r8, drv_c_accsv(i,k) + duma )
         948             :                   else
         949             :                      num_a = num_a_acc
         950             :                      drv_a = drv_a_acc
         951             :                      num_c = num_c_acc
         952             :                      drv_c = drv_c_acc
         953             :                   end if
         954             :
         955             :                   if (drv_a > 0.0_r8) then
         956             :                      if (num_a <= drv_a*voltonumbhi_amode(n)) then
         957             :                         dgncur_a(i,k,n) = dgnumhi_amode(n)
         958             :                         v2ncur_a(i,k,n) = voltonumbhi_amode(n)
         959             :                      else if (num_a >= drv_a*voltonumblo_amode(n)) then
         960             :                         dgncur_a(i,k,n) = dgnumlo_amode(n)
         961             :                         v2ncur_a(i,k,n) = voltonumblo_amode(n)
         962             :                      else
         963             :                         dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third
         964             :                         v2ncur_a(i,k,n) = num_a/drv_a
         965             :                      end if
         966             :                   else
         967             :                      dgncur_a(i,k,n) = dgnum_amode(n)
         968             :                      v2ncur_a(i,k,n) = voltonumb_amode(n)
         969             :                   end if
         970             :
         971             :                   if (drv_c > 0.0_r8) then
         972             :                      if (num_c <= drv_c*voltonumbhi_amode(n)) then
         973             :                         dgncur_c(i,k,n) = dgnumhi_amode(n)
         974             :                         v2ncur_c(i,k,n) = voltonumbhi_amode(n)
         975             :                      else if (num_c >= drv_c*voltonumblo_amode(n)) then
         976             :                         dgncur_c(i,k,n) = dgnumlo_amode(n)
         977             :                         v2ncur_c(i,k,n) = voltonumblo_amode(n)
         978             :                      else
         979             :                         dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third
         980             :                         v2ncur_c(i,k,n) = num_c/drv_c
         981             :                      end if
         982             :                   else
         983             :                      dgncur_c(i,k,n) = dgnum_amode(n)
         984             :                      v2ncur_c(i,k,n) = voltonumb_amode(n)
         985             :                   end if
         986             :
         987             :                end do
         988             :
         989             :
         990             :                !
         991             :                ! compute tendency amounts for aitken <--> accum transfer
         992             :                !
         993             :
         994             :                if ( masterproc ) then
         995             :                   if (idiagaa > 0) then
         996             :                      do j = 1, 2
         997             :                         do iq = 1, nspecfrm_renamexf(ipair)
         998             :                            do jac = 1, 2
         999             :                               if (j .eq. 1) then
        1000             :                                  if (jac .eq. 1) then
        1001             :                                     lsfrm = lspecfrma_renamexf(iq,ipair)
        1002             :                                     lstoo = lspectooa_renamexf(iq,ipair)
        1003             :                                  else
        1004             :                                     lsfrm = lspecfrmc_renamexf(iq,ipair)
        1005             :                                     lstoo = lspectooc_renamexf(iq,ipair)
        1006             :                                  end if
        1007             :                               else
        1008             :                                  if (jac .eq. 1) then
        1009             :                                     lsfrm = lspectooa_renamexf(iq,ipair)
        1010             :                                     lstoo = lspecfrma_renamexf(iq,ipair)
        1011             :                                  else
        1012             :                                     lsfrm = lspectooc_renamexf(iq,ipair)
        1013             :                                     lstoo = lspecfrmc_renamexf(iq,ipair)
        1014             :                                  end if
        1015             :                               end if
        1016             :                               write( 6, '(a,3i3,2i4)' ) 'calcsize j,iq,jac, lsfrm,lstoo',   &
        1017             :                                  j,iq,jac, lsfrm,lstoo
        1018             :                            end do
        1019             :                         end do
        1020             :                      end do
        1021             :                   end if
        1022             :                end if
        1023             :                idiagaa = -1
        1024             :
        1025             :
        1026             :                ! j=1 does aitken-->accum; j=2 does accum-->aitken
        1027             :                do  j = 1, 2
        1028             :
        1029             :                   if ((j .eq. 1 .and. ixfer_ait2acc > 0) .or. &
        1030             :                      (j .eq. 2 .and. ixfer_acc2ait > 0)) then
        1031             :
        1032             :                      jsrflx = j+2
        1033             :                      if (j .eq. 1) then
        1034             :                         xfercoef = xfercoef_vol_ait2acc
        1035             :                      else
        1036             :                         xfercoef = xfercoef_vol_acc2ait
        1037             :                      end if
        1038             :
        1039             :                      do  iq = 1, nspecfrm_renamexf(ipair)
        1040             :
        1041             :                         ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
        1042             :                         do  jac = 1, 2
        1043             :
        1044             :                            ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species
        1045             :                            ! the lspectooa_renamexf (and lspectooc_renamexf) are accum  species
        1046             :                            ! for j=1, want lsfrm=aitken species, lstoo=accum  species
        1047             :                            ! for j=2, want lsfrm=accum  species,  lstoo=aitken species
        1048             :                            if (j .eq. 1) then
        1049             :                               if (jac .eq. 1) then
        1050             :                                  lsfrm = lspecfrma_renamexf(iq,ipair)
        1051             :                                  lstoo = lspectooa_renamexf(iq,ipair)
        1052             :                               else
        1053             :                                  lsfrm = lspecfrmc_renamexf(iq,ipair)
        1054             :                                  lstoo = lspectooc_renamexf(iq,ipair)
        1055             :                               end if
        1056             :                            else
        1057             :                               if (jac .eq. 1) then
        1058             :                                  lsfrm = lspectooa_renamexf(iq,ipair)
        1059             :                                  lstoo = lspecfrma_renamexf(iq,ipair)
        1060             :                               else
        1061             :                                  lsfrm = lspectooc_renamexf(iq,ipair)
        1062             :                                  lstoo = lspecfrmc_renamexf(iq,ipair)
        1063             :                               end if
        1064             :                            end if
        1065             :
        1066             :                            if ((lsfrm > 0) .and. (lstoo > 0)) then
        1067             :                               if (jac .eq. 1) then
        1068             :                                  if (iq .eq. 1) then
        1069             :                                     xfertend = xfertend_num(j,jac)
        1070             :                                  else
        1071             :                                     xfertend = max(0.0_r8,q(i,k,lsfrm))*xfercoef
        1072             :                                  end if
        1073             :                                  dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xfertend
        1074             :                                  dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xfertend
        1075             :                               else
        1076             :                                  if (iq .eq. 1) then
        1077             :                                     xfertend = xfertend_num(j,jac)
        1078             :                                  else
        1079             :                                     fldcw => qqcw_get_field(pbuf,lsfrm,lchnk)
        1080             :                                     xfertend = max(0.0_r8,fldcw(i,k))*xfercoef
        1081             :                                  end if
        1082             :                                  dqqcwdt(i,k,lsfrm) = dqqcwdt(i,k,lsfrm) - xfertend
        1083             :                                  dqqcwdt(i,k,lstoo) = dqqcwdt(i,k,lstoo) + xfertend
        1084             :                               end if
        1085             :                               qsrflx(i,lsfrm,jsrflx,jac) = qsrflx(i,lsfrm,jsrflx,jac) - xfertend*pdel_fac
        1086             :                               qsrflx(i,lstoo,jsrflx,jac) = qsrflx(i,lstoo,jsrflx,jac) + xfertend*pdel_fac
        1087             :                            end if
        1088             :
        1089             :                         end do
        1090             :                      end do
        1091             :                   end if
        1092             :                end do
        1093             :
        1094             :             end if
        1095             :          end do
        1096             :       end do
        1097             :
        1098             :
        1099             :    end if  !  do_aitacc_transfer
        1100             :    lsfrm = -123456789   ! executable statement for debugging
        1101             :
        1102             :
        1103             :    !
        1104             :    ! apply tendencies to cloud-borne species MRs
        1105             :    !
        1106             :    do l = 1, pcnst
        1107             :       lc = l
        1108             :       if ( lc>0 .and. dotendqqcw(lc) ) then
        1109             :          fldcw=> qqcw_get_field(pbuf,l,lchnk)
        1110             :          do k = top_lev, pver
        1111             :             do i = 1, ncol
        1112             :                fldcw(i,k) = max( 0.0_r8,   &
        1113             :                   (fldcw(i,k) + dqqcwdt(i,k,lc)*deltat) )
        1114             :             end do
        1115             :          end do
        1116             :       end if
        1117             :    end do
        1118             :
        1119             :    !
        1120             :    ! do outfld calls
        1121             :    !
        1122             :
        1123             :    ! history fields for number-adjust source-sink for all modes
        1124             :    if ( .not. do_adjust ) return
        1125             :
        1126             :    do n = 1, ntot_amode
        1127             :       if (mprognum_amode(n) <= 0) cycle
        1128             :
        1129             :       do jac = 1, 2
        1130             :          if (jac == 1) then
        1131             :             l = numptr_amode(n)
        1132             :             tmpnamea = cnst_name(l)
        1133             :          else
        1134             :             l = numptrcw_amode(n)
        1135             :             tmpnamea = cnst_name_cw(l)
        1136             :          end if
        1137             :          fieldname = trim(tmpnamea) // '_sfcsiz1'
        1138             :          call outfld( fieldname, qsrflx(:,l,1,jac), pcols, lchnk)
        1139             :
        1140             :          fieldname = trim(tmpnamea) // '_sfcsiz2'
        1141             :          call outfld( fieldname, qsrflx(:,l,2,jac), pcols, lchnk)
        1142             :       end do   ! jac = ...
        1143             :
        1144             :    end do   ! n = ...
        1145             :
        1146             :
        1147             :    ! history fields for aitken-accum transfer
        1148             :    if ( .not. do_aitacc_transfer ) return
        1149             :
        1150             :    do iq = 1, nspecfrm_renamexf(ipair)
        1151             :
        1152             :       ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
        1153             :       do jac = 1, 2
        1154             :
        1155             :          ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species
        1156             :          ! the lspectooa_renamexf (and lspectooc_renamexf) are accum  species
        1157             :          if (jac .eq. 1) then
        1158             :             lsfrm = lspecfrma_renamexf(iq,ipair)
        1159             :             lstoo = lspectooa_renamexf(iq,ipair)
        1160             :          else
        1161             :             lsfrm = lspecfrmc_renamexf(iq,ipair)
        1162             :             lstoo = lspectooc_renamexf(iq,ipair)
        1163             :          end if
        1164             :          if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
        1165             :
        1166             :          if (jac .eq. 1) then
        1167             :             tmpnamea = cnst_name(lsfrm)
        1168             :             tmpnameb = cnst_name(lstoo)
        1169             :          else
        1170             :             tmpnamea = cnst_name_cw(lsfrm)
        1171             :             tmpnameb = cnst_name_cw(lstoo)
        1172             :          end if
        1173             :          if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
        1174             :
        1175             :          fieldname = trim(tmpnamea) // '_sfcsiz3'
        1176             :          call outfld( fieldname, qsrflx(:,lsfrm,3,jac), pcols, lchnk)
        1177             :
        1178             :          fieldname = trim(tmpnameb) // '_sfcsiz3'
        1179             :          call outfld( fieldname, qsrflx(:,lstoo,3,jac), pcols, lchnk)
        1180             :
        1181             :          fieldname = trim(tmpnamea) // '_sfcsiz4'
        1182             :          call outfld( fieldname, qsrflx(:,lsfrm,4,jac), pcols, lchnk)
        1183             :
        1184             :          fieldname = trim(tmpnameb) // '_sfcsiz4'
        1185             :          call outfld( fieldname, qsrflx(:,lstoo,4,jac), pcols, lchnk)
        1186             :
        1187             :       end do   ! jac = ...
        1188             :    end do   ! iq = ...
        1189             :
        1190             :    call modal_aero_calcdry(state, pbuf)
        1191             :
        1192             : #endif
        1193             :
        1194           0 : end subroutine modal_aero_calcsize_sub
        1195             :
        1196             :
        1197             : !----------------------------------------------------------------------
        1198             :
        1199             :
        1200           0 : subroutine modal_aero_calcsize_diag(state, pbuf, list_idx_in, dgnum_m, &
        1201             :                                     hygro_m, dryvol_m, dryrad_m, drymass_m, so4dryvol_m, naer_m)
        1202             :
        1203             :    !-----------------------------------------------------------------------
        1204             :    !
        1205             :    ! Calculate aerosol size distribution parameters
        1206             :    !
        1207             :    ! ***N.B.*** DGNUM for the modes in the climate list are put directly into
        1208             :    !            the physics buffer.  For diagnostic list calculations use the
        1209             :    !            optional list_idx and dgnum args.
        1210             :    !-----------------------------------------------------------------------
        1211             :
        1212             :    ! arguments
        1213             :    type(physics_state), intent(in), target :: state   ! Physics state variables
        1214             :    type(physics_buffer_desc), pointer :: pbuf(:)      ! physics buffer
        1215             :
        1216             :    integer,  optional, intent(in)   :: list_idx_in    ! diagnostic list index
        1217             :    real(r8), optional, pointer      :: dgnum_m(:,:,:) ! interstital aerosol dry number mode radius (m)
        1218             :    real(r8), optional, pointer      :: hygro_m(:,:,:)
        1219             :    real(r8), optional, pointer      :: dryvol_m(:,:,:)
        1220             :    real(r8), optional, pointer      :: dryrad_m(:,:,:)
        1221             :    real(r8), optional, pointer      :: drymass_m(:,:,:)
        1222             :    real(r8), optional, pointer      :: so4dryvol_m(:,:,:)
        1223             :    real(r8), optional, pointer      :: naer_m(:,:,:)
        1224             :
        1225             :
        1226             :    ! local
        1227             :    integer  :: i, k, l1, n
        1228             :    integer  :: lchnk, ncol
        1229             :    integer  :: list_idx, stat
        1230             :    integer  :: nmodes
        1231             :    integer  :: nspec
        1232             :
        1233           0 :    real(r8), pointer :: dgncur_a(:,:) ! (pcols,pver)
        1234             :
        1235             :
        1236             :    real(r8), parameter :: third = 1.0_r8/3.0_r8
        1237             :
        1238           0 :    real(r8), pointer :: mode_num(:,:) ! mode number mixing ratio
        1239           0 :    real(r8), pointer :: specmmr(:,:)  ! specie mmr
        1240             :    real(r8)          :: specdens      ! specie density
        1241             :
        1242             :    real(r8) :: dryvol_a(pcols,pver)   ! interstital aerosol dry volume (cm^3/mol_air)
        1243             :
        1244             :    real(r8) :: dgnum, dgnumhi, dgnumlo
        1245             :    real(r8) :: dgnyy, dgnxx           ! dgnumlo/hi of current mode
        1246             :    real(r8) :: drv_a                  ! dry volume (cm3/mol_air)
        1247             :    real(r8) :: dumfac, dummwdens      ! work variables
        1248             :    real(r8) :: num_a0                 ! initial number (#/mol_air)
        1249             :    real(r8) :: num_a                  ! final number (#/mol_air)
        1250             :    real(r8) :: voltonumbhi, voltonumblo
        1251             :    real(r8) :: v2nyy, v2nxx           ! voltonumblo/hi of current mode
        1252             :    real(r8) :: sigmag, alnsg
        1253             :    !-----------------------------------------------------------------------
        1254             :
        1255           0 :    lchnk = state%lchnk
        1256           0 :    ncol  = state%ncol
        1257             :
        1258           0 :    list_idx = 0  ! climate list by default
        1259           0 :    if (present(list_idx_in)) list_idx = list_idx_in
        1260             :
        1261           0 :    call rad_cnst_get_info(list_idx, nmodes=nmodes)
        1262             :
        1263           0 :    if (list_idx /= 0) then
        1264           0 :       if (.not. present(dgnum_m)) then
        1265             :          call endrun('modal_aero_calcsize_diag called for'// &
        1266           0 :                      'diagnostic list but dgnum_m pointer not present')
        1267             :       end if
        1268           0 :       if (.not. associated(dgnum_m)) then
        1269             :          call endrun('modal_aero_calcsize_diag called for'// &
        1270           0 :             'diagnostic list but dgnum_m not associated')
        1271             :       end if
        1272             :
        1273           0 :       if (.not. present(hygro_m)) then
        1274             :          call endrun('modal_aero_calcsize_diag called for'// &
        1275           0 :                      'diagnostic list but hygro_m pointer not present')
        1276             :       end if
        1277           0 :       if (.not. associated(hygro_m)) then
        1278             :          call endrun('modal_aero_calcsize_diag called for'// &
        1279           0 :             'diagnostic list but hygro_m not associated')
        1280             :       end if
        1281             :
        1282           0 :       if (.not. present(dryvol_m)) then
        1283             :          call endrun('modal_aero_calcsize_diag called for'// &
        1284           0 :                      'diagnostic list but dryvol_m pointer not present')
        1285             :       end if
        1286           0 :       if (.not. associated(dryvol_m)) then
        1287             :          call endrun('modal_aero_calcsize_diag called for'// &
        1288           0 :             'diagnostic list but dryvol_m not associated')
        1289             :       end if
        1290             :
        1291           0 :       if (.not. present(dryrad_m)) then
        1292             :          call endrun('modal_aero_calcsize_diag called for'// &
        1293           0 :                      'diagnostic list but dryrad_m pointer not present')
        1294             :       end if
        1295           0 :       if (.not. associated(dryrad_m)) then
        1296             :          call endrun('modal_aero_calcsize_diag called for'// &
        1297           0 :             'diagnostic list but dryrad_m not associated')
        1298             :       end if
        1299             :
        1300           0 :       if (.not. present(drymass_m)) then
        1301             :          call endrun('modal_aero_calcsize_diag called for'// &
        1302           0 :                      'diagnostic list but drymass_m pointer not present')
        1303             :       end if
        1304           0 :       if (.not. associated(drymass_m)) then
        1305             :          call endrun('modal_aero_calcsize_diag called for'// &
        1306           0 :             'diagnostic list but drymass_m not associated')
        1307             :       end if
        1308             :
        1309           0 :       if (.not. present(so4dryvol_m)) then
        1310             :          call endrun('modal_aero_calcsize_diag called for'// &
        1311           0 :                      'diagnostic list but so4dryvol_m pointer not present')
        1312             :       end if
        1313           0 :       if (.not. associated(so4dryvol_m)) then
        1314             :          call endrun('modal_aero_calcsize_diag called for'// &
        1315           0 :             'diagnostic list but so4dryvol_m not associated')
        1316             :       end if
        1317             :
        1318           0 :       if (.not. present(naer_m)) then
        1319             :          call endrun('modal_aero_calcsize_diag called for'// &
        1320           0 :                      'diagnostic list but naer_m pointer not present')
        1321             :       end if
        1322           0 :       if (.not. associated(naer_m)) then
        1323             :          call endrun('modal_aero_calcsize_diag called for'// &
        1324           0 :             'diagnostic list but naer_m not associated')
        1325             :       end if
        1326             :
        1327             :    end if
        1328             :
        1329           0 :    do n = 1, nmodes
        1330             :
        1331           0 :       if (list_idx == 0) then
        1332           0 :          call pbuf_get_field(pbuf, dgnum_idx, dgncur_a, start=(/1,1,n/), kount=(/pcols,pver,1/))
        1333             :       else
        1334           0 :          dgncur_a => dgnum_m(:,:,n)
        1335             :       end if
        1336             :
        1337             :       ! get mode properties
        1338             :       call rad_cnst_get_mode_props(list_idx, n, dgnum=dgnum, dgnumhi=dgnumhi, dgnumlo=dgnumlo, &
        1339           0 :                                    sigmag=sigmag)
        1340             :
        1341             :       ! get mode number mixing ratio
        1342           0 :       call rad_cnst_get_mode_num(list_idx, n, 'a', state, pbuf, mode_num)
        1343             :
        1344           0 :       dgncur_a(:,:) = dgnum
        1345           0 :       dryvol_a(:,:) = 0.0_r8
        1346             :
        1347             :       ! compute dry volume mixrats =
        1348             :       !      sum_over_components{ component_mass mixrat / density }
        1349           0 :       call rad_cnst_get_info(list_idx, n, nspec=nspec)
        1350           0 :       do l1 = 1, nspec
        1351             :
        1352           0 :          call rad_cnst_get_aer_mmr(list_idx, n, l1, 'a', state, pbuf, specmmr)
        1353           0 :          call rad_cnst_get_aer_props(list_idx, n, l1, density_aer=specdens)
        1354             :
        1355             :          ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
        1356           0 :          dummwdens = 1.0_r8 / specdens
        1357             :
        1358           0 :          do k=top_lev,pver
        1359           0 :             do i=1,ncol
        1360           0 :                dryvol_a(i,k) = dryvol_a(i,k)    &
        1361           0 :                   + max(0.0_r8, specmmr(i,k))*dummwdens
        1362             :             end do
        1363             :          end do
        1364             :       end do
        1365             :
        1366           0 :       alnsg  = log( sigmag )
        1367           0 :       dumfac = exp(4.5_r8*alnsg**2)*pi/6.0_r8
        1368           0 :       voltonumblo = 1._r8 / ( (pi/6._r8)*(dgnumlo**3)*exp(4.5_r8*alnsg**2) )
        1369           0 :       voltonumbhi = 1._r8 / ( (pi/6._r8)*(dgnumhi**3)*exp(4.5_r8*alnsg**2) )
        1370           0 :       v2nxx = voltonumbhi
        1371           0 :       v2nyy = voltonumblo
        1372           0 :       dgnxx = dgnumhi
        1373           0 :       dgnyy = dgnumlo
        1374             :
        1375           0 :       do k = top_lev, pver
        1376           0 :          do i = 1, ncol
        1377             :
        1378           0 :             drv_a = dryvol_a(i,k)
        1379           0 :             num_a0 = mode_num(i,k)
        1380           0 :             num_a = max( 0.0_r8, num_a0 )
        1381             :
        1382           0 :             if (drv_a > 0.0_r8) then
        1383           0 :                if (num_a <= drv_a*v2nxx) then
        1384           0 :                   dgncur_a(i,k) = dgnxx
        1385           0 :                else if (num_a >= drv_a*v2nyy) then
        1386           0 :                   dgncur_a(i,k) = dgnyy
        1387             :                else
        1388           0 :                   dgncur_a(i,k) = (drv_a/(dumfac*num_a))**third
        1389             :                end if
        1390             :             end if
        1391             :
        1392             :          end do
        1393             :       end do
        1394             :
        1395             :    end do ! nmodes
        1396             :
        1397           0 :    call modal_aero_calcdry(state, pbuf, list_idx_in, dgnum_m, hygro_m, dryvol_m, dryrad_m, drymass_m, so4dryvol_m, naer_m)
        1398             :
        1399           0 : end subroutine modal_aero_calcsize_diag
        1400             :
        1401           0 : subroutine modal_aero_calcdry(state, pbuf, list_idx_in, dgnumdry_m, hygro_m, dryvol_m, dryrad_m, drymass_m, so4dryvol_m, naer_m)
        1402             :
        1403             :    type(physics_state), target, intent(in)    :: state       ! Physics state variables
        1404             :    type(physics_buffer_desc),   pointer       :: pbuf(:)     ! physics buffer
        1405             :    integer,  optional, intent(in)             :: list_idx_in ! diagnostic list index
        1406             :    real(r8), optional,          pointer       :: dgnumdry_m(:,:,:)
        1407             :    real(r8), optional,          pointer       :: hygro_m(:,:,:)
        1408             :    real(r8), optional,          pointer       :: dryvol_m(:,:,:)
        1409             :    real(r8), optional,          pointer       :: dryrad_m(:,:,:)
        1410             :    real(r8), optional,          pointer       :: drymass_m(:,:,:)
        1411             :    real(r8), optional,          pointer       :: so4dryvol_m(:,:,:)
        1412             :    real(r8), optional,          pointer       :: naer_m(:,:,:)
        1413             :
        1414             :    real(r8), parameter :: third = 1._r8/3._r8
        1415             :    real(r8), parameter :: pi43  = pi*4.0_r8/3.0_r8
        1416             :
        1417           0 :    real(r8), pointer :: maer(:,:)        ! aerosol wet mass MR (including water) (kg/kg-air)
        1418           0 :    real(r8), pointer :: hygro(:,:,:)     ! volume-weighted mean hygroscopicity (--)
        1419           0 :    real(r8), pointer :: dryvol(:,:,:)    ! single-particle-mean dry volume (m3)
        1420           0 :    real(r8), pointer :: dryrad(:,:,:)    ! dry volume mean radius of aerosol (m)
        1421           0 :    real(r8), pointer :: drymass(:,:,:)   ! single-particle-mean dry mass  (kg)
        1422           0 :    real(r8), pointer :: so4dryvol(:,:,:) ! single-particle-mean so4 dry volume (m3)
        1423           0 :    real(r8), pointer :: naer(:,:,:)      ! aerosol number MR (bounded!) (#/kg-air)
        1424             :
        1425           0 :    real(r8), pointer :: dgncur_a(:,:,:)
        1426           0 :    real(r8), pointer :: raer(:,:)   ! aerosol species MRs (kg/kg and #/kg)
        1427             :
        1428             :    real(r8), pointer :: sulfeq(:,:,:) ! H2SO4 equilibrium mixing ratios over particles (mol/mol)
        1429             :
        1430             :    real(r8) :: dryvolmr(pcols,pver)          ! volume MR for aerosol mode (m3/kg)
        1431             :    real(r8) :: so4dryvolmr(pcols,pver)       ! volume MR for sulfate aerosol in mode (m3/kg)
        1432             :
        1433             :    real(r8) :: specdens
        1434             :    real(r8) :: spechygro, spechygro_1
        1435             :    real(r8) :: sigmag
        1436             :    real(r8) :: duma, dumb
        1437             :    real(r8) :: alnsg
        1438             :
        1439             :    real(r8) :: v2ncur_a
        1440             :    real(r8) :: drydens               ! dry particle density  (kg/m^3)
        1441             :
        1442             :    character(len=fieldname_len+3) :: fieldname
        1443             :    character(len=32) :: spectype
        1444             :
        1445             :    integer :: nmodes, lchnk, ncol, list_idx, i, k, l, m
        1446             :    integer :: nspec
        1447             :
        1448             :
        1449           0 :    lchnk = state%lchnk
        1450           0 :    ncol = state%ncol
        1451             :
        1452           0 :    list_idx = 0
        1453           0 :    if (present(list_idx_in)) then
        1454           0 :       list_idx = list_idx_in
        1455             :
        1456             :       ! check that all optional args are present
        1457           0 :       if (.not. present(dgnumdry_m)) then
        1458             :          call endrun('modal_aero_calcdry called for'// &
        1459           0 :                      'diagnostic list but required args not present')
        1460             :       end if
        1461             :
        1462             :       ! arrays for diagnostic calculations must be associated
        1463           0 :      if (.not. associated(dgnumdry_m)) then
        1464             :          call endrun('modal_aero_calcdry called for'// &
        1465           0 :                      'diagnostic list but required args not associated')
        1466             :       end if
        1467             :    end if
        1468             :
        1469             :    ! loop over all aerosol modes
        1470           0 :    call rad_cnst_get_info(list_idx, nmodes=nmodes)
        1471             :
        1472           0 :    allocate( maer(pcols,pver))
        1473             :
        1474           0 :    if (list_idx == 0) then
        1475           0 :       call pbuf_get_field(pbuf, dgnum_idx,     dgncur_a )
        1476           0 :       call pbuf_get_field(pbuf, hygro_idx,     hygro)
        1477           0 :       call pbuf_get_field(pbuf, dryvol_idx,    dryvol)
        1478           0 :       call pbuf_get_field(pbuf, dryrad_idx,    dryrad)
        1479           0 :       call pbuf_get_field(pbuf, drymass_idx,   drymass)
        1480           0 :       call pbuf_get_field(pbuf, so4dryvol_idx, so4dryvol)
        1481           0 :       call pbuf_get_field(pbuf, naer_idx,      naer)
        1482             :    else
        1483           0 :       dgncur_a    => dgnumdry_m
        1484           0 :       hygro       => hygro_m
        1485           0 :       dryvol      => dryvol_m
        1486           0 :       dryrad      => dryrad_m
        1487           0 :       drymass     => drymass_m
        1488           0 :       so4dryvol   => so4dryvol_m
        1489           0 :       naer        => naer_m
        1490             :    end if
        1491             :
        1492           0 :    hygro(:,:,:)     = 0._r8
        1493           0 :    so4dryvol(:,:,:) = 0._r8
        1494             :
        1495           0 :    do m = 1, nmodes
        1496             :
        1497           0 :       maer(:,:)      = 0._r8
        1498           0 :       dryvolmr(:,:) = 0._r8
        1499           0 :       so4dryvolmr(:,:) = 0._r8
        1500             :
        1501             :       ! get mode properties
        1502           0 :       call rad_cnst_get_mode_props(list_idx, m, sigmag=sigmag)
        1503             :
        1504             :       ! get mode info
        1505           0 :       call rad_cnst_get_info(list_idx, m, nspec=nspec)
        1506             :
        1507           0 :       do l = 1, nspec
        1508             :
        1509             :          ! get species interstitial mixing ratio ('a')
        1510           0 :          call rad_cnst_get_aer_mmr(list_idx, m, l, 'a', state, pbuf, raer)
        1511             :          call rad_cnst_get_aer_props(list_idx, m, l, density_aer=specdens, &
        1512           0 :                                      hygro_aer=spechygro, spectype=spectype)
        1513             :
        1514           0 :          if (l == 1) then
        1515             :             ! save off these values to be used as defaults
        1516           0 :             spechygro_1    = spechygro
        1517             :          end if
        1518             :
        1519           0 :          do k = top_lev, pver
        1520           0 :             do i = 1, ncol
        1521           0 :                duma          = raer(i,k)     ! kg/kg air
        1522           0 :                maer(i,k)     = maer(i,k) + duma
        1523           0 :                dumb          = duma/specdens ! m3/kg air
        1524           0 :                dryvolmr(i,k) = dryvolmr(i,k) + dumb
        1525           0 :                if (modal_strat_sulfate .and. (trim(spectype).eq.'sulfate')) then
        1526           0 :                   so4dryvolmr(i,k) = so4dryvolmr(i,k) + dumb
        1527             :                end if
        1528           0 :                hygro(i,k,m)  = hygro(i,k,m) + dumb*spechygro
        1529             :             end do
        1530             :          end do
        1531             :       end do
        1532             :
        1533           0 :       alnsg = log(sigmag)
        1534             :
        1535           0 :       do k = top_lev, pver
        1536           0 :          do i = 1, ncol
        1537             :
        1538           0 :             if (dryvolmr(i,k) > 1.0e-30_r8) then
        1539           0 :                hygro(i,k,m) = hygro(i,k,m)/dryvolmr(i,k)
        1540             :             else
        1541           0 :                hygro(i,k,m) = spechygro_1
        1542             :             end if
        1543             :
        1544             :             ! dry aerosol properties
        1545             :
        1546           0 :             v2ncur_a = 1._r8 / ( (pi/6._r8)*(dgncur_a(i,k,m)**3._r8)*exp(4.5_r8*alnsg**2._r8) )
        1547             :             ! naer = aerosol number (#/kg)
        1548           0 :             naer(i,k,m) = dryvolmr(i,k)*v2ncur_a
        1549             :
        1550             :             ! compute mean (1 particle) dry volume and mass for each mode
        1551             :             ! old coding is replaced because the new (1/v2ncur_a) is equal to
        1552             :             ! the mean particle volume
        1553             :             ! also moletomass forces maer >= 1.0e-30, so (maer/dryvolmr)
        1554             :             ! should never cause problems (but check for maer < 1.0e-31 anyway)
        1555           0 :             if (maer(i,k) .gt. 1.0e-31_r8) then
        1556           0 :                drydens = maer(i,k)/dryvolmr(i,k)        ! kg/m3 aerosol
        1557             :             else
        1558             :                drydens = 1.0_r8
        1559             :             end if
        1560           0 :             dryvol(i,k,m)   = 1.0_r8/v2ncur_a             ! m3/particle
        1561           0 :             drymass(i,k,m)  = drydens*dryvol(i,k,m)       ! kg/particle
        1562           0 :             dryrad(i,k,m)   = (dryvol(i,k,m)/pi43)**third ! m
        1563             :          end do    ! i = 1, ncol
        1564             :       end do    ! k = top_lev, pver
        1565             :
        1566             :
        1567           0 :       if (modal_strat_sulfate) then
        1568           0 :          do k = top_lev, pver
        1569           0 :             do i = 1, ncol
        1570           0 :                if (so4dryvolmr(i,k) .gt. 1.0e-31_r8) then
        1571           0 :                   so4dryvol(i,k,m) = dryvol(i,k,m)*so4dryvolmr(i,k)/dryvolmr(i,k)
        1572             :                else
        1573           0 :                   so4dryvol(i,k,m) = 0.0_r8
        1574             :                end if
        1575             :
        1576             :             end do    ! i = 1, ncol
        1577             :          end do    ! k = top_lev, pver
        1578             :
        1579             :       end if
        1580             :
        1581             :    end do    ! m = 1, nmodes
        1582             :
        1583           0 :    deallocate( maer)
        1584             :
        1585             :
        1586           0 : end subroutine modal_aero_calcdry
        1587             : !----------------------------------------------------------------------
        1588             :
        1589             : end module modal_aero_calcsize

Generated by: LCOV version 1.14