LCOV - code coverage report
Current view: | top level - chemistry/utils - modal_aero_calcsize.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 578 | 672 | 86.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 5 | 5 | 100.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        1536 : 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        1536 :   call rad_cnst_get_info(0, nmodes=nmodes)
          80             :
          81        6144 :   call pbuf_add_field('DGNUM', 'global',  dtype_r8, (/pcols, pver, nmodes/), dgnum_idx)
          82             :
          83        6144 :   call pbuf_add_field('HYGRO',     'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), hygro_idx)
          84        6144 :   call pbuf_add_field('DRYVOL',    'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), dryvol_idx)
          85        6144 :   call pbuf_add_field('DRYRAD',    'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), dryrad_idx)
          86        6144 :   call pbuf_add_field('DRYMASS',   'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), drymass_idx)
          87        6144 :   call pbuf_add_field('SO4DRYVOL', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), so4dryvol_idx)
          88        6144 :   call pbuf_add_field('NAER',      'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), naer_idx)
          89             :
          90        1536 : end subroutine modal_aero_calcsize_reg
          91             :
          92             : !===============================================================================
          93             : !===============================================================================
          94             :
          95        1536 : subroutine modal_aero_calcsize_init(pbuf2d)
          96        1536 :    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        1536 :    call phys_getopts(history_aerosol_out=history_aerosol)
         126             :
         127             :    ! init entities required for both prescribed and prognostic modes
         128             :
         129        1536 :    if (is_first_step()) then
         130             :       ! initialize fields in physics buffer
         131         768 :       call pbuf_set_field(pbuf2d, dgnum_idx, 0.0_r8)
         132             :    endif
         133             :
         134             : #ifndef MODAL_AERO
         135             :    do_adjust_default          = .false.
         136             :    do_aitacc_transfer_default = .false.
         137             : #else
         138             :    !  do_adjust_default allows adjustment to be turned on/off
         139        1536 :    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        1536 :    nait = modeptr_aitken
         145        1536 :    nacc = modeptr_accum
         146        1536 :    do_aitacc_transfer_default = .false.
         147             :    if ((modeptr_aitken > 0) .and.   &
         148        1536 :       (modeptr_accum  > 0) .and.   &
         149             :       (modeptr_aitken /= modeptr_accum)) then
         150        1536 :       do_aitacc_transfer_default = .true.
         151        1536 :       if (mprognum_amode(nait) <= 0) do_aitacc_transfer_default = .false.
         152        1536 :       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        7680 :    do n = 1, ntot_amode
         159        6144 :       if (mprognum_amode(n) <= 0) cycle
         160             :
         161       19968 :       do jac = 1, 2
         162       12288 :          if (jac == 1) then
         163        6144 :             tmpnamea = cnst_name(numptr_amode(n))
         164             :          else
         165        6144 :             tmpnamea = cnst_name_cw(numptrcw_amode(n))
         166             :          end if
         167       12288 :          unit = '#/m2/s'
         168       12288 :          fieldname = trim(tmpnamea) // '_sfcsiz1'
         169       12288 :          long_name = trim(tmpnamea) // ' calcsize number-adjust column source'
         170       12288 :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         171       12288 :          if (history_aerosol) then
         172           0 :             call add_default(fieldname, 1, ' ')
         173             :          end if
         174       12288 :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         175             :
         176       12288 :          fieldname = trim(tmpnamea) // '_sfcsiz2'
         177       12288 :          long_name = trim(tmpnamea) // ' calcsize number-adjust column sink'
         178       12288 :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         179       12288 :          if (history_aerosol) then
         180           0 :             call add_default(fieldname, 1, ' ')
         181             :          end if
         182       18432 :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         183             :       end do   ! jac = ...
         184             :    end do   ! n = ...
         185             :
         186        1536 :    if ( .not. do_aitacc_transfer_default ) return
         187             :
         188             :    ! check that renaming ipair=1 is aitken-->accum
         189        1536 :    ipair = 1
         190        1536 :    if ((modefrm_renamexf(ipair) .ne. nait) .or.   &
         191             :       (modetoo_renamexf(ipair) .ne. nacc)) then
         192             :       write( 6, '(//2a//)' )   &
         193           0 :          '*** modal_aero_calcaersize_init error -- ',   &
         194           0 :          'modefrm/too_renamexf(1) are wrong'
         195           0 :       call endrun( 'modal_aero_calcaersize_init error' )
         196             :    end if
         197             :
         198             :    ! define history fields for aitken-accum transfer
         199        9216 :    do iq = 1, nspecfrm_renamexf(ipair)
         200             :
         201             :       ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
         202       24576 :       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       15360 :          if (jac .eq. 1) then
         207        7680 :             lsfrm = lspecfrma_renamexf(iq,ipair)
         208        7680 :             lstoo = lspectooa_renamexf(iq,ipair)
         209             :          else
         210        7680 :             lsfrm = lspecfrmc_renamexf(iq,ipair)
         211        7680 :             lstoo = lspectooc_renamexf(iq,ipair)
         212             :          end if
         213       15360 :          if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
         214             :
         215       15360 :          if (jac .eq. 1) then
         216        7680 :             tmpnamea = cnst_name(lsfrm)
         217        7680 :             tmpnameb = cnst_name(lstoo)
         218             :          else
         219        7680 :             tmpnamea = cnst_name_cw(lsfrm)
         220        7680 :             tmpnameb = cnst_name_cw(lstoo)
         221             :          end if
         222             :
         223       15360 :          unit = 'kg/m2/s'
         224       15360 :          if ((tmpnamea(1:3) == 'num') .or. &
         225        3072 :             (tmpnamea(1:3) == 'NUM')) unit = '#/m2/s'
         226       15360 :          fieldname = trim(tmpnamea) // '_sfcsiz3'
         227       15360 :          long_name = trim(tmpnamea) // ' calcsize aitken-to-accum adjust column tendency'
         228       15360 :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         229       15360 :          if (history_aerosol) then
         230           0 :             call add_default(fieldname, 1, ' ')
         231             :          end if
         232       15360 :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         233             :
         234       15360 :          fieldname = trim(tmpnameb) // '_sfcsiz3'
         235       15360 :          long_name = trim(tmpnameb) // ' calcsize aitken-to-accum adjust column tendency'
         236       15360 :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         237       15360 :          if (history_aerosol) then
         238           0 :             call add_default(fieldname, 1, ' ')
         239             :          end if
         240       15360 :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         241             :
         242       15360 :          fieldname = trim(tmpnamea) // '_sfcsiz4'
         243       15360 :          long_name = trim(tmpnamea) // ' calcsize accum-to-aitken adjust column tendency'
         244       15360 :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         245       15360 :          if (history_aerosol) then
         246           0 :             call add_default(fieldname, 1, ' ')
         247             :          end if
         248       15360 :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         249             :
         250       15360 :          fieldname = trim(tmpnameb) // '_sfcsiz4'
         251       15360 :          long_name = trim(tmpnameb) // ' calcsize accum-to-aitken adjust column tendency'
         252       15360 :          call addfld( fieldname, horiz_only, 'A', unit, long_name )
         253       15360 :          if (history_aerosol) then
         254           0 :             call add_default(fieldname, 1, ' ')
         255             :          end if
         256       23040 :          if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
         257             :
         258             :       end do   ! jac = ...
         259             :    end do   ! iq = ...
         260             :
         261             : #endif
         262             :
         263        1536 : end subroutine modal_aero_calcsize_init
         264             :
         265             : !===============================================================================
         266             :
         267     1489176 : 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     1489176 :    real(r8), pointer :: t(:,:)      ! Temperature in Kelvin
         307     1489176 :    real(r8), pointer :: pmid(:,:)   ! pressure at model levels (Pa)
         308     1489176 :    real(r8), pointer :: pdel(:,:)   ! pressure thickness of levels
         309     1489176 :    real(r8), pointer :: q(:,:,:)    ! Tracer MR array
         310             :
         311     1489176 :    logical,  pointer :: dotend(:)   ! flag for doing tendency
         312     1489176 :    real(r8), pointer :: dqdt(:,:,:) ! TMR tendency array
         313             :
         314     1489176 :    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     2978352 :    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     1489176 :    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     2978352 :    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     2978352 :    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     2978352 :    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     2978352 :    real(r8) :: v2ncur_a(pcols,pver,ntot_amode)
         368     1489176 :    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     1489176 :    if (present(do_adjust_in)) then
         390           0 :       do_adjust = do_adjust_in
         391             :    else
         392     1489176 :       do_adjust = do_adjust_default
         393             :    end if
         394             :
         395     1489176 :    if (present(do_aitacc_transfer_in)) then
         396           0 :       do_aitacc_transfer = do_aitacc_transfer_in
         397             :    else
         398     1489176 :       do_aitacc_transfer = do_aitacc_transfer_default
         399             :    end if
         400             :
         401     1489176 :    lchnk = state%lchnk
         402     1489176 :    ncol  = state%ncol
         403             :
         404     1489176 :    t    => state%t
         405     1489176 :    pmid => state%pmid
         406     1489176 :    pdel => state%pdel
         407     1489176 :    q    => state%q
         408             :
         409     1489176 :    dotend => ptend%lq
         410     1489176 :    dqdt   => ptend%q
         411             :
         412     1489176 :    call pbuf_get_field(pbuf, dgnum_idx, dgncur_a)
         413             :
         414     1489176 :    dotendqqcw(:) = .false.
         415     1489176 :    dqqcwdt(:,:,:) = 0.0_r8
         416     1489176 :    qsrflx(:,:,:,:) = 0.0_r8
         417             :
         418     1489176 :    nait = modeptr_aitken
         419     1489176 :    nacc = modeptr_accum
         420             :
         421     1489176 :    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     1489176 :    tadj = 86400
         426     1489176 :    tadj = max( tadj, deltat )
         427     1489176 :    tadjinv = 1.0_r8/(tadj*(1.0_r8 + 1.0e-15_r8))
         428     1489176 :    fracadj = deltat*tadjinv
         429     1489176 :    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     7445880 :    do n = 1, ntot_amode
         441             :
         442             :
         443             :       ! initialize all parameters to the default values for the mode
         444   559930176 :       do k=top_lev,pver
         445  9256025376 :          do i=1,ncol
         446             :             !    sgcur_a(i,k,n) = sigmag_amode(n)
         447             :             !    sgcur_c(i,k,n) = sigmag_amode(n)
         448  8696095200 :             dgncur_a(i,k,n) = dgnum_amode(n)
         449  8696095200 :             dgncur_c(i,k,n) = dgnum_amode(n)
         450  8696095200 :             v2ncur_a(i,k,n) = voltonumb_amode(n)
         451  8696095200 :             v2ncur_c(i,k,n) = voltonumb_amode(n)
         452  8696095200 :             dryvol_a(i,k) = 0.0_r8
         453  9250068672 :             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    28294344 :       do l1 = 1, nspec_amode(n)
         460             :          ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
         461    22337640 :          dummwdens = 1.0_r8 / specdens_amode(l1,n)
         462    22337640 :          la = lmassptr_amode(l1,n)
         463  2099738160 :          do k=top_lev,pver
         464 34710095160 :             do i=1,ncol
         465 65220714000 :                dryvol_a(i,k) = dryvol_a(i,k)    &
         466 99908471520 :                   + max(0.0_r8,q(i,k,la))*dummwdens
         467             :             end do
         468             :          end do
         469             :
         470    22337640 :          fldcw => qqcw_get_field(pbuf,lmassptrcw_amode(l1,n),lchnk)
         471  2105694864 :          do k=top_lev,pver
         472 34710095160 :             do i=1,ncol
         473 65220714000 :                dryvol_c(i,k) = dryvol_c(i,k)    &
         474 99908471520 :                   + 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     5956704 :       lna = numptr_amode(n)
         481     5956704 :       lnc = numptrcw_amode(n)
         482     5956704 :       fldcw => qqcw_get_field(pbuf,numptrcw_amode(n),lchnk,.true.)
         483             :
         484             :
         485             :       ! go to section for appropriate number/surface diagnosed/prognosed options
         486     5956704 :       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           0 :          if (lna > 0) then
         493           0 :             dotend(lna) = .true.
         494           0 :             do k=top_lev,pver
         495           0 :                do i=1,ncol
         496           0 :                   dqdt(i,k,lna) = (dryvol_a(i,k)*voltonumb_amode(n)   &
         497           0 :                      - q(i,k,lna)) * deltatinv
         498             :                end do
         499             :             end do
         500             :          end if
         501           0 :          if (lnc > 0) then
         502           0 :             dotendqqcw(lnc) = .true.
         503           0 :             do k=top_lev,pver
         504           0 :                do i=1,ncol
         505           0 :                   dqqcwdt(i,k,lnc) = (dryvol_c(i,k)*voltonumb_amode(n)   &
         506           0 :                      - 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     5956704 :       frelaxadj = 27.0_r8
         547     5956704 :       dumfac = exp(4.5_r8*alnsg_amode(n)**2)*pi/6.0_r8
         548     5956704 :       v2nxx = voltonumbhi_amode(n)
         549     5956704 :       v2nyy = voltonumblo_amode(n)
         550     5956704 :       v2nxxrl = v2nxx/frelaxadj
         551     5956704 :       v2nyyrl = v2nyy*frelaxadj
         552     5956704 :       dgnxx = dgnumhi_amode(n)
         553     5956704 :       dgnyy = dgnumlo_amode(n)
         554     5956704 :       if ( do_aitacc_transfer ) then
         555     5956704 :          if (n == nait) v2nxx = v2nxx/1.0e6_r8
         556     5956704 :          if (n == nacc) v2nyy = v2nyy*1.0e6_r8
         557     5956704 :          v2nxxrl = v2nxx/frelaxadj   ! NEW
         558     5956704 :          v2nyyrl = v2nyy*frelaxadj   ! NEW
         559             :       end if
         560             :
         561     5956704 :       if (do_adjust) then
         562     5956704 :          dotend(lna) = .true.
         563     5956704 :          dotendqqcw(lnc) = .true.
         564             :       end if
         565             :
         566   561419352 :       do  k = top_lev, pver
         567  9256025376 :          do  i = 1, ncol
         568             :
         569  8696095200 :             drv_a = dryvol_a(i,k)
         570  8696095200 :             num_a0 = q(i,k,lna)
         571  8696095200 :             num_a = max( 0.0_r8, num_a0 )
         572  8696095200 :             drv_c = dryvol_c(i,k)
         573  8696095200 :             num_c0 = fldcw(i,k)
         574  8696095200 :             num_c = max( 0.0_r8, num_c0 )
         575             :
         576  8696095200 :             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  8696095200 :                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           0 :                   num_a = 0.0_r8
         589           0 :                   dqdt(i,k,lna) = -num_a0*deltatinv
         590           0 :                   num_c = 0.0_r8
         591           0 :                   dqqcwdt(i,k,lnc) = -num_c0*deltatinv
         592  8696095200 :                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  7255763769 :                   num_c = 0.0_r8
         596  7255763769 :                   dqqcwdt(i,k,lnc) = -num_c0*deltatinv
         597  7255763769 :                   num_a1 = num_a
         598  7255763769 :                   numbnd = max( drv_a*v2nxx, min( drv_a*v2nyy, num_a1 ) )
         599  7255763769 :                   num_a  = num_a1 + (numbnd - num_a1)*fracadj
         600  7255763769 :                   dqdt(i,k,lna) = (num_a - num_a0)*deltatinv
         601             :
         602  1440331431 :                else if (drv_a <= 0.0_r8) then
         603             :                   ! interstitial volume is zero, treat similar to above
         604           0 :                   num_a = 0.0_r8
         605           0 :                   dqdt(i,k,lna) = -num_a0*deltatinv
         606           0 :                   num_c1 = num_c
         607           0 :                   numbnd = max( drv_c*v2nxx, min( drv_c*v2nyy, num_c1 ) )
         608           0 :                   num_c  = num_c1 + (numbnd - num_c1)*fracadj
         609           0 :                   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  1440331431 :                   num_a1 = num_a
         615  1440331431 :                   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  1440331431 :                   numbnd = max( drv_a*v2nxxrl, min( drv_a*v2nyyrl, num_a1 ) )
         621  1440331431 :                   delnum_a2 = (numbnd - num_a1)*fracadj
         622  1440331431 :                   num_a2 = num_a1 + delnum_a2
         623  1440331431 :                   numbnd = max( drv_c*v2nxxrl, min( drv_c*v2nyyrl, num_c1 ) )
         624  1440331431 :                   delnum_c2 = (numbnd - num_c1)*fracadj
         625  1440331431 :                   num_c2 = num_c1 + delnum_c2
         626  1440331431 :                   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   409792488 :                         num_a1-delnum_c2 ) )
         629  1030538943 :                   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         726 :                         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  1440331431 :                   drv_t = drv_a + drv_c
         636  1440331431 :                   num_t2 = num_a2 + num_c2
         637  1440331431 :                   delnum_a3 = 0.0_r8
         638  1440331431 :                   delnum_c3 = 0.0_r8
         639  1440331431 :                   if (num_t2 < drv_t*v2nxx) then
         640   131346909 :                      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   131346909 :                      if ((num_a2 < drv_a*v2nxx) .and. (num_c2 < drv_c*v2nxx)) then
         644   111641432 :                         delnum_a3 = delnum_t3*(num_a2/num_t2)
         645   111641432 :                         delnum_c3 = delnum_t3*(num_c2/num_t2)
         646    19705477 :                      else if (num_c2 < drv_c*v2nxx) then
         647             :                         delnum_c3 = delnum_t3
         648    19647066 :                      else if (num_a2 < drv_a*v2nxx) then
         649    19647066 :                         delnum_a3 = delnum_t3
         650             :                      end if
         651  1308984522 :                   else if (num_t2 > drv_t*v2nyy) then
         652    51027259 :                      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    51027259 :                      if ((num_a2 > drv_a*v2nyy) .and. (num_c2 > drv_c*v2nyy)) then
         656    34225673 :                         delnum_a3 = delnum_t3*(num_a2/num_t2)
         657    34225673 :                         delnum_c3 = delnum_t3*(num_c2/num_t2)
         658    16801586 :                      else if (num_c2 > drv_c*v2nyy) then
         659             :                         delnum_c3 = delnum_t3
         660    15430766 :                      else if (num_a2 > drv_a*v2nyy) then
         661    15430766 :                         delnum_a3 = delnum_t3
         662             :                      end if
         663             :                   end if
         664  1440331431 :                   num_a = num_a2 + delnum_a3
         665  1440331431 :                   dqdt(i,k,lna) = (num_a - num_a0)*deltatinv
         666  1440331431 :                   num_c = num_c2 + delnum_c3
         667  1440331431 :                   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  8696095200 :             if (drv_a > 0.0_r8) then
         676  8696095200 :                if (num_a <= drv_a*v2nxx) then
         677   965830793 :                   dgncur_a(i,k,n) = dgnxx
         678   965830793 :                   v2ncur_a(i,k,n) = v2nxx
         679  7730264407 :                else if (num_a >= drv_a*v2nyy) then
         680   563150813 :                   dgncur_a(i,k,n) = dgnyy
         681   563150813 :                   v2ncur_a(i,k,n) = v2nyy
         682             :                else
         683  7167113594 :                   dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third
         684  7167113594 :                   v2ncur_a(i,k,n) = num_a/drv_a
         685             :                end if
         686             :             end if
         687  8696095200 :             pdel_fac = pdel(i,k)/gravit   ! = rho*dz
         688  8696095200 :             jac = 1
         689  8696095200 :             qsrflx(i,lna,1,jac) = qsrflx(i,lna,1,jac) + max(0.0_r8,dqdt(i,k,lna))*pdel_fac
         690  8696095200 :             qsrflx(i,lna,2,jac) = qsrflx(i,lna,2,jac) + min(0.0_r8,dqdt(i,k,lna))*pdel_fac
         691             :
         692  8696095200 :             if (drv_c > 0.0_r8) then
         693  1440331431 :                if (num_c <= drv_c*v2nxx) then
         694   520002572 :                   dgncur_c(i,k,n) = dgnumhi_amode(n)
         695   520002572 :                   v2ncur_c(i,k,n) = v2nxx
         696   920328859 :                else if (num_c >= drv_c*v2nyy) then
         697    72182613 :                   dgncur_c(i,k,n) = dgnumlo_amode(n)
         698    72182613 :                   v2ncur_c(i,k,n) = v2nyy
         699             :                else
         700   848146246 :                   dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third
         701   848146246 :                   v2ncur_c(i,k,n) = num_c/drv_c
         702             :                end if
         703             :             end if
         704  8696095200 :             jac = 2
         705  8696095200 :             qsrflx(i,lnc,1,jac) = qsrflx(i,lnc,1,jac) + max(0.0_r8,dqqcwdt(i,k,lnc))*pdel_fac
         706  8696095200 :             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  8696095200 :             if ( do_aitacc_transfer ) then
         711  8696095200 :                if (n == nait) then
         712  2174023800 :                   drv_a_aitsv(i,k) = drv_a
         713  2174023800 :                   num_a_aitsv(i,k) = num_a
         714  2174023800 :                   drv_c_aitsv(i,k) = drv_c
         715  2174023800 :                   num_c_aitsv(i,k) = num_c
         716  6522071400 :                else if (n == nacc) then
         717  2174023800 :                   drv_a_accsv(i,k) = drv_a
         718  2174023800 :                   num_a_accsv(i,k) = num_a
         719  2174023800 :                   drv_c_accsv(i,k) = drv_c
         720  2174023800 :                   num_c_accsv(i,k) = num_c
         721             :                end if
         722             :             end if
         723  8696095200 :             drv_a_sv(i,k,n) = drv_a
         724  8696095200 :             num_a_sv(i,k,n) = num_a
         725  8696095200 :             drv_c_sv(i,k,n) = drv_c
         726  9250068672 :             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     1489176 :    ixfer_ait2acc_sv(:,:) = 0
         753     1489176 :    ixfer_acc2ait_sv(:,:) = 0
         754     1489176 :    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     1489176 :       if (npair_renamexf .le. 0) then
         760           0 :          npair_renamexf = 0
         761             :          !        call modal_aero_rename_init
         762             :          if (npair_renamexf .le. 0) then
         763             :             write( 6, '(//a//)' )   &
         764           0 :                '*** modal_aero_calcaersize_sub error -- npair_renamexf <= 0'
         765           0 :             call endrun( 'modal_aero_calcaersize_sub error' )
         766             :          end if
         767             :       end if
         768             :
         769             :       ! check that renaming ipair=1 is aitken-->accum
         770     1489176 :       ipair = 1
         771     1489176 :       if ((modefrm_renamexf(ipair) .ne. nait) .or.   &
         772             :          (modetoo_renamexf(ipair) .ne. nacc)) then
         773             :          write( 6, '(//2a//)' )   &
         774           0 :             '*** modal_aero_calcaersize_sub error -- ',   &
         775           0 :             'modefrm/too_renamexf(1) are wrong'
         776           0 :          call endrun( 'modal_aero_calcaersize_sub error' )
         777             :       end if
         778             :
         779             :       ! set dotend() for species that will be transferred
         780     8935056 :       do iq = 1, nspecfrm_renamexf(ipair)
         781     7445880 :          lsfrm = lspecfrma_renamexf(iq,ipair)
         782     7445880 :          lstoo = lspectooa_renamexf(iq,ipair)
         783     7445880 :          if ((lsfrm > 0) .and. (lstoo > 0)) then
         784     7445880 :             dotend(lsfrm) = .true.
         785     7445880 :             dotend(lstoo) = .true.
         786             :          end if
         787     7445880 :          lsfrm = lspecfrmc_renamexf(iq,ipair)
         788     7445880 :          lstoo = lspectooc_renamexf(iq,ipair)
         789     8935056 :          if ((lsfrm > 0) .and. (lstoo > 0)) then
         790     7445880 :             dotendqqcw(lsfrm) = .true.
         791     7445880 :             dotendqqcw(lstoo) = .true.
         792             :          end if
         793             :       end do
         794             :
         795             :       ! identify accum species cannot be transferred to aitken mode
         796    10424232 :       noxf_acc2ait(:) = .true.
         797    10424232 :       do l1 = 1, nspec_amode(nacc)
         798     8935056 :          la = lmassptr_amode(l1,nacc)
         799    55099512 :          do iq = 1, nspecfrm_renamexf(ipair)
         800    53610336 :             if (lspectooa_renamexf(iq,ipair) == la) then
         801     5956704 :                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     1489176 :       v2nzz = sqrt(voltonumb_amode(nait)*voltonumb_amode(nacc))
         809             :
         810             :       ! loop over columns and levels
         811   139982544 :       do  k = top_lev, pver
         812  2314006344 :          do  i = 1, ncol
         813             :
         814  2174023800 :             pdel_fac = pdel(i,k)/gravit   ! = rho*dz
         815  2174023800 :             xfertend_num(:,:) = 0.0_r8
         816             :
         817             :             ! compute aitken --> accum transfer rates
         818  2174023800 :             ixfer_ait2acc = 0
         819  2174023800 :             xfercoef_num_ait2acc = 0.0_r8
         820  2174023800 :             xfercoef_vol_ait2acc = 0.0_r8
         821             :
         822  2174023800 :             drv_t = drv_a_aitsv(i,k) + drv_c_aitsv(i,k)
         823  2174023800 :             num_t = num_a_aitsv(i,k) + num_c_aitsv(i,k)
         824  2174023800 :             if (drv_t > 0.0_r8) then
         825  2174023800 :                if (num_t < drv_t*v2nzz) then
         826   121323964 :                   ixfer_ait2acc = 1
         827   121323964 :                   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   121323484 :                         (voltonumb_amode(nacc) - v2nzz)
         833             :                      xferfrac_num_ait2acc = xferfrac_vol_ait2acc*   &
         834   121323484 :                         (drv_t*voltonumb_amode(nacc)/num_t)
         835   121323484 :                      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   121323484 :                      else if ((xferfrac_num_ait2acc >= 1.0_r8) .or.   &
         840             :                         (xferfrac_vol_ait2acc >= 1.0_r8)) then
         841           0 :                         xferfrac_num_ait2acc = 1.0_r8
         842           0 :                         xferfrac_vol_ait2acc = 1.0_r8
         843             :                      end if
         844             :                   end if
         845   121323964 :                   xfercoef_num_ait2acc = xferfrac_num_ait2acc*tadjinv
         846   121323964 :                   xfercoef_vol_ait2acc = xferfrac_vol_ait2acc*tadjinv
         847   121323964 :                   xfertend_num(1,1) = num_a_aitsv(i,k)*xfercoef_num_ait2acc
         848   121323964 :                   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  2174023800 :             ixfer_acc2ait = 0
         859  2174023800 :             xfercoef_num_acc2ait = 0.0_r8
         860  2174023800 :             xfercoef_vol_acc2ait = 0.0_r8
         861             :
         862  2174023800 :             drv_t = drv_a_accsv(i,k) + drv_c_accsv(i,k)
         863  2174023800 :             num_t = num_a_accsv(i,k) + num_c_accsv(i,k)
         864  2174023800 :             drv_a_noxf = 0.0_r8
         865  2174023800 :             drv_c_noxf = 0.0_r8
         866  2174023800 :             if (drv_t > 0.0_r8) then
         867  2174023800 :                if (num_t > drv_t*v2nzz) then
         868  1432489548 :                   do l1 = 1, nspec_amode(nacc)
         869             :
         870  1432489548 :                      if ( noxf_acc2ait(l1) ) then
         871             :                         ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
         872   409282728 :                         dummwdens = 1.0_r8 / specdens_amode(l1,nacc)
         873   409282728 :                         la = lmassptr_amode(l1,nacc)
         874             :                         drv_a_noxf = drv_a_noxf    &
         875   409282728 :                            + max(0.0_r8,q(i,k,la))*dummwdens
         876   409282728 :                         lc = lmassptrcw_amode(l1,nacc)
         877             :
         878   409282728 :                         fldcw => qqcw_get_field(pbuf,lmassptrcw_amode(l1,nacc),lchnk)
         879             :                         drv_c_noxf = drv_c_noxf    &
         880   409282728 :                            + max(0.0_r8,fldcw(i,k))*dummwdens
         881             :                      end if
         882             :                   end do
         883   204641364 :                   drv_t_noxf = drv_a_noxf + drv_c_noxf
         884   204641364 :                   num_t_noxf = drv_t_noxf*voltonumblo_amode(nacc)
         885   204641364 :                   num_t0 = num_t
         886   204641364 :                   drv_t0 = drv_t
         887   204641364 :                   num_t = max( 0.0_r8, num_t - num_t_noxf )
         888   204641364 :                   drv_t = max( 0.0_r8, drv_t - drv_t_noxf )
         889             :                end if
         890             :             end if
         891             :
         892  2174023800 :             if (drv_t > 0.0_r8) then
         893  2174023800 :                if (num_t > drv_t*v2nzz) then
         894   204641364 :                   ixfer_acc2ait = 1
         895   204641364 :                   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   185428024 :                         (voltonumb_amode(nait) - v2nzz)
         901             :                      xferfrac_num_acc2ait = xferfrac_vol_acc2ait*   &
         902   185428024 :                         (drv_t*voltonumb_amode(nait)/num_t)
         903   185428024 :                      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   185428024 :                      else if ((xferfrac_num_acc2ait >= 1.0_r8) .or.   &
         908             :                         (xferfrac_vol_acc2ait >= 1.0_r8)) then
         909           0 :                         xferfrac_num_acc2ait = 1.0_r8
         910           0 :                         xferfrac_vol_acc2ait = 1.0_r8
         911             :                      end if
         912             :                   end if
         913   204641364 :                   duma = 1.0e-37_r8
         914             :                   xferfrac_num_acc2ait = xferfrac_num_acc2ait*   &
         915   204641364 :                      num_t/max( duma, num_t0 )
         916   204641364 :                   xfercoef_num_acc2ait = xferfrac_num_acc2ait*tadjinv
         917   204641364 :                   xfercoef_vol_acc2ait = xferfrac_vol_acc2ait*tadjinv
         918   204641364 :                   xfertend_num(2,1) = num_a_accsv(i,k)*xfercoef_num_acc2ait
         919   204641364 :                   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  2312517168 :             if (ixfer_ait2acc+ixfer_acc2ait > 0) then
         925   325965328 :                ixfer_ait2acc_sv(i,k) = ixfer_ait2acc
         926   325965328 :                ixfer_acc2ait_sv(i,k) = ixfer_acc2ait
         927             :
         928             :                !
         929             :                ! compute new dgncur & v2ncur for aitken & accum modes
         930             :                !
         931             :                ! currently inactive
         932   977895984 :                do n = nait, nacc, (nacc-nait)
         933   651930656 :                   if (n .eq. nait) then
         934   325965328 :                      duma = (xfertend_num(1,1) - xfertend_num(2,1))*deltat
         935   325965328 :                      num_a     = max( 0.0_r8, num_a_aitsv(i,k) - duma )
         936   325965328 :                      num_a_acc = max( 0.0_r8, num_a_accsv(i,k) + duma )
         937             :                      duma = (drv_a_aitsv(i,k)*xfercoef_vol_ait2acc -   &
         938   325965328 :                         (drv_a_accsv(i,k)-drv_a_noxf)*xfercoef_vol_acc2ait)*deltat
         939   325965328 :                      drv_a     = max( 0.0_r8, drv_a_aitsv(i,k) - duma )
         940   325965328 :                      drv_a_acc = max( 0.0_r8, drv_a_accsv(i,k) + duma )
         941   325965328 :                      duma = (xfertend_num(1,2) - xfertend_num(2,2))*deltat
         942   325965328 :                      num_c     = max( 0.0_r8, num_c_aitsv(i,k) - duma )
         943   325965328 :                      num_c_acc = max( 0.0_r8, num_c_accsv(i,k) + duma )
         944             :                      duma = (drv_c_aitsv(i,k)*xfercoef_vol_ait2acc -   &
         945   325965328 :                         (drv_c_accsv(i,k)-drv_c_noxf)*xfercoef_vol_acc2ait)*deltat
         946   325965328 :                      drv_c     = max( 0.0_r8, drv_c_aitsv(i,k) - duma )
         947   325965328 :                      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   651930656 :                   if (drv_a > 0.0_r8) then
         956   651930656 :                      if (num_a <= drv_a*voltonumbhi_amode(n)) then
         957   121735097 :                         dgncur_a(i,k,n) = dgnumhi_amode(n)
         958   121735097 :                         v2ncur_a(i,k,n) = voltonumbhi_amode(n)
         959   530195559 :                      else if (num_a >= drv_a*voltonumblo_amode(n)) then
         960   222378256 :                         dgncur_a(i,k,n) = dgnumlo_amode(n)
         961   222378256 :                         v2ncur_a(i,k,n) = voltonumblo_amode(n)
         962             :                      else
         963   307817303 :                         dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third
         964   307817303 :                         v2ncur_a(i,k,n) = num_a/drv_a
         965             :                      end if
         966             :                   else
         967           0 :                      dgncur_a(i,k,n) = dgnum_amode(n)
         968           0 :                      v2ncur_a(i,k,n) = voltonumb_amode(n)
         969             :                   end if
         970             :
         971   977895984 :                   if (drv_c > 0.0_r8) then
         972     6888742 :                      if (num_c <= drv_c*voltonumbhi_amode(n)) then
         973     4053733 :                         dgncur_c(i,k,n) = dgnumhi_amode(n)
         974     4053733 :                         v2ncur_c(i,k,n) = voltonumbhi_amode(n)
         975     2835009 :                      else if (num_c >= drv_c*voltonumblo_amode(n)) then
         976        8580 :                         dgncur_c(i,k,n) = dgnumlo_amode(n)
         977        8580 :                         v2ncur_c(i,k,n) = voltonumblo_amode(n)
         978             :                      else
         979     2826429 :                         dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third
         980     2826429 :                         v2ncur_c(i,k,n) = num_c/drv_c
         981             :                      end if
         982             :                   else
         983   645041914 :                      dgncur_c(i,k,n) = dgnum_amode(n)
         984   645041914 :                      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   325965328 :                if ( masterproc ) then
         995      369823 :                   if (idiagaa > 0) then
         996           6 :                      do j = 1, 2
         997          26 :                         do iq = 1, nspecfrm_renamexf(ipair)
         998          64 :                            do jac = 1, 2
         999          40 :                               if (j .eq. 1) then
        1000          20 :                                  if (jac .eq. 1) then
        1001          10 :                                     lsfrm = lspecfrma_renamexf(iq,ipair)
        1002          10 :                                     lstoo = lspectooa_renamexf(iq,ipair)
        1003             :                                  else
        1004          10 :                                     lsfrm = lspecfrmc_renamexf(iq,ipair)
        1005          10 :                                     lstoo = lspectooc_renamexf(iq,ipair)
        1006             :                                  end if
        1007             :                               else
        1008          20 :                                  if (jac .eq. 1) then
        1009          10 :                                     lsfrm = lspectooa_renamexf(iq,ipair)
        1010          10 :                                     lstoo = lspecfrma_renamexf(iq,ipair)
        1011             :                                  else
        1012          10 :                                     lsfrm = lspectooc_renamexf(iq,ipair)
        1013          10 :                                     lstoo = lspecfrmc_renamexf(iq,ipair)
        1014             :                                  end if
        1015             :                               end if
        1016          40 :                               write( 6, '(a,3i3,2i4)' ) 'calcsize j,iq,jac, lsfrm,lstoo',   &
        1017         100 :                                  j,iq,jac, lsfrm,lstoo
        1018             :                            end do
        1019             :                         end do
        1020             :                      end do
        1021             :                   end if
        1022             :                end if
        1023   325965328 :                idiagaa = -1
        1024             :
        1025             :
        1026             :                ! j=1 does aitken-->accum; j=2 does accum-->aitken
        1027   977895984 :                do  j = 1, 2
        1028             :
        1029   651930656 :                   if ((j .eq. 1 .and. ixfer_ait2acc > 0) .or. &
        1030   325965328 :                      (j .eq. 2 .and. ixfer_acc2ait > 0)) then
        1031             :
        1032   325965328 :                      jsrflx = j+2
        1033   325965328 :                      if (j .eq. 1) then
        1034             :                         xfercoef = xfercoef_vol_ait2acc
        1035             :                      else
        1036   204641364 :                         xfercoef = xfercoef_vol_acc2ait
        1037             :                      end if
        1038             :
        1039  1955791968 :                      do  iq = 1, nspecfrm_renamexf(ipair)
        1040             :
        1041             :                         ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
        1042  5215445248 :                         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  3259653280 :                            if (j .eq. 1) then
        1049  1213239640 :                               if (jac .eq. 1) then
        1050   606619820 :                                  lsfrm = lspecfrma_renamexf(iq,ipair)
        1051   606619820 :                                  lstoo = lspectooa_renamexf(iq,ipair)
        1052             :                               else
        1053   606619820 :                                  lsfrm = lspecfrmc_renamexf(iq,ipair)
        1054   606619820 :                                  lstoo = lspectooc_renamexf(iq,ipair)
        1055             :                               end if
        1056             :                            else
        1057  2046413640 :                               if (jac .eq. 1) then
        1058  1023206820 :                                  lsfrm = lspectooa_renamexf(iq,ipair)
        1059  1023206820 :                                  lstoo = lspecfrma_renamexf(iq,ipair)
        1060             :                               else
        1061  1023206820 :                                  lsfrm = lspectooc_renamexf(iq,ipair)
        1062  1023206820 :                                  lstoo = lspecfrmc_renamexf(iq,ipair)
        1063             :                               end if
        1064             :                            end if
        1065             :
        1066  4889479920 :                            if ((lsfrm > 0) .and. (lstoo > 0)) then
        1067  3259653280 :                               if (jac .eq. 1) then
        1068  1629826640 :                                  if (iq .eq. 1) then
        1069   325965328 :                                     xfertend = xfertend_num(j,jac)
        1070             :                                  else
        1071  1303861312 :                                     xfertend = max(0.0_r8,q(i,k,lsfrm))*xfercoef
        1072             :                                  end if
        1073  1629826640 :                                  dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xfertend
        1074  1629826640 :                                  dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xfertend
        1075             :                               else
        1076  1629826640 :                                  if (iq .eq. 1) then
        1077   325965328 :                                     xfertend = xfertend_num(j,jac)
        1078             :                                  else
        1079  1303861312 :                                     fldcw => qqcw_get_field(pbuf,lsfrm,lchnk)
        1080  1303861312 :                                     xfertend = max(0.0_r8,fldcw(i,k))*xfercoef
        1081             :                                  end if
        1082  1629826640 :                                  dqqcwdt(i,k,lsfrm) = dqqcwdt(i,k,lsfrm) - xfertend
        1083  1629826640 :                                  dqqcwdt(i,k,lstoo) = dqqcwdt(i,k,lstoo) + xfertend
        1084             :                               end if
        1085  3259653280 :                               qsrflx(i,lsfrm,jsrflx,jac) = qsrflx(i,lsfrm,jsrflx,jac) - xfertend*pdel_fac
        1086  3259653280 :                               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     1489176 :    lsfrm = -123456789   ! executable statement for debugging
        1101             :
        1102             :
        1103             :    !
        1104             :    ! apply tendencies to cloud-borne species MRs
        1105             :    !
        1106    62545392 :    do l = 1, pcnst
        1107    61056216 :       lc = l
        1108    62545392 :       if ( lc>0 .and. dotendqqcw(lc) ) then
        1109    17870112 :          fldcw=> qqcw_get_field(pbuf,l,lchnk)
        1110  1679790528 :          do k = top_lev, pver
        1111 27768076128 :             do i = 1, ncol
        1112           0 :                fldcw(i,k) = max( 0.0_r8,   &
        1113 27750206016 :                   (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     1489176 :    if ( .not. do_adjust ) return
        1125             :
        1126     7445880 :    do n = 1, ntot_amode
        1127     5956704 :       if (mprognum_amode(n) <= 0) cycle
        1128             :
        1129    19359288 :       do jac = 1, 2
        1130    11913408 :          if (jac == 1) then
        1131     5956704 :             l = numptr_amode(n)
        1132     5956704 :             tmpnamea = cnst_name(l)
        1133             :          else
        1134     5956704 :             l = numptrcw_amode(n)
        1135     5956704 :             tmpnamea = cnst_name_cw(l)
        1136             :          end if
        1137    11913408 :          fieldname = trim(tmpnamea) // '_sfcsiz1'
        1138    11913408 :          call outfld( fieldname, qsrflx(:,l,1,jac), pcols, lchnk)
        1139             :
        1140    11913408 :          fieldname = trim(tmpnamea) // '_sfcsiz2'
        1141    17870112 :          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     1489176 :    if ( .not. do_aitacc_transfer ) return
        1149             :
        1150     8935056 :    do iq = 1, nspecfrm_renamexf(ipair)
        1151             :
        1152             :       ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
        1153    23826816 :       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    14891760 :          if (jac .eq. 1) then
        1158     7445880 :             lsfrm = lspecfrma_renamexf(iq,ipair)
        1159     7445880 :             lstoo = lspectooa_renamexf(iq,ipair)
        1160             :          else
        1161     7445880 :             lsfrm = lspecfrmc_renamexf(iq,ipair)
        1162     7445880 :             lstoo = lspectooc_renamexf(iq,ipair)
        1163             :          end if
        1164    14891760 :          if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
        1165             :
        1166    14891760 :          if (jac .eq. 1) then
        1167     7445880 :             tmpnamea = cnst_name(lsfrm)
        1168     7445880 :             tmpnameb = cnst_name(lstoo)
        1169             :          else
        1170     7445880 :             tmpnamea = cnst_name_cw(lsfrm)
        1171     7445880 :             tmpnameb = cnst_name_cw(lstoo)
        1172             :          end if
        1173             :          if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
        1174             :
        1175    14891760 :          fieldname = trim(tmpnamea) // '_sfcsiz3'
        1176    14891760 :          call outfld( fieldname, qsrflx(:,lsfrm,3,jac), pcols, lchnk)
        1177             :
        1178    14891760 :          fieldname = trim(tmpnameb) // '_sfcsiz3'
        1179    14891760 :          call outfld( fieldname, qsrflx(:,lstoo,3,jac), pcols, lchnk)
        1180             :
        1181    14891760 :          fieldname = trim(tmpnamea) // '_sfcsiz4'
        1182    14891760 :          call outfld( fieldname, qsrflx(:,lsfrm,4,jac), pcols, lchnk)
        1183             :
        1184    14891760 :          fieldname = trim(tmpnameb) // '_sfcsiz4'
        1185    22337640 :          call outfld( fieldname, qsrflx(:,lstoo,4,jac), pcols, lchnk)
        1186             :
        1187             :       end do   ! jac = ...
        1188             :    end do   ! iq = ...
        1189             :
        1190     1489176 :    call modal_aero_calcdry(state, pbuf)
        1191             :
        1192             : #endif
        1193             :
        1194     1490712 : end subroutine modal_aero_calcsize_sub
        1195             :
        1196             :
        1197             : !----------------------------------------------------------------------
        1198             :
        1199             :
        1200        6192 : 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        3096 :    real(r8), pointer :: dgncur_a(:,:) ! (pcols,pver)
        1234             :
        1235             :
        1236             :    real(r8), parameter :: third = 1.0_r8/3.0_r8
        1237             :
        1238        3096 :    real(r8), pointer :: mode_num(:,:) ! mode number mixing ratio
        1239        3096 :    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        3096 :    lchnk = state%lchnk
        1256        3096 :    ncol  = state%ncol
        1257             :
        1258        3096 :    list_idx = 0  ! climate list by default
        1259           0 :    if (present(list_idx_in)) list_idx = list_idx_in
        1260             :
        1261        3096 :    call rad_cnst_get_info(list_idx, nmodes=nmodes)
        1262             :
        1263        3096 :    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       15480 :    do n = 1, nmodes
        1330             :
        1331       12384 :       if (list_idx == 0) then
        1332       49536 :          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       12384 :                                    sigmag=sigmag)
        1340             :
        1341             :       ! get mode number mixing ratio
        1342       12384 :       call rad_cnst_get_mode_num(list_idx, n, 'a', state, pbuf, mode_num)
        1343             :
        1344    19591488 :       dgncur_a(:,:) = dgnum
        1345       12384 :       dryvol_a(:,:) = 0.0_r8
        1346             :
        1347             :       ! compute dry volume mixrats =
        1348             :       !      sum_over_components{ component_mass mixrat / density }
        1349       12384 :       call rad_cnst_get_info(list_idx, n, nspec=nspec)
        1350       58824 :       do l1 = 1, nspec
        1351             :
        1352       46440 :          call rad_cnst_get_aer_mmr(list_idx, n, l1, 'a', state, pbuf, specmmr)
        1353       46440 :          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       46440 :          dummwdens = 1.0_r8 / specdens
        1357             :
        1358     4377744 :          do k=top_lev,pver
        1359    72162360 :             do i=1,ncol
        1360   135594000 :                dryvol_a(i,k) = dryvol_a(i,k)    &
        1361   207709920 :                   + max(0.0_r8, specmmr(i,k))*dummwdens
        1362             :             end do
        1363             :          end do
        1364             :       end do
        1365             :
        1366       12384 :       alnsg  = log( sigmag )
        1367       12384 :       dumfac = exp(4.5_r8*alnsg**2)*pi/6.0_r8
        1368       12384 :       voltonumblo = 1._r8 / ( (pi/6._r8)*(dgnumlo**3)*exp(4.5_r8*alnsg**2) )
        1369       12384 :       voltonumbhi = 1._r8 / ( (pi/6._r8)*(dgnumhi**3)*exp(4.5_r8*alnsg**2) )
        1370       12384 :       v2nxx = voltonumbhi
        1371       12384 :       v2nyy = voltonumblo
        1372       12384 :       dgnxx = dgnumhi
        1373       12384 :       dgnyy = dgnumlo
        1374             :
        1375     1191960 :       do k = top_lev, pver
        1376    19243296 :          do i = 1, ncol
        1377             :
        1378    18079200 :             drv_a = dryvol_a(i,k)
        1379    18079200 :             num_a0 = mode_num(i,k)
        1380    18079200 :             num_a = max( 0.0_r8, num_a0 )
        1381             :
        1382    19230912 :             if (drv_a > 0.0_r8) then
        1383    18079200 :                if (num_a <= drv_a*v2nxx) then
        1384     2716438 :                   dgncur_a(i,k) = dgnxx
        1385    15362762 :                else if (num_a >= drv_a*v2nyy) then
        1386     1760257 :                   dgncur_a(i,k) = dgnyy
        1387             :                else
        1388    13602505 :                   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       24768 :    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        3096 : end subroutine modal_aero_calcsize_diag
        1400             :
        1401     1492272 : 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     1492272 :    real(r8), pointer :: maer(:,:)        ! aerosol wet mass MR (including water) (kg/kg-air)
        1418     1492272 :    real(r8), pointer :: hygro(:,:,:)     ! volume-weighted mean hygroscopicity (--)
        1419     1492272 :    real(r8), pointer :: dryvol(:,:,:)    ! single-particle-mean dry volume (m3)
        1420     1492272 :    real(r8), pointer :: dryrad(:,:,:)    ! dry volume mean radius of aerosol (m)
        1421     1492272 :    real(r8), pointer :: drymass(:,:,:)   ! single-particle-mean dry mass  (kg)
        1422     1492272 :    real(r8), pointer :: so4dryvol(:,:,:) ! single-particle-mean so4 dry volume (m3)
        1423     1492272 :    real(r8), pointer :: naer(:,:,:)      ! aerosol number MR (bounded!) (#/kg-air)
        1424             :
        1425     1492272 :    real(r8), pointer :: dgncur_a(:,:,:)
        1426     1492272 :    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     1492272 :    lchnk = state%lchnk
        1450     1492272 :    ncol = state%ncol
        1451             :
        1452     1492272 :    list_idx = 0
        1453     1492272 :    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     1492272 :    call rad_cnst_get_info(list_idx, nmodes=nmodes)
        1471             :
        1472     1492272 :    allocate( maer(pcols,pver))
        1473             :
        1474     1492272 :    if (list_idx == 0) then
        1475     1492272 :       call pbuf_get_field(pbuf, dgnum_idx,     dgncur_a )
        1476     1492272 :       call pbuf_get_field(pbuf, hygro_idx,     hygro)
        1477     1492272 :       call pbuf_get_field(pbuf, dryvol_idx,    dryvol)
        1478     1492272 :       call pbuf_get_field(pbuf, dryrad_idx,    dryrad)
        1479     1492272 :       call pbuf_get_field(pbuf, drymass_idx,   drymass)
        1480     1492272 :       call pbuf_get_field(pbuf, so4dryvol_idx, so4dryvol)
        1481     1492272 :       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  9444589488 :    hygro(:,:,:)     = 0._r8
        1493  9444589488 :    so4dryvol(:,:,:) = 0._r8
        1494             :
        1495     7461360 :    do m = 1, nmodes
        1496             :
        1497  9443097216 :       maer(:,:)      = 0._r8
        1498     5969088 :       dryvolmr(:,:) = 0._r8
        1499     5969088 :       so4dryvolmr(:,:) = 0._r8
        1500             :
        1501             :       ! get mode properties
        1502     5969088 :       call rad_cnst_get_mode_props(list_idx, m, sigmag=sigmag)
        1503             :
        1504             :       ! get mode info
        1505     5969088 :       call rad_cnst_get_info(list_idx, m, nspec=nspec)
        1506             :
        1507    28353168 :       do l = 1, nspec
        1508             :
        1509             :          ! get species interstitial mixing ratio ('a')
        1510    22384080 :          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    22384080 :                                      hygro_aer=spechygro, spectype=spectype)
        1513             :
        1514    22384080 :          if (l == 1) then
        1515             :             ! save off these values to be used as defaults
        1516     5969088 :             spechygro_1    = spechygro
        1517             :          end if
        1518             :
        1519  2132456688 :          do k = top_lev, pver
        1520 34782257520 :             do i = 1, ncol
        1521 32678154000 :                duma          = raer(i,k)     ! kg/kg air
        1522 32678154000 :                maer(i,k)     = maer(i,k) + duma
        1523 32678154000 :                dumb          = duma/specdens ! m3/kg air
        1524 32678154000 :                dryvolmr(i,k) = dryvolmr(i,k) + dumb
        1525 32678154000 :                if (modal_strat_sulfate .and. (trim(spectype).eq.'sulfate')) then
        1526           0 :                   so4dryvolmr(i,k) = so4dryvolmr(i,k) + dumb
        1527             :                end if
        1528 34759873440 :                hygro(i,k,m)  = hygro(i,k,m) + dumb*spechygro
        1529             :             end do
        1530             :          end do
        1531             :       end do
        1532             :
        1533     5969088 :       alnsg = log(sigmag)
        1534             :
        1535   561094272 :       do k = top_lev, pver
        1536  9275268672 :          do i = 1, ncol
        1537             :
        1538  8714174400 :             if (dryvolmr(i,k) > 1.0e-30_r8) then
        1539  8604268133 :                hygro(i,k,m) = hygro(i,k,m)/dryvolmr(i,k)
        1540             :             else
        1541   109906267 :                hygro(i,k,m) = spechygro_1
        1542             :             end if
        1543             :
        1544             :             ! dry aerosol properties
        1545             :
        1546  8714174400 :             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  8714174400 :             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  8714174400 :             if (maer(i,k) .gt. 1.0e-31_r8) then
        1556  8645238371 :                drydens = maer(i,k)/dryvolmr(i,k)        ! kg/m3 aerosol
        1557             :             else
        1558             :                drydens = 1.0_r8
        1559             :             end if
        1560  8714174400 :             dryvol(i,k,m)   = 1.0_r8/v2ncur_a             ! m3/particle
        1561  8714174400 :             drymass(i,k,m)  = drydens*dryvol(i,k,m)       ! kg/particle
        1562  9269299584 :             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    13430448 :       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     1492272 :    deallocate( maer)
        1584             :
        1585             :
        1586     1492272 : end subroutine modal_aero_calcdry
        1587             : !----------------------------------------------------------------------
        1588             :
        1589             : end module modal_aero_calcsize

Generated by: LCOV version 1.14