physpkg.F90

module physpkg
  !-----------------------------------------------------------------------
  ! Purpose:
  !
  ! Provides the interface to CAM physics package
  !
  ! Revision history:
  ! Aug  2005,  E. B. Kluzek,  Creation of module from physpkg subroutine
  ! 2005-10-17  B. Eaton       Add contents of inti.F90 to phys_init().  Add
  !                            initialization of grid info in phys_state.
  ! Nov 2010    A. Gettelman   Put micro/macro physics into separate routines
  !-----------------------------------------------------------------------

  use shr_kind_mod,     only: r8 => shr_kind_r8
  use spmd_utils,       only: masterproc
  use physconst,        only: latvap, latice
  use physics_types,    only: physics_state, physics_tend, physics_state_set_grid, &
       physics_ptend, physics_tend_init, physics_update,    &
       physics_type_alloc, physics_ptend_dealloc,&
       physics_state_alloc, physics_state_dealloc, physics_tend_alloc, physics_tend_dealloc
  use phys_grid,        only: get_ncols_p
  use phys_gmean,       only: gmean_mass
  use ppgrid,           only: begchunk, endchunk, pcols, pver, pverp, psubcols
  use constituents,     only: pcnst, cnst_get_ind
  use camsrfexch,       only: cam_out_t, cam_in_t

  use cam_control_mod,  only: ideal_phys, adiabatic
  use phys_control,     only: phys_do_flux_avg, phys_getopts, waccmx_is
  use scamMod,          only: single_column, scm_crm_mode
  use flux_avg,         only: flux_avg_init
  use perf_mod
  use cam_logfile,     only: iulog
  use camsrfexch,      only: cam_export

  use phys_control,    only: use_hemco            ! Use Harmonized Emissions Component (HEMCO)

  use modal_aero_calcsize,    only: modal_aero_calcsize_init, modal_aero_calcsize_diag, modal_aero_calcsize_reg
  use modal_aero_calcsize,    only: modal_aero_calcsize_sub
  use modal_aero_wateruptake, only: modal_aero_wateruptake_init, modal_aero_wateruptake_dr, modal_aero_wateruptake_reg

  use carma_diags_mod, only: carma_diags_t

  use offline_driver,  only: offline_driver_dorun

  implicit none
  private
  save

  ! Public methods
  public phys_register ! was initindx  - register physics methods
  public phys_init   ! Public initialization method
  public phys_run1   ! First phase of the public run method
  public phys_run2   ! Second phase of the public run method
  public phys_final  ! Public finalization method

  ! Private module data

  ! Physics package options
  character(len=16) :: shallow_scheme
  character(len=16) :: macrop_scheme
  character(len=16) :: microp_scheme
  character(len=16) :: subcol_scheme
  character(len=32) :: cam_take_snapshot_before ! Physics routine to take a snapshot "before"
  character(len=32) :: cam_take_snapshot_after  ! Physics routine to take a snapshot "after"
  integer           :: cld_macmic_num_steps    ! Number of macro/micro substeps
  integer           :: cam_snapshot_before_num ! tape number for before snapshots
  integer           :: cam_snapshot_after_num  ! tape number for after snapshots
  logical           :: do_clubb_sgs
  logical           :: use_subcol_microp   ! if true, use subcolumns in microphysics
  logical           :: state_debug_checks  ! Debug physics_state.
  logical           :: clim_modal_aero     ! climate controled by prognostic or prescribed modal aerosols
  logical           :: prog_modal_aero     ! Prognostic modal aerosols present

  !  Physics buffer index
  integer ::  teout_idx          = 0

  integer ::  landm_idx          = 0
  integer ::  sgh_idx            = 0
  integer ::  sgh30_idx          = 0

  integer ::  qini_idx           = 0
  integer ::  cldliqini_idx      = 0
  integer ::  cldiceini_idx      = 0
  integer ::  totliqini_idx      = 0
  integer ::  toticeini_idx      = 0

  integer ::  prec_str_idx       = 0
  integer ::  snow_str_idx       = 0
  integer ::  prec_sed_idx       = 0
  integer ::  snow_sed_idx       = 0
  integer ::  prec_pcw_idx       = 0
  integer ::  snow_pcw_idx       = 0
  integer ::  prec_dp_idx        = 0
  integer ::  snow_dp_idx        = 0
  integer ::  prec_sh_idx        = 0
  integer ::  snow_sh_idx        = 0
  integer ::  dlfzm_idx          = 0     ! detrained convective cloud water mixing ratio.
  integer ::  ducore_idx         = 0     ! ducore index in physics buffer
  integer ::  dvcore_idx         = 0     ! dvcore index in physics buffer
  integer ::  dtcore_idx         = 0     ! dtcore index in physics buffer
  integer ::  dqcore_idx         = 0     ! dqcore index in physics buffer

!=======================================================================
contains
!=======================================================================

  subroutine phys_register
    !-----------------------------------------------------------------------
    !
    ! Purpose: Register constituents and physics buffer fields.
    !
    ! Author:    CSM Contact: M. Vertenstein, Aug. 1997
    !            B.A. Boville, Oct 2001
    !            A. Gettelman, Nov 2010 - put micro/macro physics into separate routines
    !
    !-----------------------------------------------------------------------
    use cam_abortutils,     only: endrun
    use physics_buffer,     only: pbuf_init_time, pbuf_cam_snapshot_register
    use physics_buffer,     only: pbuf_add_field, dtype_r8, pbuf_register_subcol
    use shr_kind_mod,       only: r8 => shr_kind_r8
    use constituents,       only: pcnst, cnst_add, cnst_chk_dim

    use cam_control_mod,    only: moist_physics
    use chemistry,          only: chem_register
    use mo_lightning,       only: lightning_register
    use cloud_fraction,     only: cldfrc_register
    use rk_stratiform_cam,  only: rk_stratiform_cam_register
    use microp_driver,      only: microp_driver_register
    use microp_aero,        only: microp_aero_register
    use macrop_driver,      only: macrop_driver_register
    use clubb_intr,         only: clubb_register_cam
    use conv_water,         only: conv_water_register
    use physconst,          only: mwh2o, cpwv
    use tracers,            only: tracers_register
    use check_energy,       only: check_energy_register
    use carma_intr,         only: carma_register
    use ghg_data,           only: ghg_data_register
    use vertical_diffusion, only: vd_register
    use convect_deep,       only: convect_deep_register
    use convect_shallow,    only: convect_shallow_register
    use radiation,          only: radiation_register
    use co2_cycle,          only: co2_register
    use flux_avg,           only: flux_avg_register
    use iondrag,            only: iondrag_register
    use waccmx_phys_intr,   only: waccmx_phys_ion_elec_temp_reg
    use prescribed_ozone,   only: prescribed_ozone_register
    use prescribed_volcaero,only: prescribed_volcaero_register
    use prescribed_strataero,only: prescribed_strataero_register
    use prescribed_aero,    only: prescribed_aero_register
    use prescribed_ghg,     only: prescribed_ghg_register
    use sslt_rebin,         only: sslt_rebin_register
    use aoa_tracers,        only: aoa_tracers_register
    use aircraft_emit,      only: aircraft_emit_register
    use cam_diagnostics,    only: diag_register
    use cloud_diagnostics,  only: cloud_diagnostics_register
    use cospsimulator_intr, only: cospsimulator_intr_register
    use rad_constituents,   only: rad_cnst_get_info ! Added to query if it is a modal aero sim or not
    use radheat,            only: radheat_register
    use subcol,             only: subcol_register
    use subcol_utils,       only: is_subcol_on, subcol_get_scheme
    use dyn_comp,           only: dyn_register
    use offline_driver,     only: offline_driver_reg
    use hemco_interface,    only: HCOI_Chunk_Init
    use upper_bc,           only: ubc_fixed_conc
    use surface_emissions_mod, only: surface_emissions_reg
    use elevated_emissions_mod, only: elevated_emissions_reg
    use ctem_diags_mod, only: ctem_diags_reg

    !---------------------------Local variables-----------------------------
    !
    integer  :: m        ! loop index
    integer  :: mm       ! constituent index
    integer  :: nmodes
    logical  :: has_fixed_ubc ! for upper bndy cond
    !-----------------------------------------------------------------------

    ! Get physics options
    call phys_getopts(shallow_scheme_out          = shallow_scheme, &
                      macrop_scheme_out           = macrop_scheme,   &
                      microp_scheme_out           = microp_scheme,   &
                      cld_macmic_num_steps_out    = cld_macmic_num_steps, &
                      do_clubb_sgs_out            = do_clubb_sgs,     &
                      use_subcol_microp_out       = use_subcol_microp, &
                      state_debug_checks_out      = state_debug_checks, &
                      cam_take_snapshot_before_out= cam_take_snapshot_before, &
                      cam_take_snapshot_after_out = cam_take_snapshot_after, &
                      cam_snapshot_before_num_out = cam_snapshot_before_num, &
                      cam_snapshot_after_num_out  = cam_snapshot_after_num)

    subcol_scheme = subcol_get_scheme()

    ! Initialize dyn_time_lvls
    call pbuf_init_time()

    ! Register the subcol scheme
    call subcol_register()

    ! Register water vapor.
    ! ***** N.B. ***** This must be the first call to cnst_add so that
    !                  water vapor is constituent 1.
    has_fixed_ubc = ubc_fixed_conc('Q') ! .false.
    if (moist_physics) then
       call cnst_add('Q', mwh2o, cpwv, 1.E-12_r8, mm, fixed_ubc=has_fixed_ubc, &
            longname='Specific humidity', readiv=.true., is_convtran1=.true.)
    else
       call cnst_add('Q', mwh2o, cpwv, 0.0_r8, mm, fixed_ubc=has_fixed_ubc, &
            longname='Specific humidity', readiv=.false., is_convtran1=.true.)
    end if

    ! Topography file fields.
    call pbuf_add_field('LANDM',     'global',  dtype_r8, (/pcols/),      landm_idx)
    call pbuf_add_field('SGH',       'global',  dtype_r8, (/pcols/),      sgh_idx)
    call pbuf_add_field('SGH30',     'global',  dtype_r8, (/pcols/),      sgh30_idx)

    ! Fields for physics package diagnostics
    call pbuf_add_field('QINI',      'physpkg', dtype_r8, (/pcols,pver/), qini_idx)
    call pbuf_add_field('CLDLIQINI', 'physpkg', dtype_r8, (/pcols,pver/), cldliqini_idx)
    call pbuf_add_field('CLDICEINI', 'physpkg', dtype_r8, (/pcols,pver/), cldiceini_idx)
    call pbuf_add_field('TOTLIQINI', 'physpkg', dtype_r8, (/pcols,pver/), totliqini_idx)
    call pbuf_add_field('TOTICEINI', 'physpkg', dtype_r8, (/pcols,pver/), toticeini_idx)

    ! check energy package
    call check_energy_register

    ! If using a simple physics option (e.g., held_suarez, adiabatic),
    ! the normal CAM physics parameterizations are not called.
    if (moist_physics) then

       ! register fluxes for saving across time
       if (phys_do_flux_avg()) call flux_avg_register()

       call cldfrc_register()

       ! cloud water
       if( microp_scheme == 'RK' ) then
          call rk_stratiform_cam_register()
       elseif( microp_scheme == 'MG' ) then
          if (.not. do_clubb_sgs) call macrop_driver_register()
          call microp_aero_register()
          call microp_driver_register()
       end if

       ! Register CLUBB_SGS here
       if (do_clubb_sgs) call clubb_register_cam()


       call pbuf_add_field('PREC_STR',  'physpkg',dtype_r8,(/pcols/),prec_str_idx)
       call pbuf_add_field('SNOW_STR',  'physpkg',dtype_r8,(/pcols/),snow_str_idx)
       call pbuf_add_field('PREC_PCW',  'physpkg',dtype_r8,(/pcols/),prec_pcw_idx)
       call pbuf_add_field('SNOW_PCW',  'physpkg',dtype_r8,(/pcols/),snow_pcw_idx)
       call pbuf_add_field('PREC_SED',  'physpkg',dtype_r8,(/pcols/),prec_sed_idx)
       call pbuf_add_field('SNOW_SED',  'physpkg',dtype_r8,(/pcols/),snow_sed_idx)
       if (is_subcol_on()) then
         call pbuf_register_subcol('PREC_STR', 'phys_register', prec_str_idx)
         call pbuf_register_subcol('SNOW_STR', 'phys_register', snow_str_idx)
         call pbuf_register_subcol('PREC_PCW', 'phys_register', prec_pcw_idx)
         call pbuf_register_subcol('SNOW_PCW', 'phys_register', snow_pcw_idx)
         call pbuf_register_subcol('PREC_SED', 'phys_register', prec_sed_idx)
         call pbuf_register_subcol('SNOW_SED', 'phys_register', snow_sed_idx)
       end if

    ! Who should add FRACIS?
    ! -- It does not seem that aero_intr should add it since FRACIS is used in convection
    !     even if there are no prognostic aerosols ... so do it here for now
       call pbuf_add_field('FRACIS','physpkg',dtype_r8,(/pcols,pver,pcnst/),m)

       call conv_water_register()

       ! Determine whether its a 'modal' aerosol simulation  or not
       call rad_cnst_get_info(0, nmodes=nmodes)
       clim_modal_aero = (nmodes > 0)

       if (clim_modal_aero) then
          call modal_aero_calcsize_reg()
          call modal_aero_wateruptake_reg()
       endif

       call surface_emissions_reg()
       call elevated_emissions_reg()

       ! register chemical constituents including aerosols ...
       call chem_register()

      ! add prognostic lightning flash freq pbuf fld
       call lightning_register()

       ! co2 constituents
       call co2_register()

       ! register data model ozone with pbuf
       call prescribed_volcaero_register()
       call prescribed_strataero_register()
       call prescribed_ozone_register()
       call prescribed_aero_register()
       call prescribed_ghg_register()
       call sslt_rebin_register

       ! register various data model gasses with pbuf
       call ghg_data_register()

       ! carma microphysics
       !
       call carma_register()

       ! Register iondrag variables with pbuf
       call iondrag_register()

       ! Register ionosphere variables with pbuf if mode set to ionosphere
       if( waccmx_is('ionosphere') ) then
          call waccmx_phys_ion_elec_temp_reg()
       endif

       call aircraft_emit_register()

       ! deep convection
       call convect_deep_register

       !  shallow convection
       call convect_shallow_register

       ! radiation
       call radiation_register
       call cloud_diagnostics_register
       call radheat_register

       ! COSP
       call cospsimulator_intr_register

       ! vertical diffusion
       call vd_register()
    else
       ! held_suarez/adiabatic physics option should be in simple_physics
       call endrun('phys_register: moist_physics configuration error')
    end if

    ! Register diagnostics PBUF
    call diag_register()

    ! Register age of air tracers
    call aoa_tracers_register()

    ! Register test tracers
    call tracers_register()

    call dyn_register()

    ! All tracers registered, check that the dimensions are correct
    call cnst_chk_dim()

    ! ***NOTE*** No registering constituents after the call to cnst_chk_dim.

    call offline_driver_reg()

    if (use_hemco) then
        ! initialize harmonized emissions component (HEMCO)
        call HCOI_Chunk_Init()
    endif

    ! TEM diagnostics
    call ctem_diags_reg()

    ! This needs to be last as it requires all pbuf fields to be added
    if (cam_snapshot_before_num > 0 .or. cam_snapshot_after_num > 0) then
        call pbuf_cam_snapshot_register()
    end if

  end subroutine phys_register



  !=======================================================================

  subroutine phys_inidat( cam_out, pbuf2d )
    use cam_abortutils,      only: endrun

    use physics_buffer,      only: pbuf_get_index, physics_buffer_desc, pbuf_set_field, dyn_time_lvls


    use cam_initfiles,       only: initial_file_get_id, topo_file_get_id
    use cam_grid_support,    only: cam_grid_check, cam_grid_id
    use cam_grid_support,    only: cam_grid_get_dim_names
    use pio,                 only: file_desc_t
    use ncdio_atm,           only: infld
    use dycore,              only: dycore_is
    use polar_avg,           only: polar_average
    use short_lived_species, only: initialize_short_lived_species
    use cam_control_mod,     only: aqua_planet
    use waccmx_phys_intr,    only: waccmx_phys_ion_elec_temp_inidat

    type(cam_out_t),     intent(inout) :: cam_out(begchunk:endchunk)
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
    integer :: lchnk, m, n, ncol
    type(file_desc_t), pointer :: fh_ini, fh_topo
    character(len=8) :: fieldname
    real(r8), pointer :: tptr(:,:), tptr_2(:,:), tptr3d(:,:,:), tptr3d_2(:,:,:)

    character(len=11) :: subname='phys_inidat' ! subroutine name
    integer :: tpert_idx, qpert_idx, pblh_idx

    logical :: found=.false., found2=.false.
    integer :: ierr
    character(len=8) :: dim1name, dim2name
    integer :: ixcldice, ixcldliq
    integer                   :: grid_id  ! grid ID for data mapping

    nullify(tptr,tptr_2,tptr3d,tptr3d_2)

    fh_ini  => initial_file_get_id()
    fh_topo => topo_file_get_id()

    !   dynamics variables are handled in dyn_init - here we read variables needed for physics
    !   but not dynamics

    grid_id = cam_grid_id('physgrid')
    if (.not. cam_grid_check(grid_id)) then
      call endrun(trim(subname)//': Internal error, no "physgrid" grid')
    end if
    call cam_grid_get_dim_names(grid_id, dim1name, dim2name)

    allocate(tptr(1:pcols,begchunk:endchunk))

    if (associated(fh_topo) .and. .not. aqua_planet) then
      call infld('SGH', fh_topo, dim1name, dim2name, 1, pcols, begchunk, endchunk, &
           tptr, found, gridname='physgrid')
      if(.not. found) call endrun('ERROR: SGH not found on topo file')

      call pbuf_set_field(pbuf2d, sgh_idx, tptr)

      allocate(tptr_2(1:pcols,begchunk:endchunk))
      call infld('SGH30', fh_topo, dim1name, dim2name, 1, pcols, begchunk, endchunk, &
           tptr_2, found, gridname='physgrid')
      if(found) then
        call pbuf_set_field(pbuf2d, sgh30_idx, tptr_2)
      else
        if (masterproc) write(iulog,*) 'Warning: Error reading SGH30 from topo file.'
        if (masterproc) write(iulog,*) 'The field SGH30 will be filled using data from SGH.'
        call pbuf_set_field(pbuf2d, sgh30_idx, tptr)
      end if

      deallocate(tptr_2)

      call infld('LANDM_COSLAT', fh_topo, dim1name, dim2name, 1, pcols, begchunk, endchunk, &
           tptr, found, gridname='physgrid')

      if(.not.found) call endrun(' ERROR: LANDM_COSLAT not found on topo dataset.')

      call pbuf_set_field(pbuf2d, landm_idx, tptr)

    else
      call pbuf_set_field(pbuf2d, sgh_idx, 0._r8)
      call pbuf_set_field(pbuf2d, sgh30_idx, 0._r8)
      call pbuf_set_field(pbuf2d, landm_idx, 0._r8)
    end if

    call infld('PBLH', fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, &
         tptr(:,:), found, gridname='physgrid')
    if(.not. found) then
       tptr(:,:) = 0._r8
       if (masterproc) write(iulog,*) 'PBLH initialized to 0.'
    end if
    pblh_idx = pbuf_get_index('pblh')

    call pbuf_set_field(pbuf2d, pblh_idx, tptr)

    call infld('TPERT', fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, &
         tptr(:,:), found, gridname='physgrid')
    if(.not. found) then
       tptr(:,:) = 0._r8
       if (masterproc) write(iulog,*) 'TPERT initialized to 0.'
    end if
    tpert_idx = pbuf_get_index( 'tpert')
    call pbuf_set_field(pbuf2d, tpert_idx, tptr)

    fieldname='QPERT'
    qpert_idx = pbuf_get_index( 'qpert',ierr)
    if (qpert_idx > 0) then
       call infld(fieldname, fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, &
            tptr(:,:), found, gridname='physgrid')
       if(.not. found) then
          tptr(:,:) = 0._r8
          if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.'
       end if
       call pbuf_set_field(pbuf2d, qpert_idx, tptr)
    end if

    fieldname='CUSH'
    m = pbuf_get_index('cush', ierr)
    if (m > 0) then
       call infld(fieldname, fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, &
            tptr, found, gridname='physgrid')
       if(.not.found) then
          if(masterproc) write(iulog,*) trim(fieldname), ' initialized to 1000.'
          tptr=1000._r8
       end if
       do n=1,dyn_time_lvls
          call pbuf_set_field(pbuf2d, m, tptr, start=(/1,n/), kount=(/pcols,1/))
       end do
       deallocate(tptr)
    end if

    !
    ! 3-D fields
    !

    allocate(tptr3d(pcols,pver,begchunk:endchunk))

    fieldname='CLOUD'
    m = pbuf_get_index('CLD')
    call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
         tptr3d, found, gridname='physgrid')
    if(found) then
       do n = 1, dyn_time_lvls
          call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/))
       end do
    else
       call pbuf_set_field(pbuf2d, m, 0._r8)
       if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.'
    end if

    fieldname='QCWAT'
    m = pbuf_get_index(fieldname,ierr)
    if (m > 0) then
       call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
            tptr3d, found, gridname='physgrid')
       if(.not. found) then
          call infld('Q',fh_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
               tptr3d, found, gridname='physgrid')
          if (found) then
             if (masterproc) write(iulog,*) trim(fieldname), ' initialized with Q'
             if(dycore_is('LR')) call polar_average(pver, tptr3d)
          else
             if (masterproc) write(iulog,*) trim(fieldname), ' initialized to huge()'
             tptr3d = huge(1.0_r8)
          end if
       end if
       do n = 1, dyn_time_lvls
          call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/))
       end do
    end if

    fieldname = 'ICCWAT'
    m = pbuf_get_index(fieldname, ierr)
    if (m > 0) then
       call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
          tptr3d, found, gridname='physgrid')
       if(found) then
          do n = 1, dyn_time_lvls
             call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/))
          end do
       else
          call cnst_get_ind('CLDICE', ixcldice)
          call infld('CLDICE',fh_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
             tptr3d, found, gridname='physgrid')
          if(found) then
             do n = 1, dyn_time_lvls
                call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/))
             end do
          else
             call pbuf_set_field(pbuf2d, m, 0._r8)
          end if
          if (masterproc) then
             if (found) then
                write(iulog,*) trim(fieldname), ' initialized with CLDICE'
             else
                write(iulog,*) trim(fieldname), ' initialized to 0.0'
             end if
          end if
       end if
    end if

    fieldname = 'LCWAT'
    m = pbuf_get_index(fieldname,ierr)
    if (m > 0) then
       call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
            tptr3d, found, gridname='physgrid')
       if(found) then
          do n = 1, dyn_time_lvls
             call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/))
          end do
       else
          allocate(tptr3d_2(pcols,pver,begchunk:endchunk))
          call cnst_get_ind('CLDICE', ixcldice)
          call cnst_get_ind('CLDLIQ', ixcldliq)
          call infld('CLDICE',fh_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
               tptr3d, found, gridname='physgrid')
          call infld('CLDLIQ',fh_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
               tptr3d_2, found2, gridname='physgrid')
          if(found .and. found2) then
             do lchnk = begchunk, endchunk
                ncol = get_ncols_p(lchnk)
                tptr3d(:ncol,:,lchnk)=tptr3d(:ncol,:,lchnk)+tptr3d_2(:ncol,:,lchnk)
             end do
             if (masterproc) write(iulog,*) trim(fieldname), ' initialized with CLDICE + CLDLIQ'
          else if (found) then ! Data already loaded in tptr3d
             if (masterproc) write(iulog,*) trim(fieldname), ' initialized with CLDICE only'
          else if (found2) then
             tptr3d(:,:,:)=tptr3d_2(:,:,:)
             if (masterproc) write(iulog,*) trim(fieldname), ' initialized with CLDLIQ only'
          end if

          if (found .or. found2) then
             do n = 1, dyn_time_lvls
                call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/))
             end do
             if(dycore_is('LR')) call polar_average(pver, tptr3d)
          else
             call pbuf_set_field(pbuf2d, m, 0._r8)
             if (masterproc)  write(iulog,*) trim(fieldname), ' initialized to 0.0'
          end if
          deallocate(tptr3d_2)
       end if
    end if

    deallocate(tptr3d)
    allocate(tptr3d(pcols,pver,begchunk:endchunk))

    fieldname = 'TCWAT'
    m = pbuf_get_index(fieldname,ierr)
    if (m > 0) then
       call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
            tptr3d, found, gridname='physgrid')
       if(.not.found) then
          call infld('T', fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
               tptr3d, found, gridname='physgrid')
          if (found) then
             if(dycore_is('LR')) call polar_average(pver, tptr3d)
             if (masterproc) write(iulog,*) trim(fieldname), ' initialized with T'
          else
             if (masterproc) write(iulog,*) trim(fieldname), ' initialized to huge()'
             tptr3d = huge(1._r8)
          end if
       end if
       do n = 1, dyn_time_lvls
          call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/))
       end do
    end if

    deallocate(tptr3d)
    allocate(tptr3d(pcols,pverp,begchunk:endchunk))

    fieldname = 'TKE'
    m = pbuf_get_index( 'tke')
    call infld(fieldname, fh_ini, dim1name, 'ilev', dim2name, 1, pcols, 1, pverp, begchunk, endchunk, &
         tptr3d, found, gridname='physgrid')
    if (found) then
       call pbuf_set_field(pbuf2d, m, tptr3d)
    else
       call pbuf_set_field(pbuf2d, m, 0.01_r8)
       if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.01'
    end if


    fieldname = 'KVM'
    m = pbuf_get_index('kvm')
    call infld(fieldname, fh_ini, dim1name, 'ilev', dim2name, 1, pcols, 1, pverp, begchunk, endchunk, &
         tptr3d, found, gridname='physgrid')
    if (found) then
       call pbuf_set_field(pbuf2d, m, tptr3d)
    else
       call pbuf_set_field(pbuf2d, m, 0._r8)
       if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.'
    end if


    fieldname = 'KVH'
    m = pbuf_get_index('kvh')
    call infld(fieldname, fh_ini, dim1name, 'ilev', dim2name, 1, pcols, 1, pverp, begchunk, endchunk, &
         tptr3d, found, gridname='physgrid')
    if (found) then
       call pbuf_set_field(pbuf2d, m, tptr3d)
    else
       call pbuf_set_field(pbuf2d, m, 0._r8)
       if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.'
    end if

    deallocate(tptr3d)
    allocate(tptr3d(pcols,pver,begchunk:endchunk))

    fieldname = 'CONCLD'
    m = pbuf_get_index('CONCLD',ierr)
    if (m > 0) then
       call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
            tptr3d, found, gridname='physgrid')
       if(found) then
          do n = 1, dyn_time_lvls
             call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/))
          end do
       else
          call pbuf_set_field(pbuf2d, m, 0._r8)
          if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.'
       end if

       deallocate (tptr3d)
    end if

    call initialize_short_lived_species(fh_ini, pbuf2d)

    !---------------------------------------------------------------------------------
    !  If needed, get ion and electron temperature fields from initial condition file
    !---------------------------------------------------------------------------------

    call waccmx_phys_ion_elec_temp_inidat(fh_ini,pbuf2d)

  end subroutine phys_inidat


  subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )

    !-----------------------------------------------------------------------
    !
    ! Initialization of physics package.
    !
    !-----------------------------------------------------------------------

    use physics_buffer,     only: physics_buffer_desc, pbuf_initialize, pbuf_get_index
    use physconst,          only: rair, cpair, gravit, zvir, karman
    use cam_thermo,         only: cam_thermo_init
    use ref_pres,           only: pref_edge, pref_mid

    use carma_intr,         only: carma_init
    use cam_control_mod,    only: initial_run
    use check_energy,       only: check_energy_init
    use chemistry,          only: chem_init
    use mo_lightning,       only: lightning_init
    use prescribed_ozone,   only: prescribed_ozone_init
    use prescribed_ghg,     only: prescribed_ghg_init
    use prescribed_aero,    only: prescribed_aero_init
    use aerodep_flx,        only: aerodep_flx_init
    use aircraft_emit,      only: aircraft_emit_init
    use prescribed_volcaero,only: prescribed_volcaero_init
    use prescribed_strataero,only: prescribed_strataero_init
    use cloud_fraction,     only: cldfrc_init
    use cldfrc2m,           only: cldfrc2m_init
    use co2_cycle,          only: co2_init, co2_transport
    use convect_deep,       only: convect_deep_init
    use convect_shallow,    only: convect_shallow_init
    use constituents,       only: cnst_get_ind
    use cam_diagnostics,    only: diag_init
    use gw_drag_cam,        only: gw_drag_cam_init
    use radheat,            only: radheat_init
    use radiation,          only: radiation_init
    use cloud_diagnostics,  only: cloud_diagnostics_init
    use rk_stratiform_cam,  only: rk_stratiform_cam_init
    use wv_saturation,      only: wv_sat_init
    use microp_driver,      only: microp_driver_init
    use microp_aero,        only: microp_aero_init
    use macrop_driver,      only: macrop_driver_init
    use conv_water,         only: conv_water_init
    use tracers,            only: tracers_init
    use aoa_tracers,        only: aoa_tracers_init
    use rayleigh_friction,  only: rayleigh_friction_init
    use rayleigh_friction_cam, only: rf_nl_k0, rf_nl_krange, rf_nl_tau0
    use vertical_diffusion, only: vertical_diffusion_init
    use phys_debug_util,    only: phys_debug_init
    use rad_constituents,   only: rad_cnst_init
    use aer_rad_props,      only: aer_rad_props_init
    use subcol,             only: subcol_init
    use qbo,                only: qbo_init
    use qneg_module,        only: qneg_init
    use lunar_tides,        only: lunar_tides_init
    use iondrag,            only: iondrag_init
#if ( defined OFFLINE_DYN )
    use metdata,            only: metdata_phys_init
#endif
    use epp_ionization,     only: epp_ionization_init, epp_ionization_active
    use waccmx_phys_intr,   only: waccmx_phys_ion_elec_temp_init  ! Initialization of ionosphere module (WACCM-X)
    use waccmx_phys_intr,   only: waccmx_phys_mspd_init   ! Initialization of major species diffusion module (WACCM-X)
    use clubb_intr,         only: clubb_ini_cam
    use sslt_rebin,         only: sslt_rebin_init
    use tropopause,         only: tropopause_init
    use solar_data,         only: solar_data_init
    use dadadj_cam,         only: dadadj_cam_init
    use cam_abortutils,     only: endrun
    use nudging,            only: Nudge_Model, nudging_init
    use cam_snapshot,       only: cam_snapshot_init
    use cam_history,        only: addfld, register_vector_field, add_default
    use phys_control,       only: phys_getopts
    use phys_grid_ctem,     only: phys_grid_ctem_init
    use cam_budget,         only: cam_budget_init
    use surface_emissions_mod, only: surface_emissions_init
    use elevated_emissions_mod, only: elevated_emissions_init

    use ccpp_constituent_prop_mod, only: ccpp_const_props_init
    use ctem_diags_mod, only: ctem_diags_init

    ! Input/output arguments
    type(physics_state), pointer       :: phys_state(:)
    type(physics_tend ), pointer       :: phys_tend(:)
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)

    type(cam_in_t), intent(in)         :: cam_in(begchunk:endchunk)
    type(cam_out_t),intent(inout)      :: cam_out(begchunk:endchunk)

    ! local variables
    integer :: lchnk
    integer :: ierr, ixq

    logical :: history_budget              ! output tendencies and state variables for
                                           ! temperature, water vapor, cloud
                                           ! ice, cloud liquid, U, V
    integer :: history_budget_histfile_num ! output history file number for budget fields

    ! Needed for rayleigh friction
    character(len=512) errmsg
    integer errflg

    !-----------------------------------------------------------------------

    call physics_type_alloc(phys_state, phys_tend, begchunk, endchunk, pcols)

    do lchnk = begchunk, endchunk
       call physics_state_set_grid(lchnk, phys_state(lchnk))
    end do

    !-------------------------------------------------------------------------------------------
    ! Initialize any variables in cam_thermo which are not temporally and/or spatially constant
    !-------------------------------------------------------------------------------------------
    call cam_thermo_init()

    ! Initialize debugging a physics column
    call phys_debug_init()

    call pbuf_initialize(pbuf2d)

    ! Initialize subcol scheme
    call subcol_init(pbuf2d)

    ! diag_init makes addfld calls for dynamics fields that are output from
    ! the physics decomposition
    call diag_init(pbuf2d)

    call check_energy_init()

    call tracers_init()

    ! age of air tracers
    call aoa_tracers_init()

    teout_idx = pbuf_get_index( 'TEOUT')

    ! adiabatic or ideal physics should be only used if in simple_physics
    if (adiabatic .or. ideal_phys) then
      if (adiabatic) then
        call endrun('phys_init: adiabatic configuration error')
      else
        call endrun('phys_init: ideal_phys configuration error')
      end if
    end if

    if (initial_run) then
       call phys_inidat(cam_out, pbuf2d)
    end if

    ! wv_saturation is relatively independent of everything else and
    ! low level, so init it early. Must at least do this before radiation.
    call wv_sat_init

    ! solar irradiance data modules
    call solar_data_init()

    ! Initialize rad constituents and their properties
    call rad_cnst_init()

    call radiation_init(pbuf2d)

    call aer_rad_props_init()

    ! initialize carma
    call carma_init(pbuf2d)

    if (.not. offline_driver_dorun) then

       call surface_emissions_init(pbuf2d)
       call elevated_emissions_init(pbuf2d)

       ! Prognostic chemistry.
       call chem_init(phys_state,pbuf2d)

       ! Lightning flash frq and NOx prod
       call lightning_init( pbuf2d )

       ! Prescribed tracers
       call prescribed_ozone_init()
       call prescribed_ghg_init()
       call prescribed_aero_init()
       call aerodep_flx_init()
       call aircraft_emit_init()
       call prescribed_volcaero_init()
       call prescribed_strataero_init()
    end if

    ! co2 cycle
    if (co2_transport()) then
       call co2_init()
    end if

    call gw_drag_cam_init()

    call rayleigh_friction_init(pver, rf_nl_tau0, rf_nl_krange, rf_nl_k0, masterproc, &
         iulog, errmsg, errflg)
    if (errflg /= 0) call endrun(errmsg)

    call vertical_diffusion_init(pbuf2d)

    if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
       call waccmx_phys_mspd_init ()
       ! Initialization of ionosphere module if mode set to ionosphere
       if( waccmx_is('ionosphere') ) then
          call waccmx_phys_ion_elec_temp_init(pbuf2d)
       endif
    endif

    call cloud_diagnostics_init(pbuf2d)

    call radheat_init(pref_mid)

    call convect_shallow_init(pref_edge, pbuf2d)

    call cldfrc_init()
    call cldfrc2m_init()

    call convect_deep_init(pref_edge)

    if( microp_scheme == 'RK' ) then
       call rk_stratiform_cam_init()
    elseif( microp_scheme == 'MG' ) then
       if (.not. do_clubb_sgs) call macrop_driver_init(pbuf2d)
       call microp_aero_init(phys_state,pbuf2d)
       call microp_driver_init(pbuf2d)
       call conv_water_init
    end if

    ! initiate CLUBB within CAM
    if (do_clubb_sgs) call clubb_ini_cam(pbuf2d)

    call qbo_init

    call lunar_tides_init()

    call iondrag_init(pref_mid)
    ! Geomagnetic module -- after iondrag_init
    if (epp_ionization_active) then
      call epp_ionization_init()
    endif

#if ( defined OFFLINE_DYN )
    call metdata_phys_init()
#endif
    call sslt_rebin_init()
    call tropopause_init()
    call dadadj_cam_init()

    prec_dp_idx  = pbuf_get_index('PREC_DP')
    snow_dp_idx  = pbuf_get_index('SNOW_DP')
    prec_sh_idx  = pbuf_get_index('PREC_SH')
    snow_sh_idx  = pbuf_get_index('SNOW_SH')

    dlfzm_idx = pbuf_get_index('DLFZM', ierr)

    call phys_getopts(prog_modal_aero_out=prog_modal_aero)

    ! Initialize Nudging Parameters
    !--------------------------------
    if(Nudge_Model) call nudging_init

    if (clim_modal_aero) then

       ! If climate calculations are affected by prescribed modal aerosols, the
       ! the initialization routine for the dry mode radius calculation is called
       ! here.  For prognostic MAM the initialization is called from
       ! modal_aero_initialize
       if (.not. prog_modal_aero) then
          call modal_aero_calcsize_init(pbuf2d)
       endif

       call modal_aero_wateruptake_init(pbuf2d)

    end if

    ! Initialize CAM CCPP constituent properties array
    ! for use in CCPP-ized physics schemes:
    call cnst_get_ind('Q', ixq)
    call ccpp_const_props_init(ixq)

    ! Initialize qneg3 and qneg4
    call qneg_init()

    ! Initialize phys TEM diagnostics
    call phys_grid_ctem_init()

    ! Initialize the snapshot capability
    call cam_snapshot_init(cam_in, cam_out, pbuf2d, begchunk)

    ! Initialize the budget capability
    call cam_budget_init()

    ! addfld calls for U, V tendency budget variables that are output in
    ! tphysac, tphysbc
    call addfld ( 'UTEND_DCONV', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by deep convection')
    call addfld ( 'VTEND_DCONV', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by deep convection')
    call register_vector_field ( 'UTEND_DCONV', 'VTEND_DCONV')
    call addfld ( 'UTEND_SHCONV', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by shallow convection')
    call addfld ( 'VTEND_SHCONV', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by shallow convection')
    call register_vector_field ( 'UTEND_SHCONV', 'VTEND_SHCONV')
    call addfld ( 'UTEND_MACROP', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by macrophysics')
    call addfld ( 'VTEND_MACROP', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by macrophysics')
    call register_vector_field ( 'UTEND_MACROP', 'VTEND_MACROP')
    call addfld ( 'UTEND_VDIFF', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by vert. diffus.')
    call addfld ( 'VTEND_VDIFF', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by vert. diffus.')
    call register_vector_field ( 'UTEND_VDIFF', 'VTEND_VDIFF')
    call addfld ( 'UTEND_RAYLEIGH', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by Rayleigh Fric.')
    call addfld ( 'VTEND_RAYLEIGH', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by Rayleigh Fric.')
    call register_vector_field ( 'UTEND_RAYLEIGH', 'VTEND_RAYLEIGH')
    call addfld ( 'UTEND_GWDTOT', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by all GWs')
    call addfld ( 'VTEND_GWDTOT', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by all GWs')
    call register_vector_field ( 'UTEND_GWDTOT', 'VTEND_GWDTOT')
    call addfld ( 'UTEND_QBORLX', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by QBO relaxation')
    call addfld ( 'VTEND_QBORLX', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by QBO relaxation')
    call register_vector_field ( 'UTEND_QBORLX', 'VTEND_QBORLX')
    call addfld ( 'UTEND_LUNART', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by lunar tides')
    call addfld ( 'VTEND_LUNART', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by lunar tides')
    call register_vector_field ( 'UTEND_LUNART', 'VTEND_LUNART')
    call addfld ( 'UTEND_IONDRG', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by ion drag')
    call addfld ( 'VTEND_IONDRG', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by ion drag')
    call register_vector_field ( 'UTEND_IONDRG', 'VTEND_IONDRG')
    call addfld ( 'UTEND_NDG', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by nudging')
    call addfld ( 'VTEND_NDG', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by nudging')
    call register_vector_field ( 'UTEND_NDG', 'VTEND_NDG')
    call addfld('UTEND_CORE', (/ 'lev' /), 'A', 'm/s2' , 'Zonal wind tendency due to dynamical core')
    call addfld('VTEND_CORE', (/ 'lev' /), 'A', 'm/s2' , 'Meridional wind tendency due to dynamical core')
    call register_vector_field('UTEND_CORE','VTEND_CORE')


    call phys_getopts(history_budget_out = history_budget, &
         history_budget_histfile_num_out = history_budget_histfile_num)

    if ( history_budget ) then
       call add_default ( 'UTEND_DCONV'   , history_budget_histfile_num, ' ')
       call add_default ( 'VTEND_DCONV'   , history_budget_histfile_num, ' ')
       call add_default ( 'UTEND_SHCONV'   , history_budget_histfile_num, ' ')
       call add_default ( 'VTEND_SHCONV'   , history_budget_histfile_num, ' ')
       call add_default ( 'UTEND_MACROP'   , history_budget_histfile_num, ' ')
       call add_default ( 'VTEND_MACROP'   , history_budget_histfile_num, ' ')
       call add_default ( 'UTEND_VDIFF'   , history_budget_histfile_num, ' ')
       call add_default ( 'VTEND_VDIFF'   , history_budget_histfile_num, ' ')
       call add_default ( 'UTEND_RAYLEIGH'   , history_budget_histfile_num, ' ')
       call add_default ( 'VTEND_RAYLEIGH'   , history_budget_histfile_num, ' ')
       call add_default ( 'UTEND_GWDTOT'   , history_budget_histfile_num, ' ')
       call add_default ( 'VTEND_GWDTOT'   , history_budget_histfile_num, ' ')
       call add_default ( 'UTEND_QBORLX'   , history_budget_histfile_num, ' ')
       call add_default ( 'VTEND_QBORLX'   , history_budget_histfile_num, ' ')
       call add_default ( 'UTEND_LUNART'   , history_budget_histfile_num, ' ')
       call add_default ( 'VTEND_LUNART'   , history_budget_histfile_num, ' ')
       call add_default ( 'UTEND_IONDRG'   , history_budget_histfile_num, ' ')
       call add_default ( 'VTEND_IONDRG'   , history_budget_histfile_num, ' ')
       call add_default ( 'UTEND_NDG'   , history_budget_histfile_num, ' ')
       call add_default ( 'VTEND_NDG'   , history_budget_histfile_num, ' ')
       call add_default ( 'UTEND_CORE'   , history_budget_histfile_num, ' ')
       call add_default ( 'VTEND_CORE'   , history_budget_histfile_num, ' ')
    end if

    ducore_idx = pbuf_get_index('DUCORE')
    dvcore_idx = pbuf_get_index('DVCORE')
    dtcore_idx = pbuf_get_index('DTCORE')
    dqcore_idx = pbuf_get_index('DQCORE')

    call ctem_diags_init()

  end subroutine phys_init

  !
  !-----------------------------------------------------------------------
  !

  subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d,  cam_in, cam_out)
    !-----------------------------------------------------------------------
    !
    ! Purpose:
    ! First part of atmospheric physics package before updating of surface models
    !
    !-----------------------------------------------------------------------
    use time_manager,   only: get_nstep
    use cam_diagnostics,only: diag_allocate, diag_physvar_ic
    use check_energy,   only: check_energy_gmean
    use phys_control,   only: phys_getopts
    use spmd_utils,     only: mpicom
    use physics_buffer, only: physics_buffer_desc, pbuf_get_chunk, pbuf_allocate
    use cam_history,    only: outfld, write_camiop
    use cam_abortutils, only: endrun
#if ( defined OFFLINE_DYN )
     use metdata,       only: get_met_srf1
#endif
    use ctem_diags_mod, only: ctem_diags_calc
    !
    ! Input arguments
    !
    real(r8), intent(in) :: ztodt            ! physics time step unless nstep=0
    !
    ! Input/Output arguments
    !
    type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
    type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend

    type(physics_buffer_desc), pointer, dimension(:,:) :: pbuf2d
    type(cam_in_t),                     dimension(begchunk:endchunk) :: cam_in
    type(cam_out_t),                    dimension(begchunk:endchunk) :: cam_out
    !-----------------------------------------------------------------------
    !
    !---------------------------Local workspace-----------------------------
    !
    integer :: c                                 ! indices
    integer :: nstep                             ! current timestep number
    type(physics_buffer_desc), pointer :: phys_buffer_chunk(:)

    call t_startf ('physpkg_st1')
    nstep = get_nstep()

#if ( defined OFFLINE_DYN )
    !
    ! if offline mode set SNOWH and TS for micro-phys
    !
    call get_met_srf1( cam_in )
#endif

    ! The following initialization depends on the import state (cam_in)
    ! being initialized.  This isn't true when cam_init is called, so need
    ! to postpone this initialization to here.
    if (nstep == 0 .and. phys_do_flux_avg()) call flux_avg_init(cam_in,  pbuf2d)

    ! Compute total energy of input state and previous output state
    call t_startf ('chk_en_gmean')
    call check_energy_gmean(phys_state, pbuf2d, ztodt, nstep)
    call t_stopf ('chk_en_gmean')

    call pbuf_allocate(pbuf2d, 'physpkg')
    call diag_allocate()

    ! TEM diagnostics
    call ctem_diags_calc(phys_state)

    !-----------------------------------------------------------------------
    ! Advance time information
    !-----------------------------------------------------------------------

    call phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d)

    call t_stopf ('physpkg_st1')

#ifdef TRACER_CHECK
    call gmean_mass ('before tphysbc DRY', phys_state)
#endif


    !-----------------------------------------------------------------------
    ! Tendency physics before flux coupler invocation
    !-----------------------------------------------------------------------
    !

    if (write_camiop) then
       do c=begchunk, endchunk
          call outfld('Tg',cam_in(c)%ts,pcols   ,c     )
       end do
    end if

    call t_barrierf('sync_bc_physics', mpicom)
    call t_startf ('bc_physics')
    call t_adj_detailf(+1)

!$OMP PARALLEL DO PRIVATE (C, phys_buffer_chunk)
    do c=begchunk, endchunk
      !
      ! Output physics terms to IC file
      !
      phys_buffer_chunk => pbuf_get_chunk(pbuf2d, c)

      call t_startf ('diag_physvar_ic')
      call diag_physvar_ic ( c,  phys_buffer_chunk, cam_out(c), cam_in(c) )
      call t_stopf ('diag_physvar_ic')

      call tphysbc(ztodt, phys_state(c), phys_tend(c), phys_buffer_chunk, &
                   cam_out(c), cam_in(c) )
    end do

    call t_adj_detailf(-1)
    call t_stopf ('bc_physics')

    ! Don't call the rest in CRM mode
    if(single_column.and.scm_crm_mode) return

#ifdef TRACER_CHECK
    call gmean_mass ('between DRY', phys_state)
#endif

  end subroutine phys_run1

  !
  !-----------------------------------------------------------------------
  !

  subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d,  cam_out, &
       cam_in )
    !-----------------------------------------------------------------------
    !
    ! Purpose:
    ! Second part of atmospheric physics package after updating of surface models
    !
    !-----------------------------------------------------------------------
    use physics_buffer,  only: physics_buffer_desc, pbuf_get_chunk, pbuf_deallocate, pbuf_update_tim_idx
    use mo_lightning,    only: lightning_no_prod
    use cam_diagnostics, only: diag_deallocate, diag_surf
    use carma_intr,      only: carma_accumulate_stats
    use spmd_utils,      only: mpicom
    use iop_forcing,     only: scam_use_iop_srf
#if ( defined OFFLINE_DYN )
    use metdata,         only: get_met_srf2
#endif
    use hemco_interface,  only: HCOI_Chunk_Run
    !
    ! Input arguments
    !
    real(r8), intent(in) :: ztodt                       ! physics time step unless nstep=0
    !
    ! Input/Output arguments
    !
    type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
    type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend
    type(physics_buffer_desc),pointer,  dimension(:,:)               :: pbuf2d

    type(cam_out_t),     intent(inout), dimension(begchunk:endchunk) :: cam_out
    type(cam_in_t),      intent(inout), dimension(begchunk:endchunk) :: cam_in
    !
    !-----------------------------------------------------------------------
    !---------------------------Local workspace-----------------------------
    !
    integer :: c                                 ! chunk index
    integer :: ncol                              ! number of columns
    type(physics_buffer_desc),pointer, dimension(:)     :: phys_buffer_chunk
    !
    ! If exit condition just return
    !

    if(single_column.and.scm_crm_mode) then
       call diag_deallocate()
       return
    end if
    !-----------------------------------------------------------------------
    ! if using IOP values for surface fluxes overwrite here after surface components run
    !-----------------------------------------------------------------------
    if (single_column) call scam_use_iop_srf(cam_in)


    if(use_hemco) then
        !----------------------------------------------------------
        ! run hemco (phase 2 before chemistry)
        ! only phase 2 is used currently for HEMCO-CESM
        !----------------------------------------------------------
        call HCOI_Chunk_Run(cam_in, phys_state, pbuf2d, phase=2)
    endif

    !-----------------------------------------------------------------------
    ! Tendency physics after coupler
    ! Not necessary at terminal timestep.
    !-----------------------------------------------------------------------
    !
#if ( defined OFFLINE_DYN )
    !
    ! if offline mode set SHFLX QFLX TAUX TAUY for vert diffusion
    !
    call get_met_srf2( cam_in )
#endif
    ! lightning flash freq and prod rate of NOx
    call t_startf ('lightning_no_prod')
    call lightning_no_prod( phys_state, pbuf2d, cam_in )
    call t_stopf ('lightning_no_prod')

    call t_barrierf('sync_ac_physics', mpicom)
    call t_startf ('ac_physics')
    call t_adj_detailf(+1)

!$OMP PARALLEL DO PRIVATE (C, NCOL, phys_buffer_chunk)

    do c=begchunk,endchunk
       ncol = get_ncols_p(c)
       phys_buffer_chunk => pbuf_get_chunk(pbuf2d, c)
       !
       ! surface diagnostics for history files
       !
       call t_startf('diag_surf')
       call diag_surf(cam_in(c), cam_out(c), phys_state(c), phys_buffer_chunk)
       call t_stopf('diag_surf')

       call tphysac(ztodt, cam_in(c),  &
            cam_out(c),                              &
            phys_state(c), phys_tend(c), phys_buffer_chunk)
    end do                    ! Chunk loop

    call t_adj_detailf(-1)
    call t_stopf('ac_physics')

#ifdef TRACER_CHECK
    call gmean_mass ('after tphysac FV:WET)', phys_state)
#endif

    call t_startf ('carma_accumulate_stats')
    call carma_accumulate_stats()
    call t_stopf ('carma_accumulate_stats')

    call t_startf ('physpkg_st2')
    call pbuf_deallocate(pbuf2d, 'physpkg')

    call pbuf_update_tim_idx()
    call diag_deallocate()
    call t_stopf ('physpkg_st2')

  end subroutine phys_run2

  !
  !-----------------------------------------------------------------------
  !

  subroutine phys_final( phys_state, phys_tend, pbuf2d )
    use physics_buffer, only : physics_buffer_desc, pbuf_deallocate
    use chemistry, only : chem_final
    use carma_intr, only : carma_final
    use wv_saturation, only : wv_sat_final
    use hemco_interface,  only: HCOI_Chunk_Final
    use microp_aero, only : microp_aero_final
    use phys_grid_ctem, only : phys_grid_ctem_final
    use nudging, only: Nudge_Model, nudging_final
    use ctem_diags_mod, only: ctem_diags_final

    !-----------------------------------------------------------------------
    !
    ! Purpose:
    ! Finalization of physics package
    !
    !-----------------------------------------------------------------------
    ! Input/output arguments
    type(physics_state), pointer :: phys_state(:)
    type(physics_tend ), pointer :: phys_tend(:)
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)

    if(associated(pbuf2d)) then
       call pbuf_deallocate(pbuf2d,'global')
       deallocate(pbuf2d)
    end if
    deallocate(phys_state)
    deallocate(phys_tend)
    call chem_final
    call carma_final
    call wv_sat_final
    call microp_aero_final()
    call phys_grid_ctem_final()
    if(Nudge_Model) call nudging_final()

    if(use_hemco) then
        ! cleanup hemco
        call HCOI_Chunk_Final
    endif

    call ctem_diags_final()

  end subroutine phys_final


  subroutine tphysac (ztodt,   cam_in,  &
       cam_out,  state,   tend,    pbuf)
    !-----------------------------------------------------------------------
    !
    ! Tendency physics after coupling to land, sea, and ice models.
    !
    ! Computes the following:
    !
    !   o Aerosol Emission at Surface
    !   o Source-Sink for Advected Tracers
    !   o Symmetric Turbulence Scheme - Vertical Diffusion
    !   o Rayleigh Friction
    !   o Dry Deposition of Aerosol
    !   o Enforce Charge Neutrality ( Only for WACCM )
    !   o Gravity Wave Drag
    !   o QBO Relaxation ( Only for WACCM )
    !   o Ion Drag ( Only for WACCM )
    !   o Scale Dry Mass Energy
    !-----------------------------------------------------------------------
    use physics_buffer, only: physics_buffer_desc, pbuf_set_field, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx
    use shr_kind_mod,       only: r8 => shr_kind_r8
    use chemistry,          only: chem_is_active, chem_timestep_tend, chem_emissions
    use cam_diagnostics,    only: diag_phys_tend_writeout
    use gw_drag_cam,        only: gw_drag_cam_tend
    use trb_mtn_stress_cam, only: trb_mtn_stress_tend
    use beljaars_drag_cam,  only: beljaars_drag_tend
    use vertical_diffusion, only: vertical_diffusion_tend
    use rayleigh_friction,  only: rayleigh_friction_run
    use constituents,       only: cnst_get_ind
    use physics_types,      only: physics_state, physics_tend, physics_ptend, physics_update,    &
                                  physics_dme_adjust, set_dry_to_wet, physics_state_check,       &
                                  dyn_te_idx, physics_ptend_init
    use waccmx_phys_intr,   only: waccmx_phys_mspd_tend  ! WACCM-X major diffusion
    use waccmx_phys_intr,   only: waccmx_phys_ion_elec_temp_tend ! WACCM-X
    use aoa_tracers,        only: aoa_tracers_timestep_tend
    use physconst,          only: rhoh2o, latvap,latice
    use dyn_tests_utils,    only: vc_dycore
    use aero_model,         only: aero_model_drydep
    use carma_intr,         only: carma_emission_tend, carma_timestep_tend
    use carma_flags_mod,    only: carma_do_aerosol, carma_do_emission
    use check_energy,       only: tot_energy_phys
    use check_energy,       only: check_tracers_data, check_tracers_init, check_tracers_chng
    use check_energy,       only: check_energy_cam_chng
    use time_manager,       only: get_nstep
    use cam_abortutils,     only: endrun
    use dycore,             only: dycore_is
    use cam_control_mod,    only: aqua_planet
    use mo_gas_phase_chemdr,only: map2chm
    use clybryiy_fam,       only: clybryiy_fam_set
    use charge_neutrality,  only: charge_balance
    use qbo,                only: qbo_relax
    use iondrag,            only: iondrag_calc, do_waccm_ions
    use perf_mod
    use flux_avg,           only: flux_avg_run
    use cam_history,        only: outfld
    use qneg_module,        only: qneg4
    use co2_cycle,          only: co2_cycle_set_ptend
    use nudging,            only: Nudge_Model,Nudge_ON,nudging_timestep_tend
    use cam_snapshot,       only: cam_snapshot_all_outfld_tphysac
    use cam_snapshot_common,only: cam_snapshot_ptend_outfld
    use lunar_tides,        only: lunar_tides_tend
    use cam_thermo,         only: cam_thermo_water_update
    use cam_budget,         only: thermo_budget_history
    use dyn_tests_utils,    only: vc_dycore, vc_height, vc_dry_pressure
    use air_composition,    only: cpairv, cp_or_cv_dycore
    !
    ! Arguments
    !
    real(r8), intent(in) :: ztodt                  ! Two times model timestep (2 delta-t)

    type(cam_in_t),      intent(inout) :: cam_in
    type(cam_out_t),     intent(inout) :: cam_out
    type(physics_state), intent(inout) :: state
    type(physics_tend ), intent(inout) :: tend
    type(physics_buffer_desc), pointer :: pbuf(:)


    type(check_tracers_data):: tracerint             ! tracer mass integrals and cummulative boundary fluxes

    !
    !---------------------------Local workspace-----------------------------
    !
    type(physics_ptend)     :: ptend               ! indivdual parameterization tendencies

    integer  :: nstep                  ! current timestep number
    real(r8) :: zero(pcols)            ! array of zeros

    integer :: lchnk                   ! chunk identifier
    integer :: ncol                    ! number of atmospheric columns
    integer :: i,k                     ! Longitude, level indices
    integer :: ixq

    logical :: labort                  ! abort flag

    real(r8) surfric(pcols)            ! surface friction velocity
    real(r8) obklen(pcols)             ! Obukhov length
    real(r8) :: fh2o(pcols)            ! h2o flux to balance source from methane chemistry
    real(r8) :: flx_heat(pcols)        ! Heat flux for check_energy_chng.
    real(r8) :: tmp_trac  (pcols,pver,pcnst) ! tmp space
    real(r8) :: tmp_pdel  (pcols,pver) ! tmp space
    real(r8) :: tmp_ps    (pcols)      ! tmp space
    real(r8) :: scaling(pcols,pver)
    logical  :: moist_mixing_ratio_dycore

    ! physics buffer fields for total energy and mass adjustment
    integer itim_old, ifld

    real(r8), pointer, dimension(:,:) :: cld
    real(r8), pointer, dimension(:,:) :: qini
    real(r8), pointer, dimension(:,:) :: cldliqini
    real(r8), pointer, dimension(:,:) :: cldiceini
    real(r8), pointer, dimension(:,:) :: totliqini
    real(r8), pointer, dimension(:,:) :: toticeini
    real(r8), pointer, dimension(:,:) :: dtcore
    real(r8), pointer, dimension(:,:) :: dqcore
    real(r8), pointer, dimension(:,:) :: ducore
    real(r8), pointer, dimension(:,:) :: dvcore
    real(r8), pointer, dimension(:,:) :: ast     ! relative humidity cloud fraction

    ! For aerosol budget diagnostics
    type(carma_diags_t), pointer :: carma_diags_obj

    ! For rayleigh friction CCPP calls
    character(len=512) errmsg
    integer errflg

    !-----------------------------------------------------------------------
    carma_diags_obj => carma_diags_t()
    if (.not.associated(carma_diags_obj)) then
       call endrun('tphysac: carma_diags_obj allocation failed')
    end if

    lchnk = state%lchnk
    ncol  = state%ncol

    nstep = get_nstep()
    call cnst_get_ind('Q', ixq)

    ! Adjust the surface fluxes to reduce instabilities in near sfc layer
    if (phys_do_flux_avg()) then
       call flux_avg_run(state, cam_in,  pbuf, nstep, ztodt)
    endif

    ! Validate the physics state.
    if (state_debug_checks) &
         call physics_state_check(state, name="before tphysac")

    call t_startf('tphysac_init')
    ! Associate pointers with physics buffer fields
    itim_old = pbuf_old_tim_idx()

    call pbuf_get_field(pbuf, dtcore_idx, dtcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
    call pbuf_get_field(pbuf, dqcore_idx, dqcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
    call pbuf_get_field(pbuf, ducore_idx, ducore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
    call pbuf_get_field(pbuf, dvcore_idx, dvcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )

    call pbuf_get_field(pbuf, qini_idx, qini)
    call pbuf_get_field(pbuf, cldliqini_idx, cldliqini)
    call pbuf_get_field(pbuf, cldiceini_idx, cldiceini)
    call pbuf_get_field(pbuf, totliqini_idx, totliqini)
    call pbuf_get_field(pbuf, toticeini_idx, toticeini)

    ifld = pbuf_get_index('CLD')
    call pbuf_get_field(pbuf, ifld, cld, start=(/1,1,itim_old/),kount=(/pcols,pver,1/))

    ifld = pbuf_get_index('AST')
    call pbuf_get_field(pbuf, ifld, ast, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )

    ! zero out local variables that may be written to snapshot for safety.
    fh2o(:) = 0._r8         ! used in chem_timestep_tend.
    surfric(:) = 0._r8      ! out from vertical_diffusion_tend.
    obklen(:) = 0._r8       ! out from vertical_diffusion_tend.
    flx_heat(:) = 0._r8     ! first out from gw_drag_cam.

    !
    ! accumulate fluxes into net flux array for spectral dycores
    ! jrm Include latent heat of fusion for snow
    !
    do i=1,ncol
       tend%flx_net(i) = tend%flx_net(i) + cam_in%shf(i) + (cam_out%precc(i) &
            + cam_out%precl(i))*latvap*rhoh2o &
            + (cam_out%precsc(i) + cam_out%precsl(i))*latice*rhoh2o
    end do

    ! emissions of aerosols and gas-phase chemistry constituents at surface

    if (trim(cam_take_snapshot_before) == "chem_emissions") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    call carma_diags_obj%update(cam_in, state, pbuf)

    call chem_emissions( state, cam_in, pbuf )

    call carma_diags_obj%output(state, ptend, cam_in, "CHEMEMIS", ztodt, pbuf)

    if (trim(cam_take_snapshot_after) == "chem_emissions") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    if (carma_do_emission) then
       ! carma emissions
       call carma_diags_obj%update(cam_in, state, pbuf)
       call carma_emission_tend(state, ptend, cam_in, ztodt, pbuf)
       call carma_diags_obj%output(state, ptend, cam_in, "CREMIS", ztodt, pbuf)
       call physics_update(state, ptend, ztodt, tend)
    end if

    ! get nstep and zero array for energy checker
    zero = 0._r8
    nstep = get_nstep()
    call check_tracers_init(state, tracerint)

    ! Check if latent heat flux exceeds the total moisture content of the
    ! lowest model layer, thereby creating negative moisture.

    call qneg4('TPHYSAC', lchnk, ncol, ztodt ,                                &
         state%q(1,pver,1), state%rpdel(1,pver),                              &
         cam_in%shf, cam_in%lhf, cam_in%cflx)

    call t_stopf('tphysac_init')
    !===================================================
    ! Source/sink terms for advected tracers.
    !===================================================
    call t_startf('adv_tracer_src_snk')
    ! Test tracers

    if (trim(cam_take_snapshot_before) == "aoa_tracers_timestep_tend") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if
    call aoa_tracers_timestep_tend(state, ptend, ztodt)
    if ( (trim(cam_take_snapshot_after) == "aoa_tracers_timestep_tend") .and. &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
       call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if
    call physics_update(state, ptend, ztodt, tend)
    if (trim(cam_take_snapshot_after) == "aoa_tracers_timestep_tend") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if
    call check_tracers_chng(state, tracerint, "aoa_tracers_timestep_tend", nstep, ztodt,   &
         cam_in%cflx)

    if (trim(cam_take_snapshot_before) == "co2_cycle_set_ptend") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if
    call co2_cycle_set_ptend(state, pbuf, ptend)
    if ( (trim(cam_take_snapshot_after) == "co2_cycle_set_ptend") .and.       &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
       call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if
    call physics_update(state, ptend, ztodt, tend)
    if (trim(cam_take_snapshot_after) == "co2_cycle_set_ptend") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    !===================================================
    ! Chemistry and MAM calculation
    ! MAM core aerosol conversion process is performed in the below 'chem_timestep_tend'.
    ! In addition, surface flux of aerosol species other than 'dust' and 'sea salt', and
    ! elevated emission of aerosol species are treated in 'chem_timestep_tend' before
    ! Gas chemistry and MAM core aerosol conversion.
    ! Note that surface flux is not added into the atmosphere, but elevated emission is
    ! added into the atmosphere as tendency.
    !===================================================
    if (chem_is_active()) then

       if (trim(cam_take_snapshot_before) == "chem_timestep_tend") then
          call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
       end if

       call carma_diags_obj%update(cam_in, state, pbuf)

       call chem_timestep_tend(state, ptend, cam_in, cam_out, ztodt, &
            pbuf,  fh2o=fh2o)


       if ( (trim(cam_take_snapshot_after) == "chem_timestep_tend") .and.     &
            (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
          call cam_snapshot_ptend_outfld(ptend, lchnk)
       end if

       call carma_diags_obj%output(state, ptend, cam_in, "CHEM", ztodt, pbuf)

       call physics_update(state, ptend, ztodt, tend)

       if (trim(cam_take_snapshot_after) == "chem_timestep_tend") then
          call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
       end if
       call check_energy_cam_chng(state, tend, "chem", nstep, ztodt, fh2o, zero, zero, zero)
       call check_tracers_chng(state, tracerint, "chem_timestep_tend", nstep, ztodt, &
            cam_in%cflx)
    end if
    call t_stopf('adv_tracer_src_snk')

    !===================================================
    ! Vertical diffusion/pbl calculation
    !===================================================
    call t_startf('vertical_diffusion_tend')

    !------------------------------------------
    ! Compute orographic form drag stress
    ! (Beljaars in CAM6+; TMS in CAM5)
    !
    ! Only the pbuf fields are updated in these routines (no ptend)
    ! The computed stresses are used in the PBL scheme and the vertical diffusion solver
    ! in the call to vertical_diffusion_tend below.
    !------------------------------------------
    if (trim(cam_take_snapshot_before) == "orographic_form_drag_stress") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    call trb_mtn_stress_tend(state, pbuf, cam_in)
    call beljaars_drag_tend(state, pbuf, cam_in)

    if (trim(cam_take_snapshot_after) == "orographic_form_drag_stress") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    !------------------------------------------
    ! Call vertical diffusion code (pbl, free atmosphere and molecular)
    !------------------------------------------

    if (trim(cam_take_snapshot_before) == "vertical_diffusion_section") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    call carma_diags_obj%update(cam_in, state, pbuf)

    call vertical_diffusion_tend (ztodt ,state , cam_in, &
         surfric  ,obklen   ,ptend    ,ast    ,pbuf )

    !------------------------------------------
    ! Call major diffusion for extended model
    !------------------------------------------
    if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
       call waccmx_phys_mspd_tend (ztodt    ,state    ,ptend)
    endif

    if ( (trim(cam_take_snapshot_after) == "vertical_diffusion_section") .and. &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
       call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if
    if ( ptend%lu ) then
      call outfld( 'UTEND_VDIFF', ptend%u, pcols, lchnk)
    end if
    if ( ptend%lv ) then
      call outfld( 'VTEND_VDIFF', ptend%v, pcols, lchnk)
    end if

    call carma_diags_obj%output(state, ptend, cam_in, "VDIF", ztodt, pbuf)

    call physics_update(state, ptend, ztodt, tend)

    if (trim(cam_take_snapshot_after) == "vertical_diffusion_section") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    call t_stopf ('vertical_diffusion_tend')

    !===================================================
    ! Rayleigh friction calculation
    !===================================================
    call t_startf('rayleigh_friction')
    if (trim(cam_take_snapshot_before) == "rayleigh_friction_tend") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
            fh2o, surfric, obklen, flx_heat)
    end if

    call physics_ptend_init(ptend, state%psetcols, 'rayleigh friction', ls=.true., lu=.true., lv=.true.)

    ! Initialize ptend variables to zero
    !REMOVECAM - no longer need these when CAM is retired and pcols no longer exists
    ptend%u(:,:) = 0._r8
    ptend%v(:,:) = 0._r8
    ptend%s(:,:) = 0._r8
    !REMOVECAM_END

    call rayleigh_friction_run(pver, ztodt, state%u(:ncol,:), state%v(:ncol,:), ptend%u(:ncol,:),&
         ptend%v(:ncol,:), ptend%s(:ncol,:), errmsg, errflg)
    if (errflg /= 0) call endrun(errmsg)

    if ( (trim(cam_take_snapshot_after) == "rayleigh_friction_tend") .and.      &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
       call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if
    if ( ptend%lu ) then
      call outfld( 'UTEND_RAYLEIGH', ptend%u, pcols, lchnk)
    end if
    if ( ptend%lv ) then
      call outfld( 'VTEND_RAYLEIGH', ptend%v, pcols, lchnk)
    end if
    call physics_update(state, ptend, ztodt, tend)
    if (trim(cam_take_snapshot_after) == "rayleigh_friction_tend") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
            fh2o, surfric, obklen, flx_heat)
    end if
    call t_stopf('rayleigh_friction')

    if (do_clubb_sgs) then
      call check_energy_cam_chng(state, tend, "vdiff", nstep, ztodt, zero, zero, zero, zero)
    else
      call check_energy_cam_chng(state, tend, "vdiff", nstep, ztodt, cam_in%cflx(:,1), zero, &
           zero, cam_in%shf)
    endif

    call check_tracers_chng(state, tracerint, "vdiff", nstep, ztodt, cam_in%cflx)

    !  aerosol dry deposition processes
    call t_startf('aero_drydep')

    if (trim(cam_take_snapshot_before) == "aero_model_drydep") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    call carma_diags_obj%update(cam_in, state, pbuf)

    call aero_model_drydep( state, pbuf, obklen, surfric, cam_in, ztodt, cam_out, ptend )
    if ( (trim(cam_take_snapshot_after) == "aero_model_drydep") .and.         &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
       call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if
    call carma_diags_obj%output(state, ptend, cam_in, "DRYDEPA", ztodt, pbuf)
    call physics_update(state, ptend, ztodt, tend)

   if (trim(cam_take_snapshot_after) == "aero_model_drydep") then
      call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
   end if

    call t_stopf('aero_drydep')

   ! CARMA microphysics
   !
   ! NOTE: This does both the timestep_tend for CARMA aerosols as well as doing the dry
   ! deposition for CARMA aerosols. It needs to follow vertical_diffusion_tend, so that
   ! obklen and surfric have been calculated. It needs to follow aero_model_drydep, so
   ! that cam_out%xxxdryxxx fields have already been set for CAM aerosols and cam_out
   ! can be added to for CARMA aerosols.
   if (carma_do_aerosol) then
     call t_startf('carma_timestep_tend')
     call carma_diags_obj%update(cam_in, state, pbuf)
     call carma_timestep_tend(state, cam_in, cam_out, ptend, ztodt, pbuf, obklen=obklen, ustar=surfric)
     call carma_diags_obj%output(state, ptend, cam_in, "CRTEND", ztodt, pbuf)
     call physics_update(state, ptend, ztodt, tend)

     call check_energy_cam_chng(state, tend, "carma_tend", nstep, ztodt, zero, zero, zero, zero)
     call t_stopf('carma_timestep_tend')
   end if


    !---------------------------------------------------------------------------------
    !   ... enforce charge neutrality
    !---------------------------------------------------------------------------------
    call charge_balance(state, pbuf)

    !===================================================
    ! Gravity wave drag
    !===================================================
    call t_startf('gw_tend')

    if (trim(cam_take_snapshot_before) == "gw_tend") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    call gw_drag_cam_tend(state, pbuf, ztodt, ptend, cam_in, flx_heat)

    if ( (trim(cam_take_snapshot_after) == "gw_tend") .and.                   &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
       call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if
    if ( ptend%lu ) then
      call outfld( 'UTEND_GWDTOT', ptend%u, pcols, lchnk)
    end if
    if ( ptend%lv ) then
      call outfld( 'VTEND_GWDTOT', ptend%v, pcols, lchnk)
    end if
    call physics_update(state, ptend, ztodt, tend)

    if (trim(cam_take_snapshot_after) == "gw_tend") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    ! Check energy integrals
    call check_energy_cam_chng(state, tend, "gwdrag", nstep, ztodt, zero, &
         zero, zero, flx_heat)
    call t_stopf('gw_tend')

    ! QBO relaxation

    if (trim(cam_take_snapshot_before) == "qbo_relax") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    call qbo_relax(state, pbuf, ptend)
    if ( (trim(cam_take_snapshot_after) == "qbo_relax") .and.                 &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
       call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if
    if ( ptend%lu ) then
      call outfld( 'UTEND_QBORLX', ptend%u, pcols, lchnk)
    end if
    if ( ptend%lv ) then
      call outfld( 'VTEND_QBORLX', ptend%v, pcols, lchnk)
    end if
    call physics_update(state, ptend, ztodt, tend)

    if (trim(cam_take_snapshot_after) == "qbo_relax") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    ! Check energy integrals
    call check_energy_cam_chng(state, tend, "qborelax", nstep, ztodt, zero, zero, zero, zero)

    ! Lunar tides
    call lunar_tides_tend( state, ptend )
    if ( ptend%lu ) then
      call outfld( 'UTEND_LUNART', ptend%u, pcols, lchnk)
    end if
    if ( ptend%lv ) then
      call outfld( 'VTEND_LUNART', ptend%v, pcols, lchnk)
    end if
    call physics_update(state, ptend, ztodt, tend)
    ! Check energy integrals
    call check_energy_cam_chng(state, tend, "lunar_tides", nstep, ztodt, zero, zero, zero, zero)

    ! Ion drag calculation
    call t_startf ( 'iondrag' )

    if (trim(cam_take_snapshot_before) == "iondrag_calc_section") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    if ( do_waccm_ions ) then
       call iondrag_calc( lchnk, ncol, state, ptend, pbuf,  ztodt )
    else
       call iondrag_calc( lchnk, ncol, state, ptend)
    endif
    !----------------------------------------------------------------------------
    ! Call ionosphere routines for extended model if mode is set to ionosphere
    !----------------------------------------------------------------------------
    if( waccmx_is('ionosphere') ) then
       call waccmx_phys_ion_elec_temp_tend(state, ptend, pbuf, ztodt)
    endif

    if ( (trim(cam_take_snapshot_after) == "iondrag_calc_section") .and.      &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
       call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if
    if ( ptend%lu ) then
      call outfld( 'UTEND_IONDRG', ptend%u, pcols, lchnk)
    end if
    if ( ptend%lv ) then
      call outfld( 'VTEND_IONDRG', ptend%v, pcols, lchnk)
    end if
    call physics_update(state, ptend, ztodt, tend)

    if (trim(cam_take_snapshot_after) == "iondrag_calc_section") then
       call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
    end if

    call tot_energy_phys(state, 'phAP')
    call tot_energy_phys(state, 'dyAP',vc=vc_dycore)

    !---------------------------------------------------------------------------------
    ! Enforce charge neutrality after O+ change from ionos_tend
    !---------------------------------------------------------------------------------
    if( waccmx_is('ionosphere') ) then
       call charge_balance(state, pbuf)
    endif

    ! Check energy integrals
    call check_energy_cam_chng(state, tend, "iondrag", nstep, ztodt, zero, zero, zero, zero)

    call t_stopf  ( 'iondrag' )

    ! Update Nudging values, if needed
    !----------------------------------
    if((Nudge_Model).and.(Nudge_ON)) then
      call nudging_timestep_tend(state,ptend)
      if ( ptend%lu ) then
        call outfld( 'UTEND_NDG', ptend%u, pcols, lchnk)
      end if
      if ( ptend%lv ) then
        call outfld( 'VTEND_NDG', ptend%v, pcols, lchnk)
      end if
      call physics_update(state,ptend,ztodt,tend)
      call check_energy_cam_chng(state, tend, "nudging", nstep, ztodt, zero, zero, zero, zero)
    endif

    !-------------- Energy budget checks vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv

    ! Save total energy for global fixer in next timestep
    !
    ! This call must be after the last parameterization and call to physics_update
    !
    call pbuf_set_field(pbuf, teout_idx, state%te_cur(:,dyn_te_idx), (/1,itim_old/),(/pcols,1/))

    ! FV: convert dry-type mixing ratios to moist here because physics_dme_adjust
    !     assumes moist. This is done in p_d_coupling for other dynamics. Bundy, Feb 2004.
    moist_mixing_ratio_dycore = dycore_is('LR').or. dycore_is('FV3')
    !
    ! update cp/cv for energy computation based in updated water variables
    !
    call cam_thermo_water_update(state%q(:ncol,:,:), lchnk, ncol, vc_dycore,&
         to_dry_factor=state%pdel(:ncol,:)/state%pdeldry(:ncol,:))

    ! for dry mixing ratio dycore, physics_dme_adjust is called for energy diagnostic purposes only.
    ! So, save off tracers
    if (.not.moist_mixing_ratio_dycore) then
      !
      ! for dry-mixing ratio based dycores dme_adjust takes place in the dynamical core
      !
      ! only compute dme_adjust for diagnostics purposes
      !
      if (thermo_budget_history) then
        tmp_trac(:ncol,:pver,:pcnst) = state%q(:ncol,:pver,:pcnst)
        tmp_pdel(:ncol,:pver)        = state%pdel(:ncol,:pver)
        tmp_ps(:ncol)                = state%ps(:ncol)
        if (trim(cam_take_snapshot_before) == "physics_dme_adjust") then
          call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
               fh2o, surfric, obklen, flx_heat)
        end if
        call physics_dme_adjust(state, tend, qini, totliqini, toticeini, ztodt)
        if (trim(cam_take_snapshot_after) == "physics_dme_adjust") then
          call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
               fh2o, surfric, obklen, flx_heat)
        end if
        call tot_energy_phys(state, 'phAM')
        call tot_energy_phys(state, 'dyAM', vc=vc_dycore)
        ! Restore pre-"physics_dme_adjust" tracers
        state%q(:ncol,:pver,:pcnst) = tmp_trac(:ncol,:pver,:pcnst)
        state%pdel(:ncol,:pver)     = tmp_pdel(:ncol,:pver)
        state%ps(:ncol)             = tmp_ps(:ncol)
      end if
    else
      !
      ! for moist-mixing ratio based dycores
      !
      ! Note: this operation will NOT be reverted with set_wet_to_dry after set_dry_to_wet call
      !
      call set_dry_to_wet(state, convert_cnst_type='dry')

      if (trim(cam_take_snapshot_before) == "physics_dme_adjust") then
        call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,&
                    fh2o, surfric, obklen, flx_heat)
      end if
      call physics_dme_adjust(state, tend, qini, totliqini, toticeini, ztodt)
      if (trim(cam_take_snapshot_after) == "physics_dme_adjust") then
        call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,&
             fh2o, surfric, obklen, flx_heat)
      end if

      call tot_energy_phys(state, 'phAM')
      call tot_energy_phys(state, 'dyAM', vc=vc_dycore)
    endif

    if (vc_dycore == vc_height.or.vc_dycore == vc_dry_pressure) then
      !
      ! MPAS and SE specific scaling of temperature for enforcing energy consistency
      ! (and to make sure that temperature dependent diagnostic tendencies
      !  are computed correctly; e.g. dtcore)
      !
      scaling(1:ncol,:)  = cpairv(:ncol,:,lchnk)/cp_or_cv_dycore(:ncol,:,lchnk)
      state%T(1:ncol,:)  = state%temp_ini(1:ncol,:)+&
           scaling(1:ncol,:)*(state%T(1:ncol,:)-state%temp_ini(1:ncol,:))
      tend%dtdt(:ncol,:) = scaling(:ncol,:)*tend%dtdt(:ncol,:)
      !
      ! else: do nothing for dycores with energy consistent with CAM physics
      !
    end if


    ! store T, U, and V in buffer for use in computing dynamics T-tendency in next timestep
    do k = 1,pver
       dtcore(:ncol,k) = state%t(:ncol,k)
       dqcore(:ncol,k) = state%q(:ncol,k,ixq)
       ducore(:ncol,k) = state%u(:ncol,k)
       dvcore(:ncol,k) = state%v(:ncol,k)
    end do

    !-------------- Energy budget checks ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    if (aqua_planet) then
       labort = .false.
       do i=1,ncol
          if (cam_in%ocnfrac(i) /= 1._r8) then
             labort = .true.
             if (masterproc) write(iulog,*) 'oceanfrac(',i,')=',cam_in%ocnfrac(i)
          end if
       end do
       if (labort) then
          call endrun ('TPHYSAC error: in aquaplanet mode, but grid contains non-ocean point')
       endif
    endif

    call diag_phys_tend_writeout (state, pbuf,  tend, ztodt, qini, cldliqini, cldiceini)

    call clybryiy_fam_set( ncol, lchnk, map2chm, state%q, pbuf )

    ! clean CARMA diagnostics object
    if (associated(carma_diags_obj)) then
       deallocate(carma_diags_obj)
       nullify(carma_diags_obj)
    end if

    ! output these here -- after updates by chem_timestep_tend or export_fields within the current time step
    if (associated(cam_out%nhx_nitrogen_flx)) then
       call outfld('a2x_NHXDEP', cam_out%nhx_nitrogen_flx, pcols, lchnk)
    end if
    if (associated(cam_out%noy_nitrogen_flx)) then
       call outfld('a2x_NOYDEP', cam_out%noy_nitrogen_flx, pcols, lchnk)
    end if

  end subroutine tphysac

  subroutine tphysbc (ztodt, state,  &
       tend,    pbuf,              &
       cam_out, cam_in )
    !-----------------------------------------------------------------------
    !
    ! Purpose:
    ! Evaluate and apply physical processes that are calculated BEFORE
    ! coupling to land, sea, and ice models.
    !
    ! Processes currently included are:
    !
    !  o Resetting Negative Tracers to Positive
    !  o Global Mean Total Energy Fixer
    !  o Dry Adjustment
    !  o Asymmetric Turbulence Scheme : Deep Convection & Shallow Convection
    !  o Stratiform Macro-Microphysics
    !  o Wet Scavenging of Aerosol
    !  o Radiation
    !
    ! Method:
    !
    ! Each parameterization should be implemented with this sequence of calls:
    !  1)  Call physics interface
    !  2)  Check energy
    !  3)  Call physics_update
    ! See Interface to Column Physics and Chemistry Packages
    !   http://www.ccsm.ucar.edu/models/atm-cam/docs/phys-interface/index.html
    !
    !-----------------------------------------------------------------------

    use physics_buffer,  only: physics_buffer_desc, pbuf_get_field
    use physics_buffer,  only: pbuf_get_index, pbuf_old_tim_idx
    use physics_buffer,  only: col_type_subcol, dyn_time_lvls
    use shr_kind_mod,    only: r8 => shr_kind_r8

    use dadadj_cam,      only: dadadj_tend
    use rk_stratiform_cam, only: rk_stratiform_cam_tend
    use microp_driver,   only: microp_driver_tend
    use microp_aero,     only: microp_aero_run
    use macrop_driver,   only: macrop_driver_tend
    use physics_types,   only: physics_state, physics_tend, physics_ptend, &
                               physics_update, physics_ptend_init, physics_ptend_sum, &
                               physics_state_check, physics_ptend_scale, &
                               dyn_te_idx
    use cam_diagnostics, only: diag_conv_tend_ini, diag_phys_writeout, diag_conv, diag_export, diag_state_b4_phys_write
    use cam_diagnostics, only: diag_clip_tend_writeout
    use cam_history,     only: outfld
    use physconst,       only: latvap
    use constituents,    only: pcnst, qmin, cnst_get_ind
    use air_composition, only: thermodynamic_active_species_liq_num,thermodynamic_active_species_liq_idx
    use air_composition, only: thermodynamic_active_species_ice_num,thermodynamic_active_species_ice_idx
    use convect_deep,    only: convect_deep_tend, convect_deep_tend_2, deep_scheme_does_scav_trans
    use time_manager,    only: is_first_step, get_nstep
    use convect_shallow, only: convect_shallow_tend
    use check_energy,    only: check_energy_timestep_init, check_energy_cam_chng
    use check_energy,    only: check_energy_cam_fix
    use check_energy,    only: check_tracers_data, check_tracers_init, check_tracers_chng
    use check_energy,    only: tot_energy_phys
    use dycore,          only: dycore_is
    use aero_model,      only: aero_model_wetdep
    use aero_wetdep_cam, only: wetdep_lq
    use carma_intr,      only: carma_wetdep_tend, carma_timestep_tend
    use carma_flags_mod, only: carma_do_detrain, carma_do_cldice, carma_do_cldliq,  carma_do_wetdep
    use radiation,       only: radiation_tend
    use cloud_diagnostics, only: cloud_diagnostics_calc
    use perf_mod
    use mo_gas_phase_chemdr,only: map2chm
    use clybryiy_fam,       only: clybryiy_fam_adj
    use clubb_intr,      only: clubb_tend_cam
    use sslt_rebin,      only: sslt_rebin_adv
    use tropopause,      only: tropopause_output
    use cam_abortutils,  only: endrun
    use subcol,          only: subcol_gen, subcol_ptend_avg
    use subcol_utils,    only: subcol_ptend_copy, is_subcol_on
    use qneg_module,     only: qneg3
    use subcol_SILHS,    only: subcol_SILHS_var_covar_driver, init_state_subcol
    use subcol_SILHS,    only: subcol_SILHS_fill_holes_conserv
    use subcol_SILHS,    only: subcol_SILHS_hydromet_conc_tend_lim
    use micro_pumas_cam,    only: massless_droplet_destroyer
    use cam_snapshot,    only: cam_snapshot_all_outfld_tphysbc
    use cam_snapshot_common, only: cam_snapshot_ptend_outfld
    use ssatcontrail,       only: ssatcontrail_d0
    use dyn_tests_utils, only: vc_dycore
    use surface_emissions_mod,only: surface_emissions_set
    use elevated_emissions_mod,only: elevated_emissions_set

    ! Arguments

    real(r8), intent(in) :: ztodt                          ! 2 delta t (model time increment)

    type(physics_state), intent(inout) :: state
    type(physics_tend ), intent(inout) :: tend
    type(physics_buffer_desc), pointer :: pbuf(:)

    type(cam_out_t),     intent(inout) :: cam_out
    type(cam_in_t),      intent(in)    :: cam_in


    !
    !---------------------------Local workspace-----------------------------
    !

    type(physics_ptend)   :: ptend            ! indivdual parameterization tendencies
    type(physics_ptend)   :: ptend_macp_all   ! sum of macrophysics tendencies (e.g. CLUBB) over substeps
    type(physics_state)   :: state_sc         ! state for sub-columns
    type(physics_ptend)   :: ptend_sc         ! ptend for sub-columns
    type(physics_ptend)   :: ptend_aero       ! ptend for microp_aero
    type(physics_ptend)   :: ptend_aero_sc    ! ptend for microp_aero on sub-columns
    type(physics_tend)    :: tend_sc          ! tend for sub-columns

    integer :: nstep                          ! current timestep number

    real(r8) :: net_flx(pcols)

    real(r8) :: zdu(pcols,pver)               ! detraining mass flux from deep convection
    real(r8) :: cmfmc(pcols,pverp)            ! Convective mass flux--m sub c

    real(r8) cmfcme(pcols,pver)                ! cmf condensation - evaporation

    real(r8) dlf(pcols,pver)                   ! Detraining cld H20 from shallow + deep convections
    real(r8) dlf2(pcols,pver)                  ! Detraining cld H20 from shallow convections
    real(r8) rtdt                              ! 1./ztodt

    integer lchnk                              ! chunk identifier
    integer ncol                               ! number of atmospheric columns

    integer :: i                               ! column indicex
    integer :: ixcldice, ixcldliq, ixq         ! constituent indices for cloud liquid and ice water.
    integer :: m, m_cnst
    ! for macro/micro co-substepping
    integer :: macmic_it                      ! iteration variables
    real(r8) :: cld_macmic_ztodt              ! modified timestep
    ! physics buffer fields to compute tendencies for stratiform package
    integer itim_old, ifld
    real(r8), pointer, dimension(:,:) :: cld        ! cloud fraction


    ! physics buffer fields for total energy and mass adjustment
    real(r8), pointer, dimension(:  ) :: teout
    real(r8), pointer, dimension(:,:) :: qini
    real(r8), pointer, dimension(:,:) :: cldliqini
    real(r8), pointer, dimension(:,:) :: cldiceini
    real(r8), pointer, dimension(:,:) :: totliqini
    real(r8), pointer, dimension(:,:) :: toticeini
    real(r8), pointer, dimension(:,:) :: dtcore
    real(r8), pointer, dimension(:,:) :: dqcore
    real(r8), pointer, dimension(:,:) :: ducore
    real(r8), pointer, dimension(:,:) :: dvcore

    real(r8), pointer, dimension(:,:,:) :: fracis  ! fraction of transported species that are insoluble

    real(r8), pointer :: dlfzm(:,:)                ! ZM detrained convective cloud water mixing ratio.

    ! convective precipitation variables
    real(r8),pointer :: prec_dp(:)                ! total precipitation from ZM convection
    real(r8),pointer :: snow_dp(:)                ! snow from ZM convection
    real(r8),pointer :: prec_sh(:)                ! total precipitation from Hack convection
    real(r8),pointer :: snow_sh(:)                ! snow from Hack convection

    ! carma precipitation variables
    real(r8) :: prec_sed_carma(pcols)          ! total precip from cloud sedimentation (CARMA)
    real(r8) :: snow_sed_carma(pcols)          ! snow from cloud ice sedimentation (CARMA)

    ! stratiform precipitation variables
    real(r8),pointer :: prec_str(:)    ! sfc flux of precip from stratiform (m/s)
    real(r8),pointer :: snow_str(:)     ! sfc flux of snow from stratiform   (m/s)
    real(r8),pointer :: prec_str_sc(:)  ! sfc flux of precip from stratiform (m/s) -- for subcolumns
    real(r8),pointer :: snow_str_sc(:)  ! sfc flux of snow from stratiform   (m/s) -- for subcolumns
    real(r8),pointer :: prec_pcw(:)     ! total precip from prognostic cloud scheme
    real(r8),pointer :: snow_pcw(:)     ! snow from prognostic cloud scheme
    real(r8),pointer :: prec_sed(:)     ! total precip from cloud sedimentation
    real(r8),pointer :: snow_sed(:)     ! snow from cloud ice sedimentation

    ! Local copies for substepping
    real(r8) :: prec_pcw_macmic(pcols)
    real(r8) :: snow_pcw_macmic(pcols)
    real(r8) :: prec_sed_macmic(pcols)
    real(r8) :: snow_sed_macmic(pcols)

    ! energy checking variables
    real(r8) :: zero(pcols)                    ! array of zeros
    real(r8) :: zero_sc(pcols*psubcols)        ! array of zeros
    real(r8) :: rliq(pcols)                    ! vertical integral of liquid not yet in q(ixcldliq)
    real(r8) :: rice(pcols)                    ! vertical integral of ice not yet in q(ixcldice)
    real(r8) :: rliq2(pcols)                   ! vertical integral of liquid from shallow scheme
    real(r8) :: det_s  (pcols)                 ! vertical integral of detrained static energy from ice
    real(r8) :: det_ice(pcols)                 ! vertical integral of detrained ice
    real(r8) :: flx_cnd(pcols)
    real(r8) :: flx_heat(pcols)
    type(check_tracers_data):: tracerint             ! energy integrals and cummulative boundary fluxes
    real(r8) :: zero_tracers(pcols,pcnst)

    ! For aerosol budget diagnostics
    character(len=16) :: pname      !! package name
    type(carma_diags_t), pointer :: carma_diags_obj

    !-----------------------------------------------------------------------
    carma_diags_obj => carma_diags_t()
    if (.not.associated(carma_diags_obj)) then
       call endrun('tphysbc: carma_diags_obj allocation failed')
    end if

    call t_startf('bc_init')

    zero = 0._r8
    zero_tracers(:,:) = 0._r8
    zero_sc(:) = 0._r8

    lchnk = state%lchnk
    ncol  = state%ncol

    rtdt = 1._r8/ztodt

    nstep = get_nstep()

    ! Associate pointers with physics buffer fields
    itim_old = pbuf_old_tim_idx()
    ifld = pbuf_get_index('CLD')
    call pbuf_get_field(pbuf, ifld, cld, (/1,1,itim_old/),(/pcols,pver,1/))

    call pbuf_get_field(pbuf, teout_idx, teout, (/1,itim_old/), (/pcols,1/))

    call pbuf_get_field(pbuf, qini_idx, qini)
    call pbuf_get_field(pbuf, cldliqini_idx, cldliqini)
    call pbuf_get_field(pbuf, cldiceini_idx, cldiceini)
    call pbuf_get_field(pbuf, totliqini_idx, totliqini)
    call pbuf_get_field(pbuf, toticeini_idx, toticeini)

    call pbuf_get_field(pbuf, dtcore_idx, dtcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
    call pbuf_get_field(pbuf, dqcore_idx, dqcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
    call pbuf_get_field(pbuf, ducore_idx, ducore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
    call pbuf_get_field(pbuf, dvcore_idx, dvcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )

    ifld    = pbuf_get_index('FRACIS')
    call pbuf_get_field(pbuf, ifld, fracis, start=(/1,1,1/), kount=(/pcols, pver, pcnst/)  )
    fracis (:ncol,:,1:pcnst) = 1._r8

    ! Set physics tendencies to 0
    tend %dTdt(:ncol,:pver)  = 0._r8
    tend %dudt(:ncol,:pver)  = 0._r8
    tend %dvdt(:ncol,:pver)  = 0._r8

    ! Verify state coming from the dynamics
    if (state_debug_checks) &
         call physics_state_check(state, name="before tphysbc (dycore?)")

    call clybryiy_fam_adj( ncol, lchnk, map2chm, state%q, pbuf )

    ! Since clybryiy_fam_adj operates directly on the tracers, and has no
    ! physics_update call, re-run qneg3.
    call qneg3('TPHYSBCc',lchnk  ,ncol    ,pcols   ,pver    , &
         1, pcnst, qmin  ,state%q )

    ! Validate output of clybryiy_fam_adj.
    if (state_debug_checks) &
         call physics_state_check(state, name="clybryiy_fam_adj")

    !
    ! Dump out "before physics" state
    !
    call diag_state_b4_phys_write (state)

    ! compute mass integrals of input tracers state
    call check_tracers_init(state, tracerint)

    call t_stopf('bc_init')

    ! Zero-initialize subroutine-level variables for snapshot
    cmfmc(:,:) = 0._r8
    cmfcme(:,:) = 0._r8
    zdu(:,:) = 0._r8
    rliq(:) = 0._r8
    rice(:) = 0._r8
    dlf(:,:) = 0._r8
    dlf2(:,:) = 0._r8
    rliq2(:) = 0._r8
    det_s(:) = 0._r8
    det_ice(:) = 0._r8
    net_flx(:) = 0._r8

    !===================================================
    ! Global mean total energy fixer
    !===================================================
    call t_startf('energy_fixer')

    call tot_energy_phys(state, 'phBF')
    call tot_energy_phys(state, 'dyBF',vc=vc_dycore)

    call check_energy_cam_fix(state, ptend, nstep, flx_heat)
    call physics_update(state, ptend, ztodt, tend)
    call check_energy_cam_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat)
    call outfld( 'EFIX', flx_heat    , pcols, lchnk   )

    call tot_energy_phys(state, 'phBP')
    call tot_energy_phys(state, 'dyBP',vc=vc_dycore)
    ! Save state for convective tendency calculations.
    call diag_conv_tend_ini(state, pbuf)

    call cnst_get_ind('Q', ixq)
    call cnst_get_ind('CLDLIQ', ixcldliq)
    call cnst_get_ind('CLDICE', ixcldice)
    qini     (:ncol,:pver) = state%q(:ncol,:pver,       1)
    cldliqini(:ncol,:pver) = state%q(:ncol,:pver,ixcldliq)
    cldiceini(:ncol,:pver) = state%q(:ncol,:pver,ixcldice)

    totliqini(:ncol,:pver) = 0.0_r8
    do m_cnst=1,thermodynamic_active_species_liq_num
      m = thermodynamic_active_species_liq_idx(m_cnst)
      totliqini(:ncol,:pver) = totliqini(:ncol,:pver)+state%q(:ncol,:pver,m)
    end do
    toticeini(:ncol,:pver) = 0.0_r8
    do m_cnst=1,thermodynamic_active_species_ice_num
      m = thermodynamic_active_species_ice_idx(m_cnst)
      toticeini(:ncol,:pver) = toticeini(:ncol,:pver)+state%q(:ncol,:pver,m)
    end do


    call outfld('TEOUT', teout       , pcols, lchnk   )
    call outfld('TEINP', state%te_ini(:,dyn_te_idx), pcols, lchnk   )
    call outfld('TEFIX', state%te_cur(:,dyn_te_idx), pcols, lchnk   )

    ! T, U, V tendency due to dynamics
    if( nstep > dyn_time_lvls-1 ) then
       dtcore(:ncol,:pver) = (state%t(:ncol,:pver) - dtcore(:ncol,:pver))/ztodt
       dqcore(:ncol,:pver) = (state%q(:ncol,:pver,ixq) - dqcore(:ncol,:pver))/ztodt
       ducore(:ncol,:pver) = (state%u(:ncol,:pver) - ducore(:ncol,:pver))/ztodt
       dvcore(:ncol,:pver) = (state%v(:ncol,:pver) - dvcore(:ncol,:pver))/ztodt
       call outfld( 'DTCORE', dtcore, pcols, lchnk )
       call outfld( 'DQCORE', dqcore, pcols, lchnk )
       call outfld( 'UTEND_CORE', ducore, pcols, lchnk )
       call outfld( 'VTEND_CORE', dvcore, pcols, lchnk )
    end if

    call t_stopf('energy_fixer')

    call surface_emissions_set( lchnk, ncol, pbuf )
    call elevated_emissions_set( lchnk, ncol, pbuf )

    !
    !===================================================
    ! Dry adjustment
    !===================================================
    call t_startf('dry_adjustment')

    if (trim(cam_take_snapshot_before) == "dadadj_tend") then
       call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, &
           flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
    end if

    call dadadj_tend(ztodt, state, ptend)

    if ( (trim(cam_take_snapshot_after) == "dadadj_tend") .and. &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
            call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if
    call physics_update(state, ptend, ztodt, tend)

    if (trim(cam_take_snapshot_after) == "dadadj_tend") then
       call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, &
           flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
    end if

    call t_stopf('dry_adjustment')

    !===================================================
    ! Moist convection
    !===================================================
    call t_startf('moist_convection')

    call t_startf ('convect_deep_tend')

    if (trim(cam_take_snapshot_before) == "convect_deep_tend") then
       call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, &
           flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
    end if

    call convect_deep_tend(  &
         cmfmc,      cmfcme,             &
         zdu,       &
         rliq,    rice,      &
         ztodt,   &
         state,   ptend, cam_in%landfrac, pbuf)

    if ( (trim(cam_take_snapshot_after) == "convect_deep_tend") .and. &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
            call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if

    if ( ptend%lu ) then
      call outfld( 'UTEND_DCONV', ptend%u, pcols, lchnk)
    end if
    if ( ptend%lv ) then
      call outfld( 'VTEND_DCONV', ptend%v, pcols, lchnk)
    end if
    call physics_update(state, ptend, ztodt, tend)

    if (trim(cam_take_snapshot_after) == "convect_deep_tend") then
       call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, &
           flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
    end if

    call t_stopf('convect_deep_tend')

    call pbuf_get_field(pbuf, prec_dp_idx, prec_dp )
    call pbuf_get_field(pbuf, snow_dp_idx, snow_dp )
    call pbuf_get_field(pbuf, prec_sh_idx, prec_sh )
    call pbuf_get_field(pbuf, snow_sh_idx, snow_sh )
    call pbuf_get_field(pbuf, prec_str_idx, prec_str )
    call pbuf_get_field(pbuf, snow_str_idx, snow_str )
    call pbuf_get_field(pbuf, prec_sed_idx, prec_sed )
    call pbuf_get_field(pbuf, snow_sed_idx, snow_sed )
    call pbuf_get_field(pbuf, prec_pcw_idx, prec_pcw )
    call pbuf_get_field(pbuf, snow_pcw_idx, snow_pcw )

    if (use_subcol_microp) then
      call pbuf_get_field(pbuf, prec_str_idx, prec_str_sc, col_type=col_type_subcol)
      call pbuf_get_field(pbuf, snow_str_idx, snow_str_sc, col_type=col_type_subcol)
    end if

    ! Check energy integrals, including "reserved liquid"
    flx_cnd(:ncol) = prec_dp(:ncol) + rliq(:ncol)
    snow_dp(:ncol) = snow_dp(:ncol) + rice(:ncol)
    call check_energy_cam_chng(state, tend, "convect_deep", nstep, ztodt, zero, flx_cnd, snow_dp, zero)
    snow_dp(:ncol) = snow_dp(:ncol) - rice(:ncol)

    !
    ! Call Hack (1994) convection scheme to deal with shallow/mid-level convection
    !
    call t_startf ('convect_shallow_tend')

    if (dlfzm_idx > 0) then
       call pbuf_get_field(pbuf, dlfzm_idx, dlfzm)
       dlf(:ncol,:) = dlfzm(:ncol,:)
    else
       dlf(:,:) = 0._r8
    end if

    if (trim(cam_take_snapshot_before) == "convect_shallow_tend") then
       call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, &
           flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
    end if

    call convect_shallow_tend (ztodt   , cmfmc, &
         dlf        , dlf2   ,  rliq   , rliq2, &
         state      , ptend  ,  pbuf, cam_in)
    call t_stopf ('convect_shallow_tend')

    call physics_update(state, ptend, ztodt, tend)

    if ( (trim(cam_take_snapshot_after) == "convect_shallow_tend") .and. &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
            call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if
    if ( ptend%lu ) then
      call outfld( 'UTEND_SHCONV', ptend%u, pcols, lchnk)
    end if
    if ( ptend%lv ) then
      call outfld( 'VTEND_SHCONV', ptend%v, pcols, lchnk)
    end if
    call physics_update(state, ptend, ztodt, tend)

    if (trim(cam_take_snapshot_after) == "convect_shallow_tend") then
       call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, &
           flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
    end if

    flx_cnd(:ncol) = prec_sh(:ncol) + rliq2(:ncol)
    call check_energy_cam_chng(state, tend, "convect_shallow", nstep, ztodt, zero, flx_cnd, snow_sh, zero)

    call check_tracers_chng(state, tracerint, "convect_shallow", nstep, ztodt, zero_tracers)

    call t_stopf('moist_convection')

    ! Rebin the 4-bin version of sea salt into bins for coarse and accumulation
    ! modes that correspond to the available optics data.  This is only necessary
    ! for CAM-RT.  But it's done here so that the microphysics code which is called
    ! from the stratiform interface has access to the same aerosols as the radiation
    ! code.
    call sslt_rebin_adv(pbuf,  state)

    !===================================================
    ! Calculate tendencies from CARMA bin microphysics.
    !===================================================
    !
    ! If CARMA is doing detrainment, then on output, rliq no longer represents water reserved
    ! for detrainment, but instead represents potential snow fall. The mass and number of the
    ! snow are stored in the physics buffer and will be incorporated by the MG microphysics.
    !
    ! Currently CARMA cloud microphysics is only supported with the MG microphysics.
    call t_startf('carma_timestep_tend')

    if (carma_do_cldice .or. carma_do_cldliq) then
       call carma_diags_obj%update(cam_in, state, pbuf)
       call carma_timestep_tend(state, cam_in, cam_out, ptend, ztodt, pbuf, dlf=dlf, rliq=rliq, &
            prec_str=prec_str, snow_str=snow_str, prec_sed=prec_sed_carma, snow_sed=snow_sed_carma)
       call carma_diags_obj%output(state, ptend, cam_in, "CRTEND", ztodt, pbuf)
       call physics_update(state, ptend, ztodt, tend)

       ! Before the detrainment, the reserved condensate is all liquid, but if CARMA is doing
       ! detrainment, then the reserved condensate is snow.
       if (carma_do_detrain) then
          call check_energy_cam_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str+rliq, snow_str+rliq, zero)
       else
          call check_energy_cam_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str, snow_str, zero)
       end if
    end if

    call t_stopf('carma_timestep_tend')

    if( microp_scheme == 'RK' ) then

       !===================================================
       ! Calculate stratiform tendency (sedimentation, detrain, cloud fraction and microphysics )
       !===================================================
       call t_startf('rk_stratiform_tend')

       if (trim(cam_take_snapshot_before) == "rk_stratiform_tend") then
            call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, &
                 flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
       end if

       call rk_stratiform_cam_tend(state, ptend, pbuf, ztodt, &
            cam_in%icefrac, cam_in%landfrac, cam_in%ocnfrac, &
            cam_in%snowhland, & ! sediment
            dlf, dlf2, & ! detrain
            rliq  , & ! check energy after detrain
            cmfmc,  &
            cam_in%ts,      cam_in%sst,        zdu)

       if ( (trim(cam_take_snapshot_after) == "rk_stratiform_tend") .and. &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
            call cam_snapshot_ptend_outfld(ptend, lchnk)
       end if

       call physics_update(state, ptend, ztodt, tend)

       if (trim(cam_take_snapshot_after) == "rk_stratiform_tend") then
           call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, &
               flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
       end if

       call check_energy_cam_chng(state, tend, "cldwat_tend", nstep, ztodt, zero, prec_str, snow_str, zero)

       call t_stopf('rk_stratiform_tend')

    elseif( microp_scheme == 'MG' ) then
       ! Start co-substepping of macrophysics and microphysics
       cld_macmic_ztodt = ztodt/cld_macmic_num_steps

       ! Clear precip fields that should accumulate.
       prec_sed_macmic = 0._r8
       snow_sed_macmic = 0._r8
       prec_pcw_macmic = 0._r8
       snow_pcw_macmic = 0._r8

       ! contrail parameterization
       ! see Chen et al., 2012: Global contrail coverage simulated
       !                        by CAM5 with the inventory of 2006 global aircraft emissions, JAMES
       !                        https://doi.org/10.1029/2011MS000105
       call ssatcontrail_d0(state, pbuf, ztodt, ptend)
       call physics_update(state, ptend, ztodt, tend)

       ! initialize ptend structures where macro and microphysics tendencies are
       ! accumulated over macmic substeps
       call physics_ptend_init(ptend_macp_all,state%psetcols,'macrophysics',lu=.true.,lv=.true.)

       do macmic_it = 1, cld_macmic_num_steps

          !===================================================
          ! Calculate macrophysical tendency (sedimentation, detrain, cloud fraction)
          !===================================================

          call t_startf('macrop_tend')

          ! don't call Park macrophysics if CLUBB is called
          if (macrop_scheme .ne. 'CLUBB_SGS') then

             if (trim(cam_take_snapshot_before) == "macrop_driver_tend") then
                call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, &
                     flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
             end if

             call macrop_driver_tend( &
                  state,           ptend,          cld_macmic_ztodt, &
                  cam_in%landfrac, cam_in%ocnfrac, cam_in%snowhland, & ! sediment
                  dlf,             dlf2,                             & ! detrain
                  cmfmc,                                             &
                  cam_in%ts,       cam_in%sst,     zdu,              &
                  pbuf,            det_s,          det_ice)

             ! Since we "added" the reserved liquid back in this routine, we need
             ! to account for it in the energy checker
             flx_cnd(:ncol) = -1._r8*rliq(:ncol)
             flx_heat(:ncol) = det_s(:ncol)

             ! Unfortunately, physics_update does not know what time period
             ! "tend" is supposed to cover, and therefore can't update it
             ! with substeps correctly. For now, work around this by scaling
             ! ptend down by the number of substeps, then applying it for
             ! the full time (ztodt).
             call physics_ptend_scale(ptend, 1._r8/cld_macmic_num_steps, ncol)
             if ( (trim(cam_take_snapshot_after) == "macrop_driver_tend") .and. &
                  (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
                call cam_snapshot_ptend_outfld(ptend, lchnk)
             end if
             call physics_ptend_sum(ptend,ptend_macp_all,ncol)
             call physics_update(state, ptend, ztodt, tend)

             if (trim(cam_take_snapshot_after) == "macrop_driver_tend") then
                call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, &
                     flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
             end if

             call check_energy_cam_chng(state, tend, "macrop_tend", nstep, ztodt, &
                  zero, flx_cnd(:ncol)/cld_macmic_num_steps, &
                  det_ice(:ncol)/cld_macmic_num_steps, &
                  flx_heat(:ncol)/cld_macmic_num_steps)

          else ! Calculate CLUBB macrophysics

             ! =====================================================
             !    CLUBB call (PBL, shallow convection, macrophysics)
             ! =====================================================

             if (trim(cam_take_snapshot_before) == "clubb_tend_cam") then
                call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, &
                     flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
             end if

             call clubb_tend_cam(state, ptend, pbuf, cld_macmic_ztodt,&
                cmfmc, cam_in, macmic_it, cld_macmic_num_steps, &
                dlf, det_s, det_ice)

             ! Since we "added" the reserved liquid back in this routine, we need
             ! to account for it in the energy checker
             flx_cnd(:ncol) = -1._r8*rliq(:ncol)
             flx_heat(:ncol) = cam_in%shf(:ncol) + det_s(:ncol)

             ! Unfortunately, physics_update does not know what time period
             ! "tend" is supposed to cover, and therefore can't update it
             ! with substeps correctly. For now, work around this by scaling
             ! ptend down by the number of substeps, then applying it for
             ! the full time (ztodt).
             call physics_ptend_scale(ptend, 1._r8/cld_macmic_num_steps, ncol)

             ! Update physics tendencies and copy state to state_eq, because that is
             ! input for microphysics
             if ( (trim(cam_take_snapshot_after) == "clubb_tend_cam") .and.   &
                  (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
                call cam_snapshot_ptend_outfld(ptend, lchnk)
             end if
             call physics_ptend_sum(ptend,ptend_macp_all,ncol)
             call physics_update(state, ptend, ztodt, tend)

             if (trim(cam_take_snapshot_after) == "clubb_tend_cam") then
                call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, &
                      flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
             end if

             ! Use actual qflux (not lhf/latvap) for consistency with surface fluxes and revised code
             call check_energy_cam_chng(state, tend, "clubb_tend", nstep, ztodt, &
                cam_in%cflx(:ncol,1)/cld_macmic_num_steps, &
                flx_cnd(:ncol)/cld_macmic_num_steps, &
                det_ice(:ncol)/cld_macmic_num_steps, &
                flx_heat(:ncol)/cld_macmic_num_steps)

          endif

          call t_stopf('macrop_tend')

          !===================================================
          ! Calculate cloud microphysics
          !===================================================

          if (is_subcol_on() .neqv. use_subcol_microp ) then
            call endrun("Error calculating cloud microphysics: is_subcol_on() != use_subcol_microp")
          end if

          if (is_subcol_on()) then
             ! Allocate sub-column structures.
             call physics_state_alloc(state_sc, lchnk, psubcols*pcols)
             call physics_tend_alloc(tend_sc, psubcols*pcols)

             ! Generate sub-columns using the requested scheme
             if (trim(subcol_scheme) == 'SILHS') call init_state_subcol(state, tend, state_sc, tend_sc)
             call subcol_gen(state, tend, state_sc, tend_sc, pbuf)

             !Initialize check energy for subcolumns
             call check_energy_timestep_init(state_sc, tend_sc, pbuf, col_type_subcol)
          end if

          if (trim(cam_take_snapshot_before) == "microp_section") then
             call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, &
                  flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
          end if

          call carma_diags_obj%update(cam_in, state, pbuf)

          call t_startf('microp_aero_run')
          call microp_aero_run(state, ptend_aero, cld_macmic_ztodt, pbuf)
          call t_stopf('microp_aero_run')

          call t_startf('microp_tend')

          if (use_subcol_microp) then

             if (trim(cam_take_snapshot_before) == "microp_driver_tend_subcol") then
                call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state_sc, tend_sc, cam_in, cam_out, pbuf, &
                     flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
             end if

             call microp_driver_tend(state_sc, ptend_sc, cld_macmic_ztodt, pbuf)
             ! Parameterize subcolumn effects on covariances, if enabled
             if (trim(subcol_scheme) == 'SILHS') &
                call subcol_SILHS_var_covar_driver( cld_macmic_ztodt, state_sc, ptend_sc, pbuf )

             ! Average the sub-column ptend for use in gridded update - will not contain ptend_aero
             call subcol_ptend_avg(ptend_sc, state_sc%ngrdcol, lchnk, ptend)

             ! Call the conservative hole filler.
             ! Hole filling is only necessary when using subcolumns.
             ! Note:  this needs to be called after subcol_ptend_avg but before
             !        physics_ptend_scale.
             if (trim(subcol_scheme) == 'SILHS') &
                call subcol_SILHS_fill_holes_conserv( state, cld_macmic_ztodt, &
                                                      ptend, pbuf )

             ! Destroy massless droplets - Note this routine returns with no change unless
             ! micro_do_massless_droplet_destroyer has been set to true
             call massless_droplet_destroyer( cld_macmic_ztodt, state, & ! Intent(in)
                                              ptend )                    ! Intent(inout)

             ! Limit the value of hydrometeor concentrations in order to place
             ! reasonable limits on hydrometeor drop size and keep them from
             ! becoming too large.
             ! Note:  this needs to be called after hydrometeor mixing ratio
             !        tendencies are adjusted by subcol_SILHS_fill_holes_conserv
             !        and after massless drop concentrations are removed by the
             !        subcol_SILHS_massless_droplet_destroyer, but before the
             !        call to physics_ptend_scale.
             if (trim(subcol_scheme) == 'SILHS') &
                call subcol_SILHS_hydromet_conc_tend_lim( state, cld_macmic_ztodt, ptend )

             ! Copy ptend_aero field to one dimensioned by sub-columns before summing with ptend
             call subcol_ptend_copy(ptend_aero, state_sc, ptend_aero_sc)
             call physics_ptend_sum(ptend_aero_sc, ptend_sc, state_sc%ncol)
             call physics_ptend_dealloc(ptend_aero_sc)

             ! Have to scale and apply for full timestep to get tend right
             ! (see above note for macrophysics).
             call physics_ptend_scale(ptend_sc, 1._r8/cld_macmic_num_steps, ncol)

             if ( (trim(cam_take_snapshot_after) == "microp_driver_tend_subcol") .and. &
                  (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
                call cam_snapshot_ptend_outfld(ptend, lchnk)
             end if
             call physics_update (state_sc, ptend_sc, ztodt, tend_sc)

             if (trim(cam_take_snapshot_after) == "microp_driver_tend_subcol") then
                call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state_sc, tend_sc, cam_in, cam_out, pbuf, &
                   flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
             end if

             call check_energy_cam_chng(state_sc, tend_sc, "microp_tend_subcol", &
                  nstep, ztodt, zero_sc, &
                  prec_str_sc(:state_sc%ncol)/cld_macmic_num_steps, &
                  snow_str_sc(:state_sc%ncol)/cld_macmic_num_steps, zero_sc)

             call physics_state_dealloc(state_sc)
             call physics_tend_dealloc(tend_sc)
             call physics_ptend_dealloc(ptend_sc)
          else
             call microp_driver_tend(state, ptend, cld_macmic_ztodt, pbuf)
          end if
          ! combine aero and micro tendencies for the grid
          call physics_ptend_sum(ptend_aero, ptend, ncol)
          call physics_ptend_dealloc(ptend_aero)

          ! These need to be reported before the scaling as they are based
          ! on the substep size not ztodt.
          write(pname, '(A, I2.2)') "MICROP", macmic_it
          call carma_diags_obj%output(state, ptend, cam_in, pname, ztodt/cld_macmic_num_steps, pbuf)

          ! Have to scale and apply for full timestep to get tend right
          ! (see above note for macrophysics).
          call physics_ptend_scale(ptend, 1._r8/cld_macmic_num_steps, ncol)

          call diag_clip_tend_writeout(state, ptend, ncol, lchnk, ixcldliq, ixcldice, ixq, ztodt, rtdt)

          if ( (trim(cam_take_snapshot_after) == "microp_section") .and.      &
               (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
             call cam_snapshot_ptend_outfld(ptend, lchnk)
          end if
          call physics_update (state, ptend, ztodt, tend)

          if (trim(cam_take_snapshot_after) == "microp_section") then
             call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, &
                  flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
          end if

          call check_energy_cam_chng(state, tend, "microp_tend", nstep, ztodt, &
               zero, prec_str(:ncol)/cld_macmic_num_steps, &
               snow_str(:ncol)/cld_macmic_num_steps, zero)

          call t_stopf('microp_tend')
          prec_sed_macmic(:ncol) = prec_sed_macmic(:ncol) + prec_sed(:ncol)
          snow_sed_macmic(:ncol) = snow_sed_macmic(:ncol) + snow_sed(:ncol)
          prec_pcw_macmic(:ncol) = prec_pcw_macmic(:ncol) + prec_pcw(:ncol)
          snow_pcw_macmic(:ncol) = snow_pcw_macmic(:ncol) + snow_pcw(:ncol)

       end do ! end substepping over macrophysics/microphysics

       call outfld( 'UTEND_MACROP', ptend_macp_all%u, pcols, lchnk)
       call outfld( 'VTEND_MACROP', ptend_macp_all%v, pcols, lchnk)
       call physics_ptend_dealloc(ptend_macp_all)

       prec_sed(:ncol) = prec_sed_macmic(:ncol)/cld_macmic_num_steps
       snow_sed(:ncol) = snow_sed_macmic(:ncol)/cld_macmic_num_steps
       prec_pcw(:ncol) = prec_pcw_macmic(:ncol)/cld_macmic_num_steps
       snow_pcw(:ncol) = snow_pcw_macmic(:ncol)/cld_macmic_num_steps
       prec_str(:ncol) = prec_pcw(:ncol) + prec_sed(:ncol)
       snow_str(:ncol) = snow_pcw(:ncol) + snow_sed(:ncol)

    endif

    ! Add the precipitation from CARMA to the precipitation from stratiform.
    if (carma_do_cldice .or. carma_do_cldliq) then
       prec_sed(:ncol) = prec_sed(:ncol) + prec_sed_carma(:ncol)
       snow_sed(:ncol) = snow_sed(:ncol) + snow_sed_carma(:ncol)
    end if

    if ( .not. deep_scheme_does_scav_trans() ) then

       ! -------------------------------------------------------------------------------
       ! 1. Wet Scavenging of Aerosols by Convective and Stratiform Precipitation.
       ! 2. Convective Transport of Non-Water Aerosol Species.
       !
       !  . Aerosol wet chemistry determines scavenging fractions, and transformations
       !  . Then do convective transport of all trace species except qv,ql,qi.
       !  . We needed to do the scavenging first to determine the interstitial fraction.
       ! -------------------------------------------------------------------------------

       call t_startf('aerosol_wet_processes')
       if (clim_modal_aero) then
          if (prog_modal_aero) then
             call physics_ptend_init(ptend, state%psetcols, 'aero_water_uptake', lq=wetdep_lq)
             ! Do calculations of mode radius and water uptake if:
             ! 1) modal aerosols are affecting the climate, or
             ! 2) prognostic modal aerosols are enabled
             call modal_aero_calcsize_sub(state, ptend, ztodt, pbuf)
             ! for prognostic modal aerosols the transfer of mass between aitken and accumulation
             ! modes is done in conjunction with the dry radius calculation
             call modal_aero_wateruptake_dr(state, pbuf)
             call physics_update(state, ptend, ztodt, tend)
          else
             call modal_aero_calcsize_diag(state, pbuf)
             call modal_aero_wateruptake_dr(state, pbuf)
          endif
       endif

       if (trim(cam_take_snapshot_before) == "aero_model_wetdep") then
          call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, &
                  flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
       end if

       call carma_diags_obj%update(cam_in, state, pbuf)

       call aero_model_wetdep( state, ztodt, dlf, cam_out, ptend, pbuf)
       if ( (trim(cam_take_snapshot_after) == "aero_model_wetdep") .and.      &
            (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
          call cam_snapshot_ptend_outfld(ptend, lchnk)
       end if
       call carma_diags_obj%output(state, ptend, cam_in, "WETDEPA", ztodt, pbuf)
       call physics_update(state, ptend, ztodt, tend)

       if (trim(cam_take_snapshot_after) == "aero_model_wetdep") then
          call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, &
                  flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
       end if

       if (carma_do_wetdep) then
          ! CARMA wet deposition
          !
          ! NOTE: It needs to follow aero_model_wetdep, so that cam_out%xxxwetxxx
          ! fields have already been set for CAM aerosols and cam_out can be added
          ! to for CARMA aerosols.
          call t_startf ('carma_wetdep_tend')
          call carma_diags_obj%update(cam_in, state, pbuf)
          call carma_wetdep_tend(state, ptend, ztodt, pbuf, dlf, cam_out)
          call carma_diags_obj%output(state, ptend, cam_in, "WETDEPC", ztodt, pbuf)
          call physics_update(state, ptend, ztodt, tend)
          call t_stopf ('carma_wetdep_tend')
       end if

       call t_startf ('convect_deep_tend2')
       call convect_deep_tend_2( state,   ptend,  ztodt,  pbuf )
       call physics_update(state, ptend, ztodt, tend)
       call t_stopf ('convect_deep_tend2')

       ! check tracer integrals
       call check_tracers_chng(state, tracerint, "cmfmca", nstep, ztodt,  zero_tracers)

       call t_stopf('aerosol_wet_processes')

   endif

    !===================================================
    ! Moist physical parameteriztions complete:
    ! send dynamical variables, and derived variables to history file
    !===================================================

    call t_startf('bc_history_write')
    call diag_phys_writeout(state, pbuf)
    call diag_conv(state, ztodt, pbuf)

    call t_stopf('bc_history_write')

    !===================================================
    ! Write cloud diagnostics on history file
    !===================================================

    call t_startf('bc_cld_diag_history_write')

    call cloud_diagnostics_calc(state, pbuf)

    call t_stopf('bc_cld_diag_history_write')

    !===================================================
    ! Radiation computations
    !===================================================
    call t_startf('radiation')

    if (trim(cam_take_snapshot_before) == "radiation_tend") then
       call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, &
                  flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
    end if

    call radiation_tend( &
       state, ptend, pbuf, cam_out, cam_in, net_flx)

    ! Set net flux used by spectral dycores
    do i=1,ncol
       tend%flx_net(i) = net_flx(i)
    end do

    if ( (trim(cam_take_snapshot_after) == "radiation_tend") .and.     &
         (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then
       call cam_snapshot_ptend_outfld(ptend, lchnk)
    end if
    call physics_update(state, ptend, ztodt, tend)

    if (trim(cam_take_snapshot_after) == "radiation_tend") then
       call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, &
                  flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx)
    end if

    call check_energy_cam_chng(state, tend, "radheat", nstep, ztodt, zero, zero, zero, net_flx)

    call t_stopf('radiation')

    ! Diagnose the location of the tropopause and its location to the history file(s).
    call t_startf('tropopause')
    call tropopause_output(state)
    call t_stopf('tropopause')

    ! Save atmospheric fields to force surface models
    call t_startf('cam_export')
    call cam_export (state,cam_out,pbuf)
    call t_stopf('cam_export')

    ! Write export state to history file
    call t_startf('diag_export')
    call diag_export(cam_out)
    call t_stopf('diag_export')

    ! clean CARMA diagnostics object
    if (associated(carma_diags_obj)) then
       deallocate(carma_diags_obj)
       nullify(carma_diags_obj)
    end if

  end subroutine tphysbc

subroutine phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d)
!-----------------------------------------------------------------------------------
!
! Purpose: The place for parameterizations to call per timestep initializations.
!          Generally this is used to update time interpolated fields from boundary
!          datasets.
!
!-----------------------------------------------------------------------------------
  use chemistry,           only: chem_timestep_init
  use chem_surfvals,       only: chem_surfvals_set
  use physics_types,       only: physics_state
  use physics_buffer,      only: physics_buffer_desc
  use carma_intr,          only: carma_timestep_init
  use ghg_data,            only: ghg_data_timestep_init
  use aoa_tracers,         only: aoa_tracers_timestep_init
  use vertical_diffusion,  only: vertical_diffusion_ts_init
  use radheat,             only: radheat_timestep_init
  use solar_data,          only: solar_data_advance
  use qbo,                 only: qbo_timestep_init
  use iondrag,             only: do_waccm_ions, iondrag_timestep_init
  use perf_mod

  use prescribed_ozone,    only: prescribed_ozone_adv
  use prescribed_ghg,      only: prescribed_ghg_adv
  use prescribed_aero,     only: prescribed_aero_adv
  use aerodep_flx,         only: aerodep_flx_adv
  use aircraft_emit,       only: aircraft_emit_adv
  use prescribed_volcaero, only: prescribed_volcaero_adv
  use prescribed_strataero,only: prescribed_strataero_adv
  use mo_apex,             only: mo_apex_init
  use epp_ionization,      only: epp_ionization_active
  use iop_forcing,         only: scam_use_iop_srf
  use nudging,             only: Nudge_Model, nudging_timestep_init
  use waccmx_phys_intr,    only: waccmx_phys_ion_elec_temp_timestep_init
  use phys_grid_ctem,      only: phys_grid_ctem_diags
  use surface_emissions_mod,only: surface_emissions_adv
  use elevated_emissions_mod,only: elevated_emissions_adv

  implicit none

  type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
  type(cam_in_t),      intent(inout), dimension(begchunk:endchunk) :: cam_in
  type(cam_out_t),     intent(inout), dimension(begchunk:endchunk) :: cam_out

  type(physics_buffer_desc), pointer                 :: pbuf2d(:,:)

  !-----------------------------------------------------------------------------

  if (single_column) call scam_use_iop_srf(cam_in)

  ! update geomagnetic coordinates
  if (epp_ionization_active .or. do_waccm_ions) then
     call mo_apex_init(phys_state)
  endif

  ! Chemistry surface values
  call chem_surfvals_set()
  call surface_emissions_adv(pbuf2d, phys_state)
  call elevated_emissions_adv(pbuf2d, phys_state)

  ! Solar irradiance
  call solar_data_advance()

  ! Time interpolate for chemistry.
  call chem_timestep_init(phys_state, pbuf2d)

  if( waccmx_is('ionosphere') ) then
     call waccmx_phys_ion_elec_temp_timestep_init(phys_state,pbuf2d)
  endif

  ! Prescribed tracers
  call prescribed_ozone_adv(phys_state, pbuf2d)
  call prescribed_ghg_adv(phys_state, pbuf2d)
  call prescribed_aero_adv(phys_state, pbuf2d)
  call aircraft_emit_adv(phys_state, pbuf2d)
  call prescribed_volcaero_adv(phys_state, pbuf2d)
  call prescribed_strataero_adv(phys_state, pbuf2d)

  ! prescribed aerosol deposition fluxes
  call aerodep_flx_adv(phys_state, pbuf2d, cam_out)

  ! Time interpolate data models of gasses in pbuf2d
  call ghg_data_timestep_init(pbuf2d,  phys_state)

  ! Upper atmosphere radiative processes
  call radheat_timestep_init(phys_state, pbuf2d)

  ! Time interpolate for vertical diffusion upper boundary condition
  call vertical_diffusion_ts_init(pbuf2d, phys_state)

  !----------------------------------------------------------------------
  ! update QBO data for this time step
  !----------------------------------------------------------------------
  call qbo_timestep_init

  call iondrag_timestep_init()

  call carma_timestep_init()

  ! age of air tracers
  call aoa_tracers_timestep_init(phys_state)

  ! Update Nudging values, if needed
  !----------------------------------
  if(Nudge_Model) call nudging_timestep_init(phys_state)

  ! Update TEM diagnostics
  call phys_grid_ctem_diags(phys_state)

end subroutine phys_timestep_init

end module physpkg