chem/module_cam_mam_calcsize.F

   1 #define WRF_PORT
   2 #define MODAL_AERO
   3 ! module_cam_mam_calcsize.F
   4 ! adapted from cam3 modal_aero_calcsize.F90 by r.c.easter, june 2010
   5 !
   6 ! 2010-06-24 rce notes
   7 !     > all calls to history routines (addfld, outfld) deactivated
   8 !     > diagnostics section (for testing) in modal_aero_calcsize_sub deactivated
   9
  10 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
  11
  12 ! modal_aero_calcsize.F90
  13 !----------------------------------------------------------------------
  14 !BOP
  15 !
  16 ! !MODULE: modal_aero_calcsize --- modal aerosol size calc and number adjustment
  17 !
  18 ! !INTERFACE:
  19    module modal_aero_calcsize
  20
  21 ! !USES:
  22 !  currently none
  23
  24    implicit none
  25    private
  26    save
  27
  28 !  !PUBLIC MEMBER FUNCTIONS:
  29    public modal_aero_calcsize_sub, modal_aero_calcsize_init
  30
  31 ! !PUBLIC DATA MEMBERS:
  32    logical  :: do_adjust
  33    logical  :: do_aitacc_transfer
  34
  35
  36 ! !DESCRIPTION: This module implements ...
  37 !
  38 ! !REVISION HISTORY:
  39 !
  40 !   RCE 07.04.13:  Adapted from MIRAGE2 code
  41 !
  42 !EOP
  43 !----------------------------------------------------------------------
  44 !BOC
  45 ! list private module data here
  46 !EOC
  47 !----------------------------------------------------------------------
  48
  49
  50   contains
  51
  52
  53 !----------------------------------------------------------------------
  54
  55
  56     subroutine modal_aero_calcsize_sub(            &
  57          lchnk,   ncol,    nstep,   icalcaer_flag, &
  58          loffset,                                  &
  59          aero_mmr_flag,    deltat,                 &
  60          t,       pmid,    pdel,    q,             &
  61          dqdt,    dotend,                          &
  62 #ifdef WRF_PORT
  63          qqcw,                                     &
  64 #endif
  65          dgncur_a       )
  66
  67       !-----------------------------------------------------------------------
  68       !
  69       ! Purpose:
  70       ! Calculates aerosol size distribution parameters
  71       !    mprognum_amode >  0
  72       !       calculate Dgnum from mass, number, and fixed sigmag
  73       !    mprognum_amode <= 0
  74       !       calculate number from mass, fixed Dgnum, and fixed sigmag
  75       !
  76       ! Also (optionally) adjusts prognostic number to
  77       !    be within bounds determined by mass, Dgnum bounds, and sigma bounds
  78       !
  79       ! Method:
  80       !
  81       ! Author: R. Easter
  82       !
  83       !-----------------------------------------------------------------------
  84
  85       use modal_aero_data,  only: numptrcw_amode, mprognum_amode, numptr_amode, &
  86 #ifndef WRF_PORT
  87 qqcw_get_field,
  88 #endif
  89       lmassptrcw_amode, &
  90            alnsg_amode, voltonumblo_amode, lmassptr_amode, modeptr_accum, ntot_amode, modeptr_aitken, ntot_aspectype, &
  91            dgnum_amode, lspectype_amode, specmw_amode, specdens_amode, voltonumb_amode, nspec_amode, dgnumlo_amode, &
  92            cnst_name_cw, voltonumbhi_amode, dgnumhi_amode
  93       use modal_aero_rename, only: lspectooa_renamexf, lspecfrma_renamexf, lspectooc_renamexf, lspecfrmc_renamexf, &
  94            modetoo_renamexf, nspecfrm_renamexf, npair_renamexf, modefrm_renamexf
  95
  96       use shr_kind_mod,  only:  r8 => shr_kind_r8
  97 #ifndef WRF_PORT
  98       use mo_constants,  only:  pi
  99 #endif
 100       use physconst,     only:  avogad, gravit, mwdry, rair
 101 #ifndef WRF_PORT
 102       use phys_grid,     only:  get_lat_all_p, get_lon_all_p
 103       use pmgrid,        only:  plat, plon
 104       use ppgrid,        only:  begchunk, endchunk, pcols, pver
 105       use constituents,  only:  pcnst, cnst_name
 106       use chem_mods,     only:  adv_mass, gas_pcnst
 107       use mo_tracname,   only:  solsym
 108
 109       use abortutils,    only:  endrun
 110       use cam_history,   only:  fieldname_len, outfld
 111       use spmd_utils,    only:  masterproc
 112 #else
 113       use physconst,     only: pi
 114       use module_cam_support, only: pcols, pver, pcnst =>pcnst_runtime, &
 115            endrun,fieldname_len, outfld, masterproc
 116       use constituents,  only: cnst_name
 117 #endif
 118
 119       implicit none
 120
 121       !-----------------------------------------------------------------------
 122       ! arguments
 123       integer, intent(in)  :: lchnk                ! chunk identifier
 124       integer, intent(in)  :: ncol                 ! number of columns
 125
 126       integer, intent(in)  :: nstep                ! model time-step number
 127       integer, intent(in)  :: icalcaer_flag
 128       ! identifies call location from typhysbc (1,2)
 129       integer, intent(in)  :: loffset
 130       !  integer, intent(in)  :: nsrflx               ! last dimension of qsrflx
 131
 132       logical, intent(in)  :: aero_mmr_flag        ! if .true.,  aerosol q are kg/kg-air
 133       ! if .false., aerosol q are mol/mol-air
 134
 135       real(r8), intent(in) :: deltat               ! model time-step size (s)
 136       real(r8), intent(in) :: t(pcols,pver)        ! Temperature in Kelvin
 137       real(r8), intent(in) :: pmid(pcols,pver)     ! pressure at model levels (Pa)
 138       real(r8), intent(in) :: pdel(pcols,pver)     ! pressure thickness of levels
 139       real(r8), intent(in) :: q(pcols,pver,pcnst)  ! Tracer MR array
 140
 141       logical,  intent(inout) :: dotend(pcnst)          ! flag for doing tendency
 142       real(r8), intent(inout) :: dqdt(pcols,pver,pcnst) ! TMR tendency array
 143 #ifdef WRF_PORT
 144       real(r8), target, intent(inout) :: qqcw(pcols,pver,pcnst) ! cloudborne tracer MR array
 145 #endif
 146       !  real(r8), intent(out)   :: qsrflx(pcols,pcnst,nsrflx)
 147       ! process-specific column tracer tendencies
 148       ! (1="standard" number adjust gain;
 149       !  2="standard" number adjust loss;
 150       !  3=aitken-->accum renaming; 4=accum-->aitken)
 151       real(r8), intent(inout) :: dgncur_a(pcols,pver,ntot_amode)
 152
 153       !-----------------------------------------------------------------------
 154       ! local
 155       integer  :: i, icol_diag, iduma, ipair, iq
 156       integer  :: ixfer_acc2ait, ixfer_ait2acc
 157       integer  :: ixfer_acc2ait_sv(pcols,pver), ixfer_ait2acc_sv(pcols,pver)
 158       integer  :: j, jac, jsrflx, k
 159       integer  :: l, l1, la, lc, lna, lnc, lsfrm, lstoo
 160       integer  :: n, nacc, nait
 161       integer  :: lat(pcols), lon(pcols)
 162
 163       integer, save  :: idiagaa = 1
 164
 165       logical  :: dotendqqcw(pcnst)
 166       logical  :: noxf_acc2ait(ntot_aspectype)
 167
 168       character(len=fieldname_len)   :: tmpnamea, tmpnameb
 169       character(len=fieldname_len+3) :: fieldname
 170
 171       real(r8), parameter :: third = 1.0_r8/3.0_r8
 172       real(r8), pointer :: fldcw(:,:)
 173       real(r8) :: delnum_a2, delnum_c2            !  work variables
 174       real(r8) :: delnum_a3, delnum_c3, delnum_t3 !  work variables
 175       real(r8) :: deltatinv                     ! 1/deltat
 176       real(r8) :: dgncur_c(pcols,pver,ntot_amode)
 177       real(r8) :: dgnyy, dgnxx                  ! dgnumlo/hi of current mode
 178       real(r8) :: dqqcwdt(pcols,pver,pcnst)     ! cloudborne TMR tendency array
 179       real(r8) :: drv_a, drv_c, drv_t           ! dry volume (cm3/mol_air)
 180       real(r8) :: drv_t0
 181       real(r8) :: drv_a_noxf, drv_c_noxf, drv_t_noxf
 182       real(r8) :: drv_a_acc, drv_c_acc
 183       real(r8) :: drv_a_accsv(pcols,pver), drv_c_accsv(pcols,pver)
 184       real(r8) :: drv_a_aitsv(pcols,pver), drv_c_aitsv(pcols,pver)
 185       real(r8) :: drv_a_sv(pcols,pver,ntot_amode), drv_c_sv(pcols,pver,ntot_amode)
 186       real(r8) :: dryvol_a(pcols,pver)          ! interstital aerosol dry
 187       ! volume (cm^3/mol_air)
 188       real(r8) :: dryvol_c(pcols,pver)          ! activated aerosol dry volume
 189       real(r8) :: duma, dumb, dumc, dumd        ! work variables
 190       real(r8) :: dumfac, dummwdens             ! work variables
 191       real(r8) :: frelaxadj                     ! relaxation factor applied
 192       ! to size bounds
 193       real(r8) :: fracadj                       ! deltat/tadj
 194       real(r8) :: num_a0, num_c0, num_t0        ! initial number (#/mol_air)
 195       real(r8) :: num_a1, num_c1                ! working number (#/mol_air)
 196       real(r8) :: num_a2, num_c2, num_t2        ! working number (#/mol_air)
 197       real(r8) :: num_a, num_c, num_t           ! final number (#/mol_air)
 198       real(r8) :: num_t_noxf
 199       real(r8) :: numbnd                        ! bounded number
 200       real(r8) :: num_a_acc, num_c_acc
 201       real(r8) :: num_a_accsv(pcols,pver), num_c_accsv(pcols,pver)
 202       real(r8) :: num_a_aitsv(pcols,pver), num_c_aitsv(pcols,pver)
 203       real(r8) :: num_a_sv(pcols,pver,ntot_amode), num_c_sv(pcols,pver,ntot_amode)
 204       real(r8) :: pdel_fac                      !
 205       real(r8) :: tadj                          ! adjustment time scale
 206       real(r8) :: tadjinv                       ! 1/tadj
 207       real(r8) :: v2ncur_a(pcols,pver,ntot_amode)
 208       real(r8) :: v2ncur_c(pcols,pver,ntot_amode)
 209       real(r8) :: v2nyy, v2nxx, v2nzz           ! voltonumblo/hi of current mode
 210       real(r8) :: v2nyyrl, v2nxxrl              ! relaxed voltonumblo/hi
 211       real(r8) :: xfercoef
 212       real(r8) :: xfercoef_num_acc2ait, xfercoef_vol_acc2ait
 213       real(r8) :: xfercoef_num_ait2acc, xfercoef_vol_ait2acc
 214       real(r8) :: xferfrac_num_acc2ait, xferfrac_vol_acc2ait
 215       real(r8) :: xferfrac_num_ait2acc, xferfrac_vol_ait2acc
 216       real(r8) :: xfertend, xfertend_num(2,2)
 217
 218       integer, parameter :: nsrflx = 4    ! last dimension of qsrflx
 219       real(r8) :: qsrflx(pcols,pcnst,nsrflx,2)
 220       ! process-specific column tracer tendencies
 221       ! 3rd index --
 222       !    1="standard" number adjust gain;
 223       !    2="standard" number adjust loss;
 224       !    3=aitken-->accum renaming; 4=accum-->aitken)
 225       ! 4th index --
 226       !    1="a" species; 2="c" species
 227       !-----------------------------------------------------------------------
 228
 229       dotendqqcw(:) = .false.
 230       dqqcwdt(:,:,:) = 0.0_r8
 231       qsrflx(:,:,:,:) = 0.0_r8
 232
 233       nait = modeptr_aitken
 234       nacc = modeptr_accum
 235
 236       deltatinv = 1.0/(deltat*(1.0d0 + 1.0d-15))
 237       ! tadj = adjustment time scale for number, surface when they are prognosed
 238       !           currently set to deltat
 239       tadj = deltat
 240       tadj = 86400
 241       tadj = max( tadj, deltat )
 242       tadjinv = 1.0/(tadj*(1.0d0 + 1.0d-15))
 243       fracadj = deltat*tadjinv
 244       fracadj = max( 0.0_r8, min( 1.0_r8, fracadj ) )
 245
 246
 247       !
 248       !
 249       ! the "do 40000" loop does the original (pre jan-2006)
 250       !   number adjustment, one mode at a time
 251       ! this artificially adjusts number when mean particle size is too large
 252       !   or too small
 253       !
 254       !
 255       do n = 1, ntot_amode
 256
 257
 258          ! initialize all parameters to the default values for the mode
 259          do k=1,pver
 260             do i=1,ncol
 261                !    sgcur_a(i,k,n) = sigmag_amode(n)
 262                !    sgcur_c(i,k,n) = sigmag_amode(n)
 263                dgncur_a(i,k,n) = dgnum_amode(n)
 264                dgncur_c(i,k,n) = dgnum_amode(n)
 265                v2ncur_a(i,k,n) = voltonumb_amode(n)
 266                v2ncur_c(i,k,n) = voltonumb_amode(n)
 267                dryvol_a(i,k) = 0.0_r8
 268                dryvol_c(i,k) = 0.0_r8
 269             end do
 270          end do
 271
 272          ! compute dry volume mixrats =
 273          !      sum_over_components{ component_mass mixrat / density }
 274          do l1 = 1, nspec_amode(n)
 275             if ( aero_mmr_flag ) then
 276                ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
 277                dummwdens = 1.0_r8 / specdens_amode(lspectype_amode(l1,n))
 278             else
 279                ! need qmass*dummwdens = (kmol/kmol-air) * [(kg/kmol)/(kg/m3)] = m3/mol-air
 280                dummwdens = specmw_amode(lspectype_amode(l1,n))   &
 281                     / specdens_amode(lspectype_amode(l1,n))
 282             end if
 283             la = lmassptr_amode(l1,n) - loffset
 284             do k=1,pver
 285                do i=1,ncol
 286                   dryvol_a(i,k) = dryvol_a(i,k)    &
 287                        + max(0.0_r8,q(i,k,la))*dummwdens
 288                end do
 289             end do
 290             !      lc = lmassptrcw_amode(l1,n) - loffset
 291 #ifndef WRF_PORT
 292             fldcw => qqcw_get_field(lmassptrcw_amode(l1,n),lchnk)
 293 #else
 294             fldcw => qqcw(:,:,lmassptrcw_amode(l1,n))
 295 #endif
 296             do k=1,pver
 297                do i=1,ncol
 298                   dryvol_c(i,k) = dryvol_c(i,k)    &
 299                        + max(0.0_r8,fldcw(i,k))*dummwdens
 300                end do
 301             end do
 302          end do
 303
 304          ! set "short-hand" number pointers
 305          lna = numptr_amode(n) - loffset
 306          lnc = numptrcw_amode(n) - loffset
 307 #ifndef WRF_PORT
 308          fldcw => qqcw_get_field(numptrcw_amode(n),lchnk,.true.)
 309 #else
 310          fldcw => qqcw(:,:,numptrcw_amode(n))
 311 #endif
 312
 313
 314          ! go to section for appropriate number/surface diagnosed/prognosed options
 315          if (mprognum_amode(n) <= 0) then
 316
 317             ! option 1 -- number diagnosed (fixed dgnum and sigmag)
 318             !    compute number tendencies that will bring numbers to their
 319             !    current diagnosed values
 320             !
 321             if (lna > 0) then
 322                dotend(lna) = .true.
 323                do k=1,pver
 324                   do i=1,ncol
 325                      dqdt(i,k,lna) = (dryvol_a(i,k)*voltonumb_amode(n)   &
 326                           - q(i,k,lna)) * deltatinv
 327                   end do
 328                end do
 329             end if
 330             if (lnc > 0) then
 331                dotendqqcw(lnc) = .true.
 332                do k=1,pver
 333                   do i=1,ncol
 334                      dqqcwdt(i,k,lnc) = (dryvol_c(i,k)*voltonumb_amode(n)   &
 335                           - fldcw(i,k)) * deltatinv
 336                   end do
 337                end do
 338             end if
 339          else
 340
 341
 342             !
 343             ! option 2 -- number prognosed (variable dgnum, fixed sigmag)
 344             !       Compute number tendencies to adjust numbers if they are outside
 345             !    the limits determined by current volume and dgnumlo/hi
 346             !       The interstitial and activated aerosol fractions can, at times,
 347             !    be the lower or upper tail of the "total" distribution.  Thus they
 348             !    can be expected to have a greater range of size parameters than
 349             !    what is specified for the total distribution (via dgnumlo/hi)
 350             !       When both the interstitial and activated dry volumes are positive,
 351             !    the adjustment strategy is to (1) adjust the interstitial and activated
 352             !    numbers towards relaxed bounds, then (2) adjust the total/combined
 353             !    number towards the primary bounds.
 354             !
 355             ! note
 356             !    v2nyy = voltonumblo_amode is proportional to dgnumlo**(-3),
 357             !            and produces the maximum allowed number for a given volume
 358             !    v2nxx = voltonumbhi_amode is proportional to dgnumhi**(-3),
 359             !            and produces the minimum allowed number for a given volume
 360             !    v2nxxrl and v2nyyrl are their "relaxed" equivalents.
 361             !            Setting frelaxadj=27=3**3 means that
 362             !            dgnumlo_relaxed = dgnumlo/3 and dgnumhi_relaxed = dgnumhi*3
 363             !
 364             ! if do_aitacc_transfer is .true., then
 365             !     for n=nacc, multiply v2nyy by 1.0e6 to effectively turn off the
 366             !         adjustment when number is too big (size is too small)
 367             !     for n=nait, divide   v2nxx by 1.0e6 to effectively turn off the
 368             !         adjustment when number is too small (size is too big)
 369             !OLD  however, do not change the v2nyyrl/v2nxxrl so that
 370             !OLD      the interstitial<-->activated adjustment is not changed
 371             !NEW  also change the v2nyyrl/v2nxxrl so that
 372             !NEW      the interstitial<-->activated adjustment is turned off
 373             !
 374          end if
 375          frelaxadj = 27.0
 376          dumfac = exp(4.5*alnsg_amode(n)**2)*pi/6.0
 377          v2nxx = voltonumbhi_amode(n)
 378          v2nyy = voltonumblo_amode(n)
 379          v2nxxrl = v2nxx/frelaxadj
 380          v2nyyrl = v2nyy*frelaxadj
 381          dgnxx = dgnumhi_amode(n)
 382          dgnyy = dgnumlo_amode(n)
 383          if ( do_aitacc_transfer ) then
 384             if (n == nait) v2nxx = v2nxx/1.0e6
 385             if (n == nacc) v2nyy = v2nyy*1.0e6
 386             v2nxxrl = v2nxx/frelaxadj   ! NEW
 387             v2nyyrl = v2nyy*frelaxadj   ! NEW
 388          end if
 389
 390          if (do_adjust) then
 391             dotend(lna) = .true.
 392             dotendqqcw(lnc) = .true.
 393          end if
 394
 395          do  k = 1, pver
 396             do  i = 1, ncol
 397
 398                drv_a = dryvol_a(i,k)
 399                num_a0 = q(i,k,lna)
 400                num_a = max( 0.0_r8, num_a0 )
 401                drv_c = dryvol_c(i,k)
 402                num_c0 = fldcw(i,k)
 403                num_c = max( 0.0_r8, num_c0 )
 404
 405                if ( do_adjust) then
 406
 407                   !
 408                   ! do number adjustment for interstitial and activated particles
 409                   !    adjustments that (1) make numbers non-negative or (2) make numbers
 410                   !       zero when volume is zero are applied over time-scale deltat
 411                   !    adjustments that bring numbers to within specified bounds are
 412                   !       applied over time-scale tadj
 413                   !
 414                   if ((drv_a <= 0.0_r8) .and. (drv_c <= 0.0_r8)) then
 415                      ! both interstitial and activated volumes are zero
 416                      ! adjust both numbers to zero
 417                      num_a = 0.0_r8
 418                      dqdt(i,k,lna) = -num_a0*deltatinv
 419                      num_c = 0.0_r8
 420                      dqqcwdt(i,k,lnc) = -num_c0*deltatinv
 421                   else if (drv_c <= 0.0_r8) then
 422                      ! activated volume is zero, so interstitial number/volume == total/combined
 423                      ! apply step 1 and 3, but skip the relaxed adjustment (step 2, see below)
 424                      num_c = 0.0_r8
 425                      dqqcwdt(i,k,lnc) = -num_c0*deltatinv
 426                      num_a1 = num_a
 427                      numbnd = max( drv_a*v2nxx, min( drv_a*v2nyy, num_a1 ) )
 428                      num_a  = num_a1 + (numbnd - num_a1)*fracadj
 429                      dqdt(i,k,lna) = (num_a - num_a0)*deltatinv
 430
 431                   else if (drv_a <= 0.0_r8) then
 432                      ! interstitial volume is zero, treat similar to above
 433                      num_a = 0.0_r8
 434                      dqdt(i,k,lna) = -num_a0*deltatinv
 435                      num_c1 = num_c
 436                      numbnd = max( drv_c*v2nxx, min( drv_c*v2nyy, num_c1 ) )
 437                      num_c  = num_c1 + (numbnd - num_c1)*fracadj
 438                      dqqcwdt(i,k,lnc) = (num_c - num_c0)*deltatinv
 439                   else
 440                      ! both volumes are positive
 441                      ! apply 3 adjustment steps
 442                      ! step1:  num_a,c0 --> num_a,c1 forces non-negative values
 443                      num_a1 = num_a
 444                      num_c1 = num_c
 445                      ! step2:  num_a,c1 --> num_a,c2 applies relaxed bounds to the interstitial
 446                      !    and activated number (individually)
 447                      !    if only only a or c changes, adjust the other in the opposite direction
 448                      !    as much as possible to conserve a+c
 449                      numbnd = max( drv_a*v2nxxrl, min( drv_a*v2nyyrl, num_a1 ) )
 450                      delnum_a2 = (numbnd - num_a1)*fracadj
 451                      num_a2 = num_a1 + delnum_a2
 452                      numbnd = max( drv_c*v2nxxrl, min( drv_c*v2nyyrl, num_c1 ) )
 453                      delnum_c2 = (numbnd - num_c1)*fracadj
 454                      num_c2 = num_c1 + delnum_c2
 455                      if ((delnum_a2 == 0.0) .and. (delnum_c2 /= 0.0)) then
 456                         num_a2 = max( drv_a*v2nxxrl, min( drv_a*v2nyyrl,   &
 457                              num_a1-delnum_c2 ) )
 458                      else if ((delnum_a2 /= 0.0) .and. (delnum_c2 == 0.0)) then
 459                         num_c2 = max( drv_c*v2nxxrl, min( drv_c*v2nyyrl,   &
 460                              num_c1-delnum_a2 ) )
 461                      end if
 462                      ! step3:  num_a,c2 --> num_a,c3 applies stricter bounds to the
 463                      !    combined/total number
 464                      drv_t = drv_a + drv_c
 465                      num_t2 = num_a2 + num_c2
 466                      delnum_a3 = 0.0_r8
 467                      delnum_c3 = 0.0_r8
 468                      if (num_t2 < drv_t*v2nxx) then
 469                         delnum_t3 = (drv_t*v2nxx - num_t2)*fracadj
 470                         ! if you are here then (num_a2 < drv_a*v2nxx) and/or
 471                         !                      (num_c2 < drv_c*v2nxx) must be true
 472                         if ((num_a2 < drv_a*v2nxx) .and. (num_c2 < drv_c*v2nxx)) then
 473                            delnum_a3 = delnum_t3*(num_a2/num_t2)
 474                            delnum_c3 = delnum_t3*(num_c2/num_t2)
 475                         else if (num_c2 < drv_c*v2nxx) then
 476                            delnum_c3 = delnum_t3
 477                         else if (num_a2 < drv_a*v2nxx) then
 478                            delnum_a3 = delnum_t3
 479                         end if
 480                      else if (num_t2 > drv_t*v2nyy) then
 481                         delnum_t3 = (drv_t*v2nyy - num_t2)*fracadj
 482                         ! if you are here then (num_a2 > drv_a*v2nyy) and/or
 483                         !                      (num_c2 > drv_c*v2nyy) must be true
 484                         if ((num_a2 > drv_a*v2nyy) .and. (num_c2 > drv_c*v2nyy)) then
 485                            delnum_a3 = delnum_t3*(num_a2/num_t2)
 486                            delnum_c3 = delnum_t3*(num_c2/num_t2)
 487                         else if (num_c2 > drv_c*v2nyy) then
 488                            delnum_c3 = delnum_t3
 489                         else if (num_a2 > drv_a*v2nyy) then
 490                            delnum_a3 = delnum_t3
 491                         end if
 492                      end if
 493                      num_a = num_a2 + delnum_a3
 494                      dqdt(i,k,lna) = (num_a - num_a0)*deltatinv
 495                      num_c = num_c2 + delnum_c3
 496                      dqqcwdt(i,k,lnc) = (num_c - num_c0)*deltatinv
 497                   end if
 498
 499                end if ! do_adjust
 500
 501                !
 502                ! now compute current dgn and v2n
 503                !
 504                if (drv_a > 0.0_r8) then
 505                   if (num_a <= drv_a*v2nxx) then
 506                      dgncur_a(i,k,n) = dgnxx
 507                      v2ncur_a(i,k,n) = v2nxx
 508                   else if (num_a >= drv_a*v2nyy) then
 509                      dgncur_a(i,k,n) = dgnyy
 510                      v2ncur_a(i,k,n) = v2nyy
 511                   else
 512                      dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third
 513                      v2ncur_a(i,k,n) = num_a/drv_a
 514                   end if
 515                end if
 516                pdel_fac = pdel(i,k)/gravit   ! = rho*dz
 517                jac = 1
 518                qsrflx(i,lna,1,jac) = qsrflx(i,lna,1,jac) + max(0.0_r8,dqdt(i,k,lna))*pdel_fac
 519                qsrflx(i,lna,2,jac) = qsrflx(i,lna,2,jac) + min(0.0_r8,dqdt(i,k,lna))*pdel_fac
 520
 521                if (drv_c > 0.0_r8) then
 522                   if (num_c <= drv_c*v2nxx) then
 523                      dgncur_c(i,k,n) = dgnumhi_amode(n)
 524                      v2ncur_c(i,k,n) = v2nxx
 525                   else if (num_c >= drv_c*v2nyy) then
 526                      dgncur_c(i,k,n) = dgnumlo_amode(n)
 527                      v2ncur_c(i,k,n) = v2nyy
 528                   else
 529                      dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third
 530                      v2ncur_c(i,k,n) = num_c/drv_c
 531                   end if
 532                end if
 533                jac = 2
 534                qsrflx(i,lnc,1,jac) = qsrflx(i,lnc,1,jac) + max(0.0_r8,dqqcwdt(i,k,lnc))*pdel_fac
 535                qsrflx(i,lnc,2,jac) = qsrflx(i,lnc,2,jac) + min(0.0_r8,dqqcwdt(i,k,lnc))*pdel_fac
 536
 537
 538                ! save number and dryvol for aitken <--> accum renaming
 539                if ( do_aitacc_transfer ) then
 540                   if (n == nait) then
 541                      drv_a_aitsv(i,k) = drv_a
 542                      num_a_aitsv(i,k) = num_a
 543                      drv_c_aitsv(i,k) = drv_c
 544                      num_c_aitsv(i,k) = num_c
 545                   else if (n == nacc) then
 546                      drv_a_accsv(i,k) = drv_a
 547                      num_a_accsv(i,k) = num_a
 548                      drv_c_accsv(i,k) = drv_c
 549                      num_c_accsv(i,k) = num_c
 550                   end if
 551                end if
 552                drv_a_sv(i,k,n) = drv_a
 553                num_a_sv(i,k,n) = num_a
 554                drv_c_sv(i,k,n) = drv_c
 555                num_c_sv(i,k,n) = num_c
 556
 557             end do
 558          end do
 559
 560
 561          !
 562          ! option 3 -- number and surface prognosed (variable dgnum and sigmag)
 563          !             this is not implemented
 564          !
 565       end do
 566
 567
 568       ! "do n = 1, ntot_amode"
 569
 570
 571       !
 572       !
 573       ! the following section (from here to label 49000)
 574       !    does aitken <--> accum mode transfer
 575       !
 576       ! when the aitken mode mean size is too big, the largest
 577       !    aitken particles are transferred into the accum mode
 578       !    to reduce the aitken mode mean size
 579       ! when the accum mode mean size is too small, the smallest
 580       !    accum particles are transferred into the aitken mode
 581       !    to increase the accum mode mean size
 582       !
 583       !
 584       ixfer_ait2acc_sv(:,:) = 0
 585       ixfer_acc2ait_sv(:,:) = 0
 586       if ( do_aitacc_transfer ) then
 587
 588          ! old - on time first step, npair_renamexf will be <= 0,
 589          !       in which case need to do modal_aero_rename_init
 590          ! new - init is now done through chem_init and things below it
 591          if (npair_renamexf .le. 0) then
 592             npair_renamexf = 0
 593             !        call modal_aero_rename_init
 594             if (npair_renamexf .le. 0) then
 595                write( 6, '(//a//)' )   &
 596                     '*** modal_aero_calcaersize_sub error -- npair_renamexf <= 0'
 597                call endrun( 'modal_aero_calcaersize_sub error' )
 598             end if
 599          end if
 600
 601          ! check that renaming ipair=1 is aitken-->accum
 602          ipair = 1
 603          if ((modefrm_renamexf(ipair) .ne. nait) .or.   &
 604               (modetoo_renamexf(ipair) .ne. nacc)) then
 605             write( 6, '(//2a//)' )   &
 606                  '*** modal_aero_calcaersize_sub error -- ',   &
 607                  'modefrm/too_renamexf(1) are wrong'
 608             call endrun( 'modal_aero_calcaersize_sub error' )
 609          end if
 610
 611          ! set dotend() for species that will be transferred
 612          do iq = 1, nspecfrm_renamexf(ipair)
 613             lsfrm = lspecfrma_renamexf(iq,ipair) - loffset
 614             lstoo = lspectooa_renamexf(iq,ipair) - loffset
 615             if ((lsfrm > 0) .and. (lstoo > 0)) then
 616                dotend(lsfrm) = .true.
 617                dotend(lstoo) = .true.
 618             end if
 619             lsfrm = lspecfrmc_renamexf(iq,ipair) - loffset
 620             lstoo = lspectooc_renamexf(iq,ipair) - loffset
 621             if ((lsfrm > 0) .and. (lstoo > 0)) then
 622                dotendqqcw(lsfrm) = .true.
 623                dotendqqcw(lstoo) = .true.
 624             end if
 625          end do
 626
 627          ! identify accum species cannot be transferred to aitken mode
 628          noxf_acc2ait(:) = .true.
 629          do l1 = 1, nspec_amode(nacc)
 630             la = lmassptr_amode(l1,nacc)
 631             do iq = 1, nspecfrm_renamexf(ipair)
 632                if (lspectooa_renamexf(iq,ipair) == la) then
 633                   noxf_acc2ait(l1) = .false.
 634                end if
 635             end do
 636          end do
 637
 638          ! v2nzz is voltonumb at the "geometrically-defined" mid-point
 639          ! between the aitken and accum modes
 640          v2nzz = sqrt(voltonumb_amode(nait)*voltonumb_amode(nacc))
 641
 642          ! loop over columns and levels
 643          do  k = 1, pver
 644             do  i = 1, ncol
 645
 646                pdel_fac = pdel(i,k)/gravit   ! = rho*dz
 647                xfertend_num(:,:) = 0.0_r8
 648
 649                ! compute aitken --> accum transfer rates
 650                ixfer_ait2acc = 0
 651                xfercoef_num_ait2acc = 0.0_r8
 652                xfercoef_vol_ait2acc = 0.0_r8
 653
 654                drv_t = drv_a_aitsv(i,k) + drv_c_aitsv(i,k)
 655                num_t = num_a_aitsv(i,k) + num_c_aitsv(i,k)
 656                if (drv_t > 0.0_r8) then
 657                   if (num_t < drv_t*v2nzz) then
 658                      ixfer_ait2acc = 1
 659                      if (num_t < drv_t*voltonumb_amode(nacc)) then
 660                         xferfrac_num_ait2acc = 1.0_r8
 661                         xferfrac_vol_ait2acc = 1.0_r8
 662                      else
 663                         xferfrac_vol_ait2acc = ((num_t/drv_t) - v2nzz)/   &
 664                              (voltonumb_amode(nacc) - v2nzz)
 665                         xferfrac_num_ait2acc = xferfrac_vol_ait2acc*   &
 666                              (drv_t*voltonumb_amode(nacc)/num_t)
 667                         if ((xferfrac_num_ait2acc <= 0.0_r8) .or.   &
 668                              (xferfrac_vol_ait2acc <= 0.0_r8)) then
 669                            xferfrac_num_ait2acc = 0.0_r8
 670                            xferfrac_vol_ait2acc = 0.0_r8
 671                         else if ((xferfrac_num_ait2acc >= 1.0_r8) .or.   &
 672                              (xferfrac_vol_ait2acc >= 1.0_r8)) then
 673                            xferfrac_num_ait2acc = 1.0_r8
 674                            xferfrac_vol_ait2acc = 1.0_r8
 675                         end if
 676                      end if
 677                      xfercoef_num_ait2acc = xferfrac_num_ait2acc*tadjinv
 678                      xfercoef_vol_ait2acc = xferfrac_vol_ait2acc*tadjinv
 679                      xfertend_num(1,1) = num_a_aitsv(i,k)*xfercoef_num_ait2acc
 680                      xfertend_num(1,2) = num_c_aitsv(i,k)*xfercoef_num_ait2acc
 681                   end if
 682                end if
 683
 684                ! compute accum --> aitken transfer rates
 685                ! accum may have some species (seasalt, dust, poa, lll) that are
 686                !    not in aitken mode
 687                ! so first divide the accum drv & num into not-transferred (noxf) species
 688                !    and transferred species, and use the transferred-species
 689                !    portion in what follows
 690                ixfer_acc2ait = 0
 691                xfercoef_num_acc2ait = 0.0_r8
 692                xfercoef_vol_acc2ait = 0.0_r8
 693
 694                drv_t = drv_a_accsv(i,k) + drv_c_accsv(i,k)
 695                num_t = num_a_accsv(i,k) + num_c_accsv(i,k)
 696                drv_a_noxf = 0.0_r8
 697                drv_c_noxf = 0.0_r8
 698                if (drv_t > 0.0_r8) then
 699                   if (num_t > drv_t*v2nzz) then
 700                      do l1 = 1, nspec_amode(nacc)
 701
 702                         if ( noxf_acc2ait(l1) ) then
 703                            if ( aero_mmr_flag ) then
 704                               ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
 705                               dummwdens = 1.0_r8   &
 706                                    / specdens_amode(lspectype_amode(l1,nacc))
 707                            else
 708                               ! need qmass*dummwdens = (kmol/kmol-air) * [(kg/kmol)/(kg/m3)] = m3/mol-air
 709                               dummwdens = specmw_amode(lspectype_amode(l1,nacc))   &
 710                                    / specdens_amode(lspectype_amode(l1,nacc))
 711                            end if
 712                            la = lmassptr_amode(l1,nacc) - loffset
 713                            drv_a_noxf = drv_a_noxf    &
 714                                 + max(0.0_r8,q(i,k,la))*dummwdens
 715                            lc = lmassptrcw_amode(l1,nacc) - loffset
 716 #ifndef WRF_PORT
 717                            fldcw => qqcw_get_field(lmassptrcw_amode(l1,nacc),lchnk)
 718 #else
 719                             fldcw => qqcw(:,:,lmassptrcw_amode(l1,nacc))
 720 #endif
 721                            drv_c_noxf = drv_c_noxf    &
 722                                 + max(0.0_r8,fldcw(i,k))*dummwdens
 723                         end if
 724                      end do
 725                      drv_t_noxf = drv_a_noxf + drv_c_noxf
 726                      num_t_noxf = drv_t_noxf*voltonumblo_amode(nacc)
 727                      num_t0 = num_t
 728                      drv_t0 = drv_t
 729                      num_t = max( 0.0_r8, num_t - num_t_noxf )
 730                      drv_t = max( 0.0_r8, drv_t - drv_t_noxf )
 731                   end if
 732                end if
 733
 734                if (drv_t > 0.0_r8) then
 735                   if (num_t > drv_t*v2nzz) then
 736                      ixfer_acc2ait = 1
 737                      if (num_t > drv_t*voltonumb_amode(nait)) then
 738                         xferfrac_num_acc2ait = 1.0
 739                         xferfrac_vol_acc2ait = 1.0
 740                      else
 741                         xferfrac_vol_acc2ait = ((num_t/drv_t) - v2nzz)/   &
 742                              (voltonumb_amode(nait) - v2nzz)
 743                         xferfrac_num_acc2ait = xferfrac_vol_acc2ait*   &
 744                              (drv_t*voltonumb_amode(nait)/num_t)
 745                         if ((xferfrac_num_acc2ait <= 0.0_r8) .or.   &
 746                              (xferfrac_vol_acc2ait <= 0.0_r8)) then
 747                            xferfrac_num_acc2ait = 0.0_r8
 748                            xferfrac_vol_acc2ait = 0.0_r8
 749                         else if ((xferfrac_num_acc2ait >= 1.0_r8) .or.   &
 750                              (xferfrac_vol_acc2ait >= 1.0_r8)) then
 751                            xferfrac_num_acc2ait = 1.0_r8
 752                            xferfrac_vol_acc2ait = 1.0_r8
 753                         end if
 754                      end if
 755                      duma = 1.0e-37
 756                      xferfrac_num_acc2ait = xferfrac_num_acc2ait*   &
 757                           num_t/max( duma, num_t0 )
 758                      xfercoef_num_acc2ait = xferfrac_num_acc2ait*tadjinv
 759                      xfercoef_vol_acc2ait = xferfrac_vol_acc2ait*tadjinv
 760                      xfertend_num(2,1) = num_a_accsv(i,k)*xfercoef_num_acc2ait
 761                      xfertend_num(2,2) = num_c_accsv(i,k)*xfercoef_num_acc2ait
 762                   end if
 763                end if
 764
 765                ! jump to end-of-loop if no transfer is needed at current i,k
 766                if (ixfer_ait2acc+ixfer_acc2ait > 0) then
 767                   ixfer_ait2acc_sv(i,k) = ixfer_ait2acc
 768                   ixfer_acc2ait_sv(i,k) = ixfer_acc2ait
 769
 770                   !
 771                   ! compute new dgncur & v2ncur for aitken & accum modes
 772                   !
 773                   ! currently inactive
 774                   do n = nait, nacc, (nacc-nait)
 775                      if (n .eq. nait) then
 776                         duma = (xfertend_num(1,1) - xfertend_num(2,1))*deltat
 777                         num_a     = max( 0.0_r8, num_a_aitsv(i,k) - duma )
 778                         num_a_acc = max( 0.0_r8, num_a_accsv(i,k) + duma )
 779                         duma = (drv_a_aitsv(i,k)*xfercoef_vol_ait2acc -   &
 780                              (drv_a_accsv(i,k)-drv_a_noxf)*xfercoef_vol_acc2ait)*deltat
 781                         drv_a     = max( 0.0_r8, drv_a_aitsv(i,k) - duma )
 782                         drv_a_acc = max( 0.0_r8, drv_a_accsv(i,k) + duma )
 783                         duma = (xfertend_num(1,2) - xfertend_num(2,2))*deltat
 784                         num_c     = max( 0.0_r8, num_c_aitsv(i,k) - duma )
 785                         num_c_acc = max( 0.0_r8, num_c_accsv(i,k) + duma )
 786                         duma = (drv_c_aitsv(i,k)*xfercoef_vol_ait2acc -   &
 787                              (drv_c_accsv(i,k)-drv_c_noxf)*xfercoef_vol_acc2ait)*deltat
 788                         drv_c     = max( 0.0_r8, drv_c_aitsv(i,k) - duma )
 789                         drv_c_acc = max( 0.0_r8, drv_c_accsv(i,k) + duma )
 790                      else
 791                         num_a = num_a_acc
 792                         drv_a = drv_a_acc
 793                         num_c = num_c_acc
 794                         drv_c = drv_c_acc
 795                      end if
 796
 797                      if (drv_a > 0.0_r8) then
 798                         if (num_a <= drv_a*voltonumbhi_amode(n)) then
 799                            dgncur_a(i,k,n) = dgnumhi_amode(n)
 800                            v2ncur_a(i,k,n) = voltonumbhi_amode(n)
 801                         else if (num_a >= drv_a*voltonumblo_amode(n)) then
 802                            dgncur_a(i,k,n) = dgnumlo_amode(n)
 803                            v2ncur_a(i,k,n) = voltonumblo_amode(n)
 804                         else
 805                            dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third
 806                            v2ncur_a(i,k,n) = num_a/drv_a
 807                         end if
 808                      else
 809                         dgncur_a(i,k,n) = dgnum_amode(n)
 810                         v2ncur_a(i,k,n) = voltonumb_amode(n)
 811                      end if
 812
 813                      if (drv_c > 0.0_r8) then
 814                         if (num_c <= drv_c*voltonumbhi_amode(n)) then
 815                            dgncur_c(i,k,n) = dgnumhi_amode(n)
 816                            v2ncur_c(i,k,n) = voltonumbhi_amode(n)
 817                         else if (num_c >= drv_c*voltonumblo_amode(n)) then
 818                            dgncur_c(i,k,n) = dgnumlo_amode(n)
 819                            v2ncur_c(i,k,n) = voltonumblo_amode(n)
 820                         else
 821                            dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third
 822                            v2ncur_c(i,k,n) = num_c/drv_c
 823                         end if
 824                      else
 825                         dgncur_c(i,k,n) = dgnum_amode(n)
 826                         v2ncur_c(i,k,n) = voltonumb_amode(n)
 827                      end if
 828
 829                   end do
 830
 831
 832                   !
 833                   ! compute tendency amounts for aitken <--> accum transfer
 834                   !
 835
 836                   if ( masterproc ) then
 837                      if (idiagaa > 0) then
 838                         do j = 1, 2
 839                            do iq = 1, nspecfrm_renamexf(ipair)
 840                               do jac = 1, 2
 841                                  if (j .eq. 1) then
 842                                     if (jac .eq. 1) then
 843                                        lsfrm = lspecfrma_renamexf(iq,ipair)
 844                                        lstoo = lspectooa_renamexf(iq,ipair)
 845                                     else
 846                                        lsfrm = lspecfrmc_renamexf(iq,ipair)
 847                                        lstoo = lspectooc_renamexf(iq,ipair)
 848                                     end if
 849                                  else
 850                                     if (jac .eq. 1) then
 851                                        lsfrm = lspectooa_renamexf(iq,ipair)
 852                                        lstoo = lspecfrma_renamexf(iq,ipair)
 853                                     else
 854                                        lsfrm = lspectooc_renamexf(iq,ipair)
 855                                        lstoo = lspecfrmc_renamexf(iq,ipair)
 856                                     end if
 857                                  end if
 858                                  write( 6, '(a,3i3,2i4)' ) 'calcsize j,iq,jac, lsfrm,lstoo',   &
 859                                       j,iq,jac, lsfrm,lstoo
 860                               end do
 861                            end do
 862                         end do
 863                      end if
 864                   end if
 865                   idiagaa = -1
 866
 867
 868                   ! j=1 does aitken-->accum; j=2 does accum-->aitken
 869                   do  j = 1, 2
 870
 871                      if ((j .eq. 1 .and. ixfer_ait2acc > 0) .or. &
 872                           (j .eq. 2 .and. ixfer_acc2ait > 0)) then
 873
 874                         jsrflx = j+2
 875                         if (j .eq. 1) then
 876                            xfercoef = xfercoef_vol_ait2acc
 877                         else
 878                            xfercoef = xfercoef_vol_acc2ait
 879                         end if
 880
 881                         do  iq = 1, nspecfrm_renamexf(ipair)
 882
 883                            ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
 884                            do  jac = 1, 2
 885
 886                               ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species
 887                               ! the lspectooa_renamexf (and lspectooc_renamexf) are accum  species
 888                               ! for j=1, want lsfrm=aitken species, lstoo=accum  species
 889                               ! for j=2, want lsfrm=accum  species,  lstoo=aitken species
 890                               if (j .eq. 1) then
 891                                  if (jac .eq. 1) then
 892                                     lsfrm = lspecfrma_renamexf(iq,ipair)
 893                                     lstoo = lspectooa_renamexf(iq,ipair)
 894                                  else
 895                                     lsfrm = lspecfrmc_renamexf(iq,ipair)
 896                                     lstoo = lspectooc_renamexf(iq,ipair)
 897                                  end if
 898                               else
 899                                  if (jac .eq. 1) then
 900                                     lsfrm = lspectooa_renamexf(iq,ipair)
 901                                     lstoo = lspecfrma_renamexf(iq,ipair)
 902                                  else
 903                                     lsfrm = lspectooc_renamexf(iq,ipair)
 904                                     lstoo = lspecfrmc_renamexf(iq,ipair)
 905                                  end if
 906                               end if
 907
 908                               lsfrm = lsfrm - loffset
 909
 910                               lstoo = lstoo - loffset
 911                               if ((lsfrm > 0) .and. (lstoo > 0)) then
 912                                  if (jac .eq. 1) then
 913                                     if (iq .eq. 1) then
 914                                        xfertend = xfertend_num(j,jac)
 915                                     else
 916                                        xfertend = max(0.0_r8,q(i,k,lsfrm))*xfercoef
 917                                     end if
 918                                     dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xfertend
 919                                     dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xfertend
 920                                  else
 921                                     if (iq .eq. 1) then
 922                                        xfertend = xfertend_num(j,jac)
 923                                     else
 924 #ifndef WRF_PORT
 925                                        fldcw => qqcw_get_field(lsfrm+loffset,lchnk)
 926 #else
 927                                        fldcw => qqcw(:,:,lsfrm+loffset)
 928 #endif
 929                                        xfertend = max(0.0_r8,fldcw(i,k))*xfercoef
 930                                     end if
 931                                     dqqcwdt(i,k,lsfrm) = dqqcwdt(i,k,lsfrm) - xfertend
 932                                     dqqcwdt(i,k,lstoo) = dqqcwdt(i,k,lstoo) + xfertend
 933                                  end if
 934                                  qsrflx(i,lsfrm,jsrflx,jac) = qsrflx(i,lsfrm,jsrflx,jac) - xfertend*pdel_fac
 935                                  qsrflx(i,lstoo,jsrflx,jac) = qsrflx(i,lstoo,jsrflx,jac) + xfertend*pdel_fac
 936                               end if
 937
 938                            end do
 939                         end do
 940                      end if
 941                   end do
 942
 943                end if
 944             end do
 945          end do
 946
 947
 948       end if  !  do_aitacc_transfer
 949       lsfrm = -123456789   ! executable statement for debugging
 950
 951 #ifdef MODAL_DIAGNOSTICS
 952       !+++ start of diagnostics for testing ++++++++++++++++++++++++++++++++++
 953       if (nstep > 4) call endrun( 'modal_aero_calcsize_sub TESTING HALT' )
 954
 955       call get_lat_all_p(  lchnk, ncol, lat )
 956       call get_lon_all_p(  lchnk, ncol, lon )
 957       icol_diag = -999888777
 958       do i = 1, ncol
 959          !        if ((lat(i) == plat/2) .and. (lon(i) == plon/2)) icol_diag = i
 960          if ((lat(i) ==     40) .and. (lon(i) ==     41)) icol_diag = i
 961       end do
 962       if (icol_diag <= 0) goto 59000
 963
 964       do l1 = 1, nspec_amode(nacc)
 965          la = lmassptr_amode(l1,nacc)
 966          write( *, '(a,2i5,l5,4x,a)' ) 'calcsize l1, la, noxf', &
 967               l1, la, noxf_acc2ait(l1), cnst_name(la)
 968       end do
 969
 970       i = icol_diag ; k = 8
 971       do k = 5, 21, 16
 972          write( *, '(/a,i7,f10.2)' )   &
 973               '*** modal_aero_calcsize_sub -- nstep, deltat =', nstep, deltat
 974          write( *, '( a,3i7)' )   &
 975               '*** modal_aero_calcsize_sub -- icol_diag, lat(), lon()  ', i, lat(i), lon(i)
 976          do n = 1, ntot_amode
 977             lna = numptr_amode(n) - loffset
 978             write( *, '(/a,4i11)'      ) 'mode, mprognum; i, k         ', &
 979                  n, mprognum_amode(n), i, k
 980             if ((n == nait) .or. (n == nacc)) then
 981                write( *, '(a,2i11,l11)'   ) 'ixfer_ait2acc, acc2ait, do_  ', &
 982                     ixfer_ait2acc_sv(i,k), ixfer_acc2ait_sv(i,k), do_aitacc_transfer
 983             end if
 984             write( *, '( a,1p,5e11.3)' ) 'num_a_sv, n, dndt, n+del     ',    &
 985                  num_a_sv(i,k,n), q(i,k,lna), dqdt(i,k,lna),   &
 986                  (q(i,k,lna) + deltat*dqdt(i,k,lna))
 987             duma = max( drv_a_sv(i,k,n), 1.0d-35 )
 988             dumb = max( q(i,k,lna), 1.0d-35 )
 989             duma = duma**third
 990             dumb = dumb**third
 991             if (duma < dumb*1.0e20) then
 992                dumc = ((6.0/pi)**third)*duma/dumb
 993                dumd = dumc / exp( 1.5*(alnsg_amode(n)**2) )
 994             else
 995                dumc = 1.0e20
 996                dumd = 1.0e20
 997             end if
 998             write( *, '( a,1p,5e11.3)' ) 'dgnum, dgncur, dumd/c, drv_sv',    &
 999                  dgnum_amode(n), dgncur_a(i,k,n), dumd, dumc, drv_a_sv(i,k,n)
1000             write( *, '( a )' ) 'q(lmassptr(...)) // dqdt(...) // q+del'
1001             write( *, '( 1p,10e10.2 )' )   &
1002                  ( q(i,k,lmassptr_amode(l1,n)-loffset), l1=1,nspec_amode(n) )
1003             write( *, '( 1p,10e10.2 )' )   &
1004                  ( dqdt(i,k,lmassptr_amode(l1,n)-loffset), l1=1,nspec_amode(n) )
1005             write( *, '( 1p,10e10.2 )' )   &
1006                  ( ( q(i,k,lmassptr_amode(l1,n)-loffset)   &
1007                  + deltat*dqdt(i,k,lmassptr_amode(l1,n)-loffset) ),   &
1008                  l1=1,nspec_amode(n) )
1009
1010             lnc = numptrcw_amode(n) - loffset
1011 #ifndef WRF_PORT
1012             fldcw => qqcw_get_field(numptrcw_amode(n),lchnk)
1013 #else
1014             fldcw => qqcw(:,:,numptrcw_amode(n))
1015 #endif
1016             write( *, '(/a,1p,5e11.3)' ) 'num_c_sv, n, dndt, n+del     ',    &
1017                  num_c_sv(i,k,n), fldcw(i,k), dqqcwdt(i,k,lnc),   &
1018                  (fldcw(i,k) + deltat*dqqcwdt(i,k,lnc))
1019             duma = max( drv_c_sv(i,k,n), 1.0d-35 )
1020             dumb = max( fldcw(i,k), 1.0d-35 )
1021             duma = duma**third
1022             dumb = dumb**third
1023             if (duma < dumb*1.0e20) then
1024                dumc = ((6.0/pi)**third)*duma/dumb
1025                dumd = dumc / exp( 1.5*(alnsg_amode(n)**2) )
1026             else
1027                dumc = 1.0e20
1028                dumd = 1.0e20
1029             end if
1030             write( *, '( a,1p,5e11.3)' ) 'dgnum, dgncur, dumd/c, drv_sv',    &
1031                  dgnum_amode(n), dgncur_c(i,k,n), dumd, dumc, drv_c_sv(i,k,n)
1032             write( *, '( a )' ) 'qqcw(lmassptrcw(...)) // dqqcwdt(...) // q+del'
1033             !          write( *, '( 1p,10e10.2 )' )   &
1034             !             ( qqcw(i,k,lmassptrcw_amode(l1,n)-loffset), l1=1,nspec_amode(n) )
1035             write( *, '( 1p,10e10.2 )' )   &
1036                  ( dqqcwdt(i,k,lmassptrcw_amode(l1,n)-loffset), l1=1,nspec_amode(n) )
1037             !          write( *, '( 1p,10e10.2 )' )   &
1038             !             ( ( qqcw(i,k,lmassptrcw_amode(l1,n)-loffset)   &
1039             !               + deltat*dqqcwdt(i,k,lmassptrcw_amode(l1,n)-loffset) ),   &
1040             !               l1=1,nspec_amode(n) )
1041
1042          end do
1043       end do
1044 #endif
1045       !+++++ end of diagnostics for testing ++++++++++++++++++++++++++++++++++
1046
1047       !
1048       ! apply tendencies to cloud-borne species MRs
1049       !
1050       do l = 1, pcnst
1051          lc = l - loffset
1052          if ( lc>0 .and. dotendqqcw(lc) ) then
1053 #ifndef WRF_PORT
1054             fldcw=> qqcw_get_field(l,lchnk)
1055 #else
1056             fldcw=> qqcw(:,:,l)
1057 #endif
1058             do k = 1, pver
1059                do i = 1, ncol
1060                   fldcw(i,k) = max( 0.0_r8,   &
1061                        (fldcw(i,k) + dqqcwdt(i,k,lc)*deltat) )
1062                end do
1063             end do
1064          end if
1065       end do
1066
1067       !
1068       ! do outfld calls
1069       !
1070
1071       ! history fields for number-adjust source-sink for all modes
1072       if ( .not. do_adjust ) return
1073
1074       do n = 1, ntot_amode
1075          if (mprognum_amode(n) <= 0) cycle
1076
1077          do jac = 1, 2
1078             if (jac == 1) then
1079                l = numptr_amode(n)
1080                tmpnamea = cnst_name(l)
1081             else
1082                l = numptrcw_amode(n)
1083                tmpnamea = cnst_name_cw(l)
1084             end if
1085             l = l - loffset
1086             fieldname = trim(tmpnamea) // '_sfcsiz1'
1087             call outfld( fieldname, qsrflx(:,l,1,jac), pcols, lchnk)
1088
1089             fieldname = trim(tmpnamea) // '_sfcsiz2'
1090             call outfld( fieldname, qsrflx(:,l,2,jac), pcols, lchnk)
1091          end do   ! jac = ...
1092
1093       end do   ! n = ...
1094
1095
1096       ! history fields for aitken-accum transfer
1097       if ( .not. do_aitacc_transfer ) return
1098
1099       do iq = 1, nspecfrm_renamexf(ipair)
1100
1101          ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
1102          do jac = 1, 2
1103
1104             ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species
1105             ! the lspectooa_renamexf (and lspectooc_renamexf) are accum  species
1106             if (jac .eq. 1) then
1107                lsfrm = lspecfrma_renamexf(iq,ipair)
1108                lstoo = lspectooa_renamexf(iq,ipair)
1109             else
1110                lsfrm = lspecfrmc_renamexf(iq,ipair)
1111                lstoo = lspectooc_renamexf(iq,ipair)
1112             end if
1113             if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
1114
1115             if (jac .eq. 1) then
1116                tmpnamea = cnst_name(lsfrm)
1117                tmpnameb = cnst_name(lstoo)
1118             else
1119                tmpnamea = cnst_name_cw(lsfrm)
1120                tmpnameb = cnst_name_cw(lstoo)
1121             end if
1122             lsfrm = lsfrm - loffset
1123             lstoo = lstoo - loffset
1124             if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
1125
1126             fieldname = trim(tmpnamea) // '_sfcsiz3'
1127             call outfld( fieldname, qsrflx(:,lsfrm,3,jac), pcols, lchnk)
1128
1129             fieldname = trim(tmpnameb) // '_sfcsiz3'
1130             call outfld( fieldname, qsrflx(:,lstoo,3,jac), pcols, lchnk)
1131
1132             fieldname = trim(tmpnamea) // '_sfcsiz4'
1133             call outfld( fieldname, qsrflx(:,lsfrm,4,jac), pcols, lchnk)
1134
1135             fieldname = trim(tmpnameb) // '_sfcsiz4'
1136             call outfld( fieldname, qsrflx(:,lstoo,4,jac), pcols, lchnk)
1137
1138          end do   ! jac = ...
1139       end do   ! iq = ...
1140     end subroutine modal_aero_calcsize_sub
1141
1142
1143 !----------------------------------------------------------------------
1144
1145
1146 subroutine modal_aero_calcsize_init
1147
1148 !-----------------------------------------------------------------------
1149 !
1150 ! Purpose:
1151 !    set do_adjust and do_aitken flags
1152 !    create history fields for column tendencies associated with
1153 !       modal_aero_calcsize
1154 !
1155 ! Author: R. Easter
1156 !
1157 !-----------------------------------------------------------------------
1158
1159 use modal_aero_data
1160 use modal_aero_rename
1161 #ifndef WRF_PORT
1162 use abortutils,    only :    endrun
1163 use cam_history,   only :   addfld, add_default, fieldname_len, phys_decomp
1164 use constituents,  only :  pcnst, cnst_name
1165 use spmd_utils,    only :    masterproc
1166 use phys_control,  only : phys_getopts
1167 #else
1168 use module_cam_support, only: endrun, addfld, add_default, fieldname_len, phys_decomp, &
1169      pcnst =>pcnst_runtime, masterproc
1170 use constituents,  only :  cnst_name
1171 #endif
1172
1173
1174 implicit none
1175
1176 !-----------------------------------------------------------------------
1177 ! arguments
1178
1179 !-----------------------------------------------------------------------
1180 ! local
1181    integer  :: ipair, iq
1182    integer  :: jac
1183    integer  :: lsfrm, lstoo
1184    integer  :: n, nacc, nait
1185
1186    character(len=fieldname_len)   :: tmpnamea, tmpnameb
1187    character(len=fieldname_len+3) :: fieldname
1188    character(128)                 :: long_name
1189    character(8)                   :: unit
1190
1191    logical                        :: history_aerosol      ! Output the MAM aerosol tendencies
1192
1193    !-----------------------------------------------------------------------
1194 #ifndef WRF_PORT
1195       call phys_getopts( history_aerosol_out        = history_aerosol   )
1196 #else
1197       history_aerosol = .FALSE.
1198 #endif
1199
1200
1201 !  do_adjust allows adjustment to be turned on/off
1202       do_adjust = .true.
1203
1204 !  do_aitacc_transfer allows aitken <--> accum mode transfer to be turned on/off
1205 !  *** it can only be true when aitken & accum modes are both present
1206 !      and have prognosed number and diagnosed surface/sigmag
1207       nait = modeptr_aitken
1208       nacc = modeptr_accum
1209       do_aitacc_transfer = .false.
1210       if ((modeptr_aitken > 0) .and.   &
1211           (modeptr_accum  > 0) .and.   &
1212           (modeptr_aitken /= modeptr_accum)) then
1213          do_aitacc_transfer = .true.
1214          if (mprognum_amode(nait) <= 0) do_aitacc_transfer = .false.
1215          if (mprognum_amode(nacc) <= 0) do_aitacc_transfer = .false.
1216       end if
1217 !     do_aitacc_transfer = .false.
1218
1219
1220 !  define history fields for number-adjust source-sink for all modes
1221       if ( .not. do_adjust ) return
1222       do n = 1, ntot_amode
1223          if (mprognum_amode(n) <= 0) cycle
1224
1225          do jac = 1, 2
1226             if (jac == 1) then
1227                tmpnamea = cnst_name(numptr_amode(n))
1228             else
1229                tmpnamea = cnst_name_cw(numptrcw_amode(n))
1230             end if
1231             unit = '#/m2/s'
1232             fieldname = trim(tmpnamea) // '_sfcsiz1'
1233             long_name = trim(tmpnamea) // ' calcsize number-adjust column source'
1234             call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
1235             if ( history_aerosol ) then
1236                call add_default( fieldname, 1, ' ' )
1237             endif
1238
1239             if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
1240
1241             fieldname = trim(tmpnamea) // '_sfcsiz2'
1242             long_name = trim(tmpnamea) // ' calcsize number-adjust column sink'
1243             call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
1244             if ( history_aerosol ) then
1245                call add_default( fieldname, 1, ' ' )
1246             endif
1247             if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
1248          end do   ! jac = ...
1249
1250       end do   ! n = ...
1251
1252
1253 !  define history fields for aitken-accum transfer
1254       if ( .not. do_aitacc_transfer ) return
1255
1256 ! check that renaming ipair=1 is aitken-->accum
1257       ipair = 1
1258       if ((modefrm_renamexf(ipair) .ne. nait) .or.   &
1259           (modetoo_renamexf(ipair) .ne. nacc)) then
1260          write( 6, '(//2a//)' )   &
1261             '*** modal_aero_calcaersize_init error -- ',   &
1262             'modefrm/too_renamexf(1) are wrong'
1263          call endrun( 'modal_aero_calcaersize_init error' )
1264       end if
1265
1266       do iq = 1, nspecfrm_renamexf(ipair)
1267
1268 ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
1269          do jac = 1, 2
1270
1271 ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species
1272 ! the lspectooa_renamexf (and lspectooc_renamexf) are accum  species
1273             if (jac .eq. 1) then
1274                lsfrm = lspecfrma_renamexf(iq,ipair)
1275                lstoo = lspectooa_renamexf(iq,ipair)
1276             else
1277                lsfrm = lspecfrmc_renamexf(iq,ipair)
1278                lstoo = lspectooc_renamexf(iq,ipair)
1279             end if
1280             if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
1281
1282             if (jac .eq. 1) then
1283                tmpnamea = cnst_name(lsfrm)
1284                tmpnameb = cnst_name(lstoo)
1285             else
1286                tmpnamea = cnst_name_cw(lsfrm)
1287                tmpnameb = cnst_name_cw(lstoo)
1288             end if
1289
1290             unit = 'kg/m2/s'
1291             if ((tmpnamea(1:3) == 'num') .or. &
1292                 (tmpnamea(1:3) == 'NUM')) unit = '#/m2/s'
1293             fieldname = trim(tmpnamea) // '_sfcsiz3'
1294             long_name = trim(tmpnamea) // ' calcsize aitken-to-accum adjust column tendency'
1295             call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
1296             if ( history_aerosol ) then
1297                call add_default( fieldname, 1, ' ' )
1298             endif
1299             if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
1300
1301             fieldname = trim(tmpnameb) // '_sfcsiz3'
1302             long_name = trim(tmpnameb) // ' calcsize aitken-to-accum adjust column tendency'
1303             call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
1304             if ( history_aerosol ) then
1305                call add_default( fieldname, 1, ' ' )
1306             endif
1307             if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
1308
1309             fieldname = trim(tmpnamea) // '_sfcsiz4'
1310             long_name = trim(tmpnamea) // ' calcsize accum-to-aitken adjust column tendency'
1311             call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
1312             if ( history_aerosol ) then
1313                call add_default( fieldname, 1, ' ' )
1314             endif
1315             if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
1316
1317             fieldname = trim(tmpnameb) // '_sfcsiz4'
1318             long_name = trim(tmpnameb) // ' calcsize accum-to-aitken adjust column tendency'
1319             call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
1320             if ( history_aerosol ) then
1321                call add_default( fieldname, 1, ' ' )
1322             endif
1323             if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
1324
1325          end do   ! jac = ...
1326       end do   ! iq = ...
1327
1328
1329       return
1330       end subroutine modal_aero_calcsize_init
1331
1332
1333 !----------------------------------------------------------------------
1334
1335    end module modal_aero_calcsize
1336
1337
1338