module modal_aero_calcsize ! RCE 07.04.13: Adapted from MIRAGE2 code use shr_kind_mod, only: r8 => shr_kind_r8 use spmd_utils, only: masterproc use physconst, only: pi, rhoh2o, gravit use ppgrid, only: pcols, pver use physics_types, only: physics_state, physics_ptend use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_old_tim_idx, pbuf_get_field use phys_control, only: phys_getopts use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, rad_cnst_get_aer_props, & rad_cnst_get_mode_props, rad_cnst_get_mode_num use cam_logfile, only: iulog use abortutils, only: endrun use cam_history, only: addfld, add_default, fieldname_len, phys_decomp, outfld use constituents, only: pcnst, cnst_name use ref_pres, only: top_lev => clim_modal_aero_top_lev #ifdef MODAL_AERO ! these are the variables needed for the diagnostic calculation of dry radius use modal_aero_data, only: ntot_amode, nspec_amode, & numptr_amode, & alnsg_amode, & voltonumbhi_amode, voltonumblo_amode, & dgnum_amode, dgnumhi_amode, dgnumlo_amode ! these variables are needed for the prognostic calculations to exchange mass ! between modes use modal_aero_data, only: numptrcw_amode, mprognum_amode, qqcw_get_field, lmassptrcw_amode, & lmassptr_amode, modeptr_accum, modeptr_aitken, & !ntot_aspectype, & ! lspectype_amode, & specmw_amode, specdens_amode, voltonumb_amode, & cnst_name_cw,nspec_max use modal_aero_rename, only: lspectooa_renamexf, lspecfrma_renamexf, lspectooc_renamexf, lspecfrmc_renamexf, & modetoo_renamexf, nspecfrm_renamexf, npair_renamexf, modefrm_renamexf #endif implicit none private save public modal_aero_calcsize_init, modal_aero_calcsize_sub, modal_aero_calcsize_diag logical :: do_adjust_default logical :: do_aitacc_transfer_default integer :: dgnum_idx = -1 !=============================================================================== contains !=============================================================================== subroutine modal_aero_calcsize_init() !----------------------------------------------------------------------- ! ! Purpose: ! set do_adjust_default and do_aitacc_transfer_default flags ! create history fields for column tendencies associated with ! modal_aero_calcsize ! ! Author: R. Easter ! !----------------------------------------------------------------------- ! local integer :: ipair, iq integer :: jac integer :: lsfrm, lstoo integer :: n, nacc, nait logical :: history_aerosol character(len=fieldname_len) :: tmpnamea, tmpnameb character(len=fieldname_len+3) :: fieldname character(128) :: long_name character(8) :: unit !----------------------------------------------------------------------- call phys_getopts(history_aerosol_out=history_aerosol) ! init entities required for both prescribed and prognostic modes dgnum_idx = pbuf_get_index('DGNUM') #ifndef MODAL_AERO do_adjust_default = .false. do_aitacc_transfer_default = .false. #else ! do_adjust_default allows adjustment to be turned on/off do_adjust_default = .true. ! do_aitacc_transfer_default allows aitken <--> accum mode transfer to be turned on/off ! *** it can only be true when aitken & accum modes are both present ! and have prognosed number and diagnosed surface/sigmag nait = modeptr_aitken nacc = modeptr_accum do_aitacc_transfer_default = .false. if ((modeptr_aitken > 0) .and. & (modeptr_accum > 0) .and. & (modeptr_aitken /= modeptr_accum)) then do_aitacc_transfer_default = .true. if (mprognum_amode(nait) <= 0) do_aitacc_transfer_default = .false. if (mprognum_amode(nacc) <= 0) do_aitacc_transfer_default = .false. end if if ( .not. do_adjust_default ) return ! define history fields for number-adjust source-sink for all modes do n = 1, ntot_amode if (mprognum_amode(n) <= 0) cycle do jac = 1, 2 if (jac == 1) then tmpnamea = cnst_name(numptr_amode(n)) else tmpnamea = cnst_name_cw(numptrcw_amode(n)) end if unit = '#/m2/s' fieldname = trim(tmpnamea) // '_sfcsiz1' long_name = trim(tmpnamea) // ' calcsize number-adjust column source' call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp ) if (history_aerosol) then call add_default(fieldname, 1, ' ') end if if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname fieldname = trim(tmpnamea) // '_sfcsiz2' long_name = trim(tmpnamea) // ' calcsize number-adjust column sink' call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp ) if (history_aerosol) then call add_default(fieldname, 1, ' ') end if if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname end do ! jac = ... end do ! n = ... if ( .not. do_aitacc_transfer_default ) return ! check that renaming ipair=1 is aitken-->accum ipair = 1 if ((modefrm_renamexf(ipair) .ne. nait) .or. & (modetoo_renamexf(ipair) .ne. nacc)) then write( 6, '(//2a//)' ) & '*** modal_aero_calcaersize_init error -- ', & 'modefrm/too_renamexf(1) are wrong' call endrun( 'modal_aero_calcaersize_init error' ) end if ! define history fields for aitken-accum transfer do iq = 1, nspecfrm_renamexf(ipair) ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c"); do jac = 1, 2 ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species ! the lspectooa_renamexf (and lspectooc_renamexf) are accum species if (jac .eq. 1) then lsfrm = lspecfrma_renamexf(iq,ipair) lstoo = lspectooa_renamexf(iq,ipair) else lsfrm = lspecfrmc_renamexf(iq,ipair) lstoo = lspectooc_renamexf(iq,ipair) end if if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle if (jac .eq. 1) then tmpnamea = cnst_name(lsfrm) tmpnameb = cnst_name(lstoo) else tmpnamea = cnst_name_cw(lsfrm) tmpnameb = cnst_name_cw(lstoo) end if unit = 'kg/m2/s' if ((tmpnamea(1:3) == 'num') .or. & (tmpnamea(1:3) == 'NUM')) unit = '#/m2/s' fieldname = trim(tmpnamea) // '_sfcsiz3' long_name = trim(tmpnamea) // ' calcsize aitken-to-accum adjust column tendency' call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp ) if (history_aerosol) then call add_default(fieldname, 1, ' ') end if if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname fieldname = trim(tmpnameb) // '_sfcsiz3' long_name = trim(tmpnameb) // ' calcsize aitken-to-accum adjust column tendency' call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp ) if (history_aerosol) then call add_default(fieldname, 1, ' ') end if if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname fieldname = trim(tmpnamea) // '_sfcsiz4' long_name = trim(tmpnamea) // ' calcsize accum-to-aitken adjust column tendency' call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp ) if (history_aerosol) then call add_default(fieldname, 1, ' ') end if if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname fieldname = trim(tmpnameb) // '_sfcsiz4' long_name = trim(tmpnameb) // ' calcsize accum-to-aitken adjust column tendency' call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp ) if (history_aerosol) then call add_default(fieldname, 1, ' ') end if if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname end do ! jac = ... end do ! iq = ... #endif end subroutine modal_aero_calcsize_init !=============================================================================== subroutine modal_aero_calcsize_sub(state, ptend, deltat, pbuf, do_adjust_in, & do_aitacc_transfer_in) !----------------------------------------------------------------------- ! ! Calculates aerosol size distribution parameters ! mprognum_amode > 0 ! calculate Dgnum from mass, number, and fixed sigmag ! mprognum_amode <= 0 ! calculate number from mass, fixed Dgnum, and fixed sigmag ! ! Also (optionally) adjusts prognostic number to ! be within bounds determined by mass, Dgnum bounds, and sigma bounds ! ! Author: R. Easter ! !----------------------------------------------------------------------- ! arguments type(physics_state), target, intent(in) :: state ! Physics state variables type(physics_ptend), target, intent(inout) :: ptend ! indivdual parameterization tendencies real(r8), intent(in) :: deltat ! model time-step size (s) type(physics_buffer_desc), pointer :: pbuf(:) ! physics buffer logical, optional :: do_adjust_in logical, optional :: do_aitacc_transfer_in #ifdef MODAL_AERO ! local logical :: do_adjust logical :: do_aitacc_transfer integer :: lchnk ! chunk identifier integer :: ncol ! number of columns real(r8), pointer :: t(:,:) ! Temperature in Kelvin real(r8), pointer :: pmid(:,:) ! pressure at model levels (Pa) real(r8), pointer :: pdel(:,:) ! pressure thickness of levels real(r8), pointer :: q(:,:,:) ! Tracer MR array logical, pointer :: dotend(:) ! flag for doing tendency real(r8), pointer :: dqdt(:,:,:) ! TMR tendency array real(r8), pointer :: dgncur_a(:,:,:) integer :: i, icol_diag, iduma, ipair, iq integer :: ixfer_acc2ait, ixfer_ait2acc integer :: ixfer_acc2ait_sv(pcols,pver), ixfer_ait2acc_sv(pcols,pver) integer :: j, jac, jsrflx, k integer :: l, l1, la, lc, lna, lnc, lsfrm, lstoo integer :: n, nacc, nait integer, save :: idiagaa = 1 logical :: dotendqqcw(pcnst) ! logical :: noxf_acc2ait(ntot_aspectype) logical :: noxf_acc2ait(nspec_max) character(len=fieldname_len) :: tmpnamea, tmpnameb character(len=fieldname_len+3) :: fieldname real(r8), parameter :: third = 1.0_r8/3.0_r8 real(r8), pointer :: fldcw(:,:) real(r8) :: delnum_a2, delnum_c2 ! work variables real(r8) :: delnum_a3, delnum_c3, delnum_t3 ! work variables real(r8) :: deltatinv ! 1/deltat real(r8) :: dgncur_c(pcols,pver,ntot_amode) real(r8) :: dgnyy, dgnxx ! dgnumlo/hi of current mode real(r8) :: dqqcwdt(pcols,pver,pcnst) ! cloudborne TMR tendency array real(r8) :: drv_a, drv_c, drv_t ! dry volume (cm3/mol_air) real(r8) :: drv_t0 real(r8) :: drv_a_noxf, drv_c_noxf, drv_t_noxf real(r8) :: drv_a_acc, drv_c_acc real(r8) :: drv_a_accsv(pcols,pver), drv_c_accsv(pcols,pver) real(r8) :: drv_a_aitsv(pcols,pver), drv_c_aitsv(pcols,pver) real(r8) :: drv_a_sv(pcols,pver,ntot_amode), drv_c_sv(pcols,pver,ntot_amode) real(r8) :: dryvol_a(pcols,pver) ! interstital aerosol dry ! volume (cm^3/mol_air) real(r8) :: dryvol_c(pcols,pver) ! activated aerosol dry volume real(r8) :: duma, dumb, dumc, dumd ! work variables real(r8) :: dumfac, dummwdens ! work variables real(r8) :: frelaxadj ! relaxation factor applied ! to size bounds real(r8) :: fracadj ! deltat/tadj real(r8) :: num_a0, num_c0, num_t0 ! initial number (#/mol_air) real(r8) :: num_a1, num_c1 ! working number (#/mol_air) real(r8) :: num_a2, num_c2, num_t2 ! working number (#/mol_air) real(r8) :: num_a, num_c, num_t ! final number (#/mol_air) real(r8) :: num_t_noxf real(r8) :: numbnd ! bounded number real(r8) :: num_a_acc, num_c_acc real(r8) :: num_a_accsv(pcols,pver), num_c_accsv(pcols,pver) real(r8) :: num_a_aitsv(pcols,pver), num_c_aitsv(pcols,pver) real(r8) :: num_a_sv(pcols,pver,ntot_amode), num_c_sv(pcols,pver,ntot_amode) real(r8) :: pdel_fac ! real(r8) :: tadj ! adjustment time scale real(r8) :: tadjinv ! 1/tadj real(r8) :: v2ncur_a(pcols,pver,ntot_amode) real(r8) :: v2ncur_c(pcols,pver,ntot_amode) real(r8) :: v2nyy, v2nxx, v2nzz ! voltonumblo/hi of current mode real(r8) :: v2nyyrl, v2nxxrl ! relaxed voltonumblo/hi real(r8) :: xfercoef real(r8) :: xfercoef_num_acc2ait, xfercoef_vol_acc2ait real(r8) :: xfercoef_num_ait2acc, xfercoef_vol_ait2acc real(r8) :: xferfrac_num_acc2ait, xferfrac_vol_acc2ait real(r8) :: xferfrac_num_ait2acc, xferfrac_vol_ait2acc real(r8) :: xfertend, xfertend_num(2,2) integer, parameter :: nsrflx = 4 ! last dimension of qsrflx real(r8) :: qsrflx(pcols,pcnst,nsrflx,2) ! process-specific column tracer tendencies ! 3rd index -- ! 1="standard" number adjust gain; ! 2="standard" number adjust loss; ! 3=aitken-->accum renaming; 4=accum-->aitken) ! 4th index -- ! 1="a" species; 2="c" species !----------------------------------------------------------------------- if (present(do_adjust_in)) then do_adjust = do_adjust_in else do_adjust = do_adjust_default end if if (present(do_aitacc_transfer_in)) then do_aitacc_transfer = do_aitacc_transfer_in else do_aitacc_transfer = do_aitacc_transfer_default end if lchnk = state%lchnk ncol = state%ncol t => state%t pmid => state%pmid pdel => state%pdel q => state%q dotend => ptend%lq dqdt => ptend%q call pbuf_get_field(pbuf, dgnum_idx, dgncur_a) dotendqqcw(:) = .false. dqqcwdt(:,:,:) = 0.0_r8 qsrflx(:,:,:,:) = 0.0_r8 nait = modeptr_aitken nacc = modeptr_accum deltatinv = 1.0_r8/(deltat*(1.0_r8 + 1.0e-15_r8)) ! tadj = adjustment time scale for number, surface when they are prognosed ! currently set to deltat tadj = deltat tadj = 86400 tadj = max( tadj, deltat ) tadjinv = 1.0_r8/(tadj*(1.0_r8 + 1.0e-15_r8)) fracadj = deltat*tadjinv fracadj = max( 0.0_r8, min( 1.0_r8, fracadj ) ) ! ! ! the "do 40000" loop does the original (pre jan-2006) ! number adjustment, one mode at a time ! this artificially adjusts number when mean particle size is too large ! or too small ! ! do n = 1, ntot_amode ! initialize all parameters to the default values for the mode do k=top_lev,pver do i=1,ncol ! sgcur_a(i,k,n) = sigmag_amode(n) ! sgcur_c(i,k,n) = sigmag_amode(n) dgncur_a(i,k,n) = dgnum_amode(n) dgncur_c(i,k,n) = dgnum_amode(n) v2ncur_a(i,k,n) = voltonumb_amode(n) v2ncur_c(i,k,n) = voltonumb_amode(n) dryvol_a(i,k) = 0.0_r8 dryvol_c(i,k) = 0.0_r8 end do end do ! compute dry volume mixrats = ! sum_over_components{ component_mass mixrat / density } do l1 = 1, nspec_amode(n) ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air ! dummwdens = 1.0_r8 / specdens_amode(lspectype_amode(l1,n) dummwdens = 1.0_r8 / specdens_amode(l1,n) la = lmassptr_amode(l1,n) do k=top_lev,pver do i=1,ncol dryvol_a(i,k) = dryvol_a(i,k) & + max(0.0_r8,q(i,k,la))*dummwdens end do end do fldcw => qqcw_get_field(pbuf,lmassptrcw_amode(l1,n),lchnk) do k=top_lev,pver do i=1,ncol dryvol_c(i,k) = dryvol_c(i,k) & + max(0.0_r8,fldcw(i,k))*dummwdens end do end do end do ! set "short-hand" number pointers lna = numptr_amode(n) lnc = numptrcw_amode(n) fldcw => qqcw_get_field(pbuf,numptrcw_amode(n),lchnk,.true.) ! go to section for appropriate number/surface diagnosed/prognosed options if (mprognum_amode(n) <= 0) then ! option 1 -- number diagnosed (fixed dgnum and sigmag) ! compute number tendencies that will bring numbers to their ! current diagnosed values ! if (lna > 0) then dotend(lna) = .true. do k=top_lev,pver do i=1,ncol dqdt(i,k,lna) = (dryvol_a(i,k)*voltonumb_amode(n) & - q(i,k,lna)) * deltatinv end do end do end if if (lnc > 0) then dotendqqcw(lnc) = .true. do k=top_lev,pver do i=1,ncol dqqcwdt(i,k,lnc) = (dryvol_c(i,k)*voltonumb_amode(n) & - fldcw(i,k)) * deltatinv end do end do end if else ! ! option 2 -- number prognosed (variable dgnum, fixed sigmag) ! Compute number tendencies to adjust numbers if they are outside ! the limits determined by current volume and dgnumlo/hi ! The interstitial and activated aerosol fractions can, at times, ! be the lower or upper tail of the "total" distribution. Thus they ! can be expected to have a greater range of size parameters than ! what is specified for the total distribution (via dgnumlo/hi) ! When both the interstitial and activated dry volumes are positive, ! the adjustment strategy is to (1) adjust the interstitial and activated ! numbers towards relaxed bounds, then (2) adjust the total/combined ! number towards the primary bounds. ! ! note ! v2nyy = voltonumblo_amode is proportional to dgnumlo**(-3), ! and produces the maximum allowed number for a given volume ! v2nxx = voltonumbhi_amode is proportional to dgnumhi**(-3), ! and produces the minimum allowed number for a given volume ! v2nxxrl and v2nyyrl are their "relaxed" equivalents. ! Setting frelaxadj=27=3**3 means that ! dgnumlo_relaxed = dgnumlo/3 and dgnumhi_relaxed = dgnumhi*3 ! ! if do_aitacc_transfer is .true., then ! for n=nacc, multiply v2nyy by 1.0e6 to effectively turn off the ! adjustment when number is too big (size is too small) ! for n=nait, divide v2nxx by 1.0e6 to effectively turn off the ! adjustment when number is too small (size is too big) !OLD however, do not change the v2nyyrl/v2nxxrl so that !OLD the interstitial<-->activated adjustment is not changed !NEW also change the v2nyyrl/v2nxxrl so that !NEW the interstitial<-->activated adjustment is turned off ! end if frelaxadj = 27.0_r8 dumfac = exp(4.5_r8*alnsg_amode(n)**2)*pi/6.0_r8 v2nxx = voltonumbhi_amode(n) v2nyy = voltonumblo_amode(n) v2nxxrl = v2nxx/frelaxadj v2nyyrl = v2nyy*frelaxadj dgnxx = dgnumhi_amode(n) dgnyy = dgnumlo_amode(n) if ( do_aitacc_transfer ) then if (n == nait) v2nxx = v2nxx/1.0e6_r8 if (n == nacc) v2nyy = v2nyy*1.0e6_r8 v2nxxrl = v2nxx/frelaxadj ! NEW v2nyyrl = v2nyy*frelaxadj ! NEW end if if (do_adjust) then dotend(lna) = .true. dotendqqcw(lnc) = .true. end if do k = top_lev, pver do i = 1, ncol drv_a = dryvol_a(i,k) num_a0 = q(i,k,lna) num_a = max( 0.0_r8, num_a0 ) drv_c = dryvol_c(i,k) num_c0 = fldcw(i,k) num_c = max( 0.0_r8, num_c0 ) if ( do_adjust) then ! ! do number adjustment for interstitial and activated particles ! adjustments that (1) make numbers non-negative or (2) make numbers ! zero when volume is zero are applied over time-scale deltat ! adjustments that bring numbers to within specified bounds are ! applied over time-scale tadj ! if ((drv_a <= 0.0_r8) .and. (drv_c <= 0.0_r8)) then ! both interstitial and activated volumes are zero ! adjust both numbers to zero num_a = 0.0_r8 dqdt(i,k,lna) = -num_a0*deltatinv num_c = 0.0_r8 dqqcwdt(i,k,lnc) = -num_c0*deltatinv else if (drv_c <= 0.0_r8) then ! activated volume is zero, so interstitial number/volume == total/combined ! apply step 1 and 3, but skip the relaxed adjustment (step 2, see below) num_c = 0.0_r8 dqqcwdt(i,k,lnc) = -num_c0*deltatinv num_a1 = num_a numbnd = max( drv_a*v2nxx, min( drv_a*v2nyy, num_a1 ) ) num_a = num_a1 + (numbnd - num_a1)*fracadj dqdt(i,k,lna) = (num_a - num_a0)*deltatinv else if (drv_a <= 0.0_r8) then ! interstitial volume is zero, treat similar to above num_a = 0.0_r8 dqdt(i,k,lna) = -num_a0*deltatinv num_c1 = num_c numbnd = max( drv_c*v2nxx, min( drv_c*v2nyy, num_c1 ) ) num_c = num_c1 + (numbnd - num_c1)*fracadj dqqcwdt(i,k,lnc) = (num_c - num_c0)*deltatinv else ! both volumes are positive ! apply 3 adjustment steps ! step1: num_a,c0 --> num_a,c1 forces non-negative values num_a1 = num_a num_c1 = num_c ! step2: num_a,c1 --> num_a,c2 applies relaxed bounds to the interstitial ! and activated number (individually) ! if only only a or c changes, adjust the other in the opposite direction ! as much as possible to conserve a+c numbnd = max( drv_a*v2nxxrl, min( drv_a*v2nyyrl, num_a1 ) ) delnum_a2 = (numbnd - num_a1)*fracadj num_a2 = num_a1 + delnum_a2 numbnd = max( drv_c*v2nxxrl, min( drv_c*v2nyyrl, num_c1 ) ) delnum_c2 = (numbnd - num_c1)*fracadj num_c2 = num_c1 + delnum_c2 if ((delnum_a2 == 0.0_r8) .and. (delnum_c2 /= 0.0_r8)) then num_a2 = max( drv_a*v2nxxrl, min( drv_a*v2nyyrl, & num_a1-delnum_c2 ) ) else if ((delnum_a2 /= 0.0_r8) .and. (delnum_c2 == 0.0_r8)) then num_c2 = max( drv_c*v2nxxrl, min( drv_c*v2nyyrl, & num_c1-delnum_a2 ) ) end if ! step3: num_a,c2 --> num_a,c3 applies stricter bounds to the ! combined/total number drv_t = drv_a + drv_c num_t2 = num_a2 + num_c2 delnum_a3 = 0.0_r8 delnum_c3 = 0.0_r8 if (num_t2 < drv_t*v2nxx) then delnum_t3 = (drv_t*v2nxx - num_t2)*fracadj ! if you are here then (num_a2 < drv_a*v2nxx) and/or ! (num_c2 < drv_c*v2nxx) must be true if ((num_a2 < drv_a*v2nxx) .and. (num_c2 < drv_c*v2nxx)) then delnum_a3 = delnum_t3*(num_a2/num_t2) delnum_c3 = delnum_t3*(num_c2/num_t2) else if (num_c2 < drv_c*v2nxx) then delnum_c3 = delnum_t3 else if (num_a2 < drv_a*v2nxx) then delnum_a3 = delnum_t3 end if else if (num_t2 > drv_t*v2nyy) then delnum_t3 = (drv_t*v2nyy - num_t2)*fracadj ! if you are here then (num_a2 > drv_a*v2nyy) and/or ! (num_c2 > drv_c*v2nyy) must be true if ((num_a2 > drv_a*v2nyy) .and. (num_c2 > drv_c*v2nyy)) then delnum_a3 = delnum_t3*(num_a2/num_t2) delnum_c3 = delnum_t3*(num_c2/num_t2) else if (num_c2 > drv_c*v2nyy) then delnum_c3 = delnum_t3 else if (num_a2 > drv_a*v2nyy) then delnum_a3 = delnum_t3 end if end if num_a = num_a2 + delnum_a3 dqdt(i,k,lna) = (num_a - num_a0)*deltatinv num_c = num_c2 + delnum_c3 dqqcwdt(i,k,lnc) = (num_c - num_c0)*deltatinv end if end if ! do_adjust ! ! now compute current dgn and v2n ! if (drv_a > 0.0_r8) then if (num_a <= drv_a*v2nxx) then dgncur_a(i,k,n) = dgnxx v2ncur_a(i,k,n) = v2nxx else if (num_a >= drv_a*v2nyy) then dgncur_a(i,k,n) = dgnyy v2ncur_a(i,k,n) = v2nyy else dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third v2ncur_a(i,k,n) = num_a/drv_a end if end if pdel_fac = pdel(i,k)/gravit ! = rho*dz jac = 1 qsrflx(i,lna,1,jac) = qsrflx(i,lna,1,jac) + max(0.0_r8,dqdt(i,k,lna))*pdel_fac qsrflx(i,lna,2,jac) = qsrflx(i,lna,2,jac) + min(0.0_r8,dqdt(i,k,lna))*pdel_fac if (drv_c > 0.0_r8) then if (num_c <= drv_c*v2nxx) then dgncur_c(i,k,n) = dgnumhi_amode(n) v2ncur_c(i,k,n) = v2nxx else if (num_c >= drv_c*v2nyy) then dgncur_c(i,k,n) = dgnumlo_amode(n) v2ncur_c(i,k,n) = v2nyy else dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third v2ncur_c(i,k,n) = num_c/drv_c end if end if jac = 2 qsrflx(i,lnc,1,jac) = qsrflx(i,lnc,1,jac) + max(0.0_r8,dqqcwdt(i,k,lnc))*pdel_fac qsrflx(i,lnc,2,jac) = qsrflx(i,lnc,2,jac) + min(0.0_r8,dqqcwdt(i,k,lnc))*pdel_fac ! save number and dryvol for aitken <--> accum renaming if ( do_aitacc_transfer ) then if (n == nait) then drv_a_aitsv(i,k) = drv_a num_a_aitsv(i,k) = num_a drv_c_aitsv(i,k) = drv_c num_c_aitsv(i,k) = num_c else if (n == nacc) then drv_a_accsv(i,k) = drv_a num_a_accsv(i,k) = num_a drv_c_accsv(i,k) = drv_c num_c_accsv(i,k) = num_c end if end if drv_a_sv(i,k,n) = drv_a num_a_sv(i,k,n) = num_a drv_c_sv(i,k,n) = drv_c num_c_sv(i,k,n) = num_c end do end do ! ! option 3 -- number and surface prognosed (variable dgnum and sigmag) ! this is not implemented ! end do ! do n = 1, ntot_amode ! ! ! the following section (from here to label 49000) ! does aitken <--> accum mode transfer ! ! when the aitken mode mean size is too big, the largest ! aitken particles are transferred into the accum mode ! to reduce the aitken mode mean size ! when the accum mode mean size is too small, the smallest ! accum particles are transferred into the aitken mode ! to increase the accum mode mean size ! ! ixfer_ait2acc_sv(:,:) = 0 ixfer_acc2ait_sv(:,:) = 0 if ( do_aitacc_transfer ) then ! old - on time first step, npair_renamexf will be <= 0, ! in which case need to do modal_aero_rename_init ! new - init is now done through chem_init and things below it if (npair_renamexf .le. 0) then npair_renamexf = 0 ! call modal_aero_rename_init if (npair_renamexf .le. 0) then write( 6, '(//a//)' ) & '*** modal_aero_calcaersize_sub error -- npair_renamexf <= 0' call endrun( 'modal_aero_calcaersize_sub error' ) end if end if ! check that renaming ipair=1 is aitken-->accum ipair = 1 if ((modefrm_renamexf(ipair) .ne. nait) .or. & (modetoo_renamexf(ipair) .ne. nacc)) then write( 6, '(//2a//)' ) & '*** modal_aero_calcaersize_sub error -- ', & 'modefrm/too_renamexf(1) are wrong' call endrun( 'modal_aero_calcaersize_sub error' ) end if ! set dotend() for species that will be transferred do iq = 1, nspecfrm_renamexf(ipair) lsfrm = lspecfrma_renamexf(iq,ipair) lstoo = lspectooa_renamexf(iq,ipair) if ((lsfrm > 0) .and. (lstoo > 0)) then dotend(lsfrm) = .true. dotend(lstoo) = .true. end if lsfrm = lspecfrmc_renamexf(iq,ipair) lstoo = lspectooc_renamexf(iq,ipair) if ((lsfrm > 0) .and. (lstoo > 0)) then dotendqqcw(lsfrm) = .true. dotendqqcw(lstoo) = .true. end if end do ! identify accum species cannot be transferred to aitken mode noxf_acc2ait(:) = .true. do l1 = 1, nspec_amode(nacc) la = lmassptr_amode(l1,nacc) do iq = 1, nspecfrm_renamexf(ipair) if (lspectooa_renamexf(iq,ipair) == la) then noxf_acc2ait(l1) = .false. end if end do end do ! v2nzz is voltonumb at the "geometrically-defined" mid-point ! between the aitken and accum modes v2nzz = sqrt(voltonumb_amode(nait)*voltonumb_amode(nacc)) ! loop over columns and levels do k = top_lev, pver do i = 1, ncol pdel_fac = pdel(i,k)/gravit ! = rho*dz xfertend_num(:,:) = 0.0_r8 ! compute aitken --> accum transfer rates ixfer_ait2acc = 0 xfercoef_num_ait2acc = 0.0_r8 xfercoef_vol_ait2acc = 0.0_r8 drv_t = drv_a_aitsv(i,k) + drv_c_aitsv(i,k) num_t = num_a_aitsv(i,k) + num_c_aitsv(i,k) if (drv_t > 0.0_r8) then if (num_t < drv_t*v2nzz) then ixfer_ait2acc = 1 if (num_t < drv_t*voltonumb_amode(nacc)) then xferfrac_num_ait2acc = 1.0_r8 xferfrac_vol_ait2acc = 1.0_r8 else xferfrac_vol_ait2acc = ((num_t/drv_t) - v2nzz)/ & (voltonumb_amode(nacc) - v2nzz) xferfrac_num_ait2acc = xferfrac_vol_ait2acc* & (drv_t*voltonumb_amode(nacc)/num_t) if ((xferfrac_num_ait2acc <= 0.0_r8) .or. & (xferfrac_vol_ait2acc <= 0.0_r8)) then xferfrac_num_ait2acc = 0.0_r8 xferfrac_vol_ait2acc = 0.0_r8 else if ((xferfrac_num_ait2acc >= 1.0_r8) .or. & (xferfrac_vol_ait2acc >= 1.0_r8)) then xferfrac_num_ait2acc = 1.0_r8 xferfrac_vol_ait2acc = 1.0_r8 end if end if xfercoef_num_ait2acc = xferfrac_num_ait2acc*tadjinv xfercoef_vol_ait2acc = xferfrac_vol_ait2acc*tadjinv xfertend_num(1,1) = num_a_aitsv(i,k)*xfercoef_num_ait2acc xfertend_num(1,2) = num_c_aitsv(i,k)*xfercoef_num_ait2acc end if end if ! compute accum --> aitken transfer rates ! accum may have some species (seasalt, dust, poa, lll) that are ! not in aitken mode ! so first divide the accum drv & num into not-transferred (noxf) species ! and transferred species, and use the transferred-species ! portion in what follows ixfer_acc2ait = 0 xfercoef_num_acc2ait = 0.0_r8 xfercoef_vol_acc2ait = 0.0_r8 drv_t = drv_a_accsv(i,k) + drv_c_accsv(i,k) num_t = num_a_accsv(i,k) + num_c_accsv(i,k) drv_a_noxf = 0.0_r8 drv_c_noxf = 0.0_r8 if (drv_t > 0.0_r8) then if (num_t > drv_t*v2nzz) then do l1 = 1, nspec_amode(nacc) if ( noxf_acc2ait(l1) ) then ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air ! dummwdens = 1.0_r8 / specdens_amode(lspectype_amode(l1,nacc)) dummwdens = 1.0_r8 / specdens_amode(l1,nacc) la = lmassptr_amode(l1,nacc) drv_a_noxf = drv_a_noxf & + max(0.0_r8,q(i,k,la))*dummwdens lc = lmassptrcw_amode(l1,nacc) fldcw => qqcw_get_field(pbuf,lmassptrcw_amode(l1,nacc),lchnk) drv_c_noxf = drv_c_noxf & + max(0.0_r8,fldcw(i,k))*dummwdens end if end do drv_t_noxf = drv_a_noxf + drv_c_noxf num_t_noxf = drv_t_noxf*voltonumblo_amode(nacc) num_t0 = num_t drv_t0 = drv_t num_t = max( 0.0_r8, num_t - num_t_noxf ) drv_t = max( 0.0_r8, drv_t - drv_t_noxf ) end if end if if (drv_t > 0.0_r8) then if (num_t > drv_t*v2nzz) then ixfer_acc2ait = 1 if (num_t > drv_t*voltonumb_amode(nait)) then xferfrac_num_acc2ait = 1.0_r8 xferfrac_vol_acc2ait = 1.0_r8 else xferfrac_vol_acc2ait = ((num_t/drv_t) - v2nzz)/ & (voltonumb_amode(nait) - v2nzz) xferfrac_num_acc2ait = xferfrac_vol_acc2ait* & (drv_t*voltonumb_amode(nait)/num_t) if ((xferfrac_num_acc2ait <= 0.0_r8) .or. & (xferfrac_vol_acc2ait <= 0.0_r8)) then xferfrac_num_acc2ait = 0.0_r8 xferfrac_vol_acc2ait = 0.0_r8 else if ((xferfrac_num_acc2ait >= 1.0_r8) .or. & (xferfrac_vol_acc2ait >= 1.0_r8)) then xferfrac_num_acc2ait = 1.0_r8 xferfrac_vol_acc2ait = 1.0_r8 end if end if duma = 1.0e-37_r8 xferfrac_num_acc2ait = xferfrac_num_acc2ait* & num_t/max( duma, num_t0 ) xfercoef_num_acc2ait = xferfrac_num_acc2ait*tadjinv xfercoef_vol_acc2ait = xferfrac_vol_acc2ait*tadjinv xfertend_num(2,1) = num_a_accsv(i,k)*xfercoef_num_acc2ait xfertend_num(2,2) = num_c_accsv(i,k)*xfercoef_num_acc2ait end if end if ! jump to end-of-loop if no transfer is needed at current i,k if (ixfer_ait2acc+ixfer_acc2ait > 0) then ixfer_ait2acc_sv(i,k) = ixfer_ait2acc ixfer_acc2ait_sv(i,k) = ixfer_acc2ait ! ! compute new dgncur & v2ncur for aitken & accum modes ! ! currently inactive do n = nait, nacc, (nacc-nait) if (n .eq. nait) then duma = (xfertend_num(1,1) - xfertend_num(2,1))*deltat num_a = max( 0.0_r8, num_a_aitsv(i,k) - duma ) num_a_acc = max( 0.0_r8, num_a_accsv(i,k) + duma ) duma = (drv_a_aitsv(i,k)*xfercoef_vol_ait2acc - & (drv_a_accsv(i,k)-drv_a_noxf)*xfercoef_vol_acc2ait)*deltat drv_a = max( 0.0_r8, drv_a_aitsv(i,k) - duma ) drv_a_acc = max( 0.0_r8, drv_a_accsv(i,k) + duma ) duma = (xfertend_num(1,2) - xfertend_num(2,2))*deltat num_c = max( 0.0_r8, num_c_aitsv(i,k) - duma ) num_c_acc = max( 0.0_r8, num_c_accsv(i,k) + duma ) duma = (drv_c_aitsv(i,k)*xfercoef_vol_ait2acc - & (drv_c_accsv(i,k)-drv_c_noxf)*xfercoef_vol_acc2ait)*deltat drv_c = max( 0.0_r8, drv_c_aitsv(i,k) - duma ) drv_c_acc = max( 0.0_r8, drv_c_accsv(i,k) + duma ) else num_a = num_a_acc drv_a = drv_a_acc num_c = num_c_acc drv_c = drv_c_acc end if if (drv_a > 0.0_r8) then if (num_a <= drv_a*voltonumbhi_amode(n)) then dgncur_a(i,k,n) = dgnumhi_amode(n) v2ncur_a(i,k,n) = voltonumbhi_amode(n) else if (num_a >= drv_a*voltonumblo_amode(n)) then dgncur_a(i,k,n) = dgnumlo_amode(n) v2ncur_a(i,k,n) = voltonumblo_amode(n) else dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third v2ncur_a(i,k,n) = num_a/drv_a end if else dgncur_a(i,k,n) = dgnum_amode(n) v2ncur_a(i,k,n) = voltonumb_amode(n) end if if (drv_c > 0.0_r8) then if (num_c <= drv_c*voltonumbhi_amode(n)) then dgncur_c(i,k,n) = dgnumhi_amode(n) v2ncur_c(i,k,n) = voltonumbhi_amode(n) else if (num_c >= drv_c*voltonumblo_amode(n)) then dgncur_c(i,k,n) = dgnumlo_amode(n) v2ncur_c(i,k,n) = voltonumblo_amode(n) else dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third v2ncur_c(i,k,n) = num_c/drv_c end if else dgncur_c(i,k,n) = dgnum_amode(n) v2ncur_c(i,k,n) = voltonumb_amode(n) end if end do ! ! compute tendency amounts for aitken <--> accum transfer ! if ( masterproc ) then if (idiagaa > 0) then do j = 1, 2 do iq = 1, nspecfrm_renamexf(ipair) do jac = 1, 2 if (j .eq. 1) then if (jac .eq. 1) then lsfrm = lspecfrma_renamexf(iq,ipair) lstoo = lspectooa_renamexf(iq,ipair) else lsfrm = lspecfrmc_renamexf(iq,ipair) lstoo = lspectooc_renamexf(iq,ipair) end if else if (jac .eq. 1) then lsfrm = lspectooa_renamexf(iq,ipair) lstoo = lspecfrma_renamexf(iq,ipair) else lsfrm = lspectooc_renamexf(iq,ipair) lstoo = lspecfrmc_renamexf(iq,ipair) end if end if write( 6, '(a,3i3,2i4)' ) 'calcsize j,iq,jac, lsfrm,lstoo', & j,iq,jac, lsfrm,lstoo end do end do end do end if end if idiagaa = -1 ! j=1 does aitken-->accum; j=2 does accum-->aitken do j = 1, 2 if ((j .eq. 1 .and. ixfer_ait2acc > 0) .or. & (j .eq. 2 .and. ixfer_acc2ait > 0)) then jsrflx = j+2 if (j .eq. 1) then xfercoef = xfercoef_vol_ait2acc else xfercoef = xfercoef_vol_acc2ait end if do iq = 1, nspecfrm_renamexf(ipair) ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c"); do jac = 1, 2 ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species ! the lspectooa_renamexf (and lspectooc_renamexf) are accum species ! for j=1, want lsfrm=aitken species, lstoo=accum species ! for j=2, want lsfrm=accum species, lstoo=aitken species if (j .eq. 1) then if (jac .eq. 1) then lsfrm = lspecfrma_renamexf(iq,ipair) lstoo = lspectooa_renamexf(iq,ipair) else lsfrm = lspecfrmc_renamexf(iq,ipair) lstoo = lspectooc_renamexf(iq,ipair) end if else if (jac .eq. 1) then lsfrm = lspectooa_renamexf(iq,ipair) lstoo = lspecfrma_renamexf(iq,ipair) else lsfrm = lspectooc_renamexf(iq,ipair) lstoo = lspecfrmc_renamexf(iq,ipair) end if end if if ((lsfrm > 0) .and. (lstoo > 0)) then if (jac .eq. 1) then if (iq .eq. 1) then xfertend = xfertend_num(j,jac) else xfertend = max(0.0_r8,q(i,k,lsfrm))*xfercoef end if dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xfertend dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xfertend else if (iq .eq. 1) then xfertend = xfertend_num(j,jac) else fldcw => qqcw_get_field(pbuf,lsfrm,lchnk) xfertend = max(0.0_r8,fldcw(i,k))*xfercoef end if dqqcwdt(i,k,lsfrm) = dqqcwdt(i,k,lsfrm) - xfertend dqqcwdt(i,k,lstoo) = dqqcwdt(i,k,lstoo) + xfertend end if qsrflx(i,lsfrm,jsrflx,jac) = qsrflx(i,lsfrm,jsrflx,jac) - xfertend*pdel_fac qsrflx(i,lstoo,jsrflx,jac) = qsrflx(i,lstoo,jsrflx,jac) + xfertend*pdel_fac end if end do end do end if end do end if end do end do end if ! do_aitacc_transfer lsfrm = -123456789 ! executable statement for debugging ! ! apply tendencies to cloud-borne species MRs ! do l = 1, pcnst lc = l if ( lc>0 .and. dotendqqcw(lc) ) then fldcw=> qqcw_get_field(pbuf,l,lchnk) do k = top_lev, pver do i = 1, ncol fldcw(i,k) = max( 0.0_r8, & (fldcw(i,k) + dqqcwdt(i,k,lc)*deltat) ) end do end do end if end do ! ! do outfld calls ! ! history fields for number-adjust source-sink for all modes if ( .not. do_adjust ) return do n = 1, ntot_amode if (mprognum_amode(n) <= 0) cycle do jac = 1, 2 if (jac == 1) then l = numptr_amode(n) tmpnamea = cnst_name(l) else l = numptrcw_amode(n) tmpnamea = cnst_name_cw(l) end if fieldname = trim(tmpnamea) // '_sfcsiz1' call outfld( fieldname, qsrflx(:,l,1,jac), pcols, lchnk) fieldname = trim(tmpnamea) // '_sfcsiz2' call outfld( fieldname, qsrflx(:,l,2,jac), pcols, lchnk) end do ! jac = ... end do ! n = ... ! history fields for aitken-accum transfer if ( .not. do_aitacc_transfer ) return do iq = 1, nspecfrm_renamexf(ipair) ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c"); do jac = 1, 2 ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species ! the lspectooa_renamexf (and lspectooc_renamexf) are accum species if (jac .eq. 1) then lsfrm = lspecfrma_renamexf(iq,ipair) lstoo = lspectooa_renamexf(iq,ipair) else lsfrm = lspecfrmc_renamexf(iq,ipair) lstoo = lspectooc_renamexf(iq,ipair) end if if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle if (jac .eq. 1) then tmpnamea = cnst_name(lsfrm) tmpnameb = cnst_name(lstoo) else tmpnamea = cnst_name_cw(lsfrm) tmpnameb = cnst_name_cw(lstoo) end if if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle fieldname = trim(tmpnamea) // '_sfcsiz3' call outfld( fieldname, qsrflx(:,lsfrm,3,jac), pcols, lchnk) fieldname = trim(tmpnameb) // '_sfcsiz3' call outfld( fieldname, qsrflx(:,lstoo,3,jac), pcols, lchnk) fieldname = trim(tmpnamea) // '_sfcsiz4' call outfld( fieldname, qsrflx(:,lsfrm,4,jac), pcols, lchnk) fieldname = trim(tmpnameb) // '_sfcsiz4' call outfld( fieldname, qsrflx(:,lstoo,4,jac), pcols, lchnk) end do ! jac = ... end do ! iq = ... #endif end subroutine modal_aero_calcsize_sub !---------------------------------------------------------------------- subroutine modal_aero_calcsize_diag(state, pbuf, list_idx_in, dgnum_m) !----------------------------------------------------------------------- ! ! Calculate aerosol size distribution parameters ! ! ***N.B.*** DGNUM for the modes in the climate list are put directly into ! the physics buffer. For diagnostic list calculations use the ! optional list_idx and dgnum args. !----------------------------------------------------------------------- ! arguments type(physics_state), intent(in), target :: state ! Physics state variables type(physics_buffer_desc), pointer :: pbuf(:) ! physics buffer integer, optional, intent(in) :: list_idx_in ! diagnostic list index real(r8), optional, pointer :: dgnum_m(:,:,:) ! interstital aerosol dry number mode radius (m) ! local integer :: i, k, l1, n integer :: lchnk, ncol integer :: list_idx, stat integer :: nmodes integer :: nspec real(r8), pointer :: dgncur_a(:,:) ! (pcols,pver) real(r8), parameter :: third = 1.0_r8/3.0_r8 real(r8), pointer :: mode_num(:,:) ! mode number mixing ratio real(r8), pointer :: specmmr(:,:) ! specie mmr real(r8) :: specdens ! specie density real(r8) :: dryvol_a(pcols,pver) ! interstital aerosol dry volume (cm^3/mol_air) real(r8) :: dgnum, dgnumhi, dgnumlo real(r8) :: dgnyy, dgnxx ! dgnumlo/hi of current mode real(r8) :: drv_a ! dry volume (cm3/mol_air) real(r8) :: dumfac, dummwdens ! work variables real(r8) :: num_a0 ! initial number (#/mol_air) real(r8) :: num_a ! final number (#/mol_air) real(r8) :: voltonumbhi, voltonumblo real(r8) :: v2nyy, v2nxx ! voltonumblo/hi of current mode real(r8) :: sigmag, alnsg !----------------------------------------------------------------------- lchnk = state%lchnk ncol = state%ncol list_idx = 0 ! climate list by default if (present(list_idx_in)) list_idx = list_idx_in call rad_cnst_get_info(list_idx, nmodes=nmodes) if (list_idx /= 0) then if (.not. present(dgnum_m)) then call endrun('modal_aero_calcsize_diag called for'// & 'diagnostic list but dgnum_m pointer not present') end if allocate(dgnum_m(pcols,pver,nmodes), stat=stat) if (stat > 0) then call endrun('modal_aero_calcsize_diag: allocation FAILURE: dgnum_m') end if end if do n = 1, nmodes if (list_idx == 0) then call pbuf_get_field(pbuf, dgnum_idx, dgncur_a, start=(/1,1,n/), kount=(/pcols,pver,1/)) else dgncur_a => dgnum_m(:,:,n) end if ! get mode properties call rad_cnst_get_mode_props(list_idx, n, dgnum=dgnum, dgnumhi=dgnumhi, dgnumlo=dgnumlo, & sigmag=sigmag) ! get mode number mixing ratio call rad_cnst_get_mode_num(list_idx, n, 'a', state, pbuf, mode_num) dgncur_a(:,:) = dgnum dryvol_a(:,:) = 0.0_r8 ! compute dry volume mixrats = ! sum_over_components{ component_mass mixrat / density } call rad_cnst_get_info(list_idx, n, nspec=nspec) do l1 = 1, nspec call rad_cnst_get_aer_mmr(list_idx, n, l1, 'a', state, pbuf, specmmr) call rad_cnst_get_aer_props(list_idx, n, l1, density_aer=specdens) ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air dummwdens = 1.0_r8 / specdens do k=top_lev,pver do i=1,ncol dryvol_a(i,k) = dryvol_a(i,k) & + max(0.0_r8, specmmr(i,k))*dummwdens end do end do end do alnsg = log( sigmag ) dumfac = exp(4.5_r8*alnsg**2)*pi/6.0_r8 voltonumblo = 1._r8 / ( (pi/6._r8)*(dgnumlo**3)*exp(4.5_r8*alnsg**2) ) voltonumbhi = 1._r8 / ( (pi/6._r8)*(dgnumhi**3)*exp(4.5_r8*alnsg**2) ) v2nxx = voltonumbhi v2nyy = voltonumblo dgnxx = dgnumhi dgnyy = dgnumlo do k = top_lev, pver do i = 1, ncol drv_a = dryvol_a(i,k) num_a0 = mode_num(i,k) num_a = max( 0.0_r8, num_a0 ) if (drv_a > 0.0_r8) then if (num_a <= drv_a*v2nxx) then dgncur_a(i,k) = dgnxx else if (num_a >= drv_a*v2nyy) then dgncur_a(i,k) = dgnyy else dgncur_a(i,k) = (drv_a/(dumfac*num_a))**third end if end if end do end do end do ! nmodes end subroutine modal_aero_calcsize_diag !---------------------------------------------------------------------- end module modal_aero_calcsize