module cam_diagnostics

!---------------------------------------------------------------------------------
! Module to compute a variety of diagnostics quantities for history files
!---------------------------------------------------------------------------------

use shr_kind_mod,  only: r8 => shr_kind_r8
use camsrfexch,    only: cam_in_t, cam_out_t
use physics_types, only: physics_state, physics_tend
use ppgrid,        only: pcols, pver, pverp, begchunk, endchunk
use physics_buffer, only: physics_buffer_desc, pbuf_add_field, dtype_r8, dyn_time_lvls, &
                          pbuf_get_field, pbuf_get_index, pbuf_old_tim_idx



use cam_history,   only: outfld, write_inithist, hist_fld_active
use constituents,  only: pcnst, cnst_name, cnst_longname, cnst_cam_outfld, ptendnam, dmetendnam, apcnst, bpcnst, &
                         cnst_get_ind
use chemistry,     only: chem_is
use dycore,        only: dycore_is
use phys_control,  only: phys_getopts
use wv_saturation, only: qsat, qsat_water, svp_ice
use time_manager,  only: is_first_step
use physconst,     only: cpair, rair, gravit, latvap, epsilo

use iop_data_mod,  only: single_column, wfld
use cam_abortutils,    only: endrun

implicit none
private
save

! Public interfaces

public :: &
   diag_register,      &! register pbuf space
   diag_init,          &! initialization
   diag_allocate,      &! allocate memory for module variables
   diag_deallocate,    &! deallocate memory for module variables
   diag_conv_tend_ini, &! initialize convective tendency calcs
   diag_phys_writeout, &! output diagnostics of the dynamics
   diag_phys_tend_writeout, & ! output physics tendencies
   diag_state_b4_phys_write,& ! output state before physics execution
   diag_conv,          &! output diagnostics of convective processes
   diag_surf,          &! output diagnostics of the surface
   diag_export,        &! output export state
   diag_physvar_ic,    &
   diag_readnl          ! read namelist options

logical, public :: inithist_all = .false. ! Flag to indicate set of fields to be 
                                          ! included on IC file
                                          !  .false.  include only required fields
                                          !  .true.   include required *and* optional fields

! Private data

integer :: dqcond_num                     ! number of constituents to compute convective
character(len=16) :: dcconnam(pcnst)      ! names of convection tendencies
                                          ! tendencies for
real(r8), allocatable :: dtcond(:,:,:)    ! temperature tendency due to convection
type dqcond_t
   real(r8), allocatable :: cnst(:,:,:)   ! constituent tendency due to convection
end type dqcond_t
type(dqcond_t), allocatable :: dqcond(:)

character(len=8) :: diag_cnst_conv_tend = 'q_only' ! output constituent tendencies due to convection
                                                   ! 'none', 'q_only' or 'all'

logical          :: history_amwg                   ! output the variables used by the AMWG diag package
logical          :: history_verbose                ! produce verbose history output
logical          :: history_vdiag                  ! output the variables used by the AMWG variability diag package
logical          :: history_eddy                   ! output the eddy variables
logical          :: history_budget                 ! output tendencies and state variables for CAM4
                                                   ! temperature, water vapor, cloud ice and cloud
                                                   ! liquid budgets.
integer          :: history_budget_histfile_num    ! output history file number for budget fields

!Physics buffer indices
integer  ::      qcwat_idx  = 0 
integer  ::      tcwat_idx  = 0 
integer  ::      lcwat_idx  = 0 
integer  ::      cld_idx    = 0 
integer  ::      concld_idx = 0 
integer  ::      tke_idx    = 0 
integer  ::      kvm_idx    = 0 
integer  ::      kvh_idx    = 0 
integer  ::      cush_idx   = 0 
integer  ::      t_ttend_idx = 0

integer  ::      prec_dp_idx  = 0
integer  ::      snow_dp_idx  = 0
integer  ::      prec_sh_idx  = 0
integer  ::      snow_sh_idx  = 0
integer  ::      prec_sed_idx = 0
integer  ::      snow_sed_idx = 0
integer  ::      prec_pcw_idx = 0
integer  ::      snow_pcw_idx = 0

integer  ::      wsresp_idx = 0
integer  ::      tau_est_idx = 0

integer :: tpert_idx=-1, qpert_idx=-1, pblh_idx=-1
logical :: prog_modal_aero
logical :: linearize_pbl_winds

integer,  parameter :: nbands = 152
real(r8)            :: bbands_1d(nbands+1)
real(r8), target    :: bbands(2,nbands)
real(r8), target    :: mbands(nbands)

integer,  parameter :: nhrs = 24
real(r8)            :: bhrs_1d(nhrs+1)
real(r8), target    :: bhrs(2,nhrs)  
real(r8), target    :: mhrs(nhrs)
real(r8)            :: lonb(nhrs+1)

contains

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

subroutine diag_register
   
   use cam_history_support, only: add_hist_coord
    
   ! Request physics buffer space for fields that persist across timesteps.
   call pbuf_add_field('T_TTEND', 'global', dtype_r8, (/pcols,pver,dyn_time_lvls/), t_ttend_idx)

   call add_hist_coord('pdfbands', nbands, 'Number of PDF bands',  &
        'mm/day', mbands, bounds_name='prec_bnds', bounds=bbands)

   call add_hist_coord('diurnal', nhrs, 'Diurnal cycle',  &
        'GMT hour', mhrs, bounds_name='hr_bnds', bounds=bhrs)

end subroutine diag_register

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

subroutine diag_readnl(nlfile)
  use namelist_utils,  only: find_group_name
  use units,           only: getunit, freeunit
  use mpishorthand
  use spmd_utils,      only: masterproc

  character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input

  ! Local variables
  integer :: unitn, ierr
  character(len=*), parameter :: subname = 'diag_readnl'

  namelist /cam_diag_opts/ diag_cnst_conv_tend
  !-----------------------------------------------------------------------------

  if (masterproc) then
      unitn = getunit()
      open( unitn, file=trim(nlfile), status='old' )
      call find_group_name(unitn, 'cam_diag_opts', status=ierr)
      if (ierr == 0) then
         read(unitn, cam_diag_opts, iostat=ierr)
         if (ierr /= 0) then
            call endrun(subname // ':: ERROR reading namelist')
         end if
      end if
      close(unitn)
      call freeunit(unitn)
   end if

#ifdef SPMD
   ! Broadcast namelist variables
   call mpibcast(diag_cnst_conv_tend, len(diag_cnst_conv_tend), mpichar,  0, mpicom)
#endif
   
end subroutine diag_readnl

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

subroutine diag_init()

  ! Declare the history fields for which this module contains outfld calls.

   use cam_history,        only: addfld, horiz_only, add_default
   use constituent_burden, only: constituent_burden_init
   use cam_control_mod,    only: moist_physics, ideal_phys
   use tidal_diag,         only: tidal_diag_init 

   integer :: k, m
   ! Note - this is a duplication of information in ice_constants 
   ! Cannot put in a use statement if want to swap ice models to cice4
   integer, parameter :: plevmx = 4       ! number of subsurface levels
   character(len=8), parameter :: tsnam(plevmx) = (/ 'TS1', 'TS2', 'TS3', 'TS4' /)
   integer :: ixcldice, ixcldliq ! constituent indices for cloud liquid and ice water.
   integer :: ierr

   call phys_getopts(prog_modal_aero_out = prog_modal_aero, &
                     linearize_pbl_winds_out = linearize_pbl_winds)

   ! outfld calls in diag_phys_writeout

   call addfld ('PDFPRECT',(/ 'pdfbands' /), 'A','fraction','PDF of PRECT')
   call addfld ('PDFPRECL',(/ 'pdfbands' /), 'A','fraction','PDF of PRECL')
   call addfld ('PDFPRECC',(/ 'pdfbands' /), 'A','fraction','PDF of PRECC')

   call addfld ('ADFPRECT',(/ 'pdfbands' /), 'A','mm/day','amount distribution of PRECT')
   call addfld ('ADFPRECL',(/ 'pdfbands' /), 'A','mm/day','amount distribution of PRECL')
   call addfld ('ADFPRECC',(/ 'pdfbands' /), 'A','mm/day','amount distribution of PRECC')

   call addfld ('DPRECT',(/ 'diurnal' /), 'A','mm/day','PRECT diurnal cycle ')
   call addfld ('DPRECL',(/ 'diurnal' /), 'A','mm/day','PRECL diurnal cycle ')
   call addfld ('DPRECC',(/ 'diurnal' /), 'A','mm/day','PRECC diurnal cycle ')     

   call addfld ('NSTEP',horiz_only,    'A','timestep','Model timestep')
   call addfld ('PHIS',horiz_only,    'I','m2/s2','Surface geopotential')

   call addfld ('PS',horiz_only,    'A','Pa','Surface pressure', &
    standard_name='surface_air_pressure')
   call addfld ('T',(/ 'lev' /), 'A','K','Temperature',standard_name='air_temperature')
   call addfld ('U',(/ 'lev' /), 'A','m/s','Zonal wind',standard_name='eastward_wind')
   call addfld ('V',(/ 'lev' /), 'A','m/s','Meridional wind',standard_name='northward_wind')
   call addfld (cnst_name(1),(/ 'lev' /), 'A','kg/kg',cnst_longname(1))

   ! State before physics
   call addfld ('TBP',(/ 'lev' /), 'A','K','Temperature (before physics)'       )
   call addfld (bpcnst(1) ,(/ 'lev' /), 'A','kg/kg',cnst_longname(1)//' (before physics)')
   ! State after physics
   call addfld ('TAP',(/ 'lev' /), 'A','K','Temperature (after physics)'       )
   call addfld ('UAP',(/ 'lev' /), 'A','m/s','Zonal wind (after physics)'        )
   call addfld ('VAP',(/ 'lev' /), 'A','m/s','Meridional wind (after physics)'   )
   call addfld (apcnst(1) ,(/ 'lev' /), 'A','kg/kg',cnst_longname(1)//' (after physics)')
   call addfld ('CAPE', horiz_only, 'A', 'J/kg', 'Convectively available potential energy')
   call addfld ('CIN', horiz_only, 'A', 'J/kg', 'Convective inhibition')
   
   if ( dycore_is('LR') .or. dycore_is('SE') ) then
      call addfld ('TFIX',horiz_only,    'A'     ,'K/s','T fixer (T equivalent of Energy correction)')
      call addfld ('PTTEND_RESID',(/ 'lev' /), 'A'     ,'K/s',&
                   'T-tendency due to BAB kluge at end of tphysac (diagnostic not part of T-budget)' )
   end if
   call addfld ('TTEND_TOT',(/ 'lev' /), 'A','K/s' ,'Total temperature tendency'   )
  
   ! column burdens for all constituents except water vapor
   call constituent_burden_init

   call addfld ('Z3',(/ 'lev' /), 'A','m','Geopotential Height (above sea level)',standard_name='geopotential_height')
   call addfld ('Z1000',horiz_only,    'A','m','Geopotential Z at 1000 mbar pressure surface')
   call addfld ('Z975',horiz_only,    'A','m','Geopotential Z at 975 mbar pressure surface')
   call addfld ('Z950',horiz_only,    'A','m','Geopotential Z at 950 mbar pressure surface')
   call addfld ('Z925',horiz_only,    'A','m','Geopotential Z at 925 mbar pressure surface')
   call addfld ('Z900',horiz_only,    'A','m','Geopotential Z at 900 mbar pressure surface')
   call addfld ('Z850',horiz_only,    'A','m','Geopotential Z at 850 mbar pressure surface')
   call addfld ('Z800',horiz_only,    'A','m','Geopotential Z at 800 mbar pressure surface')
   call addfld ('Z700',horiz_only,    'A','m','Geopotential Z at 700 mbar pressure surface')
   call addfld ('Z600',horiz_only,    'A','m','Geopotential Z at 600 mbar pressure surface')
   call addfld ('Z500',horiz_only,    'A','m','Geopotential Z at 500 mbar pressure surface')
   call addfld ('Z400',horiz_only,    'A','m','Geopotential Z at 400 mbar pressure surface')
   call addfld ('Z300',horiz_only,    'A','m','Geopotential Z at 300 mbar pressure surface')
   call addfld ('Z200',horiz_only,    'A','m','Geopotential Z at 200 mbar pressure surface')
   call addfld ('Z100',horiz_only,    'A','m','Geopotential Z at 100 mbar pressure surface')
   call addfld ('Z050',horiz_only,    'A','m','Geopotential Z at 50 mbar pressure surface')
   call addfld ('Z010',horiz_only,    'A','m','Geopotential Z at 10 mbar pressure surface')

   call addfld ('ZZ',(/ 'lev' /), 'A','m2','Eddy height variance' )
   call addfld ('VZ',(/ 'lev' /), 'A','m2/s','Meridional transport of geopotential energy')
   call addfld ('VT',(/ 'lev' /), 'A','K m/s   ','Meridional heat transport')
   call addfld ('VU',(/ 'lev' /), 'A','m2/s2','Meridional flux of zonal momentum' )
   call addfld ('VV',(/ 'lev' /), 'A','m2/s2','Meridional velocity squared' )

   if(prog_modal_aero) then !Only for prognostic aerosols
      call addfld ('bc_a1_2',(/ 'lev' /), 'A','kg2/kg2','bc_a1 squared')
      call addfld ('dst_a1_2',(/ 'lev' /), 'A','kg2/kg2','dst_a1 squared')
      call addfld ('dst_a3_2',(/ 'lev' /), 'A','kg2/kg2','dst_a3 squared')
      call addfld ('ncl_a1_2',(/ 'lev' /), 'A','kg2/kg2','ncl_a1 squared')
      call addfld ('ncl_a2_2',(/ 'lev' /), 'A','kg2/kg2','ncl_a2 squared')
      call addfld ('ncl_a3_2',(/ 'lev' /), 'A','kg2/kg2','ncl_a3 squared')
      call addfld ('so4_a1_2',(/ 'lev' /), 'A','kg2/kg2','so4_a1 squared')
      call addfld ('so4_a2_2',(/ 'lev' /), 'A','kg2/kg2','so4_a2 squared')
      call addfld ('so4_a3_2',(/ 'lev' /), 'A','kg2/kg2','so4_a3 squared')
      call addfld ('soa_a1_2',(/ 'lev' /), 'A','kg2/kg2','soa_a1 squared')
      call addfld ('soa_a2_2',(/ 'lev' /), 'A','kg2/kg2','soa_a2 squared')
      call addfld ('pom_a1_2',(/ 'lev' /), 'A','kg2/kg2','pom_a1 squared')
      call addfld ('Vbc_a1',(/ 'lev' /), 'A','m/skg/kg','Meridional bc_a1 transport')
      call addfld ('Vdst_a1',(/ 'lev' /), 'A','m/skg/kg','Meridional dst_a1 transport')
      call addfld ('Vdst_a3',(/ 'lev' /), 'A','m/skg/kg','Meridional dst_a3 transport')
      call addfld ('Vncl_a1',(/ 'lev' /), 'A','m/skg/kg','Meridional ncl_a1 transport')
      call addfld ('Vncl_a2',(/ 'lev' /), 'A','m/skg/kg','Meridional ncl_a2 transport')
      call addfld ('Vncl_a3',(/ 'lev' /), 'A','m/skg/kg','Meridional ncl_a3 transport')
      call addfld ('Vso4_a1',(/ 'lev' /), 'A','m/skg/kg','Meridional so4_a1 transport')
      call addfld ('Vso4_a2',(/ 'lev' /), 'A','m/skg/kg','Meridional so4_a2 transport')
      call addfld ('Vso4_a3',(/ 'lev' /), 'A','m/skg/kg','Meridional so4_a3 transport')
      call addfld ('Vsoa_a1',(/ 'lev' /), 'A','m/skg/kg','Meridional soa_a1 transport')
      call addfld ('Vsoa_a2',(/ 'lev' /), 'A','m/skg/kg','Meridional soa_a2 transport')
      call addfld ('Vpom_a1',(/ 'lev' /), 'A','m/skg/kg','Meridional pom_a1 transport')
   endif
   call addfld ('VQ',(/ 'lev' /), 'A','m/s kg/kg','Meridional water transport')
   call addfld ('QQ',(/ 'lev' /), 'A','kg2/kg2','Eddy moisture variance')
   call addfld ('OMEGAV',(/ 'lev' /) ,'A','m Pa/s2 ','Vertical flux of meridional momentum' )
   call addfld ('OMGAOMGA',(/ 'lev' /) ,'A','Pa2/s2','Vertical flux of vertical momentum' )
   call addfld ('OMEGAQ',(/ 'lev' /) ,'A','kgPa/kgs','Vertical water transport' )

   call addfld ('UU',(/ 'lev' /), 'A','m2/s2','Zonal velocity squared' )
   call addfld ('WSPEED',(/ 'lev' /), 'X','m/s','Horizontal total wind speed maximum' )
   call addfld ('WSPDSRFMX',horiz_only,    'X','m/s','Horizontal total wind speed maximum at the surface' )
   call addfld ('WSPDSRFAV',horiz_only,    'A','m/s','Horizontal total wind speed average at the surface' )

   call addfld ('OMEGA',(/ 'lev' /), 'A','Pa/s','Vertical velocity (pressure)', &
      standard_name='lagrangian_tendency_of_air_pressure')
   call addfld ('OMEGAT',(/ 'lev' /), 'A','K Pa/s  ','Vertical heat flux' )
   call addfld ('OMEGAU',(/ 'lev' /), 'A','m Pa/s2 ','Vertical flux of zonal momentum' )
   call addfld ('OMEGA1000',horiz_only,    'A','Pa/s','Vertical velocity at 1000 mbar pressure surface')
   call addfld ('OMEGA975',horiz_only,    'A','Pa/s','Vertical velocity at 975 mbar pressure surface')
   call addfld ('OMEGA950',horiz_only,    'A','Pa/s','Vertical velocity at 950 mbar pressure surface')
   call addfld ('OMEGA925',horiz_only,    'A','Pa/s','Vertical velocity at 925 mbar pressure surface')
   call addfld ('OMEGA900',horiz_only,    'A','Pa/s','Vertical velocity at 900 mbar pressure surface')
   call addfld ('OMEGA850',horiz_only,    'A','Pa/s','Vertical velocity at 850 mbar pressure surface')
   call addfld ('OMEGA800',horiz_only,    'A','Pa/s','Vertical velocity at 800 mbar pressure surface')
   call addfld ('OMEGA700',horiz_only,    'A','Pa/s','Vertical velocity at 700 mbar pressure surface')
   call addfld ('OMEGA600',horiz_only,    'A','Pa/s','Vertical velocity at 600 mbar pressure surface')
   call addfld ('OMEGA500',horiz_only,    'A','Pa/s','Vertical velocity at 500 mbar pressure surface')
   call addfld ('OMEGA400',horiz_only,    'A','Pa/s','Vertical velocity at 400 mbar pressure surface')
   call addfld ('OMEGA300',horiz_only,    'A','Pa/s','Vertical velocity at 300 mbar pressure surface')
   call addfld ('OMEGA200',horiz_only,    'A','Pa/s','Vertical velocity at 200 mbar pressure surface')
   call addfld ('OMEGA100',horiz_only,    'A','Pa/s','Vertical velocity at 100 mbar pressure surface')
   call addfld ('OMEGABOT',horiz_only,    'A','Pa/s','Lowest model level vertical velocity')

   call addfld ('RH1000',horiz_only,    'A','%','Relative humidity at 1000 mbar pressure surface')
   call addfld ('RH975',horiz_only,    'A','%','Relative humidity at 975 mbar pressure surface')
   call addfld ('RH950',horiz_only,    'A','%','Relative humidity at 950 mbar pressure surface')
   call addfld ('RH925',horiz_only,    'A','%','Relative humidity at 925 mbar pressure surface')
   call addfld ('RH900',horiz_only,    'A','%','Relative humidity at 900 mbar pressure surface')
   call addfld ('RH850',horiz_only,    'A','%','Relative humidity at 850 mbar pressure surface')
   call addfld ('RH800',horiz_only,    'A','%','Relative humidity at 800 mbar pressure surface')
   call addfld ('RH700',horiz_only,    'A','%','Relative humidity at 700 mbar pressure surface')
   call addfld ('RH600',horiz_only,    'A','%','Relative humidity at 600 mbar pressure surface')
   call addfld ('RH500',horiz_only,    'A','%','Relative humidity at 500 mbar pressure surface')
   call addfld ('RH400',horiz_only,    'A','%','Relative humidity at 400 mbar pressure surface')
   call addfld ('RH300',horiz_only,    'A','%','Relative humidity at 300 mbar pressure surface')
   call addfld ('RH200',horiz_only,    'A','%','Relative humidity at 200 mbar pressure surface')
   call addfld ('RH100',horiz_only,    'A','%','Relative humidity at 100 mbar pressure surface')
   call addfld ('RHBOT',horiz_only,    'A','%','Lowest model level relative humidity')

   call addfld ('MQ',(/ 'lev' /), 'A','kg/m2','Water vapor mass in layer')
   call addfld ('TMQ',horiz_only,    'A','kg/m2','Total (vertically integrated) precipitable water', &
   standard_name='atmosphere_mass_content_of_water_vapor')
   call addfld ('TTQ',horiz_only,   'A', 'kg/m/s','Total (vertically integrated) vapor transport')
   call addfld ('TUQ',horiz_only,    'A','kg/m/s','Total (vertically integrated) zonal water flux')
   call addfld ('TVQ',horiz_only,    'A','kg/m/s','Total (vertically integrated) meridional water flux')
   call addfld ('TUH',horiz_only,    'A','W/m',   'Total (vertically integrated) zonal MSE flux')
   call addfld ('TVH',horiz_only,    'A','W/m',   'Total (vertically integrated) meridional MSE flux')
   call addfld ('DTENDTH', horiz_only, 'A', 'W/m2',   'Dynamic Tendency of Total (vertically integrated) moist static energy')
   call addfld ('DTENDTQ', horiz_only, 'A', 'kg/m2/s','Dynamic Tendency of Total (vertically integrated) specific humidity')
   call addfld ('RELHUM',(/ 'lev' /), 'A','percent','Relative humidity', standard_name='relative_humidity')
   call addfld ('RHW',(/ 'lev' /), 'A','percent'   ,'Relative humidity with respect to liquid')
   call addfld ('RHI',(/ 'lev' /), 'A','percent'   ,'Relative humidity with respect to ice')
   call addfld ('RHCFMIP',(/ 'lev' /), 'A','percent' ,'Relative humidity with respect to water above 273 K, ice below 273 K')
   call addfld ('PSL',horiz_only,    'A','Pa','Sea level pressure', &
      standard_name='air_pressure_at_mean_sea_level')

   call addfld ('T850',horiz_only,    'A','K','Temperature at 850 mbar pressure surface')
   call addfld ('T500',horiz_only,    'A','K','Temperature at 500 mbar pressure surface')
   call addfld ('T300',horiz_only,    'A','K','Temperature at 300 mbar pressure surface')
   call addfld ('T200',horiz_only,    'A','K','Temperature at 200 mbar pressure surface')
   call addfld ('Q850',horiz_only,    'A','kg/kg','Specific Humidity at 850 mbar pressure surface')
   call addfld ('Q500',horiz_only,    'A','kg/kg','Specific Humidity at 500 mbar pressure surface')
   call addfld ('Q200',horiz_only,    'A','kg/kg','Specific Humidity at 200 mbar pressure surface')
   call addfld ('U1000',horiz_only,    'A','m/s','Zonal wind at 1000 mbar pressure surface')
   call addfld ('U975',horiz_only,    'A','m/s','Zonal wind at 975 mbar pressure surface')
   call addfld ('U950',horiz_only,    'A','m/s','Zonal wind at 950 mbar pressure surface')
   call addfld ('U925',horiz_only,    'A','m/s','Zonal wind at 925 mbar pressure surface')
   call addfld ('U900',horiz_only,    'A','m/s','Zonal wind at 900 mbar pressure surface')
   call addfld ('U850',horiz_only,    'A','m/s','Zonal wind at 850 mbar pressure surface')
   call addfld ('U800',horiz_only,    'A','m/s','Zonal wind at 800 mbar pressure surface')
   call addfld ('U700',horiz_only,    'A','m/s','Zonal wind at 700 mbar pressure surface')
   call addfld ('U600',horiz_only,    'A','m/s','Zonal wind at 600 mbar pressure surface')
   call addfld ('U500',horiz_only,    'A','m/s','Zonal wind at 500 mbar pressure surface')
   call addfld ('U400',horiz_only,    'A','m/s','Zonal wind at 400 mbar pressure surface')
   call addfld ('U300',horiz_only,    'A','m/s','Zonal wind at 300 mbar pressure surface')
   call addfld ('U250',horiz_only,    'A','m/s','Zonal wind at 250 mbar pressure surface')
   call addfld ('U200',horiz_only,    'A','m/s','Zonal wind at 200 mbar pressure surface')
   call addfld ('U100',horiz_only,    'A','m/s','Zonal wind at 100 mbar pressure surface')
   call addfld ('U010',horiz_only,    'A','m/s','Zonal wind at  10 mbar pressure surface')
   call addfld ('V1000',horiz_only,    'A','m/s','Meridional wind at 1000 mbar pressure surface')
   call addfld ('V975',horiz_only,    'A','m/s','Meridional wind at 975 mbar pressure surface')
   call addfld ('V950',horiz_only,    'A','m/s','Meridional wind at 950 mbar pressure surface')
   call addfld ('V925',horiz_only,    'A','m/s','Meridional wind at 925 mbar pressure surface')
   call addfld ('V900',horiz_only,    'A','m/s','Meridional wind at 900 mbar pressure surface')
   call addfld ('V850',horiz_only,    'A','m/s','Meridional wind at 850 mbar pressure surface')
   call addfld ('V800',horiz_only,    'A','m/s','Meridional wind at 800 mbar pressure surface')
   call addfld ('V700',horiz_only,    'A','m/s','Meridional wind at 700 mbar pressure surface')
   call addfld ('V600',horiz_only,    'A','m/s','Meridional wind at 600 mbar pressure surface')
   call addfld ('V500',horiz_only,    'A','m/s','Meridional wind at 500 mbar pressure surface')
   call addfld ('V400',horiz_only,    'A','m/s','Meridional wind at 400 mbar pressure surface')
   call addfld ('V300',horiz_only,    'A','m/s','Meridional wind at 300 mbar pressure surface')
   call addfld ('V250',horiz_only,    'A','m/s','Meridional wind at 250 mbar pressure surface')
   call addfld ('V200',horiz_only,    'A','m/s','Meridional wind at 200 mbar pressure surface')
   call addfld ('V100',horiz_only,    'A','m/s','Meridional wind at 100 mbar pressure surface')
   call addfld ('V050',horiz_only,    'A','m/s','Meridional wind at  50 mbar pressure surface')
   call addfld ('V010',horiz_only,    'A','m/s','Meridional wind at  10 mbar pressure surface')

   call addfld ('TT',(/ 'lev' /), 'A','K2','Eddy temperature variance' )

   call addfld ('UBOT',horiz_only,    'A','m/s','Lowest model level zonal wind')
   call addfld ('VBOT',horiz_only,    'A','m/s','Lowest model level meridional wind')
   call addfld ('QBOT',horiz_only,    'A','kg/kg','Lowest model level water vapor mixing ratio')
   call addfld ('ZBOT',horiz_only,    'A','m','Lowest model level height')

   call addfld ('ATMEINT',horiz_only, 'A','J/m2','Vertically integrated total atmospheric energy ')

   call addfld ('T1000',horiz_only,   'A','K','Temperature at 1000 mbar pressure surface')
   call addfld ('T975',horiz_only,   'A','K','Temperature at 975 mbar pressure surface')   
   call addfld ('T950',horiz_only,   'A','K','Temperature at 950 mbar pressure surface')   
   call addfld ('T925',horiz_only,   'A','K','Temperature at 925 mbar pressure surface')   
   call addfld ('T900',horiz_only,   'A','K','Temperature at 900 mbar pressure surface')   
   call addfld ('T800',horiz_only,   'A','K','Temperature at 800 mbar pressure surface')   
   call addfld ('T700',horiz_only,   'A','K','Temperature at 700 mbar pressure surface')
   call addfld ('T600',horiz_only,   'A','K','Temperature at 600 mbar pressure surface')
   call addfld ('T400',horiz_only,   'A','K','Temperature at 400 mbar pressure surface')
   call addfld ('T100',horiz_only,   'A','K','Temperature at 100 mbar pressure surface')
   call addfld ('T010',horiz_only,   'A','K','Temperature at 10 mbar pressure surface')
   call addfld ('Q1000',horiz_only,   'A','kg/kg','Specific Humidity at 1000 mbar pressure surface')   
   call addfld ('Q975',horiz_only,   'A','kg/kg','Specific Humidity at 975 mbar pressure surface')
   call addfld ('Q950',horiz_only,   'A','kg/kg','Specific Humidity at 950 mbar pressure surface')
   call addfld ('Q925',horiz_only,   'A','kg/kg','Specific Humidity at 925 mbar pressure surface')
   call addfld ('Q900',horiz_only,   'A','kg/kg','Specific Humidity at 900 mbar pressure surface')
   call addfld ('Q800',horiz_only,   'A','kg/kg','Specific Humidity at 800 mbar pressure surface')
   call addfld ('Q700',horiz_only,   'A','kg/kg','Specific Humidity at 700 mbar pressure surface')
   call addfld ('Q600',horiz_only,   'A','kg/kg','Specific Humidity at 600 mbar pressure surface')
   call addfld ('Q400',horiz_only,   'A','kg/kg','Specific Humidity at 400 mbar pressure surface')
   call addfld ('Q300',horiz_only,   'A','kg/kg','Specific Humidity at 300 mbar pressure surface')
   call addfld ('Q100',horiz_only,   'A','kg/kg','Specific Humidity at 100 mbar pressure surface')
   call addfld ('Q050',horiz_only,   'A','kg/kg','Specific Humidity at 050 mbar pressure surface')
   call addfld ('Q010',horiz_only,   'A','kg/kg','Specific Humidity at 010 mbar pressure surface')

   call addfld ('T7001000',horiz_only,   'A','K','Temperature difference 700 mb - 1000 mb')
   call addfld ('TH7001000',horiz_only,   'A','K','Theta difference 700 mb - 1000 mb')
   call addfld ('THE7001000',horiz_only,   'A','K','ThetaE difference 700 mb - 1000 mb')

   call addfld ('T8501000',horiz_only,   'A','K','Temperature difference 850 mb - 1000 mb')
   call addfld ('TH8501000',horiz_only,   'A','K','Theta difference 850 mb - 1000 mb')   
   call addfld ('THE8501000',horiz_only,   'A','K','ThetaE difference 850 mb - 1000 mb')
   call addfld ('T9251000',horiz_only,   'A','K','Temperature difference 925 mb - 1000 mb') 
   call addfld ('TH9251000',horiz_only,   'A','K','Theta difference 925 mb - 1000 mb')   
   call addfld ('THE9251000',horiz_only,   'A','K','ThetaE difference 925 mb - 1000 mb') 
   call addfld ('UOVERN',horiz_only,   'A','m','wind speed/brunt vaisalla frequency 800-100 mb') 


   ! Add in fields for T and U (not already included) to track Sudden Stratospheric Warming events
   ! Levels include: 250, 200, 150, 100, 50, 25, 10, 5, 2, 1, and TOP (numbers in hPa)
   call addfld ('T250',horiz_only,   'A', 'K',   'Temperature at 250 mbar pressure surface')
   call addfld ('T150',horiz_only,   'A', 'K',   'Temperature at 150 mbar pressure surface')
   call addfld ('T050',horiz_only,   'A', 'K',   'Temperature at 50 mbar pressure surface')
   call addfld ('T025',horiz_only,   'A', 'K',   'Temperature at 25 mbar pressure surface')
   call addfld ('T005',horiz_only,   'A', 'K',   'Temperature at 5 mbar pressure surface')
   call addfld ('T002',horiz_only,   'A', 'K',   'Temperature at 2 mbar pressure surface')
   call addfld ('T001',horiz_only,   'A', 'K',   'Temperature at 1 mbar pressure surface')
   call addfld ('TTOP',horiz_only,   'A', 'K',   'Highest model level temperature')
   call addfld ('U150',horiz_only,   'A', 'm/s', 'Zonal wind at 150 mbar pressure surface')
   call addfld ('U050',horiz_only,   'A', 'm/s', 'Zonal wind at 50 mbar pressure surface')
   call addfld ('U025',horiz_only,   'A', 'm/s', 'Zonal wind at 25 mbar pressure surface')
   call addfld ('U005',horiz_only,   'A', 'm/s', 'Zonal wind at 5 mbar pressure surface')
   call addfld ('U002',horiz_only,   'A', 'm/s', 'Zonal wind at 2 mbar pressure surface')
   call addfld ('U001',horiz_only,   'A', 'm/s', 'Zonal wind at 1 mbar pressure surface')
   call addfld ('UTOP',horiz_only,   'A', 'm/s', 'Highest model level zonal wind')

   call addfld ('U90M',horiz_only,    'A','m/s','Zonal wind at turbine hub height (90m above surface)')
   call addfld ('V90M',horiz_only,    'A','m/s','Meridional wind at turbine hub height (90m above surface)')

   ! YQIN  08/18/23
   call addfld ('LTS', horiz_only,   'A', 'K', 'Lower Tropspheric Stability')
   call addfld ('EIS', horiz_only,   'A', 'K', 'Estimated Inversion Strength')
   call addfld ('EIS_fxRH', horiz_only,   'A', 'K', 'Estimated Inversion Strength with surface RH = 80%')

   ! This field is added by radiation when full physics is used
   if ( ideal_phys )then
      call addfld('QRS', (/ 'lev' /), 'A', 'K/s', 'Solar heating rate')
   end if
 
   ! ----------------------------
   ! determine default variables
   ! ----------------------------
   call phys_getopts(history_amwg_out   = history_amwg    , &
                     history_verbose_out  = history_verbose  , &
                     history_vdiag_out  = history_vdiag   , &
                     history_eddy_out   = history_eddy    , &
                     history_budget_out = history_budget  , &
                     history_budget_histfile_num_out = history_budget_histfile_num)

   if (history_amwg) then
      call add_default ('PHIS    '  , 1, ' ')
      call add_default ('PS      '  , 1, ' ')
      call add_default ('T       '  , 1, ' ')
      call add_default ('U       '  , 1, ' ')
      call add_default ('V       '  , 1, ' ')
      call add_default (cnst_name(1), 1, ' ')
      call add_default ('Z3      '  , 1, ' ')
      call add_default ('OMEGA   '  , 1, ' ')
      call add_default ('VT      ', 1, ' ')
      call add_default ('VU      ', 1, ' ')
      call add_default ('VV      ', 1, ' ')
      call add_default ('VQ      ', 1, ' ')

      ! YQIN 08/18/23
      call add_default ('LTS     ', 1, ' ')
      call add_default ('EIS     ', 1, ' ')
      call add_default ('EIS_fxRH', 1, ' ')

      if(prog_modal_aero .and. history_verbose) then !Only for prognostic aerosols
         call add_default ('Vbc_a1  ', 1, ' ')
         call add_default ('Vdst_a1 ', 1, ' ')
         call add_default ('Vdst_a3 ', 1, ' ')
         call add_default ('Vncl_a1 ', 1, ' ')
         call add_default ('Vncl_a2 ', 1, ' ')
         call add_default ('Vncl_a3 ', 1, ' ')
         call add_default ('Vso4_a1 ', 1, ' ')
         call add_default ('Vso4_a2 ', 1, ' ')
         call add_default ('Vso4_a3 ', 1, ' ')
         call add_default ('Vsoa_a1 ', 1, ' ')
         call add_default ('Vsoa_a2 ', 1, ' ')
         call add_default ('Vpom_a1 ', 1, ' ')
      endif
      call add_default ('UU      ', 1, ' ')
      call add_default ('OMEGAT  ', 1, ' ')
      if(prog_modal_aero .and. history_verbose) then !Only for prognostic aerosols
         call add_default ('bc_a1_2 ', 1, ' ')
         call add_default ('dst_a1_2', 1, ' ')
         call add_default ('dst_a3_2', 1, ' ')
         call add_default ('ncl_a1_2', 1, ' ')
         call add_default ('ncl_a2_2', 1, ' ')
         call add_default ('ncl_a3_2', 1, ' ')
         call add_default ('so4_a1_2', 1, ' ')
         call add_default ('so4_a2_2', 1, ' ')
         call add_default ('so4_a3_2', 1, ' ')
         call add_default ('soa_a1_2', 1, ' ')
         call add_default ('soa_a2_2', 1, ' ')
         call add_default ('pom_a1_2', 1, ' ')
      endif
      call add_default ('TMQ     ', 1, ' ')
      call add_default ('PSL     ', 1, ' ')
      if (moist_physics) then
         call add_default ('RELHUM  ', 1, ' ')
      end if

      ! For Tier 1b global water cycle diagostics
      call add_default ('TUQ      ', 1, ' ')
      call add_default ('TVQ      ', 1, ' ')
      call add_default ('TUH      ', 1, ' ')
      call add_default ('TVH      ', 1, ' ')
      call add_default ('DTENDTQ', 1, ' ')
      call add_default ('DTENDTH', 1, ' ')
   end if
   
   if (history_vdiag) then
     call add_default ('U200', 2, ' ')
     call add_default ('V200', 2, ' ')
     call add_default ('U850', 2, ' ')
     call add_default ('U200', 3, ' ')
     call add_default ('U850', 3, ' ')
     call add_default ('OMEGA500', 3, ' ')
   end if

   if (history_eddy) then
      call add_default ('VT      ', 1, ' ')
      call add_default ('VU      ', 1, ' ')
      call add_default ('VV      ', 1, ' ')
      call add_default ('VQ      ', 1, ' ')
      call add_default ('UU      ', 1, ' ')
      call add_default ('OMEGAT  ', 1, ' ')
      call add_default ('OMEGAQ  ', 1, ' ')
      call add_default ('OMEGAU  ', 1, ' ')
      call add_default ('OMEGAV  ', 1, ' ')
   endif

   if ( history_budget ) then
      call add_default ('PHIS    '  , history_budget_histfile_num, ' ')
      call add_default ('PS      '  , history_budget_histfile_num, ' ')
      call add_default ('T       '  , history_budget_histfile_num, ' ')
      call add_default ('U       '  , history_budget_histfile_num, ' ')
      call add_default ('V       '  , history_budget_histfile_num, ' ')
      call add_default (cnst_name(1), history_budget_histfile_num, ' ')
      call add_default ('TTEND_TOT' , history_budget_histfile_num, ' ')

      ! State before physics (FV)
      call add_default ('TBP     '  , history_budget_histfile_num, ' ')
      call add_default (bpcnst(1)   , history_budget_histfile_num, ' ')
      ! State after physics (FV)
      call add_default ('TAP     '  , history_budget_histfile_num, ' ')
      call add_default ('UAP     '  , history_budget_histfile_num, ' ')
      call add_default ('VAP     '  , history_budget_histfile_num, ' ')
      call add_default (apcnst(1)   , history_budget_histfile_num, ' ')
      if ( dycore_is('LR') .or. dycore_is('SE') ) then
         call add_default ('TFIX    '    , history_budget_histfile_num, ' ')
         call add_default ('PTTEND_RESID', history_budget_histfile_num, ' ')
      end if
   end if

   ! This field is added by radiation when full physics is used
   if ( ideal_phys )then
      call add_default('QRS     ', 1, ' ')
   end if

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Exit here for adiabatic/ideal physics cases !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   if (.not. moist_physics) return


   call addfld ('PDELDRY',(/ 'lev' /), 'A','Pa','Dry pressure difference between levels')
   call addfld ('PSDRY',horiz_only,    'A','Pa','Surface pressure')

   if (chem_is('waccm_ghg') .or. chem_is('waccm_mozart') .or. chem_is('waccm_mozart_mam3')) then
      call add_default ('PS      ', 2, ' ')
      call add_default ('T       ', 2, ' ')
   end if

   ! outfld calls in diag_conv

   call cnst_get_ind('CLDLIQ', ixcldliq)
   call cnst_get_ind('CLDICE', ixcldice)
   call addfld ('DTCOND',(/ 'lev' /), 'A','K/s','T tendency - moist processes')
   call addfld ('DTCOND_24_COS',(/ 'lev' /), 'A','K/s','T tendency - moist processes 24hr. cos coeff.')
   call addfld ('DTCOND_24_SIN',(/ 'lev' /), 'A','K/s','T tendency - moist processes 24hr. sin coeff.')
   call addfld ('DTCOND_12_COS',(/ 'lev' /), 'A','K/s','T tendency - moist processes 12hr. cos coeff.')
   call addfld ('DTCOND_12_SIN',(/ 'lev' /), 'A','K/s','T tendency - moist processes 12hr. sin coeff.')

   ! determine number of constituents for which convective tendencies must be computed
   if (history_budget) then
      dqcond_num = pcnst
   else
      if (diag_cnst_conv_tend == 'none')   dqcond_num = 0
      if (diag_cnst_conv_tend == 'q_only') dqcond_num = 1
      if (diag_cnst_conv_tend == 'all')    dqcond_num = pcnst
   end if

   do m = 1, dqcond_num
      dcconnam(m) = 'DC'//cnst_name(m)
   end do

   if (diag_cnst_conv_tend == 'q_only' .or. diag_cnst_conv_tend == 'all' .or. history_budget) then
      call addfld (dcconnam(1),(/ 'lev' /),'A', 'kg/kg/s',trim(cnst_name(1))//' tendency due to moist processes')
      if ( diag_cnst_conv_tend == 'q_only' .or. diag_cnst_conv_tend == 'all' ) then
         call add_default (dcconnam(1),                           1, ' ')
      end if
      if( history_budget ) then
         call add_default (dcconnam(1), history_budget_histfile_num, ' ')
      end if
      if (diag_cnst_conv_tend == 'all' .or. history_budget) then
         do m = 2, pcnst
            call addfld (dcconnam(m),(/ 'lev' /),'A', 'kg/kg/s',trim(cnst_name(m))//' tendency due to moist processes')
            if( diag_cnst_conv_tend == 'all' ) then
               call add_default (dcconnam(m),                           1, ' ')
            end if
            if( history_budget .and. (m == ixcldliq .or. m == ixcldice) ) then
               call add_default (dcconnam(m), history_budget_histfile_num, ' ')
            end if
         end do
      end if
   end if

   call addfld ('PRECL',horiz_only,    'A','m/s','Large-scale (stable) precipitation rate (liq + ice)'                )
   call addfld ('PRECC',horiz_only,    'A','m/s','Convective precipitation rate (liq + ice)'                          )
   call addfld ('PRECT',horiz_only,    'A','m/s','Total (convective and large-scale) precipitation rate (liq + ice)'  )
   call addfld ('PREC_PCW',horiz_only,    'A','m/s','LS_pcw precipitation rate')
   call addfld ('PREC_zmc',horiz_only,    'A','m/s','CV_zmc precipitation rate')
   call addfld ('PRECTMX',horiz_only,    'X','m/s','Maximum (convective and large-scale) precipitation rate (liq+ice)'  )
   call addfld ('PRECSL',horiz_only,    'A','m/s','Large-scale (stable) snow rate (water equivalent)'                  )
   call addfld ('PRECSC',horiz_only,    'A','m/s','Convective snow rate (water equivalent)'                            )
   call addfld ('PRECCav',horiz_only,    'A','m/s','Average large-scale precipitation (liq + ice)'                      )
   call addfld ('PRECLav',horiz_only,    'A','m/s','Average convective precipitation  (liq + ice)'                      )

   if (linearize_pbl_winds) then
      call addfld ('wsresp',horiz_only,    'A','m/s/Pa','first order response of winds to stress')
      call addfld ('tau_est',horiz_only,    'A','Pa','estimated equilibrium wind stress')
   end if

   ! outfld calls in diag_surf

   call addfld ('SHFLX',horiz_only,    'A','W/m2','Surface sensible heat flux', &
   standard_name='surface_upward_sensible_heat_flux')
   call addfld ('LHFLX',horiz_only,    'A','W/m2','Surface latent heat flux', &
   standard_name = 'surface_upward_latent_heat_flux')
   call addfld ('QFLX',horiz_only,    'A','kg/m2/s','Surface water flux', &
   standard_name = 'water_evapotranspiration_flux')

   call addfld ('TAUX',horiz_only,    'A','N/m2','Zonal surface stress')
   call addfld ('TAUY',horiz_only,    'A','N/m2','Meridional surface stress')
   call addfld ('TREFHT',horiz_only,    'A','K','Reference height temperature', &
      standard_name='air_temperature')
   call addfld ('TREFHTMN',horiz_only,    'M','K','Minimum reference height temperature over output period')
   call addfld ('TREFHTMX',horiz_only,    'X','K','Maximum reference height temperature over output period')
   call addfld ('QREFHT',horiz_only,    'A','kg/kg','Reference height humidity', &
    standard_name = 'specific_humidity')
   call addfld ('U10',horiz_only,    'A','m/s','10m wind speed', &
     standard_name='wind_speed')
   call addfld ('RHREFHT',horiz_only,    'A','1','Reference height relative humidity')

   call addfld ('LANDFRAC',horiz_only,    'A','1','Fraction of sfc area covered by land')
   call addfld ('ICEFRAC',horiz_only,    'A','1','Fraction of sfc area covered by sea-ice')
   call addfld ('OCNFRAC',horiz_only,    'A','1','Fraction of sfc area covered by ocean')

   call addfld ('TREFMNAV',horiz_only,    'A','K','Average of TREFHT daily minimum')
   call addfld ('TREFMXAV',horiz_only,    'A','K','Average of TREFHT daily maximum')

   call addfld ('TS',horiz_only,    'A','K','Surface temperature (radiative)', &
       standard_name = 'surface_temperature')
   call addfld ('TSMN',horiz_only,    'M','K','Minimum surface temperature over output period')
   call addfld ('TSMX',horiz_only,    'X','K','Maximum surface temperature over output period')
   call addfld ('SNOWHLND',horiz_only,    'A','m','Water equivalent snow depth')
   call addfld ('SNOWHICE',horiz_only,    'A','m','Snow depth over ice', fill_value = 1.e30_r8)
   call addfld ('TBOT',horiz_only,    'A','K','Lowest model level temperature')

   call addfld ('ASDIR',       horiz_only,    'A',   '1','albedo: shortwave, direct')
   call addfld ('ASDIF',       horiz_only,    'A',   '1','albedo: shortwave, diffuse')
   call addfld ('ALDIR',       horiz_only,    'A',   '1','albedo: longwave, direct')
   call addfld ('ALDIF',       horiz_only,    'A',   '1','albedo: longwave, diffuse')
   call addfld ('SST',       horiz_only,    'A',     'K','sea surface temperature')

   ! defaults
   if (history_amwg) then
       call add_default ('DTCOND  ', 1, ' ')
       call add_default ('PRECL   ', 1, ' ')
       call add_default ('PRECC   ', 1, ' ')
       call add_default ('PRECSL  ', 1, ' ')
       call add_default ('PRECSC  ', 1, ' ')
       call add_default ('SHFLX   ', 1, ' ')
       call add_default ('LHFLX   ', 1, ' ')
       call add_default ('QFLX    ', 1, ' ')
       call add_default ('TAUX    ', 1, ' ')
       call add_default ('TAUY    ', 1, ' ')
       call add_default ('TREFHT  ', 1, ' ')
       call add_default ('LANDFRAC', 1, ' ')
       call add_default ('OCNFRAC ', 1, ' ')
       call add_default ('QREFHT  ', 1, ' ')
       call add_default ('U10     ', 1, ' ')
       call add_default ('ICEFRAC ', 1, ' ')
       call add_default ('TS      ', 1, ' ')
       call add_default ('TSMN    ', 1, ' ')
       call add_default ('TSMX    ', 1, ' ')
       call add_default ('SNOWHLND', 1, ' ')
       call add_default ('SNOWHICE', 1, ' ')
    endif

    call add_default ('OMEGA500', 1, ' ')
    call add_default ('TH7001000', 1, ' ')

    if (history_vdiag) then
        call add_default ('PRECT   ', 2, ' ')
        call add_default ('PRECT   ', 3, ' ')
        call add_default ('PRECT   ', 4, ' ')
    end if

   ! outfld calls in diag_phys_tend_writeout

   call addfld ('PTTEND'   ,(/ 'lev' /), 'A','K/s','T total physics tendency'                             )
   call addfld (ptendnam(       1),(/ 'lev' /), 'A',  'kg/kg/s',trim(cnst_name(       1))//' total physics tendency '      )
   call addfld (ptendnam(ixcldliq),(/ 'lev' /), 'A',  'kg/kg/s',trim(cnst_name(ixcldliq))//' total physics tendency '      )
   call addfld (ptendnam(ixcldice),(/ 'lev' /), 'A',  'kg/kg/s',trim(cnst_name(ixcldice))//' total physics tendency '      )
   if ( dycore_is('LR') )then
      call addfld (dmetendnam(       1),(/ 'lev' /), 'A','kg/kg/s', &
           trim(cnst_name(       1))//' dme adjustment tendency (FV) ')
      call addfld (dmetendnam(ixcldliq),(/ 'lev' /), 'A','kg/kg/s', &
           trim(cnst_name(ixcldliq))//' dme adjustment tendency (FV) ')
      call addfld (dmetendnam(ixcldice),(/ 'lev' /), 'A','kg/kg/s', &
           trim(cnst_name(ixcldice))//' dme adjustment tendency (FV) ')
   end if

   if ( history_budget ) then
      call add_default ('PTTEND'          , history_budget_histfile_num, ' ')
      call add_default (ptendnam(       1), history_budget_histfile_num, ' ')
      call add_default (ptendnam(ixcldliq), history_budget_histfile_num, ' ')
      call add_default (ptendnam(ixcldice), history_budget_histfile_num, ' ')
      if ( dycore_is('LR') )then
         call add_default(dmetendnam(1)       , history_budget_histfile_num, ' ')
         call add_default(dmetendnam(ixcldliq), history_budget_histfile_num, ' ')
         call add_default(dmetendnam(ixcldice), history_budget_histfile_num, ' ')
      end if
      if( history_budget_histfile_num > 1 ) then
         call add_default ('DTCOND  '         , history_budget_histfile_num, ' ')
      end if
   end if

   ! outfld calls in diag_physvar_ic

   call addfld ('QCWAT&IC',(/ 'lev' /), 'I','kg/kg','q associated with cloud water'                   )
   call addfld ('TCWAT&IC',(/ 'lev' /), 'I','kg/kg','T associated with cloud water'                   )
   call addfld ('LCWAT&IC',(/ 'lev' /), 'I','kg/kg','Cloud water (ice + liq'                          )
   call addfld ('CLOUD&IC',(/ 'lev' /), 'I','fraction','Cloud fraction'                                  )
   call addfld ('CONCLD&IC',(/ 'lev' /), 'I','fraction','Convective cloud fraction'                      )
   call addfld ('TKE&IC',(/ 'ilev' /),'I','m2/s2','Turbulent Kinetic Energy'                        )
   call addfld ('CUSH&IC',horiz_only,    'I','m','Convective Scale Height'                         )
   call addfld ('KVH&IC',(/ 'ilev' /),'I','m2/s','Vertical diffusion diffusivities (heat/moisture)')
   call addfld ('KVM&IC',(/ 'ilev' /),'I','m2/s','Vertical diffusion diffusivities (momentum)'     )
   call addfld ('PBLH&IC',horiz_only,    'I','m','PBL height'                                      )
   call addfld ('TPERT&IC',horiz_only,    'I','K','Perturbation temperature (eddies in PBL)'        )
   call addfld ('QPERT&IC',horiz_only,    'I','kg/kg','Perturbation specific humidity (eddies in PBL)'  )
   call addfld ('TBOT&IC',horiz_only,    'I','K','Lowest model level temperature'                  )


   ! Initial file - Optional fields

   if (inithist_all) then
      call add_default ('CONCLD&IC  ',0, 'I')
      call add_default ('QCWAT&IC   ',0, 'I')
      call add_default ('TCWAT&IC   ',0, 'I')
      call add_default ('LCWAT&IC   ',0, 'I')
      call add_default ('PBLH&IC    ',0, 'I')
      call add_default ('TPERT&IC   ',0, 'I')
      call add_default ('QPERT&IC   ',0, 'I')
      call add_default ('CLOUD&IC   ',0, 'I')
      call add_default ('TKE&IC     ',0, 'I')
      call add_default ('CUSH&IC    ',0, 'I')
      call add_default ('KVH&IC     ',0, 'I')
      call add_default ('KVM&IC     ',0, 'I')
      call add_default ('TBOT&IC    ',0, 'I')
   end if

   ! CAM export state 
   call addfld('a2x_BCPHIWET', horiz_only, 'A', 'kg/m2/s', 'wetdep of hydrophilic black carbon')
   call addfld('a2x_BCPHIDRY', horiz_only, 'A', 'kg/m2/s', 'drydep of hydrophilic black carbon')
   call addfld('a2x_BCPHODRY', horiz_only, 'A', 'kg/m2/s', 'drydep of hydrophobic black carbon')
   call addfld('a2x_OCPHIWET', horiz_only, 'A', 'kg/m2/s', 'wetdep of hydrophilic organic carbon')
   call addfld('a2x_OCPHIDRY', horiz_only, 'A', 'kg/m2/s', 'drydep of hydrophilic organic carbon')
   call addfld('a2x_OCPHODRY', horiz_only, 'A', 'kg/m2/s', 'drydep of hydrophobic organic carbon')
   call addfld('a2x_DSTWET1', horiz_only, 'A',  'kg/m2/s', 'wetdep of dust (bin1)')
   call addfld('a2x_DSTDRY1', horiz_only, 'A',  'kg/m2/s', 'drydep of dust (bin1)')
   call addfld('a2x_DSTWET2', horiz_only, 'A',  'kg/m2/s', 'wetdep of dust (bin2)')
   call addfld('a2x_DSTDRY2', horiz_only, 'A',  'kg/m2/s', 'drydep of dust (bin2)')
   call addfld('a2x_DSTWET3', horiz_only, 'A',  'kg/m2/s', 'wetdep of dust (bin3)')
   call addfld('a2x_DSTDRY3', horiz_only, 'A',  'kg/m2/s', 'drydep of dust (bin3)')
   call addfld('a2x_DSTWET4', horiz_only, 'A',  'kg/m2/s', 'wetdep of dust (bin4)')
   call addfld('a2x_DSTDRY4', horiz_only, 'A',  'kg/m2/s', 'drydep of dust (bin4)')

   !---------------------------------------------------------
   ! CAM history fields for CAM-DOM/CAM-CSIM 
   !---------------------------------------------------------

   ! CAM-DOM history fields
#ifdef COUP_DOM
   call addfld ('TSOCN&IC',horiz_only,    'I','m','Ocean tempertare')
   call add_default ('TSOCN&IC   ',0, 'I')
#endif

  ! CAM-CSIM history fields

  do k=1,plevmx
     call addfld (tsnam(k),horiz_only,'A','K',tsnam(k)//' subsoil temperature')
  end do
  call addfld ('SICTHK'   ,horiz_only,'A','m','Sea ice thickness')
  call addfld ('TSICE'   ,horiz_only,'A','K','Ice temperature')
  do k = 1,plevmx
     call addfld (trim(tsnam(k))//'&IC',horiz_only,'I','K',tsnam(k)//' subsoil temperature')
  end do
  call addfld ('SICTHK&IC',horiz_only,'I','m','Sea ice thickness'                      )
  call addfld ('TSICE&IC',horiz_only,'I','K','Ice temperature'                        )
  call addfld ('SNOWHICE&IC',horiz_only,'I','m','Water equivalent snow depth'            )
  call addfld ('ICEFRAC&IC',horiz_only,'I','fraction','Fraction of sfc area covered by sea-ice')
  call addfld ('TSICERAD&IC',horiz_only,'I','K','Radiatively equivalent ice temperature' )
  do k = 1,plevmx
     call add_default(trim(tsnam(k))//'&IC',0, 'I')
  end do
  call add_default ('SICTHK&IC  ',0, 'I')
  call add_default ('TSICE&IC   ',0, 'I')
  call add_default ('SNOWHICE&IC',0, 'I')
  call add_default ('ICEFRAC&IC ',0, 'I')
  if (inithist_all) then
     call add_default ('TSICERAD&IC',0, 'I')
  end if

  !---------------------------------------------------------
  ! WACCM diagnostic history fields 
  !---------------------------------------------------------

  if (chem_is('waccm_ghg') .or. chem_is('waccm_mozart') .or. chem_is('waccm_mozart_mam3')) then

    ! create history variables for fourier coefficients of the diurnal 
    ! and semidiurnal tide in T, U, V, and Z3

    call tidal_diag_init() 

  endif

  qcwat_idx  = pbuf_get_index('QCWAT',ierr)
  tcwat_idx  = pbuf_get_index('TCWAT',ierr)
  lcwat_idx  = pbuf_get_index('LCWAT',ierr)
  cld_idx    = pbuf_get_index('CLD')
  concld_idx = pbuf_get_index('CONCLD')
  
  tke_idx  = pbuf_get_index('tke')
  kvm_idx  = pbuf_get_index('kvm')
  kvh_idx  = pbuf_get_index('kvh')
  cush_idx = pbuf_get_index('cush')
  
  pblh_idx  = pbuf_get_index('pblh')
  tpert_idx = pbuf_get_index('tpert')
  qpert_idx = pbuf_get_index('qpert',ierr)

  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')
  prec_sed_idx = pbuf_get_index('PREC_SED')
  snow_sed_idx = pbuf_get_index('SNOW_SED')
  prec_pcw_idx = pbuf_get_index('PREC_PCW')
  snow_pcw_idx = pbuf_get_index('SNOW_PCW')
  pblh_idx     = pbuf_get_index('pblh')

  if (linearize_pbl_winds) then
     wsresp_idx  = pbuf_get_index('wsresp')
     tau_est_idx  = pbuf_get_index('tau_est')
  end if

  bbands_1d(1) = 0._r8
  mbands(1) = 0._r8
  do k=2, nbands
      mbands(k) = 0.03 * 1.07**(k-2)    ! smallest bin center is 0.03, 7% growth of bin sizes
  end do 
  do k=2,nbands+1
        bbands_1d(k) = sqrt(0.03**2 / 1.07) * 1.07**(k-2)
  end do
  do k=1,nbands
        bbands(1,k)=bbands_1d(k)
        bbands(2,k)=bbands_1d(k+1)
  enddo    
  bhrs_1d(1)=0._r8
  lonb(1) = -7.5_r8
  do k=2,nhrs+1
     bhrs_1d(k)=bhrs_1d(k-1)+3600._r8
     lonb(k)   =lonb(k-1)+15._r8
  end do
  lonb(1) = 352.5_r8    
  do k=1,nhrs
     bhrs(1,k)=bhrs_1d(k)
     bhrs(2,k)=bhrs_1d(k+1)
     mhrs(k)  = 0.5_r8 * (bhrs_1d(k)+bhrs_1d(k+1))
  enddo

end subroutine diag_init

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

subroutine diag_allocate()
   use infnan, only: nan, assignment(=)

   ! Allocate memory for module variables.
   ! Done at the begining of a physics step at same point as the pbuf allocate for
   ! variables with "physpkg" scope.

   ! Local variables
   character(len=*), parameter :: sub = 'diag_allocate'
   integer :: i, istat

   allocate(dtcond(pcols,pver,begchunk:endchunk), stat=istat)
   if ( istat /= 0 ) call endrun (sub//': ERROR: allocate failed')
   dtcond = nan

   if (dqcond_num > 0) then
      allocate(dqcond(dqcond_num))
      do i = 1, dqcond_num
         allocate(dqcond(i)%cnst(pcols,pver,begchunk:endchunk), stat=istat)
         if ( istat /= 0 ) call endrun (sub//': ERROR: allocate failed')
         dqcond(i)%cnst = nan
      end do
   end if

end subroutine diag_allocate

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

subroutine diag_deallocate()

! Deallocate memory for module variables.
! Done at the end of a physics step at same point as the pbuf deallocate for
! variables with "physpkg" scope.

! Local variables
   character(len=*), parameter :: sub = 'diag_deallocate'
   integer :: i, istat

   deallocate(dtcond, stat=istat)
   if ( istat /= 0 ) call endrun (sub//': ERROR: deallocate failed')

   if (dqcond_num > 0) then
      do i = 1, dqcond_num
         deallocate(dqcond(i)%cnst, stat=istat)
         if ( istat /= 0 ) call endrun (sub//': ERROR: deallocate failed')
      end do
      deallocate(dqcond, stat=istat)
      if ( istat /= 0 ) call endrun (sub//': ERROR: deallocate failed')
   end if

end subroutine diag_deallocate
!===============================================================================

subroutine diag_conv_tend_ini(state,pbuf)

! Initialize convective tendency calcs.

   
! Argument:

   type(physics_state), intent(in) :: state
   
   type(physics_buffer_desc), pointer :: pbuf(:)

! Local variables:

   integer :: i, k, m, lchnk, ncol
   real(r8), pointer, dimension(:,:) :: t_ttend

   lchnk = state%lchnk
   ncol  = state%ncol

   do k = 1, pver
      do i = 1, ncol
         dtcond(i,k,lchnk) = state%s(i,k)
      end do
   end do

   do m = 1, dqcond_num
      do k = 1, pver
         do i = 1, ncol
            dqcond(m)%cnst(i,k,lchnk) = state%q(i,k,m)
         end do
      end do
   end do

   !! initialize to pbuf T_TTEND to temperature at first timestep
   if (is_first_step()) then
      do m = 1, dyn_time_lvls
         call pbuf_get_field(pbuf, t_ttend_idx, t_ttend, start=(/1,1,m/), kount=(/pcols,pver,1/))
         t_ttend(:ncol,:) = state%t(:ncol,:)
      end do
   end if

end subroutine diag_conv_tend_ini
!===============================================================================

  subroutine diag_phys_writeout(state, psl)

!----------------------------------------------------------------------- 
! 
! Purpose: record dynamics variables on physics grid
!
!-----------------------------------------------------------------------
    use physconst,          only: gravit, rga, rair, cpair, latvap, rearth, pi, cappa
    use time_manager,       only: get_nstep
    use interpolate_data,   only: vertinterp, vertinterpz
    use constituent_burden, only: constituent_burden_comp
    use cam_control_mod,    only: moist_physics
    use co2_cycle,          only: c_i, co2_transport

    use tidal_diag,         only: tidal_diag_write
!-----------------------------------------------------------------------
!
! Arguments
!
   type(physics_state), intent(inout) :: state
   real(r8), optional , intent(out)   :: psl(pcols) 
!
!---------------------------Local workspace-----------------------------
!
    real(r8) ftem(pcols,pver) ! temporary workspace
    real(r8) ftem1(pcols,pver) ! another temporary workspace
    real(r8) ftem2(pcols,pver) ! another temporary workspace
    real(r8) ftem4(pcols,pver) ! another temporary workspace
    real(r8) ftem5(pcols,pver) ! another temporary workspace
    real(r8) psl_tmp(pcols)   ! Sea Level Pressure
    real(r8) z3(pcols,pver)   ! geo-potential height
    real(r8) p_surf(pcols)    ! data interpolated to a pressure surface
    real(r8) p_surf_t1(pcols)    ! data interpolated to a pressure surface
    real(r8) p_surf_t2(pcols)    ! data interpolated to a pressure surface
    real(r8) p_surf_q1(pcols)    ! data interpolated to a pressure surface    
    real(r8) p_surf_q2(pcols)    ! data interpolated to a pressure surface        
    real(r8) tem2(pcols,pver) ! temporary workspace
    real(r8) timestep(pcols)  ! used for outfld call
    real(r8) esl(pcols,pver)   ! saturation vapor pressures 
    real(r8) esi(pcols,pver)   ! 
    real(r8) dlon(pcols)      ! width of grid cell (meters)
    integer  plon             ! number of longitudes

    integer i, k, m, lchnk, ncol, nstep
!
!-----------------------------------------------------------------------
!
    lchnk = state%lchnk
    ncol  = state%ncol

    ! Output NSTEP for debugging
    nstep = get_nstep()
    timestep(:ncol) = nstep
    call outfld ('NSTEP   ',timestep, pcols, lchnk)

    call outfld('T       ',state%t , pcols   ,lchnk   )
    call outfld('PS      ',state%ps, pcols   ,lchnk   )
    call outfld('U       ',state%u , pcols   ,lchnk   )
    call outfld('V       ',state%v , pcols   ,lchnk   )
    do m=1,pcnst
       if ( cnst_cam_outfld(m) ) then
          call outfld(cnst_name(m),state%q(1,1,m),pcols ,lchnk )
       end if
    end do

    if (co2_transport()) then
       do m = 1,4
          call outfld(trim(cnst_name(c_i(m)))//'_BOT', state%q(1,pver,c_i(m)), pcols, lchnk)
       end do
    end if

    ! column burdens of all constituents except water vapor
    call constituent_burden_comp(state)

    if ( moist_physics) then
       call outfld('PDELDRY ',state%pdeldry, pcols, lchnk)
       call outfld('PSDRY',   state%psdry,   pcols, lchnk) 
    end if

    call outfld('PHIS    ',state%phis,    pcols,   lchnk     )



#if (defined E3SM_SCM_REPLAY )
    call outfld('phis    ',state%phis,    pcols,   lchnk     )
#endif

!
! Add height of surface to midpoint height above surface 
!
    do k = 1, pver
       z3(:ncol,k) = state%zm(:ncol,k) + state%phis(:ncol)*rga
    end do
    call outfld('Z3      ',z3,pcols,lchnk)
!           
! Output Z3 on pressure surfaces
!
    if (hist_fld_active('Z1000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 100000._r8, z3, p_surf)
       call outfld('Z1000    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z975')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 97500._r8, z3, p_surf)
       call outfld('Z975    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z950')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 95000._r8, z3, p_surf)
       call outfld('Z950    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z925')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 92500._r8, z3, p_surf)
       call outfld('Z925    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z900')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 90000._r8, z3, p_surf)
       call outfld('Z900    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z850')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 85000._r8, z3, p_surf)
       call outfld('Z850    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z800')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 80000._r8, z3, p_surf)
       call outfld('Z800    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z700')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 70000._r8, z3, p_surf)
       call outfld('Z700    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z600')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 60000._r8, z3, p_surf)
       call outfld('Z600    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z500')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 50000._r8, z3, p_surf)
       call outfld('Z500    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z400')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 40000._r8, z3, p_surf)
       call outfld('Z400    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z300')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 30000._r8, z3, p_surf)
       call outfld('Z300    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z200')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 20000._r8, z3, p_surf)
       call outfld('Z200    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z100')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 10000._r8, z3, p_surf)
       call outfld('Z100    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z050')) then
       call vertinterp(ncol, pcols, pver, state%pmid,  5000._r8, z3, p_surf)
       call outfld('Z050    ', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('Z010')) then
       call vertinterp(ncol, pcols, pver, state%pmid,  1000._r8, z3, p_surf)
       call outfld('Z010    ', p_surf, pcols, lchnk)
    end if
!
! Quadratic height fiels Z3*Z3
!
    ftem(:ncol,:) = z3(:ncol,:)*z3(:ncol,:)
    call outfld('ZZ      ',ftem,pcols,lchnk)

    ftem(:ncol,:) = z3(:ncol,:)*state%v(:ncol,:)*gravit
    call outfld('VZ      ',ftem,  pcols,lchnk)
!
! Meridional advection fields
!
    ftem(:ncol,:) = state%v(:ncol,:)*state%t(:ncol,:)
    call outfld ('VT      ',ftem    ,pcols   ,lchnk     )

    ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,1)
    call outfld ('VQ      ',ftem    ,pcols   ,lchnk     )
    if(prog_modal_aero) then !Only for prognostic aerosols
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,14)
       call outfld ('Vbc_a1  ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,15)
       call outfld ('Vdst_a1 ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,22)
       call outfld ('Vdst_a3 ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,16)
       call outfld ('Vncl_a1 ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,20)
       call outfld ('Vncl_a2 ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,23)
       call outfld ('Vncl_a3 ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,11)
       call outfld ('Vso4_a1 ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,18)
       call outfld ('Vso4_a2 ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,24)
       call outfld ('Vso4_a3 ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,13)
       call outfld ('Vsoa_a1 ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,19)
       call outfld ('Vsoa_a2 ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,12)
       call outfld ('Vpom_a1 ',ftem    ,pcols   ,lchnk     )
    endif

    ftem(:ncol,:) = state%q(:ncol,:,1)*state%q(:ncol,:,1)
    call outfld ('QQ      ',ftem    ,pcols   ,lchnk     )

    if(prog_modal_aero) then !Only for prognostic aerosols
       ftem(:ncol,:) = state%q(:ncol,:,14)*state%q(:ncol,:,14)
       call outfld ('bc_a1_2 ',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%q(:ncol,:,15)*state%q(:ncol,:,15)
       call outfld ('dst_a1_2',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%q(:ncol,:,22)*state%q(:ncol,:,22)
       call outfld ('dst_a3_2',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%q(:ncol,:,16)*state%q(:ncol,:,16)
       call outfld ('ncl_a1_2',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%q(:ncol,:,20)*state%q(:ncol,:,20)
       call outfld ('ncl_a2_2',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%q(:ncol,:,23)*state%q(:ncol,:,23)
       call outfld ('ncl_a3_2',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%q(:ncol,:,11)*state%q(:ncol,:,11)
       call outfld ('so4_a1_2',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%q(:ncol,:,18)*state%q(:ncol,:,18)
       call outfld ('so4_a2_2',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%q(:ncol,:,24)*state%q(:ncol,:,24)
       call outfld ('so4_a3_2',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%q(:ncol,:,13)*state%q(:ncol,:,13)
       call outfld ('soa_a1_2',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%q(:ncol,:,19)*state%q(:ncol,:,19)
       call outfld ('soa_a2_2',ftem    ,pcols   ,lchnk     )
       
       ftem(:ncol,:) = state%q(:ncol,:,12)*state%q(:ncol,:,12)
       call outfld ('pom_a1_2',ftem    ,pcols   ,lchnk     )
    endif

    ftem(:ncol,:) = state%v(:ncol,:)**2
    call outfld ('VV      ',ftem    ,pcols   ,lchnk     )

    ftem(:ncol,:) = state%v(:ncol,:) * state%u(:ncol,:)
    call outfld ('VU      ',ftem    ,pcols   ,lchnk     )

! zonal advection

    ftem(:ncol,:) = state%u(:ncol,:)**2
    call outfld ('UU      ',ftem    ,pcols   ,lchnk     )

! Wind speed
    ftem(:ncol,:) = sqrt( state%u(:ncol,:)**2 + state%v(:ncol,:)**2)
    call outfld ('WSPEED  ',ftem    ,pcols   ,lchnk     )
    call outfld ('WSPDSRFMX',ftem(:,pver)   ,pcols   ,lchnk     )
    call outfld ('WSPDSRFAV',ftem(:,pver)   ,pcols   ,lchnk     ) 

! Vertical velocity and advection

    if (single_column) then
       call outfld('OMEGA   ',wfld,    pcols,   lchnk     )
    else
       call outfld('OMEGA   ',state%omega,    pcols,   lchnk     )
    endif

#if (defined E3SM_SCM_REPLAY )
    call outfld('omega   ',state%omega,    pcols,   lchnk     )
#endif

    ftem(:ncol,:) = state%omega(:ncol,:)*state%t(:ncol,:)
    call outfld('OMEGAT  ',ftem,    pcols,   lchnk     )
    ftem(:ncol,:) = state%omega(:ncol,:)*state%u(:ncol,:)
    call outfld('OMEGAU  ',ftem,    pcols,   lchnk     )
    ftem(:ncol,:) = state%omega(:ncol,:)*state%v(:ncol,:)
    call outfld('OMEGAV  ',ftem,    pcols,   lchnk     )
    ftem(:ncol,:) = state%omega(:ncol,:)*state%q(:ncol,:,1)
    call outfld('OMEGAQ  ',ftem,    pcols,   lchnk     )
    ftem(:ncol,:) = state%omega(:ncol,:)*state%omega(:ncol,:)
    call outfld('OMGAOMGA',ftem,    pcols,   lchnk     )
!
! Output omega at 850 and 500 mb pressure levels
!
    if (hist_fld_active('OMEGA1000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 100000._r8, state%omega, p_surf)
       call outfld('OMEGA1000', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA975')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 97500._r8, state%omega, p_surf)
       call outfld('OMEGA975', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA950')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 95000._r8, state%omega, p_surf)
       call outfld('OMEGA950', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA925')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 92500._r8, state%omega, p_surf)
       call outfld('OMEGA925', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA900')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 90000._r8, state%omega, p_surf)
       call outfld('OMEGA900', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA850')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 85000._r8, state%omega, p_surf)
       call outfld('OMEGA850', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA800')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 80000._r8, state%omega, p_surf)
       call outfld('OMEGA800', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA700')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 70000._r8, state%omega, p_surf)
       call outfld('OMEGA700', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA600')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 60000._r8, state%omega, p_surf)
       call outfld('OMEGA600', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA500')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 50000._r8, state%omega, p_surf)
       call outfld('OMEGA500', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA400')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 40000._r8, state%omega, p_surf)
       call outfld('OMEGA400', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA300')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 30000._r8, state%omega, p_surf)
       call outfld('OMEGA300', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA200')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 20000._r8, state%omega, p_surf)
       call outfld('OMEGA200', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGA100')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 10000._r8, state%omega, p_surf)
       call outfld('OMEGA100', p_surf, pcols, lchnk)
    end if
    if (hist_fld_active('OMEGABOT')) then
       call outfld('OMEGABOT', state%omega(:,pver), pcols, lchnk)
    end if
!     
! Mass of q, by layer and vertically integrated
!
    ftem(:ncol,:) = state%q(:ncol,:,1) * state%pdel(:ncol,:) * rga
    call outfld ('MQ      ',ftem    ,pcols   ,lchnk     )

    do k=2,pver
       ftem(:ncol,1) = ftem(:ncol,1) + ftem(:ncol,k)
    end do
    call outfld ('TMQ     ',ftem, pcols   ,lchnk     )
!
! Mass of vertically integrated water vapor flux
!
    ftem4(:ncol,:) = state%u(:ncol,:)*state%q(:ncol,:,1)*state%pdel(:ncol,:)*rga
    ftem5(:ncol,:) = state%v(:ncol,:)*state%q(:ncol,:,1)*state%pdel(:ncol,:)*rga
    do k=2,pver
       ftem4(:ncol,1) = ftem4(:ncol,1) + ftem4(:ncol,k)
       ftem5(:ncol,1) = ftem5(:ncol,1) + ftem5(:ncol,k)
    end do

    ftem(:ncol,1) = sqrt( ftem4(:ncol,1)**2 + ftem5(:ncol,1)**2)

    call outfld ('TUQ     ',ftem4, pcols   ,lchnk     )
    call outfld ('TVQ     ',ftem5, pcols   ,lchnk     )
    call outfld ('TTQ     ',ftem, pcols   ,lchnk     )

!
! Mass of vertically integrated MSE flux
!
    ftem(:ncol,:) = state%u(:ncol,:)*(state%s(:ncol,:)+latvap*state%q(:ncol,:,1))*state%pdel(:ncol,:)*rga
    do k=2,pver
       ftem(:ncol,1) = ftem(:ncol,1) + ftem(:ncol,k)
    end do
    call outfld ('TUH     ',ftem, pcols   ,lchnk     )

    ftem(:ncol,:) = state%v(:ncol,:)*(state%s(:ncol,:)+latvap*state%q(:ncol,:,1))*state%pdel(:ncol,:)*rga
    do k=2,pver
       ftem(:ncol,1) = ftem(:ncol,1) + ftem(:ncol,k)
    end do
    call outfld ('TVH     ',ftem, pcols   ,lchnk     )

    if (moist_physics) then

       ! Relative humidity
       call qsat(state%t(:ncol,:), state%pmid(:ncol,:), &
            tem2(:ncol,:), ftem(:ncol,:))
       ftem(:ncol,:) = state%q(:ncol,:,1)/ftem(:ncol,:)*100._r8

       if (hist_fld_active('RELHUM')) then
          call outfld ('RELHUM  ',ftem    ,pcols   ,lchnk     )
       end if

       ! ftem defined above assuming RELHUM in one of the history tape
       if (hist_fld_active('RH1000')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 100000._r8, ftem, p_surf)
          call outfld('RH1000', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH975')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 97500._r8, ftem, p_surf)
          call outfld('RH975', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH950')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 95000._r8, ftem, p_surf)
          call outfld('RH950', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH925')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 92500._r8, ftem, p_surf)
          call outfld('RH925', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH900')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 90000._r8, ftem, p_surf)
          call outfld('RH900', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH850')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 85000._r8, ftem, p_surf)
          call outfld('RH850', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH800')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 80000._r8, ftem, p_surf)
          call outfld('RH800', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH700')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 70000._r8, ftem, p_surf)
          call outfld('RH700', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH600')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 60000._r8, ftem, p_surf)
          call outfld('RH600', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH500')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 50000._r8, ftem, p_surf)
          call outfld('RH500', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH400')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 40000._r8, ftem, p_surf)
          call outfld('RH400', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH300')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 30000._r8, ftem, p_surf)
          call outfld('RH300', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH200')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 20000._r8, ftem, p_surf)
          call outfld('RH200', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RH100')) then
          call vertinterp(ncol, pcols, pver, state%pmid, 10000._r8, ftem, p_surf)
          call outfld('RH100', p_surf, pcols, lchnk)
       end if
       if (hist_fld_active('RHBOT')) then
          call outfld('RHBOT', ftem(:,pver), pcols, lchnk)
       end if

       if (hist_fld_active('RHW') .or. hist_fld_active('RHI') .or. hist_fld_active('RHCFMIP') ) then
	  
          ! RH w.r.t liquid (water)
          call qsat_water (state%t(:ncol,:), state%pmid(:ncol,:), &
               esl(:ncol,:), ftem(:ncol,:))
          ftem(:ncol,:) = state%q(:ncol,:,1)/ftem(:ncol,:)*100._r8
          call outfld ('RHW  ',ftem    ,pcols   ,lchnk     )

          ! Convert to RHI (ice)
          do i=1,ncol
             do k=1,pver
                esi(i,k)=svp_ice(state%t(i,k))
                ftem1(i,k)=ftem(i,k)*esl(i,k)/esi(i,k)
             end do
          end do
          call outfld ('RHI  ',ftem1    ,pcols   ,lchnk     )

	  ! use temperature to decide if you populate with ftem (liquid, above 0 C) or ftem1 (ice, below 0 C)

	  ftem2(:ncol,:)=ftem(:ncol,:)

          do i=1,ncol
             do k=1,pver
		if (state%t(i,k) .gt. 273) then
                   ftem2(i,k)=ftem(i,k)  !!wrt water
 		else
                   ftem2(i,k)=ftem1(i,k) !!wrt ice
		end if
             end do
          end do
          
          call outfld ('RHCFMIP  ',ftem2    ,pcols   ,lchnk     )

       end if

    end if
!
! Sea level pressure
!
    if (present(psl) .or. hist_fld_active('PSL')) then
       call cpslec (ncol, state%pmid, state%phis, state%ps, state%t,psl_tmp, gravit, rair) 
       call outfld ('PSL     ',psl_tmp  ,pcols, lchnk     )
       if (present(psl)) then	
          psl(:ncol) = psl_tmp(:ncol)
       end if
    end if
!
! Output T,q,u,v fields on pressure surfaces
!
    if (hist_fld_active('T975')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 97500._r8, state%t, p_surf)
       call outfld('T975    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T950')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 95000._r8, state%t, p_surf)
       call outfld('T950    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T900')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 90000._r8, state%t, p_surf)
       call outfld('T900    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T850')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 85000._r8, state%t, p_surf)
       call outfld('T850    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T800')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 80000._r8, state%t, p_surf)
       call outfld('T800    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T600')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 60000._r8, state%t, p_surf)
       call outfld('T600    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T500')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 50000._r8, state%t, p_surf)
       call outfld('T500    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T400')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 40000._r8, state%t, p_surf)
       call outfld('T400    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T300')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 30000._r8, state%t, p_surf)
       call outfld('T300    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T200')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 20000._r8, state%t, p_surf)
       call outfld('T200    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T100')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 10000._r8, state%t, p_surf)
       call outfld('T100    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q975')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 97500._r8, state%q(1,1,1), p_surf)
       call outfld('Q975    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q950')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 95000._r8, state%q(1,1,1), p_surf)
       call outfld('Q950    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q900')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 90000._r8, state%q(1,1,1), p_surf)
       call outfld('Q900    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q850')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 85000._r8, state%q(1,1,1), p_surf)
       call outfld('Q850    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q800')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 80000._r8, state%q(1,1,1), p_surf)
       call outfld('Q800    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q700')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 70000._r8, state%q(1,1,1), p_surf)
       call outfld('Q700    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q600')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 60000._r8, state%q(1,1,1), p_surf)
       call outfld('Q600    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q500')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 50000._r8, state%q(1,1,1), p_surf)
       call outfld('Q500    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q400')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 40000._r8, state%q(1,1,1), p_surf)
       call outfld('Q400    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q300')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 30000._r8, state%q(1,1,1), p_surf)
       call outfld('Q300    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q200')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 20000._r8, state%q(1,1,1), p_surf)
       call outfld('Q200    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q100')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 10000._r8, state%q(1,1,1), p_surf)
       call outfld('Q100    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q050')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 5000._r8, state%q(1,1,1), p_surf)
       call outfld('Q050    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('Q010')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 1000._r8, state%q(1,1,1), p_surf)
       call outfld('Q010    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U1000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 100000._r8, state%u, p_surf)
       call outfld('U1000   ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U975')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 97500._r8, state%u, p_surf)
       call outfld('U975    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U950')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 95000._r8, state%u, p_surf)
       call outfld('U950    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U925')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 92500._r8, state%u, p_surf)
       call outfld('U925    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U900')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 90000._r8, state%u, p_surf)
       call outfld('U900    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U850')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 85000._r8, state%u, p_surf)
       call outfld('U850    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U800')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 80000._r8, state%u, p_surf)
       call outfld('U800    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U700')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 70000._r8, state%u, p_surf)
       call outfld('U700    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U600')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 60000._r8, state%u, p_surf)
       call outfld('U600    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U500')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 50000._r8, state%u, p_surf)
       call outfld('U500    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U400')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 40000._r8, state%u, p_surf)
       call outfld('U400    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U300')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 30000._r8, state%u, p_surf)
       call outfld('U300    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U250')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 25000._r8, state%u, p_surf)
       call outfld('U250    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U200')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 20000._r8, state%u, p_surf)
       call outfld('U200    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U100')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 10000._r8, state%u, p_surf)
       call outfld('U100    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U050')) then
       call vertinterp(ncol, pcols, pver, state%pmid,  5000._r8, state%u, p_surf)
       call outfld('U050    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U010')) then
       call vertinterp(ncol, pcols, pver, state%pmid,  1000._r8, state%u, p_surf)
       call outfld('U010    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V1000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 100000._r8, state%v, p_surf)
       call outfld('V1000   ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V975')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 97500._r8, state%v, p_surf)
       call outfld('V975    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V950')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 95000._r8, state%v, p_surf)
       call outfld('V950    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V925')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 92500._r8, state%v, p_surf)
       call outfld('V925    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V900')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 90000._r8, state%v, p_surf)
       call outfld('V900    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V850')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 85000._r8, state%v, p_surf)
       call outfld('V850    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V800')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 80000._r8, state%v, p_surf)
       call outfld('V800    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V700')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 70000._r8, state%v, p_surf)
       call outfld('V700    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V600')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 60000._r8, state%v, p_surf)
       call outfld('V600    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V500')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 50000._r8, state%v, p_surf)
       call outfld('V500    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V400')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 40000._r8, state%v, p_surf)
       call outfld('V400    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V300')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 30000._r8, state%v, p_surf)
       call outfld('V300    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V250')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 25000._r8, state%v, p_surf)
       call outfld('V250    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V200')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 20000._r8, state%v, p_surf)
       call outfld('V200    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V100')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 10000._r8, state%v, p_surf)
       call outfld('V100    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V050')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 5000._r8, state%v, p_surf)
       call outfld('V050    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V010')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 1000._r8, state%v, p_surf)
       call outfld('V010    ', p_surf, pcols, lchnk )
    end if

    if (hist_fld_active('U90M')) then
       call vertinterpz(ncol, pcols, pver, state%zm, 90._r8, state%u, p_surf)
       call outfld('U90M    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('V90M')) then
       call vertinterpz(ncol, pcols, pver, state%zm, 90._r8, state%v, p_surf)
       call outfld('V90M    ', p_surf, pcols, lchnk )
    endif

    ftem(:ncol,:) = state%t(:ncol,:)*state%t(:ncol,:)
    call outfld('TT      ',ftem    ,pcols   ,lchnk   )
!
! Output U, V, T, Q, P and Z at bottom level
!
    call outfld ('UBOT    ', state%u(1,pver)  ,  pcols, lchnk)
    call outfld ('VBOT    ', state%v(1,pver)  ,  pcols, lchnk)
    call outfld ('QBOT    ', state%q(1,pver,1),  pcols, lchnk)
    call outfld ('ZBOT    ', state%zm(1,pver) , pcols, lchnk)

! Total energy of the atmospheric column for atmospheric heat storage calculations

    !! temporary variable to get surface geopotential in dimensions of (ncol,pver)
    do k=1,pver
      ftem1(:ncol,k)=state%phis(:ncol)  !! surface geopotential in units (m2/s2)
    end do

    !! calculate sum of sensible, kinetic, latent, and surface geopotential energy
    !! E=CpT+PHIS+Lv*q+(0.5)*(u^2+v^2)
    ftem(:ncol,:) = (cpair*state%t(:ncol,:) +  ftem1(:ncol,:) + latvap*state%q(:ncol,:,1) + &
         0.5_r8*(state%u(:ncol,:)**2+state%v(:ncol,:)**2))*(state%pdel(:ncol,:)/gravit)
    !! vertically integrate
    do k=2,pver       
	ftem(:ncol,1) = ftem(:ncol,1) + ftem(:ncol,k)
    end do
    call outfld ('ATMEINT   ',ftem(:ncol,1)  ,pcols   ,lchnk     )

!! Boundary layer atmospheric stability, temperature, water vapor diagnostics

    if (hist_fld_active('T1000')      .or. &
        hist_fld_active('T9251000')   .or. & 
        hist_fld_active('TH9251000')  .or. &
        hist_fld_active('THE9251000') .or. &
        hist_fld_active('T8501000')   .or. &
        hist_fld_active('TH8501000')  .or. &
        hist_fld_active('THE8501000') .or. &
        hist_fld_active('T7001000')   .or. &
        hist_fld_active('TH7001000')  .or. &
        hist_fld_active('THE7001000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 100000._r8, state%t, p_surf_t1)
    end if

    if (hist_fld_active('T925')       .or. &
        hist_fld_active('T9251000')   .or. & 
        hist_fld_active('TH9251000')  .or. &
        hist_fld_active('THE9251000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 92500._r8, state%t, p_surf_t2)
    end if

    if (hist_fld_active('Q1000')      .or. &
        hist_fld_active('THE9251000') .or. &
        hist_fld_active('THE8501000') .or. &
        hist_fld_active('THE7001000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 100000._r8, state%q(1,1,1), p_surf_q1)
    end if

    if (hist_fld_active('Q925')       .or. &
        hist_fld_active('THE9251000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 92500._r8, state%q(1,1,1), p_surf_q2)
    end if

    !!! at 1000 mb and 925 mb
    if (hist_fld_active('T1000')) then
       call outfld('T1000    ', p_surf_t1, pcols, lchnk )
    end if

    if (hist_fld_active('T925')) then
       call outfld('T925    ', p_surf_t2, pcols, lchnk )
    end if

    if (hist_fld_active('Q1000')) then
       call outfld('Q1000    ', p_surf_q1, pcols, lchnk ) 
    end if

    if (hist_fld_active('Q925')) then
       call outfld('Q925    ', p_surf_q2, pcols, lchnk )  
    end if

    if (hist_fld_active('T9251000')) then
       p_surf = p_surf_t2-p_surf_t1  
       call outfld('T9251000    ', p_surf, pcols, lchnk ) 
    end if

    if (hist_fld_active('TH9251000')) then
       p_surf = (p_surf_t2*(1000.0_r8/925.0_r8)**cappa)-(p_surf_t1*(1.0_r8)**cappa)   
       call outfld('TH9251000    ', p_surf, pcols, lchnk )    
    end if

    if (hist_fld_active('THE9251000')) then
       p_surf = (p_surf_t2*(1000.0_r8/925.0_r8)**cappa)*exp((2500000.0_r8*p_surf_q2)/(1004.0_r8*p_surf_t2))- &
            (p_surf_t1*(1.0_r8)**cappa)*exp((2500000.0_r8*p_surf_q1)/(1004.0_r8*p_surf_t1))  
       call outfld('THE9251000    ', p_surf, pcols, lchnk )
    end if

    if (hist_fld_active('T8501000')  .or. &
        hist_fld_active('TH8501000') .or. &
        hist_fld_active('THE8501000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 85000._r8, state%t, p_surf_t2)
    end if

    !!! at 1000 mb and 850 mb
    if (hist_fld_active('T8501000')) then
       p_surf = p_surf_t2-p_surf_t1  
       call outfld('T8501000    ', p_surf, pcols, lchnk ) 
    end if

    if (hist_fld_active('TH8501000')) then
       p_surf = (p_surf_t2*(1000.0_r8/850.0_r8)**cappa)-(p_surf_t1*(1.0_r8)**cappa)   
       call outfld('TH8501000    ', p_surf, pcols, lchnk )   
    end if

    if (hist_fld_active('THE8501000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 85000._r8, state%q(1,1,1), p_surf_q2)
       p_surf = (p_surf_t2*(1000.0_r8/850.0_r8)**cappa)*exp((2500000.0_r8*p_surf_q2)/(1004.0_r8*p_surf_t2))- &
            (p_surf_t1*(1.0_r8)**cappa)*exp((2500000.0_r8*p_surf_q1)/(1004.0_r8*p_surf_t1))  
       call outfld('THE8501000    ', p_surf, pcols, lchnk ) 
    end if

    if (hist_fld_active('T7001000')  .or. &
        hist_fld_active('TH7001000') .or. &
        hist_fld_active('T700') .or. &
        hist_fld_active('THE7001000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 70000._r8, state%t, p_surf_t2)
    end if

   !!! at 700 mb
    if (hist_fld_active('T700')) then
       call outfld('T700    ', p_surf_t2, pcols, lchnk )
    end if

    !!! at 1000 mb and 700 mb
    if (hist_fld_active('T7001000')) then
       p_surf = p_surf_t2-p_surf_t1
       call outfld('T7001000    ', p_surf, pcols, lchnk )
    end if

    if (hist_fld_active('TH7001000')) then
       p_surf = (p_surf_t2*(1000.0_r8/700.0_r8)**cappa)-(p_surf_t1*(1.0_r8)**cappa)
       call outfld('TH7001000    ', p_surf, pcols, lchnk )
    end if

    if (hist_fld_active('THE7001000')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 70000._r8, state%q(1,1,1), p_surf_q2)
       p_surf = (p_surf_t2*(1000.0_r8/700.0_r8)**cappa)*exp((2500000.0_r8*p_surf_q2)/(1004.0_r8*p_surf_t2))- &
            (p_surf_t1*(1.0_r8)**cappa)*exp((2500000.0_r8*p_surf_q1)/(1004.0_r8*p_surf_t1))
       call outfld('THE7001000    ', p_surf, pcols, lchnk )
    end if

    if (hist_fld_active('T250')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 25000._r8, state%t, p_surf)
       call outfld('T250           ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T150')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 15000._r8, state%t, p_surf)
       call outfld('T150           ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T050')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 5000._r8, state%t, p_surf)
       call outfld('T050           ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T025')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 2500._r8, state%t, p_surf)
       call outfld('T025           ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T010')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 1000._r8, state%t, p_surf)
       call outfld('T010           ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T005')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 500._r8, state%t, p_surf)
       call outfld('T005           ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T002')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 200._r8, state%t, p_surf)
       call outfld('T002           ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('T001')) then
       call vertinterp(ncol, pcols, pver, state%pmid, 100._r8, state%t, p_surf)
       call outfld('T001           ', p_surf, pcols, lchnk )
    end if
    call outfld ('TTOP    ', state%t(:,1)  ,  pcols, lchnk)
    if (hist_fld_active('U150')) then
       call vertinterp(ncol, pcols, pver, state%pmid,  15000._r8, state%u, p_surf)
       call outfld('U150    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U050')) then
       call vertinterp(ncol, pcols, pver, state%pmid,  5000._r8, state%u, p_surf)
       call outfld('U050    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U025')) then
       call vertinterp(ncol, pcols, pver, state%pmid,  2500._r8, state%u, p_surf)
       call outfld('U025    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U005')) then
       call vertinterp(ncol, pcols, pver, state%pmid,  500._r8, state%u, p_surf)
       call outfld('U005    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U002')) then
       call vertinterp(ncol, pcols, pver, state%pmid,  200._r8, state%u, p_surf)
       call outfld('U002    ', p_surf, pcols, lchnk )
    end if
    if (hist_fld_active('U001')) then
       call vertinterp(ncol, pcols, pver, state%pmid,  100._r8, state%u, p_surf)
       call outfld('U001    ', p_surf, pcols, lchnk )
    end if
    call outfld ('UTOP    ', state%u(:,1)  ,  pcols, lchnk)


  !---------------------------------------------------------
  ! WACCM tidal diagnostics 
  !---------------------------------------------------------

  if (chem_is('waccm_ghg') .or. chem_is('waccm_mozart') .or. chem_is('waccm_mozart_mam3')) then

    call tidal_diag_write(state) 

  endif

    return
  end subroutine diag_phys_writeout
!===============================================================================

subroutine diag_conv(state, ztodt, pbuf)

!----------------------------------------------------------------------- 
! 
! Output diagnostics associated with all convective processes.
!
!-----------------------------------------------------------------------
   use physconst,     only: cpair
   use tidal_diag,    only: get_tidal_coeffs
   use time_manager,  only: get_curr_date

! Arguments:

   real(r8),            intent(in) :: ztodt   ! timestep for computing physics tendencies
   type(physics_state), intent(in) :: state
   type(physics_buffer_desc), pointer :: pbuf(:)

! 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
   real(r8), pointer :: prec_sed(:)                ! total precipitation   from ZM convection
   real(r8), pointer :: snow_sed(:)                ! snow from ZM   convection
   real(r8), pointer :: prec_pcw(:)                ! total precipitation   from Hack convection
   real(r8), pointer :: snow_pcw(:)                ! snow from Hack   convection
   real(r8), pointer :: pblh(:)                    ! PBLH depth

   real(r8), pointer :: wsresp(:)                  ! first order response of winds to stress
   real(r8), pointer :: tau_est(:)                 ! estimated equilibrium stress

! Local variables:
   
   integer :: i, k, m, lchnk, ncol

   real(r8) :: rtdt

   real(r8):: precc(pcols)                ! convective precip rate
   real(r8):: precl(pcols)                ! stratiform precip rate
   real(r8):: snowc(pcols)                ! convective snow rate
   real(r8):: snowl(pcols)                ! stratiform snow rate
   real(r8):: prect(pcols)                ! total (conv+large scale) precip rate
   real(r8):: cape(pcols)                 ! CAPE
   real(r8):: cin(pcols)                  ! CIN
   integer :: pblt(pcols)                 ! Indicee of PBL top
   real(r8) :: dcoef(4)                   ! for tidal component of T tend

   real(r8):: prect1(pcols)
   real(r8):: precl1(pcols)
   real(r8):: precc1(pcols)
   real(r8):: adfprecl(pcols,nbands)
   real(r8):: adfprecc(pcols,nbands)
   real(r8):: adfprect(pcols,nbands)
   real(r8):: pdfprecl(pcols,nbands)
   real(r8):: pdfprecc(pcols,nbands)
   real(r8):: pdfprect(pcols,nbands)
   real(r8):: dprecl(pcols,nhrs)
   real(r8):: dprecc(pcols,nhrs)
   real(r8):: dprect(pcols,nhrs)
   real(r8):: tod1
   integer  :: year, month
   integer  :: day              ! day of month
   integer  :: tod              ! time of day (seconds past 0Z) 
   integer  :: itim_old

   lchnk = state%lchnk
   ncol  = state%ncol

   rtdt = 1._r8/ztodt

   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_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)
   call pbuf_get_field(pbuf, pblh_idx, pblh)
   if (linearize_pbl_winds) then
      call pbuf_get_field(pbuf, wsresp_idx, wsresp)
      call pbuf_get_field(pbuf, tau_est_idx, tau_est)
   end if

! Precipitation rates (multi-process)
   precc(:ncol) = prec_dp(:ncol)  + prec_sh(:ncol)
   precl(:ncol) = prec_sed(:ncol) + prec_pcw(:ncol)
   snowc(:ncol) = snow_dp(:ncol)  + snow_sh(:ncol)
   snowl(:ncol) = snow_sed(:ncol) + snow_pcw(:ncol)
   prect(:ncol) = precc(:ncol)    + precl(:ncol)

   call outfld('PRECC   ', precc, pcols, lchnk )
   call outfld('PRECL   ', precl, pcols, lchnk )
   call outfld('PREC_PCW', prec_pcw,pcols   ,lchnk )
   call outfld('PREC_zmc', prec_dp ,pcols   ,lchnk )
   call outfld('PRECSC  ', snowc, pcols, lchnk )
   call outfld('PRECSL  ', snowl, pcols, lchnk )
   call outfld('PRECT   ', prect, pcols, lchnk )
   call outfld('PRECTMX ', prect, pcols, lchnk )

   call outfld('PRECLav ', precl, pcols, lchnk )
   call outfld('PRECCav ', precc, pcols, lchnk )

   if (linearize_pbl_winds) then
      call outfld('wsresp', wsresp, pcols, lchnk)
      call outfld('tau_est', tau_est, pcols, lchnk)
   end if

#if ( defined E3SM_SCM_REPLAY )
   call outfld('Prec   ' , prect, pcols, lchnk )
#endif

   ! Total convection tendencies.

   do k = 1, pver
      do i = 1, ncol
         dtcond(i,k,lchnk) = (state%s(i,k) - dtcond(i,k,lchnk))*rtdt / cpair
      end do
   end do
   call outfld('DTCOND  ', dtcond(:,:,lchnk), pcols, lchnk)

   ! output tidal coefficients
   call get_tidal_coeffs( dcoef )
   call outfld( 'DTCOND_24_SIN', dtcond(:ncol,:,lchnk)*dcoef(1), ncol, lchnk )
   call outfld( 'DTCOND_24_COS', dtcond(:ncol,:,lchnk)*dcoef(2), ncol, lchnk )
   call outfld( 'DTCOND_12_SIN', dtcond(:ncol,:,lchnk)*dcoef(3), ncol, lchnk )
   call outfld( 'DTCOND_12_COS', dtcond(:ncol,:,lchnk)*dcoef(4), ncol, lchnk )

   do m = 1, dqcond_num
      if ( cnst_cam_outfld(m) ) then
         do k = 1, pver
            do i = 1, ncol
               dqcond(m)%cnst(i,k,lchnk) = (state%q(i,k,m) - dqcond(m)%cnst(i,k,lchnk))*rtdt
            end do
         end do
         call outfld(dcconnam(m), dqcond(m)%cnst(:,:,lchnk), pcols, lchnk)
      end if
   end do

   prect1  = prect*1.e3_r8*86400._r8
   precl1  = precl*1.e3_r8*86400._r8
   precc1  = precc*1.e3_r8*86400._r8
   pdfprecl= 0._r8
   pdfprecc= 0._r8
   pdfprect= 0._r8
   adfprecl= 0._r8
   adfprecc= 0._r8
   adfprect= 0._r8
   dprecl= 0._r8
   dprecc= 0._r8
   dprect= 0._r8

   call get_curr_date(year, month, day, tod)
   tod1=dble(tod)

   do i=1, ncol
     do k=1,nbands
        if (prect1(i).gt.bbands_1d(k).and.prect1(i).le.bbands_1d(k+1)) then
         pdfprect(i,k)=1._r8
         adfprect(i,k)=prect1(i)
        endif
        if (precl1(i).gt.bbands_1d(k).and.precl1(i).le.bbands_1d(k+1)) then
         pdfprecl(i,k)=1._r8
         adfprecl(i,k)=precl1(i)
        endif
        if (precc1(i).gt.bbands_1d(k).and.precc1(i).le.bbands_1d(k+1)) then
         pdfprecc(i,k)=1._r8
         adfprecc(i,k)=precc1(i)
        endif
     end do

     do k=1, nhrs
      if (tod1 >= bhrs_1d(k) .and. tod1 < bhrs_1d(k+1)) then 
         dprect(i,k)=prect1(i)
         dprecl(i,k)=precl1(i)
         dprecc(i,k)=precc1(i)
      endif
     end do
   end do

   call outfld('PDFPRECT', pdfprect * (1./0.07), pcols, lchnk )  ! bin width 0.07 
   call outfld('PDFPRECL', pdfprecl * (1./0.07), pcols, lchnk )
   call outfld('PDFPRECC', pdfprecc * (1./0.07), pcols, lchnk )
   call outfld('ADFPRECT', adfprect * (1./0.07), pcols, lchnk )  ! bin width 0.07 
   call outfld('ADFPRECL', adfprecl * (1./0.07), pcols, lchnk ) 
   call outfld('ADFPRECC', adfprecc * (1./0.07), pcols, lchnk )  
   call outfld('DPRECT', dprect * nhrs, pcols, lchnk )
   call outfld('DPRECL', dprecl * nhrs, pcols, lchnk )
   call outfld('DPRECC', dprecc * nhrs, pcols, lchnk )
   
   ! Add CAPE and CIN calculation here
   ! Note that this routine needs the input pressures to be in units of hPa 
   call diag_CAPEandCIN(ncol,state%q(:pcols,:pver,1),state%t,&
          0.01_r8*state%pmid,0.01_r8*state%pint,state%zm,state%zi,pblh,&
          state%phis,cape,cin)
          
   call outfld('CAPE', cape, pcols, lchnk )
   call outfld('CIN', cin, pcols, lchnk )

end subroutine diag_conv

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

! YQIN 08/18/23 diagnose EIS
!subroutine diag_surf (cam_in, cam_out, ps, trefmxav, trefmnav )
subroutine diag_surf(cam_in, cam_out, ps, pmid, phis, t, zm, trefmxav, trefmnav)

!----------------------------------------------------------------------- 
! 
! Purpose: record surface diagnostics
!
!-----------------------------------------------------------------------

   use time_manager,     only: is_end_curr_day
   use co2_cycle,        only: c_i, co2_transport
   use constituents,     only: sflxnam
   ! YQIN
   use physconst,          only: gravit, rga, rair, cpair, latvap, rearth, pi, cappa
   use interpolate_data,   only: vertinterp

!-----------------------------------------------------------------------
!
! Input arguments
!
    type(cam_in_t),  intent(in) :: cam_in
    type(cam_out_t), intent(in) :: cam_out

    real(r8), intent(inout) :: trefmnav(pcols) ! daily minimum tref  
    real(r8), intent(inout) :: trefmxav(pcols) ! daily maximum tref

    real(r8), intent(in)    :: ps(pcols)       ! Surface pressure.
    ! YQIN 
    real(r8), intent(in)    :: pmid(pcols)
    real(r8), intent(in)    :: phis(pcols)
    real(r8), intent(in)    :: t(pcols,pver)
    real(r8), intent(in)    :: zm(pcols,pver)

!
!---------------------------Local workspace-----------------------------
!
    integer :: i, k, m      ! indexes
    integer :: lchnk        ! chunk identifier
    integer :: ncol         ! longitude dimension
    real(r8) tem2(pcols)    ! temporary workspace
    real(r8) ftem(pcols)    ! temporary workspace

    ! YQIN
    real(r8) :: z3(pcols,pver)
    real(r8) :: TS(pcols), SLP(pcols), T700(pcols), Z700(pcols) 
    real(r8) :: LTS(pcols), EIS(pcols)
    real(r8) :: EIS_fxRH(pcols)
    real(r8) :: rh_ref(pcols) ! RH at the reference height

!
!-----------------------------------------------------------------------
!
    lchnk = cam_in%lchnk
    ncol  = cam_in%ncol

    call outfld('SHFLX',    cam_in%shf,       pcols, lchnk)
    call outfld('LHFLX',    cam_in%lhf,       pcols, lchnk)
    call outfld('QFLX',     cam_in%cflx(1,1), pcols, lchnk)

    call outfld('TAUX',     cam_in%wsx,       pcols, lchnk)
    call outfld('TAUY',     cam_in%wsy,       pcols, lchnk)
    call outfld('TREFHT  ', cam_in%tref,      pcols, lchnk)
    call outfld('TREFHTMX', cam_in%tref,      pcols, lchnk)
    call outfld('TREFHTMN', cam_in%tref,      pcols, lchnk)
    call outfld('QREFHT',   cam_in%qref,      pcols, lchnk)
    call outfld('U10',      cam_in%u10,       pcols, lchnk)
! 
! Calculate and output reference height RH (RHREFHT)

   call qsat(cam_in%tref(:ncol), ps(:ncol), tem2(:ncol), ftem(:ncol))
       ftem(:ncol) = cam_in%qref(:ncol)/ftem(:ncol)*100._r8

      
    call outfld('RHREFHT',   ftem,      pcols, lchnk)

    ! YQIN 08/18/23  EIS diagnostics 
    !! get reference height RH 
    rh_ref(:ncol) = ftem(:ncol)/100._r8 ! no unit

    !! get Z700 
    do k=1,pver
        z3(:ncol,k) = zm(:ncol,k) + phis(:ncol)*rga
    end do 
    call vertinterp(ncol, pcols, pver, pmid, 70000._r8, z3, Z700)

    !! get T700 
    call vertinterp(ncol, pcols, pver, pmid, 70000._r8, t, T700)

    !! get SLP 
    call cpslec (ncol, pmid, phis, ps, t, SLP, gravit, rair) 

    !! calculate EIS
    call calc_EIS(ncol, pcols, cam_in%tref, rh_ref, SLP, T700, Z700, LTS, EIS, EIS_fxRH)

    call outfld('LTS', LTS, pcols, lchnk)
    call outfld('EIS', EIS, pcols, lchnk)
    call outfld('EIS_fxRH', EIS_fxRH, pcols, lchnk) 


#if (defined E3SM_SCM_REPLAY )
    call outfld('shflx   ',cam_in%shf,   pcols,   lchnk)
    call outfld('lhflx   ',cam_in%lhf,   pcols,   lchnk)
    call outfld('trefht  ',cam_in%tref,  pcols,   lchnk)
    call outfld('Tg', cam_in%ts, pcols, lchnk)
    call outfld('Tsair',cam_in%ts, pcols, lchnk)
#endif
!
! Ouput ocn and ice fractions
!
    call outfld('LANDFRAC', cam_in%landfrac, pcols, lchnk)
    call outfld('ICEFRAC',  cam_in%icefrac,  pcols, lchnk)
    call outfld('OCNFRAC',  cam_in%ocnfrac,  pcols, lchnk)
!
! Compute daily minimum and maximum of TREF
!
    do i = 1,ncol
       trefmxav(i) = max(cam_in%tref(i),trefmxav(i))
       trefmnav(i) = min(cam_in%tref(i),trefmnav(i))
    end do
    if (is_end_curr_day()) then
       call outfld('TREFMXAV', trefmxav,pcols,   lchnk     )
       call outfld('TREFMNAV', trefmnav,pcols,   lchnk     )
       trefmxav(:ncol) = -1.0e36_r8
       trefmnav(:ncol) =  1.0e36_r8
    endif

    call outfld('TBOT',     cam_out%tbot,     pcols, lchnk)
    call outfld('TS',       cam_in%ts,        pcols, lchnk)
    call outfld('TSMN',     cam_in%ts,        pcols, lchnk)
    call outfld('TSMX',     cam_in%ts,        pcols, lchnk)
    call outfld('SNOWHLND', cam_in%snowhland, pcols, lchnk)
    call outfld('SNOWHICE', cam_in%snowhice,  pcols, lchnk)
    call outfld('ASDIR',    cam_in%asdir,     pcols, lchnk)
    call outfld('ASDIF',    cam_in%asdif,     pcols, lchnk)
    call outfld('ALDIR',    cam_in%aldir,     pcols, lchnk)
    call outfld('ALDIF',    cam_in%aldif,     pcols, lchnk)
    call outfld('SST',      cam_in%sst,       pcols, lchnk)

    if (co2_transport()) then
       do m = 1,4
          call outfld(sflxnam(c_i(m)), cam_in%cflx(:,c_i(m)), pcols, lchnk)
       end do
    end if

end subroutine diag_surf

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

subroutine diag_export(cam_out)

!----------------------------------------------------------------------- 
! 
! Purpose: Write export state to history file
!
!-----------------------------------------------------------------------

   ! arguments
   type(cam_out_t), intent(inout) :: cam_out

   ! Local variables:
   integer :: lchnk        ! chunk identifier
   logical :: atm_dep_flux ! true ==> sending deposition fluxes to coupler.
                           ! Otherwise, set them to zero.
   !-----------------------------------------------------------------------

   lchnk = cam_out%lchnk

   call phys_getopts(atm_dep_flux_out=atm_dep_flux)

   if (.not. atm_dep_flux) then
      ! set the fluxes to zero before outfld and sending them to the
      ! coupler
      cam_out%bcphiwet = 0.0_r8
      cam_out%bcphidry = 0.0_r8
      cam_out%bcphodry = 0.0_r8
      cam_out%ocphiwet = 0.0_r8
      cam_out%ocphidry = 0.0_r8
      cam_out%ocphodry = 0.0_r8
      cam_out%dstwet1  = 0.0_r8
      cam_out%dstdry1  = 0.0_r8
      cam_out%dstwet2  = 0.0_r8
      cam_out%dstdry2  = 0.0_r8
      cam_out%dstwet3  = 0.0_r8
      cam_out%dstdry3  = 0.0_r8
      cam_out%dstwet4  = 0.0_r8
      cam_out%dstdry4  = 0.0_r8
   end if

   call outfld('a2x_BCPHIWET', cam_out%bcphiwet, pcols, lchnk)
   call outfld('a2x_BCPHIDRY', cam_out%bcphidry, pcols, lchnk)
   call outfld('a2x_BCPHODRY', cam_out%bcphodry, pcols, lchnk)
   call outfld('a2x_OCPHIWET', cam_out%ocphiwet, pcols, lchnk)
   call outfld('a2x_OCPHIDRY', cam_out%ocphidry, pcols, lchnk)
   call outfld('a2x_OCPHODRY', cam_out%ocphodry, pcols, lchnk)
   call outfld('a2x_DSTWET1',  cam_out%dstwet1,  pcols, lchnk)
   call outfld('a2x_DSTDRY1',  cam_out%dstdry1,  pcols, lchnk)
   call outfld('a2x_DSTWET2',  cam_out%dstwet2,  pcols, lchnk)
   call outfld('a2x_DSTDRY2',  cam_out%dstdry2,  pcols, lchnk)
   call outfld('a2x_DSTWET3',  cam_out%dstwet3,  pcols, lchnk)
   call outfld('a2x_DSTDRY3',  cam_out%dstdry3,  pcols, lchnk)
   call outfld('a2x_DSTWET4',  cam_out%dstwet4,  pcols, lchnk)
   call outfld('a2x_DSTDRY4',  cam_out%dstdry4,  pcols, lchnk)

end subroutine diag_export

!#######################################################################

   subroutine diag_physvar_ic (lchnk,  pbuf, cam_out, cam_in)
!
!---------------------------------------------
!
! Purpose: record physics variables on IC file
!
!---------------------------------------------
!

!
! Arguments
!
   integer       , intent(in) :: lchnk  ! chunk identifier
   type(physics_buffer_desc), pointer :: pbuf(:)
   
   type(cam_out_t), intent(inout) :: cam_out
   type(cam_in_t),  intent(inout) :: cam_in 
!
!---------------------------Local workspace-----------------------------
!
   integer  :: k                 ! indices
   integer  :: itim_old          ! indices

   real(r8), pointer, dimension(:,:) :: cwat_var
   real(r8), pointer, dimension(:,:) :: conv_var_3d
   real(r8), pointer, dimension(:  ) :: conv_var_2d
   real(r8), pointer :: tpert(:), pblh(:), qpert(:)
!
!-----------------------------------------------------------------------
!
   if( write_inithist() ) then

      !
      ! Associate pointers with physics buffer fields
      !
      itim_old = pbuf_old_tim_idx()

      if (qcwat_idx > 0) then
         call pbuf_get_field(pbuf, qcwat_idx, cwat_var, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         call outfld('QCWAT&IC   ',cwat_var, pcols,lchnk)
      end if

      if (tcwat_idx > 0) then
         call pbuf_get_field(pbuf, tcwat_idx,  cwat_var, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         call outfld('TCWAT&IC   ',cwat_var, pcols,lchnk)
      end if

      if (lcwat_idx > 0) then
         call pbuf_get_field(pbuf, lcwat_idx,  cwat_var, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         call outfld('LCWAT&IC   ',cwat_var, pcols,lchnk)
      end if

      call pbuf_get_field(pbuf, cld_idx,    cwat_var, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
      call outfld('CLOUD&IC   ',cwat_var, pcols,lchnk)

      call pbuf_get_field(pbuf, concld_idx, cwat_var, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
      call outfld('CONCLD&IC   ',cwat_var, pcols,lchnk)

      call pbuf_get_field(pbuf, tke_idx, conv_var_3d)
      call outfld('TKE&IC    ',conv_var_3d, pcols,lchnk)

      call pbuf_get_field(pbuf, kvm_idx,  conv_var_3d)
      call outfld('KVM&IC    ',conv_var_3d, pcols,lchnk)

      call pbuf_get_field(pbuf, kvh_idx,  conv_var_3d)
      call outfld('KVH&IC    ',conv_var_3d, pcols,lchnk)
 
      call pbuf_get_field(pbuf, cush_idx, conv_var_2d ,(/1,itim_old/),  (/pcols,1/))
      call outfld('CUSH&IC   ',conv_var_2d, pcols,lchnk)

      if (qpert_idx > 0) then 
         call pbuf_get_field(pbuf, qpert_idx, qpert)
         call outfld('QPERT&IC   ', qpert, pcols, lchnk)
      end if

      call pbuf_get_field(pbuf, pblh_idx,  pblh)
      call outfld('PBLH&IC    ', pblh,  pcols, lchnk)

      call pbuf_get_field(pbuf, tpert_idx, tpert)
      call outfld('TPERT&IC   ', tpert, pcols, lchnk)

      ! The following is only needed for cam-csim
      call outfld('TBOT&IC    ', cam_out%tbot, pcols, lchnk)
   end if

   end subroutine diag_physvar_ic


!#######################################################################

subroutine diag_phys_tend_writeout(state, pbuf,  tend, ztodt, tmp_q, tmp_cldliq, tmp_cldice, &
                                   tmp_t, qini, cldliqini, cldiceini)

   !---------------------------------------------------------------
   !
   ! Purpose:  Dump physics tendencies for moisture and temperature
   !
   !---------------------------------------------------------------

   use check_energy,    only: check_energy_get_integrals
   use physconst,       only: cpair
   
   ! Arguments

   type(physics_state), intent(in   ) :: state 
   
   type(physics_buffer_desc), pointer :: pbuf(:)
   type(physics_tend ), intent(in   ) :: tend
   real(r8)           , intent(in   ) :: ztodt                  ! physics timestep
   real(r8)           , intent(inout) :: tmp_q     (pcols,pver) ! As input, holds pre-adjusted tracers (FV)
   real(r8)           , intent(inout) :: tmp_cldliq(pcols,pver) ! As input, holds pre-adjusted tracers (FV)
   real(r8)           , intent(inout) :: tmp_cldice(pcols,pver) ! As input, holds pre-adjusted tracers (FV)
   real(r8)           , intent(inout) :: tmp_t     (pcols,pver) ! holds last physics_updated T (FV)
   real(r8)           , intent(in   ) :: qini      (pcols,pver) ! tracer fields at beginning of physics
   real(r8)           , intent(in   ) :: cldliqini (pcols,pver) ! tracer fields at beginning of physics
   real(r8)           , intent(in   ) :: cldiceini (pcols,pver) ! tracer fields at beginning of physics

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

   integer  :: m      ! constituent index
   integer  :: lchnk  ! chunk index
   integer  :: ncol   ! number of columns in chunk
   real(r8) :: ftem2(pcols     ) ! Temporary workspace for outfld variables
   real(r8) :: ftem3(pcols,pver) ! Temporary workspace for outfld variables
   real(r8) :: rtdt
   real(r8) :: heat_glob         ! global energy integral (FV only)
   integer  :: ixcldice, ixcldliq! constituent indices for cloud liquid and ice water.
   ! CAM pointers to get variables from the physics buffer
   real(r8), pointer, dimension(:,:) :: t_ttend  
   integer  :: itim_old

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

   lchnk = state%lchnk
   ncol  = state%ncol
   rtdt  = 1._r8/ztodt
   call cnst_get_ind('CLDLIQ', ixcldliq)
   call cnst_get_ind('CLDICE', ixcldice)

   ! Dump out post-physics state (FV only)

   if (dycore_is('LR') .or. dycore_is('SE')) then
      tmp_t(:ncol,:pver) = (tmp_t(:ncol,:pver) - state%t(:ncol,:pver))/ztodt
      call outfld('PTTEND_RESID', tmp_t, pcols, lchnk   )
   end if
   call outfld('TAP', state%t, pcols, lchnk   )
   call outfld('UAP', state%u, pcols, lchnk   )
   call outfld('VAP', state%v, pcols, lchnk   )

   if ( cnst_cam_outfld(       1) ) call outfld (apcnst(       1), state%q(1,1,       1), pcols, lchnk)
   if ( cnst_cam_outfld(ixcldliq) ) call outfld (apcnst(ixcldliq), state%q(1,1,ixcldliq), pcols, lchnk)
   if ( cnst_cam_outfld(ixcldice) ) call outfld (apcnst(ixcldice), state%q(1,1,ixcldice), pcols, lchnk)

   ! T-tendency due to FV Energy fixer (remove from total physics tendency diagnostic)

   if (dycore_is('LR') .or. dycore_is('SE')) then
      call check_energy_get_integrals( heat_glob_out=heat_glob )
      ftem2(:ncol)  = heat_glob/cpair
      call outfld('TFIX', ftem2, pcols, lchnk   )
      ftem3(:ncol,:pver)  = tend%dtdt(:ncol,:pver) - heat_glob/cpair
   else
      ftem3(:ncol,:pver)  = tend%dtdt(:ncol,:pver)
   end if

   ! Total physics tendency for Temperature

   call outfld('PTTEND',ftem3, pcols, lchnk )

   ! Tendency for dry mass adjustment of q (valid for FV only)

   if (dycore_is('LR')) then
      tmp_q     (:ncol,:pver) = (state%q(:ncol,:pver,       1) - tmp_q     (:ncol,:pver))*rtdt
      tmp_cldliq(:ncol,:pver) = (state%q(:ncol,:pver,ixcldliq) - tmp_cldliq(:ncol,:pver))*rtdt
      tmp_cldice(:ncol,:pver) = (state%q(:ncol,:pver,ixcldice) - tmp_cldice(:ncol,:pver))*rtdt
      if ( cnst_cam_outfld(       1) ) call outfld (dmetendnam(       1), tmp_q     , pcols, lchnk)
      if ( cnst_cam_outfld(ixcldliq) ) call outfld (dmetendnam(ixcldliq), tmp_cldliq, pcols, lchnk)
      if ( cnst_cam_outfld(ixcldice) ) call outfld (dmetendnam(ixcldice), tmp_cldice, pcols, lchnk)
   end if

   ! Total physics tendency for moisture and other tracers

   if ( cnst_cam_outfld(       1) ) then
      ftem3(:ncol,:pver) = (state%q(:ncol,:pver,       1) - qini     (:ncol,:pver) )*rtdt
      call outfld (ptendnam(       1), ftem3, pcols, lchnk)
   end if
   if ( cnst_cam_outfld(ixcldliq) ) then
      ftem3(:ncol,:pver) = (state%q(:ncol,:pver,ixcldliq) - cldliqini(:ncol,:pver) )*rtdt
      call outfld (ptendnam(ixcldliq), ftem3, pcols, lchnk)
   end if
   if ( cnst_cam_outfld(ixcldice) ) then
      ftem3(:ncol,:pver) = (state%q(:ncol,:pver,ixcldice) - cldiceini(:ncol,:pver) )*rtdt
      call outfld (ptendnam(ixcldice), ftem3, pcols, lchnk)
   end if

   ! Total (physics+dynamics, everything!) tendency for Temperature

   !! get temperature stored in physics buffer
   itim_old = pbuf_old_tim_idx()
   call pbuf_get_field(pbuf, t_ttend_idx, t_ttend, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))

   !! calculate and outfld the total temperature tendency
   ftem3(:ncol,:) = (state%t(:ncol,:) - t_ttend(:ncol,:))/ztodt
   call outfld('TTEND_TOT', ftem3, pcols, lchnk)

   !! update physics buffer with this time-step's temperature
   t_ttend(:ncol,:) = state%t(:ncol,:)

end subroutine diag_phys_tend_writeout

!#######################################################################

   subroutine diag_state_b4_phys_write (state)
!
!---------------------------------------------------------------
!
! Purpose:  Dump state just prior to executing physics
!
!---------------------------------------------------------------
!
! Arguments
!
   type(physics_state), intent(in) :: state 
!
!---------------------------Local workspace-----------------------------
!
   integer :: ixcldice, ixcldliq ! constituent indices for cloud liquid and ice water.
   integer :: lchnk              ! chunk index
!
!-----------------------------------------------------------------------
!
   lchnk = state%lchnk

   call cnst_get_ind('CLDLIQ', ixcldliq)
   call cnst_get_ind('CLDICE', ixcldice)
   call outfld('TBP', state%t, pcols, lchnk   )
   if ( cnst_cam_outfld(       1) ) call outfld (bpcnst(       1), state%q(1,1,       1), pcols, lchnk)
   if ( cnst_cam_outfld(ixcldliq) ) call outfld (bpcnst(ixcldliq), state%q(1,1,ixcldliq), pcols, lchnk)
   if ( cnst_cam_outfld(ixcldice) ) call outfld (bpcnst(ixcldice), state%q(1,1,ixcldice), pcols, lchnk)

   end subroutine diag_state_b4_phys_write

!#######################################################################

   subroutine diag_CAPEandCIN(ncol, &
                     q, t, p, pf, zm, zi, pblh, phis, &
                     cape, cin)

   !-----------------------------------------------------------------------
   !
   ! Purpose:
   ! Compute CAPE and CIN for diagnostic output only
   !
   ! Method:
   ! Originated from the CAM3.5 subroutine "buoyan" that computed CAPE
   !  in the ZM convection scheme, with modifications.  Also extended
   !  to compute and output CIN
   !
   ! Author:
   ! Original:          P. Rasch, April 1996
   ! Modifications:     P. Bogenschutz, 2020
   !
   !-----------------------------------------------------------------------
      implicit none
   !-----------------------------------------------------------------------

   !  INPUTS arguments
      integer, intent(in) :: ncol             ! number of atmospheric columns

      real(r8), intent(in) :: q(pcols,pver)   ! spec. humidity [kg/kg]
      real(r8), intent(in) :: t(pcols,pver)   ! temperature [K]
      real(r8), intent(in) :: p(pcols,pver)   ! pressure [hPa]
      real(r8), intent(in) :: pf(pcols,pverp) ! pressure at interfaces [hPa]
      real(r8), intent(in) :: zm(pcols,pver)  ! midpoint height [m]
      real(r8), intent(in) :: zi(pcols,pverp) ! interface height [m]
      real(r8), intent(in) :: pblh(pcols)     ! PBL height [m]
      real(r8), intent(in) :: phis(pcols)     ! Surface geopotential [m]

   !  OUTPUT arguments
      real(r8), intent(out) :: cape(pcols)    ! convective aval. pot. energy [J/kg]
      real(r8), intent(out) :: cin(pcols)     ! covective inhibition [J/kg]

   !
   !--------------------------Local Variables------------------------------
   !

      real(r8) tv(pcols,pver)       ! virtual temperature [K]
      real(r8) tpv(pcols,pver)      ! virtual temperature of parcel [K]
      real(r8) buoy(pcols,pver)     ! Buoyancy for CAPE calculations [K]
      real(r8) neg_buoy(pcols,pver) ! "Negative" buoyancy for CIN calculations [K]

      ! height with respect to sea level [m]
      real(r8) :: z_sl(pcols,pver)

      real(r8) :: tp(pcols,pver)   ! parcel temperature [K]
      real(r8) :: qstp(pcols,pver) ! saturation mixing ratio of parcel [kg/kg]
      real(r8) :: tl(pcols)        ! parcel temperature at lcl [K]

      integer lcl(pcols)        ! level of lifting condensation level
      integer lfc(pcols)        ! level of free convection
      integer lel(pcols)        ! level of equilibrium level
      integer lon(pcols)        ! level of onset of deep convection
      integer mx(pcols)         ! level of max moist static energy
      integer pblt(pcols)       ! integer of PBL height

      ! Miscellanous 2d arrays needed for diagnostic
      real(r8) :: a1(pcols), a2(pcols), estp(pcols), pl(pcols)
      real(r8) :: plexp(pcols), hmax(pcols), hmn(pcols), y(pcols)
      integer :: knt(pcols)

      real(r8) :: e, tiedke_add, rgravit

      ! Indicees
      integer :: i, k, n

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

      ! Optional argument used in CAM/EAM to increase buoyancy
      !  according to Tiedke.  For diagnostics purposes, set to zero.
      tiedke_add = 0.0_r8

      rgravit = 1._r8/gravit

   !  Compute midpoint height above sea level, needed for MSE computation
      do k = 1,pver
         do i = 1, ncol
            z_sl(i,k) = zm(i,k) + phis(i)*rgravit
         end do
      end do

   !  Find indicee of PBL height, using mid point and interface heights
   !    from surface.
      do k = pver -1, 1, -1
        do i = 1, ncol
          if (abs(zm(i,k)-pblh(i)) < (zi(i,k)-zi(i,k+1))*0.5_r8) pblt(i) = k
        end do
      end do

   !  Initialize variables
      do i = 1,ncol
         cin(i) = 0._r8
         lon(i) = pver
         knt(i) = 0
         lel(i) = pver
         mx(i) = lon(i)
         cape(i) = 0._r8
         hmax(i) = 0._r8
      end do

      tp(:ncol,:) = t(:ncol,:)
      qstp(:ncol,:) = q(:ncol,:)

   ! Initialize virtual temperature and buoyancy
      tv(:ncol,:) = t(:ncol,:) *(1._r8+1.608_r8*q(:ncol,:))/ (1._r8+q(:ncol,:))
      tpv(:ncol,:) = tv(:ncol,:)
      buoy(:ncol,:) = 0._r8

   ! set "launching" level(mx) to be at maximum moist static energy.
   !   search for this level stops at planetary boundary layer top.
      do k = pver,1,-1
         do i = 1,ncol
            hmn(i) = cpair*t(i,k) + gravit*z_sl(i,k) + latvap*q(i,k)
            if (k >= pblt(i) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then
               hmax(i) = hmn(i)
               mx(i) = k
            end if
         end do
      end do

   !  Following computation to compute the Temperature and pressure of the
   !   LCL following Bolton 1980.
      do i = 1,ncol
         lcl(i) = mx(i)
         e = p(i,mx(i))*q(i,mx(i))/ (epsilo+q(i,mx(i)))
         tl(i) = 2840._r8/ (3.5_r8*log(t(i,mx(i)))-log(e)-4.805_r8) + 55._r8
         if (tl(i) < t(i,mx(i))) then
            plexp(i) = (1._r8/ (0.2854_r8* (1._r8-0.28_r8*q(i,mx(i)))))
            pl(i) = p(i,mx(i))* (tl(i)/t(i,mx(i)))**plexp(i)
         else
            tl(i) = t(i,mx(i))
            pl(i) = p(i,mx(i))
         end if
      end do

   ! Find the index of the LCL
      do k = pver,2,-1
         do i = 1,ncol
            if (k <= mx(i) .and. (p(i,k) > pl(i) .and. p(i,k-1) <= pl(i))) then
               lcl(i) = k - 1
            end if
         end do
      end do

   ! initialize parcel properties in sub-cloud layer below lcl.
      do k = pver,1,-1
         do i=1,ncol
             if (k > lcl(i)) then
               tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
               qstp(i,k) = q(i,mx(i))
               tp(i,k) = t(i,mx(i))* (p(i,k)/p(i,mx(i)))**(0.2854_r8* (1._r8-0.28_r8*q(i,mx(i))))
               tpv(i,k) = tp(i,k)*(1._r8+1.608_r8*q(i,mx(i)))/ (1._r8+q(i,mx(i)))
               buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
               neg_buoy(i,k) = tv(i,k) - tpv(i,k) + tiedke_add
            end if
         end do
      end do

   ! define parcel properties at lcl (i.e. level immediately above pl).
      do k = pver,1,-1
         do i=1,ncol
            if (k == lcl(i)) then
               tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
               qstp(i,k) = q(i,mx(i))
               tp(i,k) = tl(i)* (p(i,k)/pl(i))**(0.2854_r8* (1._r8-0.28_r8*qstp(i,k)))
               call qsat_hPa(tp(i,k), p(i,k), estp(i), qstp(i,k))
               a1(i) = cpair / latvap + qstp(i,k) * (1._r8+ qstp(i,k) / epsilo) * latvap * epsilo / &
                       (rair * tp(i,k) ** 2)
               a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/epsilo*qstp(i,k))* &
                       (1._r8+qstp(i,k)/epsilo)*epsilo**2*latvap*latvap/ &
                       (rair**2*tp(i,k)**4)-qstp(i,k)* &
                       (1._r8+qstp(i,k)/epsilo)*2._r8*epsilo*latvap/ &
                       (rair*tp(i,k)**3))
               a1(i) = 1._r8/a1(i)
               a2(i) = -a2(i)*a1(i)**3
               y(i) = q(i,mx(i)) - qstp(i,k)
               tp(i,k) = tp(i,k) + a1(i)*y(i) + a2(i)*y(i)**2
               call qsat_hPa(tp(i,k), p(i,k), estp(i), qstp(i,k))
               tpv(i,k) =  tp(i,k) * (1._r8+1.608_r8*qstp(i,k)) / (1._r8+q(i,mx(i)))
               buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
               neg_buoy(i,k) = tv(i,k) - tpv(i,k) + tiedke_add
            end if
         end do
      end do

   ! main buoyancy calculation.
      do k = pver - 1,1,-1
         do i=1,ncol
            if (k < lcl(i)) then
               tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
               qstp(i,k) = qstp(i,k+1)
               tp(i,k) = tp(i,k+1)* (p(i,k)/p(i,k+1))**(0.2854_r8* (1._r8-0.28_r8*qstp(i,k)))
               call qsat_hPa(tp(i,k), p(i,k), estp(i), qstp(i,k))
               a1(i) = cpair/latvap + qstp(i,k)* (1._r8+qstp(i,k)/epsilo)*latvap*epsilo/ (rair*tp(i,k)**2)
               a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/epsilo*qstp(i,k))* &
                       (1._r8+qstp(i,k)/epsilo)*epsilo**2*latvap*latvap/ &
                       (rair**2*tp(i,k)**4)-qstp(i,k)* &
                       (1._r8+qstp(i,k)/epsilo)*2._r8*epsilo*latvap/ &
                       (rair*tp(i,k)**3))
               a1(i) = 1._r8/a1(i)
               a2(i) = -a2(i)*a1(i)**3
               y(i) = qstp(i,k+1) - qstp(i,k)
               tp(i,k) = tp(i,k) + a1(i)*y(i) + a2(i)*y(i)**2
               call qsat_hPa(tp(i,k), p(i,k), estp(i), qstp(i,k))
               tpv(i,k) = tp(i,k) * (1._r8+1.608_r8*qstp(i,k))/(1._r8+q(i,mx(i)))
               buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
               neg_buoy(i,k) = tv(i,k) - tpv(i,k) + tiedke_add
            end if
         end do
      end do

   ! determine equilibrim level (EL)
      do k = 2,pver
         do i = 1,ncol
            ! must reside above the LCL
            if (k < lcl(i)) then
               ! should be level where buoyancy becomes negative
               if (buoy(i,k+1) > 0._r8 .and. buoy(i,k) <= 0._r8) then
                  lel(i) = k
               end if
            end if
         end do
      end do

   ! initialize the level of free convection (LFC) to be LCL
      do i = 1, ncol
        lfc(i) = lcl(i)
      end do

   ! determine the LFC
      do k = 2, pver
         do i = 1, ncol
           ! must reside above LCL and below EL
           if (k < lcl(i) .and. k > lel(i)) then
              ! should be the level where buoyancy becomes positive
              if (buoy(i,k+1) < 0._r8 .and. buoy(i,k) >= 0._r8) then
                 lfc(i) = k
              end if
           end if
         end do
      end do

   ! calculate convective available potential energy (CAPE), buoyancy
   !  integrated from LFC to the EL
      do k = 1,pver
         do i = 1,ncol
            if (k <= lfc(i) .and. k > lel(i)) then
               cape(i) = cape(i) + rair*buoy(i,k)*log(pf(i,k+1)/pf(i,k))
            end if
         end do
      end do

   ! Compute CIN based on information of levels computed above, which
   !  is the negative buoyancy integrated from surface to the LFC
      do k = 1, pver
        do i = 1, ncol
           if (k > lfc(i)) then
             cin(i) = cin(i) + rair*neg_buoy(i,k) * log(pf(i,k+1)/pf(i,k))
           endif
        end do
      end do

   ! put lower bound on cape and cin for diagnostic purposes.
      do i = 1,ncol
         cape(i) = max(cape(i), 0._r8)
         cin(i) = max(cin(i), 0._r8)
      end do
   !
      return
   end subroutine diag_CAPEandCIN

!#######################################################################

   ! Wrapper for qsat_water that does translation between Pa and hPa
   ! qsat_water uses Pa internally, so get it right, need to pass in Pa.
   ! Afterward, set es back to hPa.
   elemental subroutine qsat_hPa(t, p, es, qm)
     use wv_saturation, only: qsat_water

     ! Inputs
     real(r8), intent(in) :: t    ! Temperature (K)
     real(r8), intent(in) :: p    ! Pressure (hPa)
     ! Outputs
     real(r8), intent(out) :: es  ! Saturation vapor pressure (hPa)
     real(r8), intent(out) :: qm  ! Saturation mass mixing ratio
                                  ! (vapor mass over dry mass, kg/kg)

     call qsat_water(t, p*100._r8, es, qm)

     es = es*0.01_r8

   end subroutine qsat_hPa

! ########################################
! YQIN 08/18/23
subroutine calc_EIS(ncol, ncold, TS, rh_ref, SLP, T700, z700, LTS, EIS, EIS_fxRH)

integer , intent(in) :: ncol
integer , intent(in) :: ncold
real(r8), intent(in) :: TS(ncold)   ! surface temperature (K)
real(r8), intent(in) :: rh_ref(ncold) ! RH at the reference height (%)
real(r8), intent(in) :: SLP(ncold)  ! sea level pressure (Pa)
real(r8), intent(in) :: T700(ncold) ! temperature at 700 hPa (K)
real(r8), intent(in) :: z700(ncold) ! height at 700 hPa (m)

real(r8), intent(out) :: LTS(ncold) ! Lower tropspheric stability (LTS)
real(r8), intent(out) :: EIS(ncold) ! estimated inversion strength (EIS)
real(r8), intent(out) :: EIS_fxRH(ncold) ! estimated inversion strength (EIS)


real(r8) :: es
real(r8) :: e, Td, p_LCL, T_LCL, LCL
real(r8) :: e_fxRH, Td_fxRH, p_LCL_fxRH, T_LCL_fxRH, LCL_fxRH
real(r8) :: theta700, theta_slp
real(r8) :: T_bar, es_bar, qs_bar, gamma_m

integer :: i

do i=1,ncol
    !calculate the LCL height
    es=6.11_r8*exp(((2.5_r8*10._r8**6)/461._r8)*((1._r8/273.15_r8)-(1._r8/TS(i))))  !Clausius-Clapeyron Equation (Equation 4.23 from Curry and Webster)

    ! used varied RH 
    e=es*rh_ref(i)  ! use calculated RH at the reference height
    Td=((1/273.15_r8)-(461._r8/(2.5_r8*10._r8**6))*log(e/6.11_r8))**(-1) !form of Clausius-Clapeyron Equation
    p_LCL=SLP(i)*(((TS(i)-Td)/223.15_r8)+1._r8)**(-3.5_r8)                 !Equation (12) from Georgakakos and Bras (1984)
    T_LCL=TS(i)*(((TS(i)-Td)/223.15_r8)+1._r8)**(-1)                    !Equation (13) from Georgakakos and Bras (1984)
    LCL=((287._r8*(TS(i)+T_LCL)/2._r8)/9.8_r8)*log(SLP(i)/p_LCL)         !Hypsometric equation

    ! use fixed RH 
    e_fxRH=es*80._r8/100._r8                                            !assume RH is 80%, as in WB06
    Td_fxRH=((1/273.15_r8)-(461._r8/(2.5_r8*10._r8**6))*log(e_fxRH/6.11_r8))**(-1) !form of Clausius-Clapeyron Equation
    p_LCL_fxRH=SLP(i)*(((TS(i)-Td_fxRH)/223.15_r8)+1._r8)**(-3.5_r8)                 !Equation (12) from Georgakakos and Bras (1984)
    T_LCL_fxRH=TS(i)*(((TS(i)-Td_fxRH)/223.15_r8)+1._r8)**(-1)                    !Equation (13) from Georgakakos and Bras (1984)
    LCL_fxRH=((287._r8*(TS(i)+T_LCL_fxRH)/2._r8)/9.8_r8)*log(SLP(i)/p_LCL_fxRH)         !Hypsometric equation

    !calculate LTS
    theta700=T700(i)*(1000._r8/700._r8)**(287._r8/1004._r8)
    theta_slp=TS(i)*(1000._r8/(SLP(i)*0.01_r8))**(287._r8/1004._r8)
    LTS(i)=theta700-theta_slp

    !calculate the moist adaiabtic lapse rate at 850 hPa
    T_bar=(TS(i)+T700(i))/2._r8 !approximation of the temperature at 850 hPa, following WB06
    es_bar=6.11_r8*exp(((2.5_r8*10._r8**6)/461._r8)*((1._r8/273.15_r8)-(1._r8/T_bar))) !Clausius-Clapeyron Equation (Equation 4.23 from Curry and Webster)
    qs_bar=0.622_r8*es_bar/(850._r8-es_bar) !Equation 4.37 of Curry and Webster
    gamma_m=(9.8_r8/1004._r8)*(1._r8-((1._r8+(((2.5_r8*10._r8**6)*qs_bar)/(287._r8*T_bar)))/(1._r8+((((2.5_r8*10._r8**6)**2)*qs_bar)/(1004._r8*461._r8*(T_bar**2)))))) !Equation (5) of WB06

    EIS_fxRH(i)=LTS(i)-gamma_m*(z700(i)-LCL_fxRH) !Equation (4) of WB06
    EIS(i)=LTS(i)-gamma_m*(z700(i)-LCL) !Equation (4) of WB06

end do ! i

end subroutine calc_EIS
! ########################################
   
end module cam_diagnostics
