#include <misc.h>
#include <params.h>

subroutine tphysac (ztodt,   pblh,    qpert,   tpert,  shf,  & 1,48
                    taux,    tauy,    cflx,    sgh,    lhf,  &
                    landfrac,snowh,   tref,    precc,  precl,&
                    precsc, precsl,   state,   tend,    pbuf,&
                    ocnfrac, fsds, icefrac, fv, ram1 )
! Purpose: 
! Tendency physics after coupling to land, sea, and ice models.
! Computes the following:
!   o Radon surface flux and decay (optional)
!   o Vertical diffusion and planetary boundary layer
!   o Multiple gravity wave drag
! Method: 
! <Describe the algorithm(s) used in the routine.> 
! <Also include any applicable external references.> 
! Author: CCM1, CMS Contact: J. Truesdale
   use shr_kind_mod, only: r8 => shr_kind_r8
   use ppgrid,             only: pcols, pver
   use chemistry,          only: trace_gas, chem_timestep_tend
   use gw_drag,            only: gw_intr
   use vertical_diffusion, only: vd_intr
   use physics_types,      only: physics_state, physics_tend, physics_ptend, physics_update,    &
                                 physics_ptend_init, physics_dme_adjust
   use phys_buffer,        only: pbuf_size_max, pbuf_fld, pbuf_old_tim_idx, pbuf_get_fld_idx
   use constituents,       only: ppcnst, qmin
   use tracers,       only: tracers_timestep_tend, set_state_pdry, set_dry_to_wet
   use physconst,          only: zvir, gravit, rhoh2o, latvap,latice, cpair, rair
   use aerosol_intr,       only: aerosol_emis_intr, aerosol_drydep_intr, aerosol_srcsnk_intr
   use check_energy,    only: check_energy_chng
   use check_energy,       only: check_tracers_data, check_tracers_init, check_tracers_chng
   use time_manager,    only: get_nstep
   use abortutils,      only: endrun
   use dycore,          only: dycore_is

   implicit none

#include <comctl.h>
! Arguments
   real(r8), intent(in) :: ztodt                  ! Two times model timestep (2 delta-t)
   real(r8), intent(in) :: landfrac(pcols)        ! Land fraction
   real(r8), intent(in) :: icefrac(pcols)         ! ice fraction
   real(r8), intent(in) :: ocnfrac(pcols)         ! Land fraction
   real(r8), intent(in) :: snowh(pcols)           ! snow depth (liquid water equivalent)
   real(r8), intent(in) :: fv(pcols)              ! for dry deposition velocities for dust from land model
   real(r8), intent(in) :: ram1(pcols)            ! for dry deposition velocities for dust from land model
   real(r8), intent(in) :: tref(pcols)            ! 2m air temperature
   real(r8), intent(in) :: precc(pcols)           ! convective precipitation
   real(r8), intent(in) :: precl(pcols)           ! large-scale precipitation
   real(r8), intent(in) :: fsds(pcols)            ! down solar flux
   real(r8), intent(in) :: precsc(pcols)           ! convective snow
   real(r8), intent(in) :: precsl(pcols)           ! large-scale snow
   real(r8), intent(out) :: pblh(pcols)           ! Planetary boundary layer height
   real(r8), intent(out) :: qpert(pcols,ppcnst)   ! Moisture/constit. perturbation (PBL)
   real(r8), intent(out) :: tpert(pcols)          ! Temperature perturbation (PBL)
   real(r8), intent(inout) :: shf(pcols)          ! Sensible heat flux (w/m^2)
   real(r8), intent(in) :: taux(pcols)            ! X surface stress (zonal)
   real(r8), intent(in) :: tauy(pcols)            ! Y surface stress (meridional)
   real(r8), intent(inout) :: cflx(pcols,ppcnst)  ! Surface constituent flux (kg/m^2/s)
   real(r8), intent(in) :: sgh(pcols)             ! Std. deviation of orography for gwd
   real(r8), intent(inout) :: lhf(pcols)          ! Latent heat flux (w/m^2)

   type(physics_state), intent(inout) :: state
   type(physics_tend ), intent(inout) :: tend
   type(pbuf_fld),      intent(inout), dimension(pbuf_size_max) :: pbuf

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

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

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

   integer :: lchnk                                ! chunk identifier
   integer :: ncol                                 ! number of atmospheric columns
   integer i,k,m                 ! Longitude, level indices
   integer :: yr, mon, day, tod       ! components of a date

   logical :: labort                            ! abort flag

   real(r8) tvm(pcols,pver)           ! virtual temperature
   real(r8) prect(pcols)              ! total precipitation
   real(r8) surfric(pcols)              ! surface friction velocity
   real(r8) obklen(pcols)             ! Obukhov length
   real(r8) :: fh2o(pcols)            ! h2o flux to balance source from methane chemistry

! physics buffer fields for total energy and mass adjustment
   integer itim, ifld
   real(r8), pointer, dimension(:  ) :: teout
   real(r8), pointer, dimension(:,:) :: qini
   real(r8), pointer, dimension(:,:) :: tini
   real(r8), pointer, dimension(:,:) :: cld

   lchnk = state%lchnk
   ncol  = state%ncol

! Associate pointers with physics buffer fields
   itim = pbuf_old_tim_idx()
   ifld = pbuf_get_fld_idx('TEOUT')
   teout => pbuf(ifld)%fld_ptr(1,1:pcols,1,lchnk,itim)
   ifld = pbuf_get_fld_idx('QINI')
   qini  => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
   ifld = pbuf_get_fld_idx('TINI')
   tini  => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
   ifld = pbuf_get_fld_idx('CLD')
   cld => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,itim)
! accumulate fluxes into net flux array
! jrm Include latent heat of fusion for snow
   do i=1,ncol
      tend%flx_net(i) = tend%flx_net(i) + shf(i) + (precc(i) + precl(i))*latvap*rhoh2o &
           + (precsc(i) + precsl(i))*latice*rhoh2o
   end do

! Initialize parameterization tendency structure

   call physics_ptend_init(ptend)

! emission of aerosols at surface
   call aerosol_emis_intr (state, ptend, cflx, ztodt)
   call physics_update (state, tend, ptend, ztodt)

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

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

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

! Source/sink terms for advected tracers.
! Test tracers

      call tracers_timestep_tend(state, ptend, cflx, landfrac, ztodt)      
      call physics_update (state, tend, ptend, ztodt)
      call check_tracers_chng(state, tracerint, "tracers_timestep_tend", nstep, ztodt, cflx)

! Advected greenhouse trace gases:

   if (trace_gas) then
      call chem_timestep_tend(state, ptend, cflx, ztodt, fh2o)
      call physics_update (state, tend, ptend, ztodt)
! Check energy integrals
      call check_energy_chng(state, tend, "chem", nstep, ztodt, fh2o, zero, zero, zero)
! Check tracer integrals
      call check_tracers_chng(state, tracerint, "chem_timestep_tend", nstep, ztodt, cflx)
   end if

! Vertical diffusion/pbl calculation
! Call vertical diffusion code (pbl, free atmosphere and molecular)

   call vd_intr (ztodt    ,state    ,taux     ,tauy     , shf    ,&
                 cflx     ,pblh     ,tpert    ,qpert    , surfric  ,&
                 obklen   ,ptend    ,cld      ,ocnfrac  , landfrac, sgh )

   call physics_update (state, tend, ptend, ztodt)
! Check energy integrals
   call check_energy_chng(state, tend, "vdiff", nstep, ztodt, cflx(:,1), zero, zero, shf)
   call check_tracers_chng(state, tracerint, "vdiff", nstep, ztodt, cflx)

!  aerosol dry deposition processes
   call aerosol_drydep_intr (state, ptend, cflx(:,1), ztodt, &
       fsds, obklen, tref, surfric, prect, snowh, pblh, shf, landfrac, &
       icefrac, ocnfrac, fv, ram1)

   call physics_update (state, tend, ptend, ztodt)

! Gravity wave drag

   call gw_intr (state   ,sgh     ,pblh    ,ztodt   , ptend , landfrac)
   call physics_update (state, tend, ptend, ztodt)
! Check energy integrals
   call check_energy_chng(state, tend, "gwdrag", nstep, ztodt, zero, zero, zero, zero)

! conversion of aerosols from one category to another by non-wet processes
   call aerosol_srcsnk_intr (state, ptend, ztodt, ocnfrac)
   call physics_update (state, tend, ptend, ztodt)

!-------------- Energy budget checks vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
   teout(:ncol) = state%te_cur(:ncol)

!*** BAB's FV heating kludge *** apply the heating as temperature tendency.
!*** BAB's FV heating kludge *** modify the temperature in the state structure
   state%t(:ncol,:pver) = tini(:ncol,:pver) + ztodt*tend%dtdt(:ncol,:pver)

! FV: convert dry-type mixing ratios to moist here because physics_dme_adjust
!     assumes moist. This is done in p_d_coupling for other dynamics. Bundy, Feb 2004.
   if ( dycore_is('LR') ) call set_dry_to_wet(state)    ! Physics had dry, dynamics wants moist

! Scale dry mass and energy (does nothing if dycore is EUL or SLD)
   call physics_dme_adjust(state, tend, qini, ztodt)
!!!   call check_energy_chng(state, tend, "drymass", nstep, ztodt, zero, zero, zero, zero)
!-------------- Energy budget checks ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

   if (aqua_planet) then
      labort = .false.
      do i=1,ncol
         if (ocnfrac(i) /= 1.) labort = .true.
      end do
      if (labort) then
         call endrun ('TPHYSAC error:  grid contains non-ocean point')

end subroutine tphysac