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


subroutine tphysbc (ztodt,   pblht,   tpert,   ts,      sst,     & 1,120
                    qpert,   precl,   precc,   precsl,  precsc,  &
                    asdir,   asdif,   aldir,   aldif,   snowh,   &
                    qrs,     qrl,     flwds,   fsns,    fsnt,    &
                    flns,    flnt,    lwup,    srfrad,  sols,    &
                    soll,    solsd,   solld,   state,   tend,    &
                    pbuf,    prcsnw,  fsds ,   landm,   landfrac,&
		    ocnfrac, icefrac  )
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Tendency physics BEFORE coupling to land, sea, and ice models.
! 
! Method: 
! Call physics subroutines and compute the following:
!     o cloud calculations (cloud fraction, emissivity, etc.)
!     o radiation calculations
! Pass surface fields for separate surface flux calculations
! Dump appropriate fields to history file.
! 
! Author: CCM1, CMS Contact: J. Truesdale
! 
!-----------------------------------------------------------------------

   use shr_kind_mod,    only: r8 => shr_kind_r8
   use ppgrid
   use phys_grid,       only: get_rlat_all_p, get_rlon_all_p, get_lon_all_p, get_lat_all_p
   use phys_buffer,     only: pbuf_size_max, pbuf_fld, pbuf_old_tim_idx, pbuf_get_fld_idx
   use cldcond,         only: cldcond_tend, cldcond_zmconv_detrain, cldcond_sediment
   use param_cldoptics, only: param_cldoptics_calc
   use physics_types,   only: physics_state, physics_tend, physics_ptend, physics_update, physics_ptend_init
   use diagnostics,     only: diag_dynvar
   use history,         only: outfld
   use physconst,       only: gravit, latvap, cpair, tmelt, cappa, zvir, rair, rga
   use radheat,         only: radheat_net
   use constituents,    only: pcnst, pnats, ppcnst, qmin
   use constituents,    only: dcconnam, cnst_get_ind
   use zm_conv,         only: zm_conv_evap, zm_convr
   use time_manager,    only: is_first_step, get_nstep, get_curr_calday
   use moistconvection, only: cmfmca
   use check_energy,    only: check_energy_chng, check_energy_fix
   use check_energy,    only: check_tracers_data, check_tracers_init, check_tracers_chng
   use dycore,          only: dycore_is
   use cloudsimulator,  only: doisccp, cloudsimulator_run
   use aerosol_intr, only: aerosol_wet_intr

   implicit none

#include <comctl.h>
!
! Arguments
!
   real(r8), intent(in) :: ztodt                          ! 2 delta t (model time increment)
   real(r8), intent(in) :: ts(pcols)                      ! surface temperature
   real(r8), intent(in) :: sst(pcols)                     ! sea surface temperature
   real(r8), intent(inout) :: pblht(pcols)                ! Planetary boundary layer height
   real(r8), intent(inout) :: tpert(pcols)                ! Thermal temperature excess
   real(r8), intent(inout) :: qpert(pcols,ppcnst)         ! Thermal humidity & constituent excess
   real(r8), intent(in) :: asdir(pcols)                  ! Albedo: shortwave, direct
   real(r8), intent(in) :: asdif(pcols)                  ! Albedo: shortwave, diffuse
   real(r8), intent(in) :: aldir(pcols)                  ! Albedo: longwave, direct
   real(r8), intent(in) :: aldif(pcols)                  ! Albedo: longwave, diffuse
   real(r8), intent(in) :: snowh(pcols)                  ! Snow depth (liquid water equivalent)
   real(r8), intent(inout) :: qrs(pcols,pver)            ! Shortwave heating rate
   real(r8), intent(inout) :: qrl(pcols,pver)            ! Longwave  heating rate
   real(r8), intent(inout) :: flwds(pcols)               ! Surface longwave down flux
   real(r8), intent(inout) :: fsns(pcols)                   ! Surface solar absorbed flux
   real(r8), intent(inout) :: fsnt(pcols)                   ! Net column abs solar flux at model top
   real(r8), intent(inout) :: flns(pcols)                   ! Srf longwave cooling (up-down) flux
   real(r8), intent(inout) :: flnt(pcols)                   ! Net outgoing lw flux at model top
   real(r8), intent(in) :: lwup(pcols)                    ! Surface longwave up flux
   real(r8), intent(out) :: srfrad(pcols)                 ! Net surface radiative flux (watts/m**2)
   real(r8), intent(inout) :: sols(pcols)                   ! Direct beam solar rad. onto srf (sw)
   real(r8), intent(inout) :: soll(pcols)                   ! Direct beam solar rad. onto srf (lw)
   real(r8), intent(inout) :: solsd(pcols)                  ! Diffuse solar radiation onto srf (sw)
   real(r8), intent(inout) :: solld(pcols)                  ! Diffuse solar radiation onto srf (lw)
   real(r8), intent(out) :: precl(pcols)                  ! Large-scale precipitation rate
   real(r8), intent(out) :: precc(pcols)                  ! Convective-scale preciptn rate
   real(r8), intent(out) :: precsl(pcols)                 ! L.S. snowfall rate
   real(r8), intent(out) :: precsc(pcols)                 ! C.S. snowfall rate
   real(r8), intent(out) :: prcsnw(pcols)                 ! snowfall rate (precsl + precsc)
   real(r8), intent(out) :: fsds(pcols)                   ! Surface solar down flux
   real(r8), intent(in) :: landm(pcols)                   ! land fraction ramp
   real(r8), intent(in) :: landfrac(pcols)                ! land fraction
   real(r8), intent(in) :: ocnfrac(pcols)                ! land fraction
   real(r8), intent(in) :: icefrac(pcols)                ! land fraction

   type(physics_state), intent(inout) :: state
   type(physics_tend ), intent(inout) :: tend
   type(pbuf_fld),      intent(inout), dimension(pbuf_size_max) :: pbuf
!
!---------------------------Local workspace-----------------------------
!
   real(r8) :: rhdfda(pcols,pver)            ! dRh/dcloud, old 
   real(r8) :: rhu00 (pcols,pver)            ! Rh threshold for cloud, old

   type(physics_ptend)   :: ptend                  ! indivdual parameterization tendencies

   integer :: nstep                          ! current timestep number
   integer      lat(pcols)                   ! current latitudes(indices)
   integer      lon(pcols)                   ! current longtitudes(indices)

   real(r8) :: calday                        ! current calendar day
   real(r8) :: clat(pcols)                   ! current latitudes(radians)
   real(r8) :: clon(pcols)                   ! current longitudes(radians)

   real(r8) :: zdu(pcols,pver)               ! detraining mass flux from deep convection
   real(r8) :: ftem(pcols,pver)              ! Temporary workspace for outfld variables
   real(r8) :: zmrprd(pcols,pver)            ! rain production in ZM convection
   real(r8) :: cmfmc(pcols,pverp)            ! Convective mass flux--m sub c
   real(r8) :: cmfsl(pcols,pver)             ! Moist convection lw stat energy flux
   real(r8) :: cmflq(pcols,pver)             ! Moist convection total water flux
   real(r8) :: dtcond(pcols,pver)            ! dT/dt due to moist processes
   real(r8) :: dqcond(pcols,pver,ppcnst)     ! dq/dt due to moist processes

   real(r8) cldst(pcols,pver)
   real(r8) cltot(pcols)                      ! Diagnostic total cloud cover
   real(r8) cllow(pcols)                      !       "     low  cloud cover
   real(r8) clmed(pcols)                      !       "     mid  cloud cover
   real(r8) clhgh(pcols)                      !       "     hgh  cloud cover
   real(r8) cmfcme(pcols,pver)                ! cmf condensation - evaporation
   real(r8) cmfdqr2(pcols,pver)               ! dq/dt due to moist convective rainout
   real(r8) cmfmc2(pcols,pverp)               ! Moist convection cloud mass flux
   real(r8) cmfsl2(pcols,pver)                ! Moist convection lw stat energy flux
   real(r8) cmflq2(pcols,pver)                ! Moist convection total water flux
   real(r8) cnt(pcols)                        ! Top level of convective activity
   real(r8) cnb(pcols)                        ! Lowest level of convective activity
   real(r8) cnt2(pcols)                       ! Top level of convective activity
   real(r8) cnb2(pcols)                       ! Bottom level of convective activity
   real(r8) concld(pcols,pver)             
   real(r8) coszrs(pcols)                     ! Cosine solar zenith angle
   real(r8) dlf(pcols,pver)                   ! Detraining cld H20 from convection
   real(r8) pflx(pcols,pverp)                 ! Conv rain flux thru out btm of lev
   real(r8) prect(pcols)                      ! total (conv+large scale) precip rate
   real(r8) dlf2(pcols,pver)                   ! dq/dt due to rainout terms
   real(r8) qpert2(pcols,ppcnst)              ! Perturbation q
   real(r8) rtdt                              ! 1./ztodt
   real(r8) tpert2(pcols)                     ! Perturbation T
   real(r8) pmxrgn(pcols,pverp)               ! Maximum values of pressure for each
!                                             !    maximally overlapped region.
!                                             !    0->pmxrgn(i,1) is range of pressure for
!                                             !    1st region,pmxrgn(i,1)->pmxrgn(i,2) for
!                                             !    2nd region, etc
   integer lchnk                              ! chunk identifier
   integer ncol                               ! number of atmospheric columns

   integer nmxrgn(pcols)                      ! Number of maximally overlapped regions
   integer  i,k,m                             ! Longitude, level, constituent indices
   integer :: ixcldice, ixcldliq              ! constituent indices for cloud liquid and ice water.
                                           
!  real(r8) engt                              ! Thermal   energy integral
!  real(r8) engk                              ! Kinetic   energy integral
!  real(r8) engp                              ! Potential energy integral
   real(r8) rel(pcols,pver)                   ! Liquid cloud particle effective radius
   real(r8) rei(pcols,pver)                   ! Ice effective drop size (microns)
   real(r8) emis(pcols,pver)                  ! Cloud longwave emissivity
   real(r8) clc(pcols)                        ! Total convective cloud (cloud scheme)
   real(r8) :: cicewp(pcols,pver)             ! in-cloud cloud ice water path
   real(r8) :: cliqwp(pcols,pver)             ! in-cloud cloud liquid water path
!
   real(r8) dellow(pcols)                     ! delta p for bottom three levels of model
   real(r8) tavg(pcols)                       ! mass weighted average temperature for 

! physics buffer fields to compute tendencies for cloud condensation package
   integer itim, ifld
   real(r8), pointer, dimension(:,:) :: qcwat, tcwat, lcwat, cld

! physics buffer fields for total energy and mass adjustment
   real(r8), pointer, dimension(:  ) :: teout
   real(r8), pointer, dimension(:,:) :: qini
   real(r8), pointer, dimension(:,:) :: tini
!                                          
! Used for OUTFLD only                     
!                                          
   real(r8) icwmr1(pcols,pver)                ! in cloud water mixing ration for zhang scheme
   real(r8) icwmr2(pcols,pver)                ! in cloud water mixing ration for hack scheme
   real(r8) fracis(pcols,pver,ppcnst)         ! fraction of transported species that are insoluble
   real(r8) timestep(pcols)
!
!     Variables for doing deep convective transport outside of zm_convr
!
   real(r8) mu2(pcols,pver)
   real(r8) eu2(pcols,pver)
   real(r8) du2(pcols,pver)
   real(r8) md2(pcols,pver)
   real(r8) ed2(pcols,pver)
   real(r8) dp(pcols,pver)
   real(r8) dpdry(pcols,pver)
   real(r8) dsubcld(pcols)
   real(r8) conicw(pcols,pver)
   real(r8) cmfdqrt(pcols,pver)               ! dq/dt due to moist convective rainout

! stratiform precipitation variables
   real(r8) :: prec_pcw(pcols)                ! total precip from prognostic cloud scheme
   real(r8) :: snow_pcw(pcols)                ! snow from prognostic cloud scheme
   real(r8) :: prec_sed(pcols)                ! total precip from cloud sedimentation
   real(r8) :: snow_sed(pcols)                ! snow from cloud ice sedimentation

! convective precipitation variables
   real(r8) :: prec_zmc(pcols)                ! total precipitation from ZM convection
   real(r8) :: snow_zmc(pcols)                ! snow from ZM convection
   real(r8) :: prec_cmf(pcols)                ! total precipitation from Hack convection
   real(r8) :: snow_cmf(pcols)                ! snow from Hack convection

! energy checking variables
   real(r8) :: zero(pcols)                    ! array of zeros
   real(r8) :: rliq(pcols)                    ! vertical integral of liquid not yet in q(ixcldliq)
   real(r8) :: rliq2(pcols)                   ! vertical integral of liquid from shallow scheme
   real(r8) :: flx_cnd(pcols)
   real(r8) :: flx_heat(pcols)
   logical  :: conserve_energy = .true.       ! flag to carry (QRS,QRL)*dp across time steps
   type(check_tracers_data):: tracerint             ! energy integrals and cummulative boundary fluxes

   integer jt(pcols)
   integer maxg(pcols)
   integer ideep(pcols)
   integer lengath
   real(r8) cldc(pcols,pver)
   real(r8) nevapr(pcols,pver)
   real(r8) qme(pcols,pver)
   real(r8) prain(pcols,pver)
   real(r8) cflx(pcols,ppcnst)
!-----------------------------------------------------------------------
   zero = 0.
!
   lchnk = state%lchnk
   ncol  = state%ncol

   rtdt = 1./ztodt

   nstep = get_nstep()
   calday = get_curr_calday()

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

!  dry surface pressure
   call outfld ('PSDRY',  state%psdry, pcols, lchnk)

! Associate pointers with physics buffer fields
   itim = pbuf_old_tim_idx()
   ifld = pbuf_get_fld_idx('QCWAT')
   qcwat => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,itim)
   ifld = pbuf_get_fld_idx('TCWAT')
   tcwat => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,itim)
   ifld = pbuf_get_fld_idx('LCWAT')
   lcwat => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,itim)
   ifld = pbuf_get_fld_idx('CLD')
   cld => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,itim)

   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)

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

   call physics_ptend_init (ptend) ! Initialize parameterization tendency structure
!
! Make sure that input tracers are all positive (probably unnecessary)
!
   call qneg3('TPHYSBCb',lchnk  ,ncol    ,pcols   ,pver    , &
              ppcnst,qmin  ,state%q )
!
! Setup q and t accumulation fields
!
   dqcond(:ncol,:,:) = state%q(:ncol,:,:)
   dtcond(:ncol,:)   = state%s(:ncol,:)

   fracis (:ncol,:,1:ppcnst) = 1.

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

!===================================================
! Global mean total energy fixer
!===================================================
   !*** BAB's FV heating kludge *** save the initial temperature
   tini(:ncol,:pver) = state%t(:ncol,:pver)
   if (dycore_is('LR')) then
      call check_energy_fix(state, ptend, nstep, flx_heat)
      call physics_update(state, tend, ptend, ztodt)
      call check_energy_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat)
   end if
   qini(:ncol,:pver) = state%q(:ncol,:pver,1)

   call outfld('TEOUT', teout       , pcols, lchnk   )
   call outfld('TEINP', state%te_ini, pcols, lchnk   )
   call outfld('TEFIX', state%te_cur, pcols, lchnk   )
!
!===================================================
! Dry adjustment
!===================================================

! Copy state info for input to dadadj
! This is a kludge, so that dadadj does not have to be correctly reformulated in dry static energy

   ptend%s(:ncol,:pver)   = state%t(:ncol,:pver)
   ptend%q(:ncol,:pver,1) = state%q(:ncol,:pver,1)

   call dadadj (lchnk, ncol, state%pmid,  state%pint,  state%pdel,  &
                ptend%s, ptend%q(1,1,1))
   ptend%name  = 'dadadj'
   ptend%ls    = .TRUE.
   ptend%lq(1) = .TRUE.
   ptend%s(:ncol,:)   = (ptend%s(:ncol,:)   - state%t(:ncol,:)  )/ztodt * cpair
   ptend%q(:ncol,:,1) = (ptend%q(:ncol,:,1) - state%q(:ncol,:,1))/ztodt
   call physics_update (state, tend, ptend, ztodt)
!
!===================================================
! Moist convection
!===================================================
!
! Since the PBL doesn't pass constituent perturbations, they
! are zeroed here for input to the moist convection routine
!
   qpert(:ncol,2:ppcnst) = 0.0
!
! Begin with Zhang-McFarlane (1996) convection parameterization
!
   call t_startf ('zm_convr')
   call zm_convr( lchnk,    ncol, &
                  state%t,   state%q,    prec_zmc,   cnt,     cnb,      &
                  pblht,   state%zm, state%phis,    state%zi,   ptend%q(:,:,1),     &
                  ptend%s, state%pmid,   state%pint,  state%pdel,       &
                  .5*ztodt,cmfmc,    cmfcme,             &
                  tpert,   dlf,      pflx,    zdu,     zmrprd,   &
                  mu2,      md2,     du2,     eu2,     ed2,      &
                  dp,       dsubcld, jt,      maxg,    ideep,    &
                  lengath, icwmr1,   rliq    )
   ptend%name  = 'zm_convr'
   ptend%ls    = .TRUE.
   ptend%lq(1) = .TRUE.
   cmfsl (:ncol,:) = 0. ! This is not returned from zm, hence it is zeroed.
   cmflq (:ncol,:) = 0. ! This is not returned from zm, hence it is zeroed.

   ftem(:ncol,:pver) = ptend%s(:ncol,:pver)/cpair
   call outfld('ZMDT    ',ftem           ,pcols   ,lchnk   )
   call outfld('ZMDQ    ',ptend%q(1,1,1) ,pcols   ,lchnk   )
   call t_stopf('zm_convr')

   call physics_update(state, tend, ptend, ztodt)
!
! Determine the phase of the precipitation produced and add latent heat of fusion
! Evaporate some of the precip directly into the environment (Sundqvist)
   call zm_conv_evap(state, ptend, zmrprd, cld, ztodt, prec_zmc, snow_zmc, .false.)
   call physics_update(state, tend, ptend, ztodt)
! Check energy integrals, including "reserved liquid"
   flx_cnd(:ncol) = prec_zmc(:ncol) + rliq(:ncol)
   call check_energy_chng(state, tend, "zm_evap", nstep, ztodt, zero, flx_cnd, snow_zmc, zero)

! Transport cloud water and ice only
!
   call cnst_get_ind('CLDLIQ', ixcldliq)
   call cnst_get_ind('CLDICE', ixcldice)
   ptend%name = 'convtran1'
   ptend%lq(ixcldice) = .true.
   ptend%lq(ixcldliq) = .true.
   call convtran (lchnk,                                        &
                  ptend%lq,state%q, ppcnst,  mu2,     md2,   &
                  du2,     eu2,     ed2,     dp,      dsubcld,  &
                  jt,      maxg,    ideep,   1,       lengath,  &
                  nstep,   fracis,  ptend%q, dpdry  )
   call physics_update (state, tend, ptend, ztodt)
!
! Convert mass flux from reported mb/s to kg/m^2/s
!
   cmfmc(:ncol,:pver) = cmfmc(:ncol,:pver) * 100./gravit
!
! Call Hack (1994) convection scheme to deal with shallow/mid-level convection
!
   call t_startf('cmfmca')
   tpert2(:ncol  ) =0.
   qpert2(:ncol,:) = qpert(:ncol,:)  ! BAB Why is this not zero, if tpert2=0???
   call cmfmca (lchnk,   ncol, &
                nstep,   ztodt,   state%pmid,  state%pdel,   &
                state%rpdel,   state%zm,      tpert2,  qpert2,  state%phis,     &
                pblht,   state%t,   state%q,   ptend%s,   ptend%q,      &
                cmfmc2,  cmfdqr2, cmfsl2,  cmflq2,  prec_cmf,   &
                dlf2,     cnt2,    cnb2,    icwmr2   , rliq2, & 
                state%pmiddry, state%pdeldry, state%rpdeldry)
   ptend%name  = 'cmfmca'
   ptend%ls    = .TRUE.
   ptend%lq(:) = .TRUE.

! Add shallow cloud water detrainment to cloud water detrained from ZM
   dlf(:ncol,:pver) = dlf(:ncol,:pver) + dlf2(:ncol,:pver)
   rliq(:ncol) = rliq(:ncol) + rliq2(:ncol)
   
   ftem(:ncol,:pver) = ptend%s(:ncol,:pver)/cpair
   call outfld('CMFDT   ',ftem          ,pcols   ,lchnk   )
   call outfld('CMFDQ   ',ptend%q(1,1,1),pcols   ,lchnk   )
   call t_stopf('cmfmca')
   call physics_update (state, tend, ptend, ztodt)
!
! Determine the phase of the precipitation produced and add latent heat of fusion
   call zm_conv_evap(state, ptend, cmfdqr2, cld, ztodt, prec_cmf, snow_cmf, .true.)
   call physics_update(state, tend, ptend, ztodt)
   flx_cnd(:ncol) = prec_cmf(:ncol) + rliq2(:ncol)
   call check_energy_chng(state, tend, "hk_evap", nstep, ztodt, zero, flx_cnd, snow_cmf, zero)
!
! Merge shallow/mid-level output with prior results from Zhang-McFarlane
!
   do i=1,ncol
      if (cnt2(i) < cnt(i)) cnt(i) = cnt2(i)
      if (cnb2(i) > cnb(i)) cnb(i) = cnb2(i)
   end do
!
   cmfmc (:ncol,:pver) = cmfmc (:ncol,:pver) + cmfmc2 (:ncol,:pver)
   cmfsl (:ncol,:pver) = cmfsl (:ncol,:pver) + cmfsl2 (:ncol,:pver)
   cmflq (:ncol,:pver) = cmflq (:ncol,:pver) + cmflq2 (:ncol,:pver)
   call outfld('CMFMC' , cmfmc  , pcols, lchnk)
!  output new partition of cloud condensate variables, as well as precipitation 
   call outfld('QC      ',dlf2           ,pcols   ,lchnk   )
   call outfld('PRECSH  ',prec_cmf      ,pcols   ,lchnk       )
   call outfld('CMFDQR', cmfdqr2, pcols, lchnk)
   call outfld('CMFSL' , cmfsl  , pcols, lchnk)
   call outfld('CMFLQ' , cmflq  , pcols, lchnk)
   call outfld('DQP'   , dlf2    , pcols, lchnk)

! Allow the cloud liquid drops and ice particles to sediment
! Occurs before adding convectively detrained cloud water, because the phase of the
! of the detrained water is unknown.
   call t_startf('cldwat_sediment')
   call cldcond_sediment(state, ptend, ztodt,cld, icefrac, landfrac, ocnfrac, prec_sed, &
                         snow_sed, landm, snowh)
   call physics_update(state, tend, ptend, ztodt)
   call t_stopf('cldwat_sediment')

! check energy integrals
   call check_energy_chng(state, tend, "cldwat_sediment", nstep, ztodt, zero, prec_sed, snow_sed, zero)

! Put the detraining cloud water from convection into the cloud and environment. 

   call cldcond_zmconv_detrain(dlf, cld, state, ptend)
   call physics_update(state, tend, ptend, ztodt)

! check energy integrals, reserved liquid has now been used
   flx_cnd(:ncol) = -rliq(:ncol)
   call check_energy_chng(state, tend, "cldwat_detrain", nstep, ztodt, zero, flx_cnd, zero, zero)
!
! cloud fraction after transport and convection,
! derive the relationship between rh and cld from 
! the employed cloud scheme
!
   call cldnrh(lchnk,   ncol,                                &
               state%pmid,    state%t,   state%q(1,1,1),   state%omega, &
               cnt,     cnb,     cld,    clc,     state%pdel,   &
               cmfmc,   cmfmc2, landfrac,snowh,   concld,  cldst,    &
               ts,      sst, state%pint(1,pverp),       zdu,  ocnfrac, &
               rhdfda,   rhu00 , state%phis)
   call outfld('CONCLD  ',concld, pcols,lchnk)
   call outfld('CLDST   ',cldst,  pcols,lchnk)
   call outfld('CNVCLD  ',clc,    pcols,lchnk)

! cloud water and ice parameterizations
   call t_startf('cldwat_tend')
   call cldcond_tend(state, ptend, ztodt, &
       tcwat, qcwat, lcwat, prec_pcw, snow_pcw, icefrac, rhdfda, rhu00, cld, nevapr, prain, qme, snowh)
   call physics_update (state, tend, ptend, ztodt)
   call t_stopf('cldwat_tend')

! check energy integrals
   call check_energy_chng(state, tend, "cldwat_tend", nstep, ztodt, zero, prec_pcw, snow_pcw, zero)

! Save off q and t after cloud water
   do k=1,pver
      qcwat(:ncol,k) = state%q(:ncol,k,1)
      tcwat(:ncol,k) = state%t(:ncol,k)
      lcwat(:ncol,k) = state%q(:ncol,k,ixcldice) + state%q(:ncol,k,ixcldliq)
   end do
!
!  aerosol wet chemistry determines scavenging fractions, and transformations
   call get_lat_all_p(lchnk, ncol, lat)
   call get_lon_all_p(lchnk, ncol, lon)
   call get_rlat_all_p(lchnk, ncol, clat)
   conicw(:ncol,:) = icwmr1(:ncol,:) + icwmr2(:ncol,:)
   cmfdqrt(:ncol,:) = zmrprd(:ncol,:) + cmfdqr2(:ncol,:)
   call aerosol_wet_intr (state, ptend, cflx, nstep, ztodt, lat, clat, lon,&
        qme, prain, &
       nevapr, cldc, cld, fracis, calday, cmfdqrt, conicw)
   call physics_update (state, tend, ptend, ztodt)

!
!     Convective transport of all trace species except water vapor and
!     cloud liquid and ice done here because we need to do the scavenging first
!     to determine the interstitial fraction.
!
   ptend%name  = 'convtran2'
   ptend%lq(:) = .true.
   ptend%lq(ixcldice) = .false.
   ptend%lq(ixcldliq) = .false.
   dpdry = 0
   do i = 1,lengath
      dpdry(i,:) = state%pdeldry(ideep(i),:)/100.
   end do
   call convtran (lchnk,                                           &
                  ptend%lq,state%q, ppcnst,     mu2,     md2,      &
                  du2,     eu2,     ed2,        dp,      dsubcld,  &
                  jt,      maxg,    ideep,      1,       lengath,  &
                  nstep,   fracis,  ptend%q,    dpdry)

   call physics_update (state, tend, ptend, ztodt)

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

!
! Compute rates of temperature and constituent change due to moist processes
!
   dtcond(:ncol,:) = (state%s(:ncol,:) - dtcond(:ncol,:))*rtdt / cpair
   dqcond(:ncol,:,:) = (state%q(:ncol,:,:) - dqcond(:ncol,:,:))*rtdt
   call outfld('DTCOND  ',dtcond, pcols   ,lchnk   )
   do m=1,ppcnst
      call outfld(dcconnam(m),dqcond(1,1,m),pcols   ,lchnk )
   end do

! Compute total convective and stratiform precipitation and snow rates
   do i=1,ncol
      precc (i) = prec_zmc(i) + prec_cmf(i)
      precl (i) = prec_sed(i) + prec_pcw(i)
      precsc(i) = snow_zmc(i) + snow_cmf(i)
      precsl(i) = snow_sed(i) + snow_pcw(i)
! jrm These checks should not be necessary if they exist in the parameterizations
      if(precc(i).lt.0.) precc(i)=0.
      if(precl(i).lt.0.) precl(i)=0.
      if(precsc(i).lt.0.) precsc(i)=0.
      if(precsl(i).lt.0.) precsl(i)=0.
      if(precsc(i).gt.precc(i)) precsc(i)=precc(i)
      if(precsl(i).gt.precl(i)) precsl(i)=precl(i)
! end jrm
   end do
   prcsnw(:ncol) = precsc(:ncol) + precsl(:ncol)   ! total snowfall rate: needed by slab ocean model
!
!===================================================
! Moist physical parameteriztions complete: 
! send dynamical variables, and derived variables to history file
!===================================================
!
   call diag_dynvar (lchnk, ncol, state)
!
!===================================================
! Radiation computations
!===================================================
!
! Cosine solar zenith angle for current time step
!
   call get_rlat_all_p(lchnk, ncol, clat)
   call get_rlon_all_p(lchnk, ncol, clon)
   call zenith (calday, clat, clon, coszrs, ncol)

   if (dosw .or. dolw) then

! Compute cloud water/ice paths and optical properties for input to radiation
      call t_startf('cldoptics')
      call param_cldoptics_calc(state, cld, landfrac, landm,icefrac, &
                                cicewp, cliqwp, emis, rel, rei, pmxrgn, nmxrgn, snowh)
      call t_stopf('cldoptics')
!
! Complete radiation calculations
!
      call t_startf ('radctl')
      call radctl (lchnk, ncol, lwup, emis, state%pmid,             &
                   state%pint, state%lnpmid, state%lnpint, state%t, state%q,   &
                   cld, cicewp, cliqwp, coszrs, asdir, asdif,               &
                   aldir, aldif, pmxrgn, nmxrgn, fsns, fsnt    ,flns    ,flnt    , &
                   qrs, qrl, flwds, rel, rei,                       &
                   sols, soll, solsd, solld,                  &
                   landfrac, state%zm, state, fsds)
      call t_stopf ('radctl')
!
! Cloud cover diagnostics
! radctl can change pmxrgn and nmxrgn so cldsav needs to follow 
! radctl.
!
      call cldsav (lchnk, ncol, cld, state%pmid, cltot, &
                   cllow, clmed, clhgh, nmxrgn, pmxrgn)
!
! Dump cloud field information to history tape buffer (diagnostics)
!
      call outfld('CLDTOT  ',cltot  ,pcols,lchnk)
      call outfld('CLDLOW  ',cllow  ,pcols,lchnk)
      call outfld('CLDMED  ',clmed  ,pcols,lchnk)
      call outfld('CLDHGH  ',clhgh  ,pcols,lchnk)
      call outfld('CLOUD   ',cld    ,pcols,lchnk)
      if (doisccp) then
         call cloudsimulator_run(state, ts, concld, cld, cliqwp, &
                                 cicewp, rel, rei, emis, coszrs  )
      end if
   else

! convert radiative heating rates from Q*dp to Q for energy conservation
      if (conserve_energy) then
         do k =1 , pver
            do i = 1, ncol
               qrs(i,k) = qrs(i,k)/state%pdel(i,k)
               qrl(i,k) = qrl(i,k)/state%pdel(i,k)
            end do
         end do
      end if
         
   end if

!
! Compute net flux
! Since fsns, fsnt, flns, and flnt are in the buffer, array values will be carried across
! timesteps when the radiation code is not invoked.
!
   do i=1,ncol
      tend%flx_net(i) = fsnt(i) - fsns(i) - flnt(i) + flns(i)
   end do
!
! Compute net radiative heating
!
   call radheat_net (state, ptend, qrl, qrs)
!
! Add radiation tendencies to cummulative model tendencies and update profiles
!
   call physics_update(state, tend, ptend, ztodt)

! check energy integrals
   call check_energy_chng(state, tend, "radheat", nstep, ztodt, zero, zero, zero, tend%flx_net)
!
! Compute net surface radiative flux for use by surface temperature code.
! Note that units have already been converted to mks in RADCTL.  Since
! fsns and flwds are in the buffer, array values will be carried across
! timesteps when the radiation code is not invoked.
!
   srfrad(:ncol) = fsns(:ncol) + flwds(:ncol)
   call outfld('SRFRAD  ',srfrad,pcols,lchnk)
!
! Save atmospheric fields to force surface models
!
   call srfxfer (lchnk, ncol, state%ps, state%u(1,pver), state%v(1,pver),    &
                 state%t(1,pver), state%q(1,pver,1), state%exner(1,pver), state%zm(1,pver), &
                    state%pmid,      &
                 state%rpdel(1,pver))

!---------------------------------------------------------------------------------------
! Save history variables. These should move to the appropriate parameterization interface
!---------------------------------------------------------------------------------------

   call outfld('PRECL   ',precl   ,pcols   ,lchnk       )
   call outfld('PRECC   ',precc   ,pcols   ,lchnk       )
   call outfld('PRECSL  ',precsl  ,pcols   ,lchnk       )
   call outfld('PRECSC  ',precsc  ,pcols   ,lchnk       )
   
   prect(:ncol) = precc(:ncol) + precl(:ncol)
   call outfld('PRECT   ',prect   ,pcols   ,lchnk       )
   call outfld('PRECTMX ',prect   ,pcols   ,lchnk       )

#if ( defined COUP_CSM )
   call outfld('PRECLav ',precl   ,pcols   ,lchnk   )
   call outfld('PRECCav ',precc   ,pcols   ,lchnk   )
#endif

#if ( defined BFB_CAM_SCAM_IOP )
   call outfld('Prec   ',prect   ,pcols   ,lchnk       )
#endif
!     
! Compute heating rate for dtheta/dt 
!
   do k=1,pver
      do i=1,ncol
         ftem(i,k) = (qrs(i,k) + qrl(i,k))/cpair * (1.e5/state%pmid(i,k))**cappa
      end do
   end do
   call outfld('HR      ',ftem    ,pcols   ,lchnk   )

! convert radiative heating rates to Q*dp for energy conservation
   if (conserve_energy) then
      do k =1 , pver
         do i = 1, ncol
            qrs(i,k) = qrs(i,k)*state%pdel(i,k)
            qrl(i,k) = qrl(i,k)*state%pdel(i,k)
         end do
      end do
   end if

   return
end subroutine tphysbc