module nlmsas_ad_mod
!$$$   module documentation block
!                .      .    .                                       .
! module:    nlmsas_ad_mod   module wrapper around subroutine nlmsas_ad
!   prgmmr: parrish          org: np22                date: 2013-01-26
!
! abstract: This module has been added as a wrapper around subroutine nlmsas_ad
!            to eliminate type mismatch compile errors when using the debug
!            compile option on WCOSS.
!
! program history log:
!   2012-01-26  parrish
!
! subroutines included:
!  nlmsas_ad
!  nlmsas_ad_1_1_
!  nlmsas_ad_im_ix_

  implicit none

! set default to private
  private
! set subroutines to public
  public :: nlmsas_ad

  interface nlmsas_ad
     module procedure nlmsas_ad_1_1_
     module procedure nlmsas_ad_im_ix_
  end interface

contains

subroutine nlmsas_ad_1_1_(im,ix,km,jcap,delt,del_,sl_,rcs_,&
     slimsk_,xkt2_,ncloud,ps_, &
     t0_,q0_,cwm0_,u0_,v0_,dot0_, &
     t1_,q1_,cwm1_,u1_,v1_,rn1_, &
     cldwrk_,kbot_,ktop_,jmin_,kuo_,kb_, &
     t0_ad_,q0_ad_,cwm0_ad_,u0_ad_,v0_ad_,dot0_ad_, &
     t1_ad_,q1_ad_,cwm1_ad_,u1_ad_,v1_ad_,rn1_ad_, &
     adjoint)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    nlmsas_ad_1_1_
!   prgmmr: parrish          org: np22                date: 2013-01-26
!
! abstract:  interface for nlmsas_ad, where im=1,ix=1, and calling routine has
!              no dimension index corresponding to im and ix.
!
! program history log:
!   2013-01-26  parrish - initial documentation
!
!   input argument list:
!     im,ix,km,jcap,ncloud,adjoint,delt
!     del_(km),sl_(km),rcs_,slimsk_,xkt2_,ps_,t0_(km),q0_(km),cwm0_(km),u0_(km),v0_(km),dot0_(km)
!     t1_(km),q1_(km),cwm1_(km),u1_(km),v1_(km),rn1_,cldwrk_,kbot_,ktop_,jmin_,kuo_,kb_
!     t0_ad_(km),q0_ad_(km),cwm0_ad_(km),u0_ad_(km),v0_ad_(km),dot0_ad_(km)
!     t1_ad_(km),q1_ad_(km),cwm1_ad_(km),u1_ad_(km),v1_ad_(km),rn1_ad_
!
!   output argument list:
!     t1_(km),q1_(km),cwm1_(km),u1_(km),v1_(km),rn1_,cldwrk_,kbot_,ktop_,jmin_,kuo_,kb_
!     t0_ad_(km),q0_ad_(km),cwm0_ad_(km),u0_ad_(km),v0_ad_(km),dot0_ad_(km)
!     t1_ad_(km),q1_ad_(km),cwm1_ad_(km),u1_ad_(km),v1_ad_(km),rn1_ad_

  use kinds, only: r_kind,i_kind
  implicit none

  integer(i_kind), intent(in   ) :: im,ix,km,jcap,ncloud
  logical,         intent(in   ) :: adjoint
  real(r_kind),    intent(in   ) :: delt

  real(r_kind),    intent(in   ) :: del_(km),sl_(km),rcs_,slimsk_
  real(r_kind),    intent(in   ) :: xkt2_,ps_,t0_(km),q0_(km),cwm0_(km),u0_(km),v0_(km)
  real(r_kind),    intent(in   ) :: dot0_(km)
  real(r_kind),    intent(inout) :: t1_(km),q1_(km),cwm1_(km)
  real(r_kind),    intent(inout) :: u1_(km),v1_(km),rn1_,cldwrk_
  integer(i_kind), intent(inout) :: kbot_,ktop_,jmin_,kuo_
  integer(i_kind), intent(inout) :: kb_
  real(r_kind),    intent(inout) :: t0_ad_(km),q0_ad_(km),cwm0_ad_(km)
  real(r_kind),    intent(inout) :: u0_ad_(km),v0_ad_(km),dot0_ad_(km)
  real(r_kind),    intent(inout) :: t1_ad_(km),q1_ad_(km),cwm1_ad_(km),u1_ad_(km),v1_ad_(km)
  real(r_kind),    intent(inout) :: rn1_ad_


  real(r_kind)     :: del(km,ix),sl(km,ix),rcs(ix),slimsk(ix)
  real(r_kind)     :: xkt2(ix),ps(ix),t0(km,ix),q0(km,ix),cwm0(km,ix),u0(km,ix),v0(km,ix)
  real(r_kind)     :: dot0(km,ix)
  real(r_kind)     :: t1(km,ix),q1(km,ix),cwm1(km,ix)
  real(r_kind)     :: u1(km,ix),v1(km,ix),rn1(ix),cldwrk(ix)
  integer(i_kind)  :: kbot(ix),ktop(ix),jmin(ix),kuo(ix)
  integer(i_kind)  :: kb(ix)
  real(r_kind)     :: t0_ad(km,ix),q0_ad(km,ix),cwm0_ad(km,ix)
  real(r_kind)     :: u0_ad(km,ix),v0_ad(km,ix),dot0_ad(km,ix)
  real(r_kind)     :: t1_ad(km,ix),q1_ad(km,ix),cwm1_ad(km,ix),u1_ad(km,ix),v1_ad(km,ix)
  real(r_kind)     :: rn1_ad(ix)

  integer(i_kind) k

  if( im /= 1 .or. ix /= 1 ) then
     write(6,*)' NLMSAS_AD_1_1_, IM,IX=',IM,IX,' -- BOTH MUST BE 1.  PROGRAM FAILS'
     stop
  end if

  rcs(1)=rcs_
  slimsk(1)=slimsk_
  xkt2(1)=xkt2_
  ps(1)=ps_
  rn1(1)=rn1_
  cldwrk(1)=cldwrk_
  kbot(1)=kbot_
  ktop(1)=ktop_
  jmin(1)=jmin_
  kuo(1)=kuo_
  kb(1)=kb_
  rn1_ad(1)=rn1_ad_
  do k=1,km
     del(k,1)=del_(k)
     sl(k,1)=sl_(k)
     t0(k,1)=t0_(k)
     q0(k,1)=q0_(k)
     cwm0(k,1)=cwm0_(k)
     u0(k,1)=u0_(k)
     v0(k,1)=v0_(k)
     dot0(k,1)=dot0_(k)
     t1(k,1)=t1_(k)
     q1(k,1)=q1_(k)
     cwm1(k,1)=cwm1_(k)
     u1(k,1)=u1_(k)
     v1(k,1)=v1_(k)
     t0_ad(k,1)=t0_ad_(k)
     q0_ad(k,1)=q0_ad_(k)
     cwm0_ad(k,1)=cwm0_ad_(k)
     u0_ad(k,1)=u0_ad_(k)
     v0_ad(k,1)=v0_ad_(k)
     dot0_ad(k,1)=dot0_ad_(k)
     t1_ad(k,1)=t1_ad_(k)
     q1_ad(k,1)=q1_ad_(k)
     cwm1_ad(k,1)=cwm1_ad_(k)
     u1_ad(k,1)=u1_ad_(k)
     v1_ad(k,1)=v1_ad_(k)
  end do

  call nlmsas_ad_im_ix_(im,ix,km,jcap,delt,del,sl,rcs,&
     slimsk,xkt2,ncloud,ps, &
     t0,q0,cwm0,u0,v0,dot0, &
     t1,q1,cwm1,u1,v1,rn1, &
     cldwrk,kbot,ktop,jmin,kuo,kb, &
     t0_ad,q0_ad,cwm0_ad,u0_ad,v0_ad,dot0_ad, &
     t1_ad,q1_ad,cwm1_ad,u1_ad,v1_ad,rn1_ad, &
     adjoint)

  rn1_=rn1(1)
  cldwrk_=cldwrk(1)
  kbot_=kbot(1)
  ktop_=ktop(1)
  jmin_=jmin(1)
  kuo_=kuo(1)
  kb_=kb(1)
  rn1_ad_=rn1_ad(1)
  do k=1,km
     t1_(k)=t1(k,1)
     q1_(k)=q1(k,1)
     cwm1_(k)=cwm1(k,1)
     u1_(k)=u1(k,1)
     v1_(k)=v1(k,1)
     t0_ad_(k)=t0_ad(k,1)
     q0_ad_(k)=q0_ad(k,1)
     cwm0_ad_(k)=cwm0_ad(k,1)
     u0_ad_(k)=u0_ad(k,1)
     v0_ad_(k)=v0_ad(k,1)
     dot0_ad_(k)=dot0_ad(k,1)
     t1_ad_(k)=t1_ad(k,1)
     q1_ad_(k)=q1_ad(k,1)
     cwm1_ad_(k)=cwm1_ad(k,1)
     u1_ad_(k)=u1_ad(k,1)
     v1_ad_(k)=v1_ad(k,1)
  end do

end subroutine nlmsas_ad_1_1_

subroutine nlmsas_ad_im_ix_(im,ix,km,jcap,delt,del,sl,rcs,&
     slimsk,xkt2,ncloud,ps, &
     t0,q0,cwm0,u0,v0,dot0, &
     t1,q1,cwm1,u1,v1,rn1, &
     cldwrk,kbot,ktop,jmin,kuo,kb, &
     t0_ad,q0_ad,cwm0_ad,u0_ad,v0_ad,dot0_ad, &
     t1_ad,q1_ad,cwm1_ad,u1_ad,v1_ad,rn1_ad, &
     adjoint)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    nlmsas_ad    forward and adjoint model for GFS SASCNV
!     prgmmr:    treadon     org: np23                date: 2003-12-18
!
! abstract:  This routine computes convective heating and moistening using a 
!   one cloud type Arakawa-Schubert convection scheme originally developed
!   by George Grell.  The scheme includes updraft and downdraft effects.
!   The closure is the cloud work function. both updraft and downdraft
!   are assumed to be saturated and the heating and moistening are
!   accomplished by the compensating environment. The name, SASCNV, comes
!   from "simplified arakawa-schubert convection parameterization".
!
!   This routine also includes the adjoint of SASCNV
!
!
! program history log:
!   1992-03-01  Hua-Lu Pan   - initial SASCNV (forward model only)
!   2004-06-14  Russ Treadon - reformat documentation
!   2004-10-28  treadon - replace parameters "tiny" and "qsmall"
!                         with global constant "tiny_r_kind"
!   2006-01-06  treadon - fix cnvflg bug
!   2006-04-12  treadon - change del and sl from 1d to 2d arrays
!   2008-04-29  safford - rm unused vars
!   2008-10-29  min-jeong kim - make consistent with global_fcst
!   2013-01-26  parrish - module added as a wrapper around subroutine nlmsas_ad
!                            to eliminate type mismatch compile errors when using the debug
!                            compile option on WCOSS.
!
!  input argument list:
!     im       - integer number of points
!     ix       - leading dimension of qn,tn,q1,t0,spd
!     km       - integer number of levels
!     jcap     - integer spectral truncation
!     delt     - time step in seconds
!     del      - sigma layer thickness
!     sl       - sigma values
!     rcs      - 1/sin(lat)**2
!     slimsk   - sea (=0), land (=1), ice (=2) flag
!     xkt2     - random number used to select cloud top
!     ncloud   - parameter to turn on cloud water adjustment
!     ps       - surface pressure in kilopascals (cb)
!     t0       - temperature
!     q0       - specific humidity
!     cwm0     - cloud condensate mixing ratio
!     u0       - zonal wind component
!     v0       - meridional wind component
!     dot0     - vertical velocity
!     t1_ad    - temperature perturbation
!     q1_ad    - specific humidity perturbation
!     cwm1_ad  - cloud condensate mixing ratio perturbation
!     u1_ad    - zonal wind perturbation
!     v1_ad    - meridional wind perturbation
!     rn1_ad   - convective rain rate perturbation
!     adjoint  - logical flag (.false.=forward model only, .true.=forward and ajoint)
!
!   output argument list:
!     t1      - temperature adjusted by convective processes
!     q1      - specific humidity adjusted by convective processes
!     cwm1    - cloud condensate mixing ratio adjusted by convective processes
!     u1      - zonal wind adjusted by convective processes
!     v1      - meridional wind adjusted by convective processes
!     rn1     - convective rain rate
!     cldwrk  - cloud work function
!     kbot    - integer model level of cloud base 
!     ktop    - integer model level of cloud top
!     jmin    - integer model level of downdraft origin
!     kuo     - bit flag indicating presence (=1) of convection (0=no convection)
!     kb      - integer model level of updraft origin
!     t0_ad   - partial derivative of convective rain rate with respect to t
!     q0_ad   - partial derivative of convective rain rate with respect to q
!     cwm0_ad - partial derivative of convective rain rate with respect to cwm
!     u0_ad   - partial derivative of convective rain rate with respect to u
!     v0_ad   - partial derivative of convective rain rate with respect to v
!     dot0_ad - partial derivative of convective rain rate with respect to vertical velocity
!     
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only: c0,tiny_r_kind,el2orc,factor1,factor2,pcpeff3,h300,&
       elocp,pcpeff2,pcpeff0,pcpeff1,epsm1,delta,eps,zero,one_tenth,one,two,three,cp,hvap,&
       half,rd,grav,r1000,r3600
  implicit none  

  integer(i_kind), intent(in   ) :: im,ix,km,jcap,ncloud
  logical adjoint
  real(r_kind),    intent(in   ) :: delt
  real(r_kind),    intent(in   ) :: del(km,ix),sl(km,ix),rcs(ix),slimsk(ix)
  real(r_kind),    intent(in   ) :: xkt2(ix),ps(ix),t0(km,ix),q0(km,ix),cwm0(km,ix),u0(km,ix),v0(km,ix)
  real(r_kind),    intent(in   ) :: dot0(km,ix)
  real(r_kind),    intent(inout) :: t1(km,ix),q1(km,ix),cwm1(km,ix)
  real(r_kind),    intent(inout) :: u1(km,ix),v1(km,ix),rn1(ix),cldwrk(ix)
  integer(i_kind), intent(inout) :: kbot(ix),ktop(ix),jmin(ix),kuo(ix)
  integer(i_kind), intent(inout) :: kb(ix)
  real(r_kind),    intent(inout) :: t0_ad(km,ix),q0_ad(km,ix),cwm0_ad(km,ix)
  real(r_kind),    intent(inout) :: u0_ad(km,ix),v0_ad(km,ix),dot0_ad(km,ix)
  real(r_kind),    intent(inout) :: t1_ad(km,ix),q1_ad(km,ix),cwm1_ad(km,ix),u1_ad(km,ix),v1_ad(km,ix)
  real(r_kind),    intent(inout) :: rn1_ad(ix)

! Local variables and arrays
  logical,dimension(ix):: cnvflg,dwnflg,dwnflg2,flg
  logical,dimension(km,ix):: flgk

  integer(i_kind) kbcon(ix),lmin(ix),ktcon(ix),kbdtr(ix)
 
  real(r_kind) term,r1200
  real(r_kind) to2(km,ix),qo2(km,ix)
  real(r_kind) to3(km,ix),qo3(km,ix)
  real(r_kind) p(km,ix),pdot(ix),acrtfct(ix),acrtfct1(ix)
  real(r_kind) acrtfct0(ix)
  real(r_kind) to(km,ix),qo(km,ix),uo(km,ix),vo(km,ix)
  real(r_kind) qeso(km,ix)
  real(r_kind) tvo(km,ix),dbyo(km,ix)
  real(r_kind) zo(km,ix)
  real(r_kind) heo(km,ix),heso(km,ix)
  real(r_kind) qrcd(km,ix),dellah(km,ix),dellaq(km,ix)
  real(r_kind) hcko(km,ix)
  real(r_kind) qcko(km,ix),eta(km,ix)
  real(r_kind) etad(km,ix)
  real(r_kind) qrcdo(km,ix),dtconv0(ix)
  real(r_kind) pwo(km,ix),pwdo(km,ix),dtconv(ix),acrt(ix)
  real(r_kind) psfc(ix),hmax(ix)
  real(r_kind) hkbo(ix),qkbo(ix),pbcdif(ix)
  real(r_kind) hmin(ix),pwavo(ix)
  real(r_kind) aa1(ix),vshear(ix)
  real(r_kind) edt(ix)
  real(r_kind) edto(ix),pwevo(ix)
  real(r_kind) hcdo(ix),qcdo(ix)
  real(r_kind) xhkb(ix),xqkb(ix),xpwav(ix),xpwev(ix),xhcd(ix)
  real(r_kind) xaa0(ix),f(ix),xk(ix),xmb(ix)
  real(r_kind) edtx(ix),xqcd(ix)
  real(r_kind) hsbar(ix),xmbmax(ix),xlamb(ix),xlamd(ix)
  real(r_kind) xlamdet(ix),xlamdet0(ix),xlamdet1(ix)
  real(r_kind) pcrit(15), acritt(15), acrit(15)
  real(r_kind) mbdt
  real(r_kind) rntot(ix)
  
  integer(i_kind) kt2(ix)
      
!new
  real(r_kind) qol(km,ix),uol(km,ix),vol(km,ix),tol(km,ix)
  real(r_kind) heol(km,ix),hesol(km,ix),qesol(km,ix),qol0(km,ix)
  real(r_kind) xqo(km,ix),xto(km,ix),xtvo(km,ix),xzo(km,ix),xqo0(km,ix)
  real(r_kind) xheo(km,ix),xheso(km,ix),xqeso(km,ix),xqrcd(km,ix)
  real(r_kind) xhcko(km,ix),xqcko(km,ix),xtol(km,ix),xqol(km,ix)
  real(r_kind) xqesol(km,ix),aa10(ix),xqcko0(km,ix),xqrch0(km,ix)
  real(r_kind) dpscl(ix),tdpscl(ix),dpcld(ix),tdpcld(ix)
  real(r_kind) tcwfup(ix),tcwfdn(ix),tdetrn(ix),tdpdwn(ix)
  real(r_kind) eta0(km,ix),hcko0(km,ix),hckod(km,ix),hcko1(km,ix)
  real(r_kind) qcko0(km,ix),dqk(km,ix),edt0(ix),qcko00(km,ix)
  real(r_kind) xtemp0(km,ix),xmb0(ix),rn11(ix)
  real(r_kind) qeso2(km,ix),rn0k(km,ix),delqevk(km,ix)
  real(r_kind) edto1(ix),edto10(ix),edtx1(ix),edtx10(ix)
!new

!ajm
  real(r_kind) to_ad(km,ix),qo_ad(km,ix),qeso_ad(km,ix),tvo_ad(km,ix)
  real(r_kind) zo_ad(km,ix),heo_ad(km,ix),heso_ad(km,ix)
  real(r_kind) tol_ad(km,ix),qol_ad(km,ix),hesol_ad(km,ix),heol_ad(km,ix)
  real(r_kind) qesol_ad(km,ix),hkbo_ad(ix),xlamb_ad(ix),eta_ad(km,ix)
  real(r_kind) hcko_ad(km,ix),qcko_ad(km,ix),xlamdet_ad(ix)
  real(r_kind) hcko1_ad(km,ix),dbyo_ad(km,ix),hckod_ad(km,ix)
  real(r_kind) pwo_ad(km,ix),aa1_ad(ix),pwavo_ad(ix),tcwfup_ad(ix)
  real(r_kind) tcwfdn_ad(ix),vshear_ad(ix),edt_ad(ix),edto_ad(ix)
  real(r_kind) edtx_ad(ix),uo_ad(km,ix),vo_ad(km,ix),uol_ad(km,ix)
  real(r_kind) vol_ad(km,ix),xlamd_ad(ix),etad_ad(km,ix),qrcdo_ad(km,ix)
  real(r_kind) pwdo_ad(km,ix),pwevo_ad(ix),qcdo_ad(ix),hcdo_ad(ix)
  real(r_kind) dellah_ad(km,ix),dellaq_ad(km,ix),xto_ad(km,ix)
  real(r_kind) xqo_ad(km,ix),xqeso_ad(km,ix),xtvo_ad(km,ix)
  real(r_kind) xzo_ad(km,ix),xtol_ad(km,ix),xqol_ad(km,ix)
  real(r_kind) xqesol_ad(km,ix),xheo_ad(km,ix),xheso_ad(km,ix)
  real(r_kind) xhcko_ad(km,ix),xqcko_ad(km,ix),xaa0_ad(ix),xpwav_ad(ix)
  real(r_kind) xpwev_ad(ix),xqrcd_ad(km,ix),xhcd_ad(ix),pdot_ad(ix)
  real(r_kind) acrtfct1_ad(ix),acrtfct_ad(ix),dtconv_ad(ix)
  real(r_kind) f_ad(ix),xk_ad(ix),xmb_ad(ix),to2_ad(km,ix)
  real(r_kind) qo2_ad(km,ix),qeso2_ad(km,ix),rntot_ad(ix)
  real(r_kind) to3_ad(km,ix),qo3_ad(km,ix),edto1_ad(ix),edtx1_ad(ix)
!ajm

!new
  integer(i_kind) ktcon0(ix)
  real(r_kind) sumz(km,ix),sumh(km,ix),cwmo_ad(km,ix)
  real(r_kind) ucko(km,ix),vcko(km,ix),etau(km,ix),etau0(km,ix)
  real(r_kind) ucko0(km,ix),vcko0(km,ix),ucdo(ix),vcdo(ix)
  real(r_kind) uckod(km,ix),vckod(km,ix),dellalk_ad(km,ix)
  real(r_kind) ukbo(ix),vkbo(ix),qlko_ktcon(ix),dellal(ix),dellal1(ix)
  real(r_kind) dellalk(km,ix),cwmo2_ad(km,ix),dellal1_ad(ix)
  real(r_kind) dellau(km,ix),dellav(km,ix),uo2(km,ix),vo2(km,ix)
  real(r_kind) cwmo(km,ix),cwmo2(km,ix)
  real(r_kind) qo0(km,ix),qo0_ad(km,ix),dellal_ad(ix),xqol0_ad(km,ix)
  real(r_kind) uo2_ad(km,ix),vo2_ad(km,ix),xqol0(km,ix),xqo0_ad(km,ix)
  real(r_kind) ucko_ad(km,ix),vcko_ad(km,ix),ukbo_ad(ix),vkbo_ad(ix)
  real(r_kind) ucdo_ad(ix),vcdo_ad(ix),ukcod_ad(km,ix),vckod_ad(km,ix)
  real(r_kind) dellau_ad(km,ix),dellav_ad(km,ix),qlko_ktcon_ad(ix)
  real(r_kind) etau_ad(km,ix),xlamdet1_ad(ix),xlamdet0_ad(ix)
  real(r_kind) sumz_ad(km,ix),sumh_ad(km,ix),qol0_ad(km,ix)
  
  integer(i_kind) k,indx,jmn,i
  integer(i_kind),dimension(ix):: kbmax,kbm,kmax
  
  real(r_kind) rterm,fjcap
  real(r_kind) etah_ad,esl_ad,dqsdp_ad,xpw_ad,xqrch_ad,qlk_ad,dz_ad,dz1_ad
  real(r_kind) dvv1_ad,dvu1_ad,dvq1_ad,dv1_ad,pprime_ad,es_ad,desdt_ad,qs_ad
  real(r_kind) termq_ad,term1_ad,term2_ad,term3_ad,term4_ad,dv3v_ad,dv2v_ad
  real(r_kind) dv1v_ad,dv3u_ad,dv2u_ad,dv1u_ad,dv1q_ad
  real(r_kind) dv2_ad,dv3_ad,detad_ad,detau_ad,termu_ad,termv_ad
  real(r_kind) shear_ad,rain_ad,dellat_ad,ratio_ad,dt_ad,dqs_ad,es0_ad
  real(r_kind) evef_ad,qevap_ad,delqev_ad,rn0_ad,term5_ad,dhh_ad,term_ad
  real(r_kind) dv2q_ad,e1_ad,temp_ad,qrch_ad,dv3q_ad,onemf_ad,xtemp_ad
  real(r_kind) factor_ad,dqsdt_ad,xqc_ad,vol2_ad,uol2_ad,tem1_ad,qc_ad
  real(r_kind) qcond_ad,fuv_ad,onemfu_ad,xdby_ad,xpwd_ad,dq_ad,rfact_ad
  real(r_kind) heol2_ad,dh_ad,dg_ad,term6_ad,tem2_ad,gamma_ad
  
  real(r_kind) dv1,dv2,dv3,dv1q,dv2q,dv3q,dv1u,dv2u,dv3u,dv1v,dv2v,dv3v,dp
  real(r_kind) detau,detad,term1,term2,term3,term4,termq,shear,termu,termv
  real(r_kind) betas,beta,betal,xtemp,aup,adw,rain,dvu1,dvq1,dhh,dg,dh,gamma
  real(r_kind) factor,onemf,xdby,xqrch,dq,qc,qlk,qrch,temp,etah,dellat,edtmax
  real(r_kind) edtmaxl,edtmaxs,dz1,rfact,term6,qevap0,qevap1,delq2,qcond,term5
  real(r_kind) w1l,w2l,w3l,w4l,w1s,w2s,w3s,w4s,w1,w2,w3,w4,ratio,xpwd,evef,rn0
  real(r_kind) xpw,xqc,dqs,heol2,fuv,onemfu,uol2,vol2,dvv1,fixed,fkm,dtmin,dtmax
  real(r_kind) dt2,alphal,alphas,adt,es0,po,esl,alpha,pdetrn,pdpdwn,xlambu
  real(r_kind) dz,qevap,desdt,dqsdt,dqsdp,dt,pprime,qs,tem1,tem2,e1,dif,val,es
  real(r_kind) delqev,dlnsig,evfact
  
  save pcrit, acritt
  data pcrit/850._r_kind,800._r_kind,750._r_kind,700._r_kind,&
       650._r_kind,600._r_kind,550._r_kind,500._r_kind,&
       450._r_kind,400._r_kind,350._r_kind,h300,&
       250._r_kind,200._r_kind,150._r_kind/
  data acritt/.0633_r_kind,.0445_r_kind,.0553_r_kind,.0664_r_kind,&
       .075_r_kind,.1082_r_kind,.1521_r_kind,.2216_r_kind,&
       .3151_r_kind,.3677_r_kind,.41_r_kind,.5255_r_kind,&
       .7663_r_kind,1.1686_r_kind,1.6851_r_kind/

!-----------------------------------------------------------------------
  r1200=1200._r_kind
  do i = 1,ix
     do k = 1,km
        cwmo_ad(k,i)=zero
        cwmo2_ad(k,i)=zero
        dellalk_ad(k,i)=zero
        t0_ad(k,i) = zero
        q0_ad(k,i) = zero
        u0_ad(k,i) = zero
        v0_ad(k,i) = zero
        cwm0_ad(k,i)=zero
        dot0_ad(k,i) = zero
        tvo_ad(k,i) = zero
        qeso_ad(k,i) = zero
        qo_ad(k,i) = zero
        to_ad(k,i) = zero
        uo_ad(k,i) = zero
        vo_ad(k,i) = zero
        zo_ad(k,i) = zero
        heso_ad(k,i) = zero
        heo_ad(k,i) = zero
        tol_ad(k,i) = zero
        qol_ad(k,i) = zero
        uol_ad(k,i) = zero
        vol_ad(k,i) = zero
        heol_ad(k,i) = zero
        hesol_ad(k,i) = zero
        qesol_ad(k,i) = zero
        eta_ad(k,i) = zero
        hcko_ad(k,i) = zero
        qcko_ad(k,i) = zero
        hcko1_ad(k,i) = zero
        dbyo_ad(k,i) = zero
        hckod_ad(k,i) = zero
        qcko_ad(k,i) = zero
        pwo_ad(k,i) = zero
        etad_ad(k,i) = zero
        qrcdo_ad(k,i) = zero
        pwdo_ad(k,i) = zero
        dellah_ad(k,i) = zero
        dellaq_ad(k,i) = zero
        xto_ad(k,i) = zero
        xqo_ad(k,i) = zero
        xqeso_ad(k,i) = zero
        xtvo_ad(k,i) = zero
        xzo_ad(k,i) = zero
        xtol_ad(k,i) = zero
        xqol_ad(k,i) = zero
        xqesol_ad(k,i) = zero
        xheo_ad(k,i) = zero
        xheso_ad(k,i) = zero
        xhcko_ad(k,i) = zero
        xqcko_ad(k,i) = zero
        xqrcd_ad(k,i) = zero
        to2_ad(k,i) = zero
        qo2_ad(k,i) = zero
        qeso2_ad(k,i) = zero
        to3_ad(k,i) = zero
        qo3_ad(k,i) = zero
        qo0_ad(k,i)=zero
        xqol0_ad(k,i)=zero
        uo2_ad(k,i)=zero
        vo2_ad(k,i)=zero
        xqo0_ad(k,i)=zero
        ucko_ad(k,i)=zero
        vcko_ad(k,i)=zero
        ukcod_ad(k,i)=zero
        vckod_ad(k,i)=zero
        dellau_ad(k,i)=zero
        dellav_ad(k,i)=zero
        etau_ad(k,i)=zero
        sumz_ad(k,i)=zero
        sumh_ad(k,i)=zero
        qol0_ad(k,i)=zero
!
        p(k,i) = zero
        to(k,i) = zero
        qo(k,i) = zero
        to2(k,i) = zero
        qo2(k,i) = zero
        qeso(k,i) = zero
        qeso2(k,i) = zero
        to3(k,i) = zero
        qo3(k,i) = zero
        tvo(k,i) = zero
        dbyo(k,i) = zero
        zo(k,i) = zero
        heo(k,i) = zero
        heso(k,i) = zero
        xqo(k,i) = zero
        xqeso(k,i) = zero
        xto(k,i) = zero
        xtvo(k,i) = zero
        xzo(k,i) = zero
        xheo(k,i) = zero
        xheso(k,i) = zero
        tol(k,i) = zero
        qol(k,i) = zero
        qesol(k,i) = zero
        heol(k,i) = zero
        hesol(k,i) = zero
        uol(k,i) = zero
        vol(k,i) = zero
        xtol(k,i) = zero
        xqol(k,i) = zero
        xqesol(k,i) = zero
        qrcd(k,i) = zero
        dellah(k,i) = zero
        dellaq(k,i) = zero
        hcko(k,i) = zero
        qcko(k,i) = zero
        eta(k,i) = zero
        hcko0(k,i) = zero
        qcko0(k,i) = zero
        qcko00(k,i)= zero
        dqk(k,i) = zero
        eta0(k,i) = zero
        etad(k,i) = zero
        xhcko(k,i) = zero
        xqcko(k,i) = zero
        xqcko0(k,i) = zero
        xtemp0(k,i) = zero
        xqrch0(k,i) = zero
        qrcdo(k,i) = zero
        pwo(k,i) = zero
        pwdo(k,i) = zero
        hcko1(k,i) = zero
        hckod(k,i) = zero
        xqrcd(k,i) = zero
        rn0k(k,i) = zero
        delqevk(k,i) = zero
        flgk(k,i) = .true.
!new
        sumz(k,i) = zero
        sumh(k,i) = zero
        ucko(k,i) = zero
        vcko(k,i) = zero
        etau(k,i) = zero
        uckod(k,i)= zero
        vckod(k,i)=zero
        hckod(k,i)=zero
        
     end do
  end do
  do i = 1,ix
     xlamdet1_ad(i)=zero
     xlamdet0_ad(i)=zero
     ukbo_ad(i)=zero
     vkbo_ad(i)=zero
     ucdo_ad(i)=zero
     vcdo_ad(i)=zero
     qlko_ktcon_ad(i)=zero
     dellal_ad(i)=zero
     dellal1_ad(i)=zero
     hkbo_ad(i) = zero
     xlamb_ad(i) = zero
     xlamdet_ad(i) = zero
     aa1_ad(i) = zero
     pwavo_ad(i) = zero
     tcwfup_ad(i) = zero
     tcwfdn_ad(i) = zero
     vshear_ad(i) = zero
     edt_ad(i) = zero
     edto_ad(i) = zero
     edtx_ad(i) = zero
     edto1_ad(i) = zero
     edtx1_ad(i) = zero
     xlamd_ad(i) = zero
     pwevo_ad(i) = zero
     qcdo_ad(i) = zero
     hcdo_ad(i) = zero
     xaa0_ad(i) = zero
     xpwav_ad(i) = zero
     xpwev_ad(i) = zero
     xhcd_ad(i) = zero
     pdot_ad(i) = zero
     acrtfct1_ad(i) = zero
     acrtfct_ad(i) = zero
     dtconv_ad(i) = zero
     f_ad(i) = zero
     xk_ad(i) = zero
     xmb_ad(i) = zero
     rntot_ad(i) = zero
     
     dpscl(i) = zero
     dpcld(i) = zero
     tdpscl(i) = zero
     tdpcld(i) = zero
     tdetrn(i) = zero
     tcwfup(i) = zero
     tcwfdn(i) = zero
     acrtfct(i) = zero
     acrtfct0(i) = zero
     acrtfct1(i) = zero
     acrt(i) = zero
     psfc(i) = zero
     hmax(i) = zero
     kb(i) = 0
     hkbo(i) = zero
     qkbo(i) = zero
     kbcon(i) = 0
     pbcdif(i) = zero
     hmin(i) = zero
     lmin(i) = 0
     jmin(i) = 0
     pwavo(i) = zero
     aa1(i) = zero
     vshear(i) = zero
     edt(i) = zero
     edt0(i) = zero
     edto(i) = zero
     edto1(i) = zero
     edto10(i) = zero
     pwevo(i) = zero
     hcdo(i) = zero
     qcdo(i) = zero
     xhkb(i) = zero
     xqkb(i) = zero
     xpwav(i) = zero
     xpwev(i) = zero
     xhcd(i) = zero
     xaa0(i) = zero
     f(i) = zero
     xk(i) = zero
     xmb(i) = zero
     xmb0(i) = zero
     ktcon(i) = 0
     edtx(i) = zero
     edtx1(i) = zero
     edtx10(i) = zero
     xqcd(i) = zero
     hsbar(i) = zero
     xmbmax(i) = zero
     xlamb(i) = zero
     xlamdet(i) = zero
     xlamdet0(i)=zero
     xlamdet1(i)=zero
     xlamd(i) = zero
     kbdtr(i) = 0
     pdot(i) = zero
     cldwrk(i) = zero
     kbot(i) = km+1
     ktop(i) = 0
     kuo(i) = 0
     rn1(i) = zero
     aa10(i) = zero
     rn11(i) = zero
     qlko_ktcon(i)=zero
  end do


! Initialize arrays
  do i = 1,im
     cnvflg(i)  = .true.
     dwnflg(i)  = .true.
     dwnflg2(i) = .true.
     flg(i)     = .true.
     dtconv(i)  = r3600
     dtconv0(i) = dtconv(i)
     xmbmax(i)  = one_tenth
  enddo
  do k = 1, 15
     acrit(k) = acritt(k) * (975._r_kind - pcrit(k))
  enddo
  dt2   = delt
  dtmin = max(dt2,r1200)
  dtmax = max(dt2,r3600)

! Model tunable parameters are all here
  mbdt    = 10._r_kind
  edtmaxl = 0.3_r_kind
  edtmaxs = 0.3_r_kind
  alphal  = half
  alphas  = half
  betal   = 0.15_r_kind
  betas   = 0.15_r_kind
  betal   = 0.05_r_kind
  betas   = 0.05_r_kind
  evfact = 0.3_r_kind
  pdpdwn  = zero
  pdetrn  = 200._r_kind
  xlambu  = 1.e-4_r_kind
  fjcap = (float(jcap) / 126._r_kind) ** 2
  val   =  one
  fjcap = max(fjcap,val)
  fkm = (float(km) / 28._r_kind) ** 2
  fkm = max(fkm,one)
  w1l = -8.e-3_r_kind
  w2l = -4.e-2_r_kind
  w3l = -5.e-3_r_kind
  w4l = -5.e-4_r_kind
  w1s = -2.e-4_r_kind
  w2s = -2.e-3_r_kind
  w3s = -1.e-3_r_kind
  w4s = -2.e-5_r_kind

!  Define top layer for search of the downdraft originating layer
!  and the maximum thetae for updraft

  do i=1,im
     kbmax(i) = km
     kbm(i)   = km
     kmax(i)  = km
     do k = 1, km
        if (sl(k,i)  >  0.45_r_kind) kbmax(i) = k + 1
        if (sl(k,i)  >  0.7_r_kind)  kbm(i)   = k + 1
        if (sl(k,i)  >  0.04_r_kind) kmax(i)  = k + 1
     enddo
  end do
! kbm=10, kbmax=14, kmax=25 (=23 if sl(k) > 06)


! Convert surface pressure to mb from cb
  do i = 1, im
     psfc(i) = ps(i) * 10._r_kind
  enddo
  do i = 1,im
     do k = 1,km
        p(k,i)  = psfc(i) * sl(k,i)
        to(k,i) = t0(k,i)
!       qo(k,i) = q0(k,i)
        if (q0(k,i)>zero) then
           qo(k,i) = q0(k,i)
        else
           qo(k,i) = tiny_r_kind
        endif
        uo(k,i) = u0(k,i)
        vo(k,i) = v0(k,i)
        cwmo(k,i)= cwm0(k,i)
        
        qo0(k,i) = q0(k,i)
        
     enddo
  enddo

! Column variables
! p is pressure of the layer (mb)
! t is temperature at t-dt (k)..tn
! q is specific humidity at t-dt (kg/kg)..qn
! to is temperature at t+dt (k)... this is after advection and turbulan
! qo is specific humidity at t+dt (kg/kg)..q1
  do i = 1,im
     do k = 1, kmax(i)
        call fpvsx_ad(to(k,i),es0,adt,es0_ad,.false.)
        es = 10.0_r_kind*es0
        qeso(k,i) = eps*es / (p(k,i) + epsm1*es)
        
!new code bounds qeso to be >= 1.e-8, qo>=1.e-10

        tvo(k,i)  = to(k,i) + delta * to(k,i) * qo(k,i)
     enddo
  enddo

! Hydrostatic height assume zero terr
  do i = 1, im
     dlnsig = log(sl(1,i))
     zo(1,i) = zero - dlnsig * rd / grav * tvo(1,i)
  enddo

!round
!     reverse order of this loop leads to slightly different 
!     final results.  change order of loop in sascnv2 to
!     restore similar results.
  do i = 1,im
     do k = 2, kmax(i)
        dlnsig  = log(sl(k,i) / sl(k-1,i))
        term1   = dlnsig * rd / grav
        term2   = half * (tvo(k,i) + tvo(k-1,i))
        zo(k,i) = zo(k-1,i) - term1*term2
     enddo
  enddo

! Compute moist static energy
  do i = 1,im
     do k = 1, kmax(i)
        heo(k,i)  = grav * zo(k,i) + cp * to(k,i) + hvap * qo(k,i)
        heso(k,i) = grav * zo(k,i) + cp * to(k,i) + hvap * qeso(k,i)
     enddo
  enddo

! Determine level with largest moist static energy
! this is the level where updraft starts
  do i = 1, im
     hmax(i) = heo(1,i)
     kb(i) = 1
  enddo
  do i = 1,im
     do k = 2, kbm(i)
        if (heo(k,i) > hmax(i)) then
           kb(i) = k
           hmax(i) = heo(k,i)
        endif
     enddo
  enddo

! Search for downdraft originating level above theta-e minimum
!      do i = 1, im
!         hmin(i) = heso(1,i)
!         lmin(i) = kbmax(i)
!         jmin(i) = kbmax(i)
!      enddo
!      do i = 1,im
!         do k = 2, kbmax(i)
!            if (heso(k,i) < hmin(i)) then
!               lmin(i) = k + 1
!               hmin(i) = heso(k,i)
!            endif
!         enddo
!      enddo
  do i = 1,im
     do k = 1, kmax(i) - 1
        dz = half * (zo(k+1,i) - zo(k,i))
        dp = half * (p(k+1,i) - p(k,i))
        call fpvsx_ad(to(k+1,i),es0,adt,es0_ad,.false.)
        es = 10._r_kind*es0
        pprime = p(k+1,i) + epsm1 * es
        qs     = eps * es / pprime
        dqsdp  = - qs / pprime
        desdt = es * (factor1 / to(k+1,i) + factor2 / (to(k+1,i)**2))
        dqsdt = qs * p(k+1,i) * desdt / (es * pprime)
        gamma = el2orc * qeso(k+1,i) / (to(k+1,i)**2)
        dt    = (grav * dz + hvap * dqsdp * dp) / (cp * (one + gamma))
        dq    = dqsdt * dt + dqsdp * dp
        tol(k,i) = to(k+1,i) + dt
        qol0(k,i) = qo(k+1,i) + dq
        if (qol0(k,i)>zero) then
           qol(k,i) = qol0(k,i)
        else
           qol(k,i) = tiny_r_kind   
        endif
     end do
  end do

  do i = 1,im
     do k = 1,kmax(i)-1
        po    = half * (p(k,i) + p(k+1,i))
        call fpvsx_ad(tol(k,i),es0,adt,es0_ad,.false.)
        esl   = 10.0_r_kind*es0
        qesol(k,i) = eps*esl / (po + epsm1*esl)
        heol(k,i)  = half * grav * (zo(k,i) + zo(k+1,i)) + &
             cp * tol(k,i) + hvap * qol(k,i)
        hesol(k,i) = half * grav * (zo(k,i) + zo(k+1,i)) + &
             cp * tol(k,i) + hvap * qesol(k,i)
        uol(k,i)   = half * (uo(k,i) + uo(k+1,i))
        vol(k,i)   = half * (vo(k,i) + vo(k+1,i))
     enddo
  enddo
  do i = 1, im
     k = kmax(i)
     heol(k,i) = heo(k,i)
     hesol(k,i) = heso(k,i)
  enddo

! Look for convective cloud base as the level of free convection
  do i = 1, im
     indx     = kb(i)
     hkbo(i)  = heol(indx,i)
     ukbo(i)  = uol(indx,i)
     vkbo(i)  = vol(indx,i)
     kbcon(i) = kmax(i)
     flg(i)   = .true.
  enddo
  do i = 1,im
     do k = 1, kbmax(i)
        if (flg(i) .and. k > kb(i)) then
           hsbar(i) = hesol(k,i)
           if (hkbo(i) > hsbar(i)) then
              flg(i)   = .false.
              kbcon(i) = k
           endif
        endif
     enddo
  enddo
  do i = 1, im
     dpscl(i) = -p(kbcon(i),i) + p(kb(i),i)
!         if (abs(dpscl(i)) > tiny_r_kind) then
!            tdpscl(i) = funsqr(150.,dpscl(i),1)
!         else
!            tdpscl(i) = one
!         endif
     tdpscl(i) = one
     pdot(i) = 10._r_kind* dot0(kbcon(i),i)
     if (dpscl(i) > 150._r_kind) tdpscl(i) = zero
     if (kbcon(i)==kmax(i)) tdpscl(i) = zero
  enddo

! Found LFC, can define rest of variables
! determine entrainment rate between kb and kbcon
  do i = 1, im
     alpha = alphas
     if (nint(slimsk(i)) == 1) alpha = alphal
     if (kb(i)==1) then
        dz = half * (zo(kbcon(i),i) + zo(kbcon(i)-1,i)) - zo(1,i)
     else
        dz = half * (zo(kbcon(i),i) + zo(kbcon(i)-1,i)) &
             - half * (zo(kb(i),i) + zo(kb(i)-1,i))
     endif
     if (kbcon(i)/=kb(i)) then
        xlamb(i) = -log(alpha) / dz
     else
        xlamb(i) = zero
     endif
  enddo

! Determine updraft mass flux
  do i = 1,im
     do k = 1, kmax(i)
        eta(k,i) = one
        etau(k,i)= one
     enddo
  enddo
  do i = 1,im
     do k = kbmax(i), 2, -1
        if (k < kbcon(i) .and. k >= kb(i)) then
           dz       = half * (zo(k+1,i) - zo(k-1,i))
           eta(k,i) = eta(k+1,i) * exp(-xlamb(i) * dz)
           etau(k,i)= eta(k,i)
        endif
     enddo
  enddo
  do i = 1, im
     if (kb(i)==1 .and. kbcon(i) > 1) then
        dz = half * (zo(2,i) - zo(1,i))
        eta(1,i) = eta(2,i) * exp(-xlamb(i) * dz)
        etau(1,i)= eta(1,i)
     endif
  enddo

! Work up updraft cloud properties
  do i = 1, im
     indx = kb(i)
     hcko(indx,i) = hkbo(i)
     ucko(indx,i) = ukbo(i)
     vcko(indx,i) = vkbo(i)
  enddo

! Cloud property below cloud base is modified by the entrainment proces
  do i = 1,im
     do k = 2, kmax(i) - 1
        if (k > kb(i) .and. k <= kbcon(i)) then
           factor = eta(k-1,i) / eta(k,i)
           onemf = one - factor
           hcko(k,i) = factor * hcko(k-1,i) + onemf * &
                half * (heol(k,i) + heol(k+1,i))
           ucko(k,i) = factor * ucko(k-1,i) + onemf * &
                half * (uol(k,i) + uol(k+1,i))
           vcko(k,i) = factor * vcko(k-1,i) + onemf * &
                half * (vol(k,i) + vol(k+1,i))
        endif
        if (k > kbcon(i)) then
           hcko(k,i) = hcko(kbcon(i),i)
           ucko(k,i) = ucko(kbcon(i),i)
           vcko(k,i) = vcko(kbcon(i),i)
        endif
     enddo
  enddo

! Load eta into new work array to preserve previous value
  do i = 1,im
     do k = 1,km
        eta0(k,i) = eta(k,i)
        hcko0(k,i) = hcko(k,i)
        ucko0(k,i) = ucko(k,i)
        vcko0(k,i) = vcko(k,i)
        etau0(k,i) = etau(k,i)
     end do
  end do

! Determine cloud top
  do i = 1, im
     flg(i) = .true.
     ktcon(i) = 1
  enddo
  do i = 1,im
     do k = 2, kmax(i)
        dif = hcko(k,i) - hesol(k,i)
        if (dif < zero .and. flg(i) .and. k > kbcon(i)) then
           ktcon(i) = k
           flg(i) = .false.
        endif
     enddo
  enddo


! Make sure that jmin is within the cloud
!new
! Search for downdraft origin level above theta-e minimum
  do i=1,im
     hmin(i) = heol(kbcon(i),i)
     lmin(i) = kbmax(i)
     jmin(i) = kbmax(i)
     if(cnvflg(i)) then
        do k = kbcon(i), kbmax(i)
           if(heol(k,i) < hmin(i)) then
              lmin(i) = k + 1
              hmin(i) = heol(k,i)
           endif
        enddo
     endif
  end do
  
! Make sure jmin is within the cloud      
  do i = 1, im
     jmin(i) = min(lmin(i),ktcon(i)-1)
     jmin(i) = max(jmin(i),kbcon(i)+1)
     if (ktcon(i) < kbcon(i)) jmin(i) = 1
  enddo

  do i = 1, im
     dpcld(i)  = p(kbcon(i),i) - p(ktcon(i),i)
!        tdpcld(i) = funsqr(dpcld(i),150.,1)
     tdpcld(i) = one
     if (dpcld(i) < 150._r_kind) tdpcld(i) = zero
  enddo

!new
!  KTCON is model level of deepest cloud in ensemble.
!  Now randomly select cloud from ensemble of clouds
!  ranging from one level above the downdraft origin to
!  to level of neutral buoyancy (ktcon)

  do i=1,im
     do k=2,kmax(i)-1
        if (k > jmin(i) .and. k <= ktcon(i)) then
           sumz(k,i) = sumz(k-1,i) + half * (zo(k+1,i) - zo(k-1,i))
           sumh(k,i) = sumh(k-1,i) + half * (zo(k+1,i) - zo(k-1,i)) &
                * heol(k,i)
        endif
     end do
  end do

! Select cloud from ensemble
  do i=1,im
     kt2(i) = nint(xkt2(i)*float(ktcon(i)-jmin(i))-half) + jmin(i) + 1
     tem1 = hcko(jmin(i),i) - hesol(kt2(i),i)
     tem2 = sumz(kt2(i),i) * hesol(kt2(i),i) - sumh(kt2(i),i)
     if (abs(tem2)  >  0.000001_r_kind) then
        xlamdet0(i) = tem1 / tem2
     else
        xlamdet0(i) = zero
        cnvflg(i) = .false.
     endif
     if (xlamdet0(i) >= zero) then
        xlamdet1(i) = xlamdet0(i)
     else
        xlamdet1(i) = zero
     endif
     if (sumz(kt2(i),i)  >  0.000001_r_kind) then
        term = 2.3_r_kind/sumz(kt2(i),i)
        if (xlamdet1(i) < term) then
           xlamdet(i) = xlamdet1(i)
        else
           xlamdet(i) = term
        endif
     else
        xlamdet(i) = xlamdet1(i)
     endif
  end do

! Apply logic to determine if convection is possible
  do i=1,im
     if (kt2(i) >= ktcon(i)) dwnflg(i) = .false.
     if(xlamdet(i) <= 1.e-30_r_kind .or. &
          hcko(jmin(i),i)-hesol(kt2(i),i) <= 1.e-30_r_kind) &
          dwnflg(i) = .false.
     do k = jmin(i), kt2(i)
        if(dwnflg(i) .and. &
             heol(k,i) > hesol(kt2(i),i)) dwnflg(i)=.false.
     enddo
  end do
         
         


! Detraining cloud
  do i = 1, im
!        if (abs(dpcld(i)) > tiny_r_kind) then
!           tdetrn(i) = funsqr(pdetrn,dpcld(i),1)
!        else
!           tdetrn(i) = one
!        endif
     tdetrn(i) = one
     if (dpcld(i) < pdpdwn) dwnflg2(i)=.false.
!        if (abs(pdpdwn) > tiny_r_kind) then
!           tdpdwn(i) = funsqr(dpcld(i),pdpdwn,1)
!        else
!           tdpdwn(i) = one
!        endif
     tdpdwn(i) = one
     if (dwnflg2(i) .and. jmin(i) <= kbcon(i)) then
        dwnflg2(i) = .false.
        tdetrn(i)  = zero
        tdpdwn(i)  = zero
     endif
     
  enddo

! Cloud property above cloud top is modified by the detrainment process
  do i=1,im
     if (dwnflg(i)) then
        do k=2,kmax(i)-1
           if (k > jmin(i) .and. k <= kt2(i)) then
              dz = half * (zo(k+1,i) - zo(k-1,i))
!   To simplify matters, we will take the linear approach here
              eta(k,i)  = eta(k-1,i)  * (one + xlamdet(i) * dz)
              etau(k,i) = etau(k-1,i) * (one +(xlamdet(i)+xlambu)* dz)
           endif
        end do
     else
        do k=2,kmax(i)-1
           if (k > jmin(i) .and. k <= ktcon(i)) then
              dz = half * (zo(k+1,i) - zo(k-1,i))
              etau(k,i) = etau(k-1,i) * (one + xlambu * dz)
           endif
        end do
     endif
     
  enddo
  do i=1,im
     if (dwnflg(i)) then
        ktcon0(i) = ktcon(i)
        ktcon(i) = kt2(i)
     else
        ktcon0(i) = ktcon(i)
     endif
  end do

! Load hcko at cloud base into hckod
  do i = 1,im
     if (dwnflg(i)) then
        indx = kbcon(i)
        hckod(indx,i) = hcko(indx,i)
        uckod(indx,i) = ucko(indx,i)
        vckod(indx,i) = vcko(indx,i)
     endif
  end do

! Adjustment to hckod by detraining cloud
  do i = 1,im
     do k = 2,kmax(i)-1
        if (k > kbcon(i) .and. k <= ktcon(i)) then
           factor = eta(k-1,i) / eta(k,i)
           onemf  = one - factor
           fuv    = etau(k-1,i)/etau(k,i)
           onemfu = one - fuv
           heol2  = half*(heol(k,i) + heol(k+1,i))
           hckod(k,i) = factor*hckod(k-1,i) + onemf*heol2
           
           uol2       = half*(uol(k,i)+uol(k+1,i))
           uckod(k,i) = fuv * uckod(k-1,i) + onemfu*uol2
           
           vol2       = half*(vol(k,i)+vol(k+1,i))
           vckod(k,i) = fuv * vckod(k-1,i) + onemfu*vol2
           
        endif
     end do
  end do

! Load hckod and hcko into hcko1 and calculate dbyo
  do i=1,im
     do k=2,kmax(i)-1
        if (k > kb(i).and.k <= kbcon(i)) then
           hcko1(k,i) = hcko(k,i)
        elseif (k > kbcon(i).and.k <= ktcon(i)) then
           if (.not.dwnflg(i)) then
              hcko1(k,i) = hcko(k,i)
           else
              hcko1(k,i) = hckod(k,i)
           endif
        elseif (k > ktcon(i)) then
           hcko1(k,i)  = hcko(k,i)
        endif
        if (k > kb(i)) then
           dbyo(k,i) = hcko1(k,i) - hesol(k,i)
        endif
     end do
  end do


! Compute cloud moisture property and precipitation
  do i = 1, im
     indx = kb(i)
     qcko(indx,i) = qol(indx,i)
  enddo
  do k = 1,km
     do i = 1,im
        qcko0(k,i) = qcko(k,i)
     end do
  end do
  do i = 1,im
     do k = 1, kmax(i)
        if (k > kb(i) .and. k < ktcon(i)) then
           factor = eta(k-1,i) / eta(k,i)
           onemf  = one - factor
           temp   = factor * qcko(k-1,i) + onemf * &
                half * (qol(k,i) + qol(k+1,i))
           qcko0(k,i) = temp
           gamma  = el2orc * qesol(k,i) / (tol(k,i)**2)
           qrch   = qesol(k,i) &
                + gamma * dbyo(k,i) / (hvap * (one + gamma))
           dq = eta(k,i) * (temp - qrch)
           dqk(k,i) = dq

!          Below LFC check if there is excess moisture to release latent heat
           if (dq > zero) then
              dz   = half * (zo(k+1,i) - zo(k-1,i))
              dz1  = (zo(k,i) - zo(k-1,i))
              etah = half * (eta(k,i) + eta(k-1,i))
              qlk  = dq / (eta(k,i) + etah * c0 * dz)
              aa1(i)    = aa1(i) - dz1 * grav * qlk
              pwo(k,i)  = etah * c0 * dz * qlk
              pwavo(i)  = pwavo(i) + pwo(k,i)
              qc        = qlk + qrch
              qcko(k,i) = qc
           else
              qcko(k,i) = temp   
           endif
        endif
     enddo
  enddo

! This section is for cloud liquid water
  if (ncloud > 0) then
     do i=1,im
        k = ktcon(i)
        if (cnvflg(i)) then
           gamma = el2orc * qesol(k,i) / (tol(k,i)**2)
           qrch = qesol(k,i) &
                + gamma * dbyo(k,i) / (hvap * (one + gamma))
           dq = qcko(k-1,i) - qrch
           qcko00(k-1,i) = qcko(k-1,i)

!          Check if there is excess moisture to release latent heat
           if(dq > zero) then
              qlko_ktcon(i) = dq
              qcko(k-1,i) = qrch
           endif
        endif
     end do
  endif
               

! Calculate cloud work function at t+dt
  do i = 1,im
     do k = 1, kmax(i)
        if (k > kbcon(i) .and. k <= ktcon(i)) then
           dz1    = zo(k,i) - zo(k-1,i)
           gamma  = el2orc * qesol(k-1,i) / (tol(k-1,i)**2)
           rfact  =  one + delta * cp * gamma &
                * tol(k-1,i) / hvap
           term1 = dz1*grav/cp
           term2 = one/tol(k-1,i)
           term3 = dbyo(k-1,i)
           term4 = one/(one+gamma)
           term5 = rfact
           aa1(i) = aa1(i) + term1*term2*term3*term4*term5
           
           dqs    = qesol(k-1,i) - qol(k-1,i)
           if (dqs > zero) then
              aa1(i) = aa1(i)+ dz1 * grav * delta * dqs
           endif
        endif
     enddo
  enddo

! Set cloud mass flux threshold function based on updaft cwf.
  do i = 1,im
     term    = 1.e3_r_kind*aa1(i)
     aa10(i) = aa1(i)
     if (abs(term)  <  100._r_kind) then
        tcwfup(i) = ftanh(term)
     else
        if (term  <  -100._r_kind) then
           tcwfup(i) = zero
        else
           tcwfup(i) = one
        endif
     endif
  end do


! ----Downdraft calculations follow----

! Determine downdraft strength in terms of windshear
  do i = 1,im
     do k = 1, kmax(i)
        if (k >= kb(i) .and. k <= ktcon(i)) then
           termu = (uol(k+1,i) - uol(k,i))**2
           termv = (vol(k+1,i) - vol(k,i))**2
           if (termu+termv > tiny_r_kind) then
              shear = rcs(i) * sqrt(termu+termv)
           else
              shear = zero
           endif
           vshear(i)= vshear(i) + shear
        endif
     enddo
  enddo
  do i = 1, im
     dz = zo(ktcon(i),i)-zo(kb(i),i)
     if (abs(dz) > tiny_r_kind) then
        term   = 1.e3_r_kind * vshear(i) / dz
        e1     = pcpeff0 + pcpeff1*term + pcpeff2*term**2 + pcpeff3*term**3
        edt(i) = one_tenth
     else
        edt(i) = zero
     endif
     edt0(i) = edt(i)
  end do

! Impose upper and lower bounds on precipitation efficiency
  do i = 1,im
     if (edt(i) > 0.9_r_kind) then
        edt(i) = 0.9_r_kind
     elseif (edt(i) < zero) then
        edt(i) = zero
     endif
     edto(i) = edt(i)
     edtx(i) = edt(i)
  enddo

! Determine detrainment rate between 1 and kbdtr
  do i = 1, im
     kbdtr(i) = kbcon(i)
     beta = betas
     if (nint(slimsk(i)) == 1) beta = betal
     kbdtr(i) = kbcon(i)
     kbdtr(i) = max(kbdtr(i),1)
     xlamd(i) = zero
     if (kbdtr(i) > 1) then
        dz = half * zo(kbdtr(i),i) + half * zo(kbdtr(i)-1,i) &
             - zo(1,i)
        xlamd(i) = log(beta) / dz
     endif
  enddo

! Determine downdraft mass flux
  do i = 1,im
     do k = 1, kmax(i)
        etad(k,i) = one
     enddo
  enddo
  do i = 1,im
     do k = kbmax(i), 2, -1
        if (k < kbdtr(i)) then
           dz = half * (zo(k+1,i) - zo(k-1,i))
           etad(k,i) = etad(k+1,i) * exp(xlamd(i) * dz)
        endif
     enddo
  enddo
  k = 1
  do i = 1, im
     if (kbdtr(i) > 1) then
        dz = half * (zo(2,i) - zo(1,i))
        etad(k,i) = etad(k+1,i) * exp(xlamd(i) * dz)
     endif
  enddo

! Downdraft moisture properties
  do i = 1,im
     do k = 1,km
        qrcdo(k,i) = zero
        pwdo(k,i) = zero
     end do
  end do
  do i = 1, im
     pwevo(i) = zero
     jmn     = jmin(i)
     hcdo(i) = heol(jmn,i)
     qrcdo(jmn,i) = qesol(jmn,i)
     ucdo(i) = uol(jmn,i)
     vcdo(i) = vol(jmn,i)
  enddo
  do i = 1,im
     do k = kmax(i)-1, 1, -1
        if (k < jmin(i)) then
           dq    = qesol(k,i)
           dt    = tol(k,i)
           gamma = el2orc * dq / dt**2
           dh    = hcdo(i) - hesol(k,i)
           term5 = (one/hvap) * gamma
           term6 = one/(one+gamma)
           qrcdo(k,i) = dq + term5*term6*dh
           detad = etad(k+1,i) - etad(k,i)
           if (k < jmin(i)-1) then
              term1 = etad(k+1,i)*qcdo(i)
              term2 = etad(k,i) * qrcdo(k,i)
              term3 = detad*half*(qrcdo(k,i)+qrcdo(k+1,i))
              pwdo(k,i) = term1 - term2 - term3
           endif
           pwevo(i) = pwevo(i) + pwdo(k,i)
        endif
     enddo
  enddo
  do i = 1,im
     qcdo(i) = qrcdo(1,i)
  end do

! Final downdraft strength dependent on precip
! efficiency (edt), normalized condensate (pwav), and
! evaporate (pwev)
  do i = 1, im
     edtmax = edtmaxl
     if (nint(slimsk(i)) == 0) edtmax = edtmaxs
     if (dwnflg2(i)) then
        if (pwevo(i) < zero) then
           edto1(i) = -edto(i) * pwavo(i) / pwevo(i)
           edto10(i) = edto1(i)
           if (edto1(i) > edtmax) then
              edto1(i) = edtmax
           endif
        else
           edto1(i) = zero
        endif
     else
        edto1(i)   = zero
     endif
  enddo


! Downdraft cloudwork functions
!round reverse order of loop --> get slightly different result
  do i = 1,im
     do k = kmax(i)-1, 1, -1
        if (dwnflg2(i) .and. k < jmin(i)) then
           gamma = el2orc * qesol(k+1,i) / tol(k+1,i)**2
           dhh   = hcdo(i)
           dt    = tol(k+1,i)
           dg    = gamma
           dh    = hesol(k+1,i)
           dz    = zo(k,i) - zo(k+1,i)
           
           term1 = edto1(i)
           term2 = dz
           term3 = grav/(cp*dt)
           term4 = dhh-dh
           term5 = one/(one+dg)
           term6 = one + delta*cp*dg*dt/hvap
           aa1(i) = aa1(i) &
                + term1*term2*term3*term4*term5*term6
           
           dqs   = qesol(k+1,i)-qol(k+1,i)
           if (dqs > zero) then
              aa1(i)= aa1(i) + edto1(i)*dz*grav*delta*dqs
           endif
        endif
     enddo
  enddo

! Compute threshold function based on final cloud work function
  do i = 1,im
     term = 1.e3_r_kind*aa1(i)
     if (abs(term)  <  100._r_kind) then
        tcwfdn(i) = ftanh(term)
     else
        if (term  <  -100._r_kind) then
           tcwfdn(i) = zero
        else
           tcwfdn(i) = one
        endif
     endif
  end do


! What would the change be, that a cloud with unit mass
! will do to the environment?
  do i = 1,im
     do k = 1, kmax(i)
        dellah(k,i) = zero
        dellaq(k,i) = zero
        dellau(k,i) = zero
        dellav(k,i) = zero
     enddo
  enddo
  do i = 1, im
     if (ktcon(i)/=1) then
        dp = 100.0_r_kind * psfc(i) * del(1,i)
        dellah(1,i) = edto1(i) * etad(1,i) * (hcdo(i) &
             - heol(1,i)) * grav / dp
        dellaq(1,i) = edto1(i) * etad(1,i) * (qcdo(i) &
             - qol(1,i)) * grav / dp
        dellau(1,i) = edto1(i) * etad(1,i) * (ucdo(i) &
             - uol(1,i))  * grav / dp
        dellav(1,i) = edto1(i) * etad(1,i) * (vcdo(i) &
             - vol(1,i))  * grav / dp
     endif
  enddo

! Changed due to subsidence and entrainment
!round reverse order of loop yields different results
  do i = 1,im
     do k = 2, kmax(i)-1
        if ( k < ktcon(i) .and. ktcon(i)/=1 ) then
           aup = one
           if (k <= kb(i)) aup = zero
           adw = one
           if (k > jmin(i)) adw = zero
           dv1   = heol(k,i)
           dv2   = half * (heol(k,i) + heol(k+1,i))
           dv3   = heol(k-1,i)
           dv1q  = qol(k,i)
           dv2q  = half * (qol(k,i) + qol(k+1,i))
           dv3q  = qol(k-1,i)
           
           dv1u  = uol(k,i)
           dv2u  = half * (uol(k,i) + uol(k+1,i))
           dv3u  = uol(k-1,i)
           
           dv1v  = vol(k,i)
           dv2v  = half * (vol(k,i) + vol(k+1,i))
           dv3v  = vol(k-1,i)
           
           
           dp    = 100.0_r_kind * psfc(i) * del(k,i)
           detau = eta(k,i) - eta(k-1,i)
           detad = etad(k,i) - etad(k-1,i)
           
           term1 = aup*eta(k,i) - adw*edto1(i)*etad(k,i)
           term2 = aup*detau
           term3 = aup*eta(k-1,i) - adw*edto1(i)*etad(k-1,i)
           term4 = adw*edto1(i)*detad
           termq = half*(qrcdo(k,i)+qrcdo(k-1,i))
           
           dellah(k,i) = dellah(k,i) + (grav/dp) * &
                (term1*dv1 - term3*dv3 - term2*dv2 + term4*hcdo(i))
           dellaq(k,i) = dellaq(k,i) + (grav/dp) * &
                ( term1*dv1q - term3*dv3q - term2*dv2q + term4*termq )
           
           dellau(k,i) = dellau(k,i) + (grav/dp) * &
                (term1*dv1u - term3*dv3u - term2*dv2u + term4*ucdo(i))
           dellav(k,i) = dellav(k,i) + (grav/dp) * &
                (term1*dv1v - term3*dv3v - term2*dv2v + term4*vcdo(i))
           
        endif
     enddo
  enddo

! Cloud top
  do i = 1, im
     if (ktcon(i)/=1) then
        indx = ktcon(i)
        dp   = 100.0_r_kind * psfc(i) * del(indx,i)
        dv1  = heol(indx-1,i)
        dellah(indx,i) = eta(indx-1,i) * &
             (hcko1(indx-1,i) - dv1) * grav / dp
        dvq1 = qol(indx-1,i)
        dellaq(indx,i) = eta(indx-1,i) * &
             (qcko(indx-1,i) - dvq1) * grav / dp
        
        dvu1 = uol(indx-1,i)
        dellau(indx,i) = eta(indx-1,i) * &
             (ucko(indx-1,i) - dvu1) * grav / dp
        dvv1 = vol(indx-1,i)
        dellav(indx,i) = eta(indx-1,i) * &
             (vcko(indx-1,i) - dvv1) * grav / dp
!cloud liquid water
        dellal(i) = eta(indx-1,i) * qlko_ktcon(i) *grav/dp
        
     endif
  enddo

! Final changed variable per unit mass flux
  do i = 1,im
     do k = 1, kmax(i)
        if ( k <= ktcon(i) .and. ktcon(i)/=1 ) then
           xqo0(k,i) = dellaq(k,i) * mbdt + qo0(k,i)
           if (xqo0(k,i)>zero) then
              xqo(k,i) = xqo0(k,i)
           else
              xqo(k,i) = tiny_r_kind
           endif
           dellat   = (dellah(k,i) - hvap * dellaq(k,i)) / cp
           xto(k,i) = dellat * mbdt + to(k,i)
        else
           xqo(k,i) = qo0(k,i)
           xto(k,i) = to(k,i)
        endif
     enddo
  enddo


! The above changed environment is now used to calulate the
! effect the arbitrary cloud (with unit mass flux) would have
! on the stability, which then is used to calculate the real 
! mass flux, necessary to keep this change in balance with !
! the large-scale destabilization.
!
! Environmental conditions again, first heights
  do i = 1,im
     do k = 1, kmax(i)
        call fpvsx_ad(xto(k,i),es0,adt,es0_ad,.false.)
        es = 10.0_r_kind*es0
        xqeso(k,i) = eps*es / (p(k,i) + epsm1*es)
        xtvo(k,i)  = xto(k,i) + delta * xto(k,i) * xqo(k,i)
     enddo
  enddo

! Hydrostatic height assume zero terr
  do i = 1, im
     dlnsig = log(sl(1,i))
     xzo(1,i) = zero - dlnsig * rd / grav * xtvo(1,i)
  enddo
  do i = 1,im
     do k = 2, kmax(i)
        dlnsig = log(sl(k,i) / sl(k-1,i))
        term1 = dlnsig * rd / grav
        term2 = half * (xtvo(k,i) + xtvo(k-1,i))
        xzo(k,i) = xzo(k-1,i) - term1*term2
     enddo
  enddo

! Moist static energy
  do i = 1,im
     do k = 1, kmax(i) - 1
        dz     = half * (xzo(k+1,i) - xzo(k,i))
        dp     = half * (p(k+1,i) - p(k,i))
        call fpvsx_ad(xto(k+1,i),es0,adt,es0_ad,.false.)
        es     = 10.0_r_kind*es0
        pprime = p(k+1,i) + epsm1 * es
        qs     = eps * es / pprime
        dqsdp  = - qs / pprime
        desdt = es * (factor1 / xto(k+1,i) + factor2 / (xto(k+1,i)**2))
        dqsdt = qs * p(k+1,i) * desdt / (es * pprime)
        gamma = el2orc * xqeso(k+1,i) / (xto(k+1,i)**2)
        dt    = (grav * dz + hvap * dqsdp * dp) / (cp * (one + gamma))
        dq    = dqsdt * dt + dqsdp * dp
        xtol(k,i) = xto(k+1,i) + dt
        xqol0(k,i) = xqo(k+1,i) + dq
        if (xqol0(k,i)>zero) then
           xqol(k,i) = xqol0(k,i)
        else
           xqol(k,i) = tiny_r_kind
        endif
     end do
  end do
  do i = 1,im
     do k = 1,kmax(i)-1
        po    = half * (p(k,i) + p(k+1,i))
        call fpvsx_ad(xtol(k,i),es0,adt,es0_ad,.false.)
        esl   = 10.0_r_kind*es0
        xqesol(k,i) = eps*esl / (po + epsm1*esl)
        xheo(k,i)   = half * grav * (xzo(k,i) + xzo(k+1,i)) + &
             cp * xtol(k,i) + hvap * xqol(k,i)
        xheso(k,i)  = half * grav * (xzo(k,i) + xzo(k+1,i)) + &
             cp * xtol(k,i) + hvap * xqesol(k,i)
     enddo
  enddo
  do i = 1, im
     k = kmax(i)
     xheo(k,i)  = grav*xzo(k,i) + cp*xtol(k,i) + hvap*xqol(k,i)
     xheso(k,i) = grav*xzo(k,i) + cp*xtol(k,i) + hvap*xqesol(k,i)
  enddo
  do i = 1, im
     indx = kb(i)
     xhcko(indx,i) = xheo(indx,i)
     xqcko(indx,i) = xqol(indx,i)
  enddo


! Static control

! Moisture and cloud work functions
  do i = 1,im
     do k = 2, kmax(i) - 1
        if (k > kb(i) .and. k <= ktcon(i)) then
           factor = eta(k-1,i) / eta(k,i)
           onemf  = one - factor
           xhcko(k,i) = factor * xhcko(k-1,i) + onemf * &
                half * (xheo(k,i) + xheo(k+1,i))
        endif
     enddo
  enddo

  do i = 1,im
     xaa0(i) = zero
     xpwav(i) = zero
  end do
  do i = 1,im
     do k = 2,kmax(i)-1
        if (k > kb(i) .and. k < ktcon(i)) then

!          xtemp calculation
           factor = eta(k-1,i)/eta(k,i)
           onemf  = one - factor
           termq  = half*(xqol(k,i) + xqol(k+1,i))
           xtemp  = factor*xqcko(k-1,i) + onemf*termq
           xtemp0(k,i) = xtemp

!          xqrch calculation
           xdby = xhcko(k,i) - xheso(k,i)
           if (xdby > zero) then
              gamma = el2orc*xqesol(k,i)/xtol(k,i)**2
              xqrch = xqesol(k,i) + gamma*xdby/(hvap*(one+gamma))
           else
              xqrch = xqesol(k,i)
           endif
           xqrch0(k,i) = xqrch

!          dq calculation
           dq = eta(k,i)*(xtemp-xqrch)

!          actions dependent upon sign of dq
           if (dq > zero) then
              dz   = half * (xzo(k+1,i) - xzo(k-1,i))
              dz1  = xzo(k,i) - xzo(k-1,i)
              etah = half * (eta(k,i) + eta(k-1,i))
              qlk  = dq / (eta(k,i) + etah * c0 * dz)
              xqc  = qlk + xqrch
              xqcko(k,i) = xqc
              xaa0(i)  = xaa0(i) - dz1*grav*qlk
              xpw      = etah * c0 * dz * qlk
              xpwav(i) = xpwav(i) + xpw
           else
              xqcko(k,i) = xtemp
           endif
!
        endif
     end do
  end do
!
  do i = 1,im
     do k = 2,kmax(i)-1
        if (k > kbcon(i) .and. k <= ktcon(i)) then
           dz1   = xzo(k,i) - xzo(k-1,i)
           gamma = el2orc * xqesol(k-1,i) / (xtol(k-1,i)**2)
           rfact =  one + delta * cp * gamma &
                * xtol(k-1,i) / hvap
           xdby  = xhcko(k-1,i) - xheso(k-1,i)
           term1 = dz1 * (grav/cp)
           term2 = one/xtol(k-1,i)
           term3 = xdby
           term4 = one/(one+gamma)
           term5 = rfact
           xaa0(i) = xaa0(i) &
                + term1*term2*term3*term4*term5
           dqs = xqesol(k-1,i) - xqol(k-1,i)
           if (dqs > zero) then
              xaa0(i) = xaa0(i)+ dz1 * grav * delta * dqs
           endif
        endif
     enddo
  enddo

! Downdraft calculations
! Downdraft moisture properties
  do i = 1, im
     xpwev(i) = zero
  enddo
  do i = 1, im
     if (dwnflg2(i)) then
        jmn     = jmin(i)
        xhcd(i) = xheo(jmn,i)
        xqrcd(jmn,i) = xqesol(jmn,i)
     endif
  enddo
  do i = 1,im
     do k = kmax(i)-1, 1, -1
        if (dwnflg2(i) .and. k < jmin(i)) then
           dq = xqesol(k,i)
           dt = xtol(k,i)
           gamma = el2orc*dq/dt**2
           dh = xhcd(i) - xheso(k,i)
           term5 = (one/hvap) * gamma
           term6 = one/(one+gamma)
           
           xqrcd(k,i) = dq + term5*term6*dh
           
           detad = etad(k+1,i) - etad(k,i)
           term1 = etad(k+1,i)*xqrcd(k+1,i)
           term2 = etad(k,i)*xqrcd(k,i)
           term3 = detad * half * (xqrcd(k,i) + xqrcd(k+1,i))
           xpwd = term1 - term2 - term3
           xpwev(i) = xpwev(i) + xpwd
        endif
     enddo
  enddo

  do i = 1, im
     edtmax = edtmaxl
     if (nint(slimsk(i)) == 0) edtmax = edtmaxs
     if (dwnflg2(i)) then
        if (xpwev(i) >= zero) then
           edtx1(i) = zero
        else
           edtx1(i) = -edtx(i) * xpwav(i) / xpwev(i)
           edtx10(i) = edtx1(i)
           if (edtx1(i) > edtmax) then
              edtx1(i) = edtmax
           endif
        endif
     else
        edtx1(i) = zero
     endif
  enddo

! Downdraft cloudwork functions
!round reverse order of loop yields slightly difference results
!      difference only in q and rn, not t
  do i = 1,im
     do k = kmax(i)-1, 1, -1
        if (dwnflg2(i) .and. k < jmin(i)) then
           gamma   = el2orc * xqesol(k+1,i) / xtol(k+1,i)**2
           dhh     = xhcd(i)
           dt      = xtol(k+1,i)
           dg      = gamma
           dh      = xheso(k+1,i)
           dz      = xzo(k,i) - xzo(k+1,i)
           
           term1 = edtx1(i)
           term2 = dz
           term3 = grav/(cp*dt)
           term4 = dhh-dh
           term5 = one/(one+dg)
           term6 = one + delta*cp*dg*dt/hvap
           dqs = xqesol(k+1,i)-xqol(k+1,i)
           
           xaa0(i) = xaa0(i) &
                + term1*term2*term3*term4*term5*term6
           
           if (dqs > zero) then
              xaa0(i) = xaa0(i) + edtx1(i)*dz*grav*delta*dqs
           endif
        endif
     enddo
  enddo

! Calculate critical cloud work function
  do i = 1, im
     acrt(i) = zero
     if (p(ktcon(i),i) < pcrit(15))then
        acrt(i) = acrit(15)*(975.0_r_kind-p(ktcon(i),i)) &
             /(975.0_r_kind-pcrit(15))
     elseif (p(ktcon(i),i) > pcrit(1))then
        acrt(i) = acrit(1)
     else
        k =  int((850.0_r_kind - p(ktcon(i),i))/50.0_r_kind) + 2
        k = min(k,15)
        k = max(k,2)
        acrt(i) = acrit(k)+(acrit(k-1)-acrit(k))* &
             (p(ktcon(i),i)-pcrit(k))/(pcrit(k-1)-pcrit(k))
     endif
  enddo
  do i = 1, im
     acrtfct(i) = one
     if (nint(slimsk(i)) == 1) then
        w1 = w1l
        w2 = w2l
        w3 = w3l
        w4 = w4l
     else
        w1 = w1s
        w2 = w2s
        w3 = w3s
        w4 = w4s
     endif
     if (pdot(i) <= w4) then
        acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
     elseif (pdot(i) >= -w4) then
        acrtfct(i) = -one*(pdot(i) + w4) / (w4 - w3)
     else
        acrtfct(i) = zero
     endif
     acrtfct0(i) = acrtfct(i)

     if (acrtfct(i) < -one) then
        acrtfct(i) = -one
     endif
     if (acrtfct(i) > one) then
        acrtfct(i) = one
     endif

     acrtfct1(i) = one - acrtfct(i)

     dtconv(i)  = dt2 + (1800.0_r_kind - dt2) * &
          (pdot(i) - w2) / (w1 - w2)
     dtconv0(i) = dtconv(i)
     if (dtconv(i) < dtmin) then
        dtconv(i) = dtmin
     endif
     if (dtconv(i) > dtmax) then
        dtconv(i) = dtmax
     endif
  enddo

! Large scale forcing
  do i= 1, im
     f(i)  = (aa1(i) - acrt(i) * acrtfct1(i)) / dtconv(i)
     xk(i) = (xaa0(i) - aa1(i)) / mbdt

!    Kernel, cloud base mass flux
     if (abs(xk(i)) > tiny_r_kind) then
        ratio = -f(i)/xk(i)
     else
        ratio = zero
     endif
     if ( (f(i) > zero) .and. (xk(i) < zero) ) then
        fixed = tdpscl(i)*tdpcld(i)*tdetrn(i)*tdpdwn(i)
        xmb(i)  = fixed * ratio*tcwfup(i)*tcwfdn(i)
     else
        xmb(i) = zero
     endif
     xmb0(i) = xmb(i)
     if (xmb(i) > xmbmax(i)) then
        xmb(i) = xmbmax(i)
     endif
  enddo


! Feedback: simply the changes from the cloud with unit mass flux
!           multiplied by  the mass flux necessary to keep the
!           equilibrium with the larger-scale.
!round - reverse order of loop, different q and rn
  do i = 1,im
     do k = 1,km
        dellal1(i)   = zero
        dellalk(k,i) = zero
        dellalk(k,i) = dellal(i)
        if (k <= ktcon(i) .and. ktcon(i)/=1) then
           dellat     = (dellah(k,i)-hvap*dellaq(k,i))/cp
           to2(k,i)   = to(k,i) + dellat*xmb(i)*dt2
           qo2(k,i)   = qo0(k,i) + dellaq(k,i)*xmb(i)*dt2
           uo2(k,i)   = uo(k,i) + dellau(k,i)*xmb(i)*dt2
           vo2(k,i)   = vo(k,i) + dellav(k,i)*xmb(i)*dt2
           call fpvsx_ad(to2(k,i),es0,adt,es0_ad,.false.)
           es         = 10.0_r_kind*es0
           qeso2(k,i) = eps*es / (p(k,i)+epsm1*es)
           if (k==ktcon(i) .and. ncloud > 0 .and. cnvflg(i)) then
              dp = 100.0_r_kind * psfc(i) * del(k,i)
              cwmo2(k,i) = cwmo(k,i) + dellalk(k,i) * xmb(i) * dt2
              dellal1(i) = dellalk(k,i)*xmb(i)*dp/grav
           else
              cwmo2(k,i) = cwmo(k,i)
              dellal1(i) = dellalk(k,i)
           endif
        else
           to2(k,i) = to(k,i)
           qo2(k,i) = qo(k,i)
           uo2(k,i) = uo(k,i)
           vo2(k,i) = vo(k,i)
           cwmo2(k,i)= cwmo(k,i)
           call fpvsx_ad(to2(k,i),es0,adt,es0_ad,.false.)
           es         = 10._r_kind*es0
           qeso2(k,i) = eps*es / (p(k,i)+epsm1*es)
           dellal1(i) = dellalk(k,i)
        endif
        dellal(i) = dellal1(i)
        dellal1(i)= zero
     enddo
  enddo

  do i = 1, im
     rntot(i)  = zero
  enddo
!round different results when reverse loop
  do i = 1,im
     do k = kmax(i), 1, -1
        if (k <= ktcon(i) .and. ktcon(i)/=1) then
           aup = one
           if (k <= kb(i)) aup = zero
           adw = one
           if (k > jmin(i)) adw = zero
           rain =  aup * pwo(k,i) + adw * edto1(i) * pwdo(k,i)
           rntot(i) = rntot(i) + rain * xmb(i) * .001_r_kind * dt2
        endif
     enddo
  enddo

!round very slight difference in q and r
  do i = 1,im
     rn0    = zero
     delqev = zero
     flg(i) = .true.
!
     do k = km,1,-1
        qevap = zero
        qcond = zero
        delq2 = zero

!       Accumulate precipitation from cloud level
        if (k <= ktcon(i) .and. ktcon(i)/=1) then
           aup = one
           if (k <= kb(i)) aup = zero
           adw = one
           if (k > jmin(i)) adw = zero
           rain = aup*pwo(k,i) + adw*edto1(i)*pwdo(k,i)
           rn0  = rn0 + rain*xmb(i)*.001_r_kind*dt2
        endif

!       Save layer specific quantities for use in adjoint
        rn0k(k,i) = rn0
        flgk(k,i) = flg(i)
        delqevk(k,i) = delqev

!       Check if any falling precipitation will evaporate
        if (flg(i) .and. k <= ktcon(i)) then
!               if (nint(slimsk(i)) == 1) then
!                  evef = 0.07_r_kind
!               else
!                  evef = edt(i)*evfact
!               endif
           evef = edt(i)*evfact
           term1 = evef
           term2 = qo2(k,i) - qeso2(k,i)
           term3 = el2orc*qeso2(k,i)/to2(k,i)**2
           term4 = one / (one+term3)
           qcond = term1*term2*term4
           dp = 100.0_r_kind * psfc(i) * del(k,i)
!
!          Compute amount of precipitation subject to evaporation
           if (rn0  >  zero .and. qcond <  zero) then
              qevap = -qcond * (one-exp(-.32_r_kind*sqrt(dt2*rn0)))
              
              if ( qevap  >  rn0*r1000*grav/dp ) then
                 qevap = rn0*r1000*grav/dp
              endif
              
              delq2 = delqev + .001_r_kind * qevap * dp / grav
              if (delq2 > rntot(i)) then
                 qevap  = r1000* grav * (rntot(i)-delqev) / dp
                 flg(i) = .false.
              endif
           endif

!          Adjust t and q profiles plus net rain in response
!          to evaporation of falling precipitation.
           if (rn0  >  zero .and. qevap >  zero) then
              qo3(k,i) = qo2(k,i) + qevap
              to3(k,i) = to2(k,i) - elocp * qevap
              rn0      = rn0  - .001_r_kind * qevap * dp / grav
              
              delqev   = delqev + .001_r_kind*dp*qevap/grav

!          Otherwise, no adjustment to t and q profiles
           else
              to3(k,i) = to2(k,i)
              qo3(k,i) = qo2(k,i)
           endif

!       Above cloud top (or not cloud present). 
!       Therefore, no changes to t and q.
        else
           to3(k,i) = to2(k,i)
           qo3(k,i) = qo2(k,i)
        endif

     enddo

!    Save net surface precipitation.
     rn1(i) = rn0
  enddo


! Load output arrays.
  do i = 1,im

!    If total precipitation is negative, restore initial profiles.
     rn11(i) = rn1(i)
     if (rn1(i) < zero) then
        rn1(i) = zero
        do k = 1,km
           t1(k,i) = to(k,i)
           q1(k,i) = qo0(k,i)
        end do

!    Otherwise, load modified t,q profile to output arrays
     else
        do k = 1,km
           t1(k,i) = to3(k,i)
           q1(k,i) = qo3(k,i)
        end do
     endif

!    Load various diagostic output arrays
     ktop(i)   = ktcon(i)
     kbot(i)   = kbcon(i)
     kuo(i)    = 1
     cldwrk(i) = aa1(i)
  end do


! Transfer input u,v,cwm to ouptut arrays
  do i = 1,im
     do k=1,km
        u1(k,i)  = uo2(k,i)
        v1(k,i)  = vo2(k,i)
        cwm1(k,i)= cwmo2(k,i)
     end do
  end do
  
  if (.not.adjoint) return

!
!
!
!=====================================================================
!
!
! Begin adjoint code here

! Adjoint of u,v,cwm transfer
  do i = 1,im
     do k=1,km
        uo2_ad(k,i)   = uo2_ad(k,i) + u1_ad(k,i)
        vo2_ad(k,i)   = vo2_ad(k,i) + v1_ad(k,i)
        cwmo2_ad(k,i) = cwmo2_ad(k,i) + cwm1_ad(k,i)
        u1_ad(k,i) = zero
        v1_ad(k,i) = zero
        cwm1_ad(k,i) = zero
     end do
  end do

!
! Adjoint of loading of output arrays
  do i = 1,im

!    zero or negative precip case.
     if (rn11(i) < zero) then
        rn1_ad(i) = zero
        do k = km,1,-1
           qo0_ad(k,i)= qo0_ad(k,i)+ q1_ad(k,i)
           to_ad(k,i) = to_ad(k,i) + t1_ad(k,i)
           q1_ad(k,i) = zero
           t1_ad(k,i) = zero
        end do

!    precipitation greater than zero case
     else
        do k = km,1,-1
           qo3_ad(k,i) = qo3_ad(k,i) + q1_ad(k,i)
           to3_ad(k,i) = to3_ad(k,i) + t1_ad(k,i)
           q1_ad(k,i) = zero
           t1_ad(k,i) = zero
        end do
     endif
  end do


! Adjoint of final precipitation calculation (including
! evaporation of falling precipitation)
  do i = 1,im

!    Zero local adjoint variables
     rn0_ad    = zero
     delqev_ad = zero

!    Adjoint of rn0 to rn1 transfer
     rn0_ad = rn0_ad + rn1_ad(i)
!
     do k = 1,km
!
!       Restore saved values from nlm
        rn0    = rn0k(k,i)
        flg(i) = flgk(k,i)
        delqev = delqevk(k,i)
        qevap  = zero
        qcond  = zero
        delq2  = zero

!       Zero local adjoint variables
        qevap_ad = zero
        qcond_ad = zero
!
!       Adjoint of code for evaporation of falling precipitation
        if (flg(i) .and. k <= ktcon(i)) then
!
!           Redo nlm
!           if (nint(slimsk(i)) == 1) then
!              evef = 0.07_r_kind
!           else
!              evef = edt(i)*evfact
!           endif
           evef = edt(i)*evfact
           term1 = evef
           term2 = qo2(k,i) - qeso2(k,i)
           term3 = el2orc*qeso2(k,i)/to2(k,i)**2
           term4 = one / (one+term3)
           qcond = term1*term2*term4
           dp = 100._r_kind * psfc(i) * del(k,i)

           if (rn0 >  zero .and. qcond < zero) then
              qevap  = -qcond * (one-exp(-.32_r_kind*sqrt(dt2*rn0)))
              qevap0 = qevap
              if ( qevap  >  (rn0*r1000*grav/dp) ) then
                 qevap = rn0*r1000*grav/dp
              endif
              delq2 = delqev + 0.001_r_kind * qevap * dp / grav
              qevap1 = qevap
              if (delq2 > rntot(i)) then
                 qevap  = r1000* grav * (rntot(i)-delqev) / dp
              endif
           endif
!
!          Adjoint of t,q,rn0 adjustments due to evaporation
           if ( rn0 > zero .and. qevap > zero ) then
              qevap_ad = qevap_ad + 0.001_r_kind*delqev_ad*dp/grav

              qevap_ad = qevap_ad - 0.001_r_kind*rn0_ad*dp/grav

              qevap_ad = qevap_ad - elocp*to3_ad(k,i)
              to2_ad(k,i) = to2_ad(k,i) + to3_ad(k,i)

              qevap_ad = qevap_ad + qo3_ad(k,i)
              qo2_ad(k,i) = qo2_ad(k,i) + qo3_ad(k,i)
!
!             Otherwise adjoint of no adjustment case
           else
              qo2_ad(k,i) = qo2_ad(k,i) + qo3_ad(k,i)
              to2_ad(k,i) = to2_ad(k,i) + to3_ad(k,i)
           endif
!
!          Adjoint of qevap calculation
           if (rn0 > zero .and. qcond < zero) then
              if (delq2 > rntot(i)) then
                 delqev_ad = delqev_ad - r1000*grav*qevap_ad/dp
                 rntot_ad(i) = rntot_ad(i) + r1000*grav*qevap_ad/dp
              endif
                  
              if (qevap0  >  rn0*r1000*grav/dp) then
                 rn0_ad = rn0_ad + qevap_ad*r1000*grav/dp
              endif

              rn0_ad = rn0_ad &
                   + qcond * (exp(-0.32_r_kind*sqrt(dt2*rn0)) * &
                   (-0.32_r_kind*half/sqrt(dt2*rn0)) * dt2*qevap_ad)
              qcond_ad = qcond_ad &
                   - qevap_ad * (one-exp(-0.32_r_kind*sqrt(dt2*rn0)))
           endif
!     
!          Adjoint of qcond calculation
           term4_ad = zero
           term3_ad = zero
           term2_ad = zero
           term1_ad = zero
           evef_ad  = zero

           term4_ad = term4_ad + term1*term2*qcond_ad
           term2_ad = term2_ad + term1*term4*qcond_ad
           term1_ad = term1_ad + term2*term4*qcond_ad

           term3_ad = term3_ad - one/(one+term3)**2 * term4_ad
               
           to2_ad(k,i) = to2_ad(k,i) &
                - two*el2orc*qeso2(k,i)/to2(k,i)**3 * term3_ad
           qeso2_ad(k,i) = qeso2_ad(k,i) &
                + el2orc*term3_ad/to2(k,i)**2

           qeso2_ad(k,i) = qeso2_ad(k,i) - term2_ad
           qo2_ad(k,i) = qo2_ad(k,i) + term2_ad

           evef_ad = evef_ad + term1_ad
!           if (nint(slimsk(i)) == 1) then
!              evef_ad = zero
!           else
!              edt_ad(i) = edt_ad(i) + evef_ad*evfact
!           endif
              edt_ad(i) = edt_ad(i) + evef_ad*evfact
!
!              
!       Adjoint of "2" to "3" transfer.  This is the case
!       when above cloud top or no cloud present.
        else
           to2_ad(k,i) = to2_ad(k,i) + to3_ad(k,i)
           qo2_ad(k,i) = qo2_ad(k,i) + qo3_ad(k,i)
        endif
!
!       Adjoint of layer rain calculation
        if (k <= ktcon(i) .and. ktcon(i)/=1) then
!
!          Redo nlm
           aup = one
           if (k <= kb(i)) aup = zero
           adw = one
           if (k > jmin(i)) adw = zero
           rain = aup*pwo(k,i) + adw*edto1(i)*pwdo(k,i)
!
!          Adjoint of tlm code.
           rain_ad   = zero
           xmb_ad(i) = xmb_ad(i) + rain*rn0_ad*0.001_r_kind*dt2
           rain_ad   = rain_ad + rn0_ad*xmb(i)*0.001_r_kind*dt2

           pwdo_ad(k,i) = pwdo_ad(k,i) + adw*edto1(i)*rain_ad
           edto1_ad(i)   = edto1_ad(i) + adw*rain_ad*pwdo(k,i)
           pwo_ad(k,i)  = pwo_ad(k,i) + aup*rain_ad
        endif
!
!       Zero local adjoint variables
        qevap_ad = zero
        qcond_ad = zero
!
!    End of loop over vertical coordinate
     end do
!
!    Zero local adjoint variables
     rn0_ad    = zero
     delqev_ad = zero
!
! End of loop over profiles
  end do

!
!
! Adjoint of intial total precipitation calculation
  do i = 1,im
     do k = 1,kmax(i)
        if (k <= ktcon(i) .and. ktcon(i)/=1) then
!
!          Redo nlm
           aup = one
           if (k <= kb(i)) aup = zero
           adw = one
           if (k > jmin(i)) adw = zero
           rain =  AUP * PWO(k,i) + ADW * EDTO1(I) * PWDO(k,i)
!
!          Adjoint of tlm code
           rain_ad   = zero
           xmb_ad(i) = xmb_ad(i) + rain*0.001_r_kind*dt2 * rntot_ad(i)
           rain_ad   = rain_ad + 0.001_r_kind*xmb(i)*dt2 * rntot_ad(i)
               
           pwdo_ad(k,i) = pwdo_ad(k,i) + adw*edto1(i)*rain_ad
           edto1_ad(i)   = edto1_ad(i) + adw*pwdo(k,i)*rain_ad
           pwo_ad(k,i)  = pwo_ad(k,i) + aup*rain_ad
        endif
     end do
  end do
!
! Adjoint of to2 and qo2 calculation (prior to evaportation)
  do i = 1,im
     do k = km,1,-1
        dellal1_ad(i) = zero
        dellal1_ad(i) = dellal1_ad(i) + dellal_ad(i)
        if (k <= ktcon(i) .and. ktcon(i)/=1) then
!
!          Redo nlm
           dellat = (dellah(k,i)-hvap*dellaq(k,i))/cp
           call fpvsx_ad(to2(k,i),es0,adt,es0_ad,.false.)
           es = 10.0_r_kind*es0
!
!          Adjoint of tlm code
           dellat_ad = zero
           es_ad     = zero
           es0_ad    = zero

           if (k==ktcon(i) .and. ncloud > 0 .and. cnvflg(i)) then
              dp = 100.0_r_kind*psfc(i)*del(k,i)

              dellalk_ad(k,i) = dellalk_ad(k,i) + &
                   dellal1_ad(i)*xmb(i)*dp/grav
              xmb_ad(i) = xmb_ad(i) + &
                   dellalk(k,i)*dellal1_ad(i)*dp/grav

              dellalk_ad(k,i) = dellalk_ad(k,i) + &
                   cwmo2_ad(k,i)*xmb(i)*dt2
              xmb_ad(i) = xmb_ad(i) + &
                   dellalk(k,i)*cwmo2_ad(k,i)*dt2
              cwmo_ad(k,i)  = cwmo_ad(k,i) + cwmo2_ad(k,i)

              cwmo2_ad(k,i) = zero
           else
              cwmo_ad(k,i) = cwmo_ad(k,i) + cwmo2_ad(k,i)
              dellalk_ad(k,i) = dellalk_ad(k,i) + dellal1_ad(i)

              cwmo2_ad(k,i) = zero
           endif
!
           es_ad = es_ad - eps*es/(p(k,i)+epsm1*es)**2 * &
                epsm1*qeso2_ad(k,i)
           es_ad = es_ad + eps*qeso2_ad(k,i)/(p(k,i)+epsm1*es)

           es0_ad = es0_ad + 10*es_ad

           adt=zero
           call fpvsx_ad(to2(k,i),es0,adt,es0_ad,adjoint)
           to2_ad(k,i) = to2_ad(k,i) + adt
               

           vo_ad(k,i) = vo_ad(k,i) + vo2_ad(k,i)
           xmb_ad(i)  = xmb_ad(i)  + &
                vo2_ad(k,i)*dellav(k,i)*dt2
           dellav_ad(k,i) = dellav_ad(k,i) + &
                vo2_ad(k,i)*xmb(i)*dt2

           uo_ad(k,i) = uo_ad(k,i) + uo2_ad(k,i)
           xmb_ad(i)  = xmb_ad(i)  + &
                uo2_ad(k,i)*dellau(k,i)*dt2
           dellau_ad(k,i) = dellau_ad(k,i) + &
                uo2_ad(k,i)*xmb(i)*dt2

           qo_ad(k,i) = qo_ad(k,i) + qo2_ad(k,i)
           xmb_ad(i)  = xmb_ad(i)  + &
                qo2_ad(k,i)*dellaq(k,i)*dt2
           dellaq_ad(k,i) = dellaq_ad(k,i) + &
                qo2_ad(k,i)*xmb(i)*dt2

           to_ad(k,i) = to_ad(k,i) + to2_ad(k,i)
           xmb_ad(i)  = xmb_ad(i)  + to2_ad(k,i)*dellat*dt2
           dellat_ad  = dellat_ad  + to2_ad(k,i)*xmb(i)*dt2

           dellaq_ad(k,i) = dellaq_ad(k,i) - hvap*dellat_ad/cp
           dellah_ad(k,i) = dellah_ad(k,i) + dellat_ad/cp

        else
!
!          Redo nlm
           call fpvsx_ad(to2(k,i),es0,adt,es0_ad,.false.)
           es = 10.0_r_kind*es0

!
!          Adjoint of tlm code
           es_ad = zero
           es0_ad= zero

           dellalk_ad(k,i) = dellalk_ad(k,i) + dellal1_ad(i)

           es_ad = es_ad - eps*es/(p(k,i)+epsm1*es)**2 * &
                epsm1*qeso2_ad(k,i)
           es_ad = es_ad + eps*qeso2_ad(k,i)/(p(k,i)+epsm1*es)

           es0_ad = es0_ad + 10.0_r_kind*es_ad
           adt = zero
           call fpvsx_ad(to2(k,i),es0,adt,es0_ad,adjoint)
           to2_ad(k,i) = to2_ad(k,i) + adt

           cwmo_ad(k,i) = cwmo_ad(k,i) + cwmo2_ad(k,i)
           vo_ad(k,i)   = vo_ad(k,i) + vo2_ad(k,i)
           uo_ad(k,i)   = uo_ad(k,i) + uo2_ad(k,i)
           qo_ad(k,i)   = qo_ad(k,i) + qo2_ad(k,i)
           to_ad(k,i)   = to_ad(k,i) + to2_ad(k,i)

        endif
        dellal_ad(i)    = dellal_ad(i) + dellalk_ad(k,i)
        dellalk_ad(k,i) = zero
        dellal1_ad(i)   = zero
!
     end do
  end do
!
! Adjoint of kernel (cloud mass flux) calculation
  do i = 1,im
!
!    Redo nlm
     if (abs(xk(i)) > tiny_r_kind) then
        ratio = -f(i)/xk(i)
     else
        ratio = zero
     endif
     fixed = tdpscl(i)*tdpcld(i)*tdetrn(i)*tdpdwn(i)
!
!    Adjoint of tlm
     if (xmb0(i) > xmbmax(i)) then
!       do nothing since tlm = zero
        xmb_ad(i) = zero
     endif

     ratio_ad = zero
     if ( (f(i) > zero) .and. (xk(i) < zero) ) then
        tcwfdn_ad(i) = tcwfdn_ad(i) &
             + fixed*ratio*tcwfup(i)*xmb_ad(i)
        tcwfup_ad(i) = tcwfup_ad(i) &
             + fixed*ratio*xmb_ad(i)*tcwfdn(i)
        ratio_ad = ratio_ad + fixed*xmb_ad(i)*tcwfup(i)*tcwfdn(i)
     else
!       do nothing since tlm = zero
        xmb_ad(i) = zero
     endif

     if (abs(xk(i)) > tiny_r_kind) then
        xk_ad(i) = xk_ad(i) + f(i)/xk(i)**2 * ratio_ad
        f_ad(i)  = f_ad(i) - ratio_ad/xk(i)
     else
!       do nothing since tlm = zero
        ratio_ad = zero
     endif
!
!    Adjoint of xk tlm
     aa1_ad(i)  = aa1_ad(i)  - xk_ad(i)/mbdt
     xaa0_ad(i) = xaa0_ad(i) + xk_ad(i)/mbdt
!
!    Adjoint of f tlm
     dtconv_ad(i) = dtconv_ad(i) &
          - (aa1(i) - acrt(i)*acrtfct1(i))/dtconv(i)**2 &
          * f_ad(i)
     acrtfct1_ad(i) = acrtfct1_ad(i) &
          - acrt(i)*f_ad(i) / dtconv(i)
     aa1_ad(i) = aa1_ad(i) + f_ad(i) / dtconv(i)
!
  end do
!
! Adjoint of acrtfct1 and dtconv calculations
!
  do i = 1,im
!
!    Redo nlm 
     if (nint(slimsk(i)) == 1) then
        w1 = w1l
        w2 = w2l
        w3 = w3l
        w4 = w4l
     else
        w1 = w1s
        w2 = w2s
        w3 = w3s
        w4 = w4s
     endif
!
!    Adjoint of dtconv calculation
     if (dtconv0(i) > dtmax) then
!       do nothing since tlm = zero
        dtconv_ad(i) = zero
     endif
     if (dtconv0(i) < dtmin) then
!       do nothing since tlm = zero
        dtconv_ad(i) = zero
     endif
     pdot_ad(i) = pdot_ad(i) + (1800.0_r_kind-dt2) &
          * dtconv_ad(i)/(w1-w2)
!
!    Adjoint of acrtfct calculation
     acrtfct_ad(i) = acrtfct_ad(i) - acrtfct1_ad(i)
     if (acrtfct0(i) > one) then
!       do nothing since tlm = zero
        acrtfct_ad(i) = zero
     endif
     if (acrtfct0(i) < -one) then
!       do nothing since tlm = zero
        acrtfct_ad(i) = zero
     endif
     if (pdot(i) <= w4) then
        pdot_ad(i) = pdot_ad(i) + acrtfct_ad(i)/(w3-w4)
     elseif (pdot(i) >= -w4) then
        pdot_ad(i) = pdot_ad(i) - acrtfct_ad(i)/(w4-w3)
     else
!       do nothing since tlm = zero
        acrtfct_ad(i) = zero
     endif
!
  end do
!
! Adjoint of xaa0 downdraft cloudwork function calculation
  do i = 1,im
     do k = 1,kmax(i)-1
        if (dwnflg2(i) .and. k < jmin(i)) then
!
!          Redo nlm calculations
           GAMMA = EL2ORC * XQESOL(k+1,i) / XTOL(k+1,i)**2
           DHH   = XHCD(I)
           DT    = XTOL(k+1,i)
           DG    = GAMMA
           DH    = XHESO(k+1,i)
           dz    = xzo(k,i) - xzo(k+1,i)
!
           term1 = edtx1(i)
           term2 = dz
           term3 = grav/(cp*dt)
           term4 = (dhh-dh)
           term5 = one/(one+dg)
           term6 = one + delta*cp*dg*dt/hvap
           dqs = xqesol(k+1,i) - xqol(k+1,i)

!
!          Adjoint of tlm code
           dqs_ad = zero
           dz_ad  = zero
           dt_ad  = zero
           dg_ad  = zero
           dh_ad  = zero
           dhh_ad = zero
           if (dqs > zero) then
              dqs_ad = dqs_ad + edtx1(i)*dz*grav*delta  * xaa0_ad(i)
              dz_ad  = dz_ad  + edtx1(i)*grav*delta*dqs * xaa0_ad(i)
              edtx1_ad(i) = edtx1_ad(i) + dz*grav*delta*dqs &
                   * xaa0_ad(i)
           endif
!     
           term1_ad = term2*term3*term4*term5*term6*xaa0_ad(i)
           term2_ad = term1*term3*term4*term5*term6*xaa0_ad(i)
           term3_ad = term1*term2*term4*term5*term6*xaa0_ad(i)
           term4_ad = term1*term2*term3*term5*term6*xaa0_ad(i)
           term5_ad = term1*term2*term3*term4*term6*xaa0_ad(i)
           term6_ad = term1*term2*term3*term4*term5*xaa0_ad(i)

           xqol_ad(k+1,i)   = xqol_ad(k+1,i)   - dqs_ad
           xqesol_ad(k+1,i) = xqesol_ad(k+1,i) + dqs_ad

!     
           dt_ad  = (delta*cp/hvap) * dg*term6_ad
           dg_ad  = (delta*cp/hvap) * term6_ad*dt
           dg_ad  = dg_ad  - one/(one+dg)**2 * term5_ad
           dh_ad  = dh_ad  - term4_ad
           dhh_ad = dhh_ad + term4_ad
           dt_ad  = dt_ad  - grav/(cp*dt*dt) * term3_ad
           dz_ad  = dz_ad  + term2_ad
           edtx1_ad(i) = edtx1_ad(i) + term1_ad 

           xzo_ad(k,i)    = xzo_ad(k,i)   + dz_ad
           xzo_ad(k+1,i)  = xzo_ad(k+1,i) - dz_ad
           xheso_ad(k+1,i)= xheso_ad(k+1,i) + dh_ad
           gamma_ad       = dg_ad
           xtol_ad(k+1,i) = xtol_ad(k+1,i) + dt_ad
           xhcd_ad(i)     = xhcd_ad(i)   + dhh_ad
           xtol_ad(k+1,i) = xtol_ad(k+1,i) &
                - two*el2orc*xqesol(k+1,i)/xtol(k+1,i)**3 &
                * gamma_ad
           xqesol_ad(k+1,i) = xqesol_ad(k+1,i) + &
                el2orc*gamma_ad/xtol(k+1,i)**2
        endif
     end do
  end do
!
! Adjoint of final edtx1 calculation
  do i = 1,im
     edtmax = edtmaxl
     if (nint(slimsk(i)) == 0) edtmax = edtmaxs
     if (dwnflg2(i)) then
        if (xpwev(i) >= zero) then
!          do nothing since tlm is zero
           edtx1_ad(i) = zero
        else
           if (edtx10(i) > edtmax) then
!             do nothing since tlm is zero
              edtx1_ad(i) = zero
           else
              xpwev_ad(i) = xpwev_ad(i) &
                   + edtx(i)*xpwav(i)/xpwev(i)**2 * edtx1_ad(i)
              xpwav_ad(i) = xpwav_ad(i) &
                   - edtx(i)*edtx1_ad(i)/xpwev(i)
              edtx_ad(i) = edtx_ad(i) &
                   - edtx1_ad(i)*xpwav(i)/xpwev(i)
           endif
        endif
     else
!       do nothing since tlm is zero             
        edtx1_ad(i) = zero
     endif
  enddo
!
! Adjoint of downdraft xpwev calculation
  do i = 1,im
     do k = 1,kmax(i)-1
        if (dwnflg2(i) .and. k < jmin(i)) then
!
!          Redo nlm
           dq    = xqesol(k,i)
           dt    = xtol(k,i)
           gamma = el2orc*dq/dt**2
           dh = xhcd(i) - xheso(k,i)
           term5 = (one/hvap) * gamma
           term6 = one/(one+gamma)
           detad = etad(k+1,i) - etad(k,i)
           term1 = etad(k+1,i)*xqrcd(k+1,i)
           term2 = etad(k,i)*xqrcd(k,i)
           term3 = detad * half * (xqrcd(k,i) + xqrcd(k+1,i))
!
!          Adjoint of tlm code
           xpwd_ad  = zero
           detad_ad = zero
           dh_ad    = zero
           gamma_ad = zero
           dt_ad    = zero
           dq_ad    = zero
           term1_ad = zero
           term2_ad = zero
           term3_ad = zero
           term5_ad = zero
           term6_ad = zero

           xpwd_ad = xpwd_ad + xpwev_ad(i)

           term1_ad = term1_ad + xpwd_ad
           term2_ad = term2_ad - xpwd_ad
           term3_ad = term3_ad - xpwd_ad

           xqrcd_ad(k+1,i) = xqrcd_ad(k+1,i) + etad(k+1,i)*term1_ad
           etad_ad(k+1,i)  = etad_ad(k+1,i)  + term1_ad*xqrcd(k+1,i)

           xqrcd_ad(k,i)   = xqrcd_ad(k,i) + etad(k,i)*term2_ad
           etad_ad(k,i)    = etad_ad(k,i)  + term2_ad*xqrcd(k,i)

           detad_ad        = detad_ad + term3_ad*half*(xqrcd(k,i)+xqrcd(k+1,i))
           xqrcd_ad(k,i)   = xqrcd_ad(k,i) + detad*half*term3_ad
           xqrcd_ad(k+1,i) = xqrcd_ad(k+1,i) + detad*half*term3_ad

           etad_ad(k,i)    = etad_ad(k,i) - detad_ad
           etad_ad(k+1,i)  = etad_ad(k+1,i) + detad_ad

           dh_ad    = dh_ad    + term5*term6*xqrcd_ad(k,i)
           term6_ad = term6_ad + term5*dh*xqrcd_ad(k,i)
           term5_ad = term5_ad + term6*dh*xqrcd_ad(k,i)
           dq_ad = dq_ad + xqrcd_ad(k,i)

           gamma_ad = gamma_ad -one/(one+gamma)**2 * term6_ad
           gamma_ad = gamma_ad + (one/hvap)*term5_ad

           xheso_ad(k,i) = xheso_ad(k,i) - dh_ad
           xhcd_ad(i)  = xhcd_ad(i) + dh_ad

           dt_ad = dt_ad - two*el2orc*dq/dt**3 * gamma_ad
           dq_ad = dq_ad + el2orc*gamma_ad/dt**2

           xtol_ad(k,i) = xtol_ad(k,i) + dt_ad
           xqesol_ad(k,i) = xqesol_ad(k,i) + dq_ad

        endif
     end do
  end do
!
! Adjoint of xheo,xqesol --> xhcd,xqrcd transfers
  do i = 1,im
     if (dwnflg2(i)) then
        jmn = jmin(i)
        xqesol_ad(jmn,i) = xqesol_ad(jmn,i) + xqrcd_ad(jmn,i)
        xheo_ad(jmn,i)   = xheo_ad(jmn,i) + xhcd_ad(i)
     endif
  end do
!
! Adjoint of xaa0 calculations
!
  do i = 1,im
     do k = kmax(i)-1,2,-1
        if (k > kbcon(i).and.k <= ktcon(i)) then
!
!          Recalculate nlm 
           dz1   = xzo(k,i) - xzo(k-1,i)
           gamma = el2orc*xqesol(k-1,i)/xtol(k-1,i)**2
           rfact =  one + delta*cp*gamma*xtol(k-1,i)/hvap
           xdby  = xhcko(k-1,i) - xheso(k-1,i)
           term1 = dz1 * (grav/cp)
           term2 = one/xtol(k-1,i)
           term3 = xdby
           term4 = one/(one+gamma)
           term5 = rfact
           dqs   = xqesol(k-1,i) - xqol(k-1,i)
!
!          Zero local adjoint variables
           dqs_ad = zero
           dz1_ad = zero
           term1_ad = zero
           term2_ad = zero
           term3_ad = zero
           term4_ad = zero
           term5_ad = zero
           rfact_ad = zero
           gamma_ad = zero
           xdby_ad  = zero
!
!          Adjoint of forward model tlm
           if (dqs > zero) then
              dqs_ad   = dqs_ad + grav*delta * dz1*xaa0_ad(i)
              dz1_ad   = dz1_ad + grav*delta * dqs*xaa0_ad(i)
           endif
           xqol_ad(k-1,i)   = xqol_ad(k-1,i)   - dqs_ad
           xqesol_ad(k-1,i) = xqesol_ad(k-1,i) + dqs_ad

           term5_ad = term5_ad + term1*term2*term3*term4*xaa0_ad(i)
           term4_ad = term4_ad + term1*term2*term3*xaa0_ad(i)*term5
           term3_ad = term3_ad + term1*term2*xaa0_ad(i)*term4*term5
           term2_ad = term2_ad + term1*xaa0_ad(i)*term3*term4*term5
           term1_ad = term1_ad + xaa0_ad(i)*term2*term3*term4*term5

           rfact_ad = rfact_ad + term5_ad
           gamma_ad = gamma_ad - one/(one+gamma)**2 * term4_ad
           xdby_ad  = xdby_ad + term3_ad
           xtol_ad(k-1,i) = xtol_ad(k-1,i) &
                - one/(xtol(k-1,i)**2) * term2_ad
           dz1_ad   = dz1_ad + term1_ad * (grav/cp)
           
           xheso_ad(k-1,i) = xheso_ad(k-1,i) - xdby_ad
           xhcko_ad(k-1,i) = xhcko_ad(k-1,i) + xdby_ad

           xtol_ad(k-1,i) = xtol_ad(k-1,i) &
                + (delta*cp/hvap)*gamma*rfact_ad
           gamma_ad = gamma_ad &
                + (delta*cp/hvap)*rfact_ad*xtol(k-1,i)

           xtol_ad(k-1,i) = xtol_ad(k-1,i) &
                - two*el2orc*xqesol(k-1,i)/(xtol(k-1,i)**3) &
                * gamma_ad
           xqesol_ad(k-1,i) = xqesol_ad(k-1,i) &
                + el2orc*gamma_ad/xtol(k-1,i)**2

           xzo_ad(k-1,i) = xzo_ad(k-1,i) - dz1_ad
           xzo_ad(k,i   )   = xzo_ad(k,i)      + dz1_ad

        endif
     end do
  end do
!
! Adjoint of xpwav calculation and initial xaa0 calculation
!
  do i = 1,im
     do k = kmax(i)-1,2,-1
        if (k > kb(i) .and. k < ktcon(i)) then
!
!          Redo nlm
           factor = eta(k-1,i)/eta(k,i)
           onemf  = one - factor
           termq  = half*(xqol(k,i) + xqol(k+1,i))
           gamma  = el2orc*xqesol(k,i)/xtol(k,i)**2
           xdby   = xhcko(k,i) - xheso(k,i)
           xtemp  = xtemp0(k,i)
           xqrch  = xqrch0(k,i)
           dq     = eta(k,i)*(xtemp-xqrch)
!
!          Adjoint code
           xtemp_ad  = zero
           termq_ad  = zero
           onemf_ad  = zero
           factor_ad = zero
           xqrch_ad  = zero
           gamma_ad  = zero
           xdby_ad   = zero
           xqc_ad    = zero
           qlk_ad    = zero
           xpw_ad    = zero
           dz_ad     = zero
           etah_ad   = zero
           dz1_ad    = zero
           dq_ad     = zero
!
!          Adjoint of actions dependent on sign of dq
           if (dq > zero) then
!
!             Redo nlm calculations
              dz   = half * (xzo(k+1,i) - xzo(k-1,i))
              dz1  = xzo(k,i) - xzo(k-1,i)
              etah = half * (eta(k,i) + eta(k-1,i))
              qlk  = dq / (eta(k,i) + etah * c0 * dz)
!
!             Adjoint code
              xpw_ad   = xpw_ad  + xpwav_ad(i)
              qlk_ad   = qlk_ad  + etah*c0*dz*xpw_ad
              dz_ad    = dz_ad   + etah*c0*xpw_ad*qlk
              etah_ad  = etah_ad + xpw_ad*c0*dz*qlk

              qlk_ad   = qlk_ad  - dz1*grav*xaa0_ad(i)
              dz1_ad   = dz1_ad  - xaa0_ad(i)*grav*qlk

              xqc_ad   = xqc_ad + xqcko_ad(k,i)
              xqrch_ad = xqrch_ad + xqc_ad
              qlk_ad   = qlk_ad + xqc_ad

              dz_ad    = dz_ad &
                   - (dq/(eta(k,i)+etah*c0*dz)**2) &
                   * etah*c0*qlk_ad
              etah_ad  = etah_ad &
                   - (dq/(eta(k,i)+etah*c0*dz)**2) &
                   * qlk_ad*c0*dz
              eta_ad(k,i) = eta_ad(k,i) &
                   - (dq/(eta(k,i)+etah*c0*dz)**2) &
                   * qlk_ad
              dq_ad = dq_ad + qlk_ad/(eta(k,i)+etah*c0*dz)

              eta_ad(k-1,i) = eta_ad(k-1,i) + half*etah_ad
              eta_ad(k,i)   = eta_ad(k,i)      + half*etah_ad
              xzo_ad(k-1,i) = xzo_ad(k-1,i) - dz1_ad
              xzo_ad(k,i)   = xzo_ad(k,i)      + dz1_ad
              xzo_ad(k-1,i) = xzo_ad(k-1,i) - half*dz_ad
              xzo_ad(k+1,i) = xzo_ad(k+1,i) + half*dz_ad
           else
              xtemp_ad = xtemp_ad + xqcko_ad(k,i)
           endif
!
!          Adjoint of dq calculation
           xqrch_ad = xqrch_ad - eta(k,i)*dq_ad
           xtemp_ad = xtemp_ad + eta(k,i)*dq_ad
           eta_ad(k,i) = eta_ad(k,i) + dq_ad*(xtemp-xqrch)
!
!          Adjoint of xqrch calculation
           xdby = xhcko(k,i) - xheso(k,i)
           if (xdby > zero) then
              xqesol_ad(k,i) = xqesol_ad(k,i) + xqrch_ad
              gamma_ad = gamma_ad &
                   - gamma*xdby/((hvap*(one+gamma))**2) &
                   * hvap*xqrch_ad
              xdby_ad  = xdby_ad + gamma*xqrch_ad/(hvap*(one+gamma))
              gamma_ad = gamma_ad + xqrch_ad*xdby/(hvap*(one+gamma))
              xtol_ad(k,i) = xtol_ad(k,i) &
                   - two*el2orc*xqesol(k,i)/xtol(k,i)**3 &
                   * gamma_ad
              xqesol_ad(k,i) = xqesol_ad(k,i) &
                   + el2orc*gamma_ad/xtol(k,i)**2
           else
              xqesol_ad(k,i) = xqesol_ad(k,i) + xqrch_ad
           endif
           xheso_ad(k,i) = xheso_ad(k,i) - xdby_ad
           xhcko_ad(k,i) = xhcko_ad(k,i) + xdby_ad
!
!          Adjoint of xtemp calculation
           termq_ad = termq_ad + onemf*xtemp_ad
           onemf_ad = onemf_ad + xtemp_ad*termq
           xqcko_ad(k-1,i) = xqcko_ad(k-1,i) + factor*xtemp_ad
           factor_ad = factor_ad + xtemp_ad*xqcko(k-1,i)

           xqol_ad(k+1,i) = xqol_ad(k+1,i) + half*termq_ad
           xqol_ad(k,i)   = xqol_ad(k,i)      + half*termq_ad

           factor_ad = factor_ad - onemf_ad

           eta_ad(k,i) = eta_ad(k,i) &
                - eta(k-1,i)/eta(k,i)**2 * factor_ad
           eta_ad(k-1,i) = eta_ad(k-1,i) + factor_ad/eta(k,i)
!
        endif
     end do
  end do
!
! Adjoint of xhcko calculation
  do i = 1,im
     do k = kmax(i)-1,2,-1
        if (k > kb(i).and.k <= ktcon(i)) then
!
!          Recalculate nlm
           factor = eta(k-1,i)/eta(k,i)
           onemf  = one - factor
!
!          Adjoint of tlm code
           onemf_ad  = zero
           factor_ad = zero

           xheo_ad(k+1,i) = xheo_ad(k+1,i) + onemf*half*xhcko_ad(k,i)
           xheo_ad(k,i)   = xheo_ad(k,i)   + onemf*half*xhcko_ad(k,i)
           onemf_ad = onemf_ad + &
                xhcko_ad(k,i) * half*(xheo(k,i)+xheo(k+1,i))
           xhcko_ad(k-1,i) = xhcko_ad(k-1,i) + factor*xhcko_ad(k,i)
           factor_ad = factor_ad + xhcko_ad(k,i)*xhcko(k-1,i)

           factor_ad = factor_ad - onemf_ad

           eta_ad(k,i) = eta_ad(k,i) - &
                (eta(k-1,i)/eta(k,i)**2) * factor_ad
           eta_ad(k-1,i) = eta_ad(k-1,i) + factor_ad/eta(k,i)
        endif
     end do
  end do
!
! Adjoint of xheo,xqol --> xhcko,xqcko transfer
  do i = 1,im
     indx = kb(i)
     xheo_ad(indx,i) = xheo_ad(indx,i) + xhcko_ad(indx,i)
     xqol_ad(indx,i) = xqol_ad(indx,i) + xqcko_ad(indx,i)
  end do
!
! Adjoint of xheo and xheso calculation at k=kmax(i)
  do i = 1,im
     k=kmax(i)
     xqesol_ad(k,i) = xqesol_ad(k,i) + hvap*xheso_ad(k,i)
     xtol_ad(k,i)   = xtol_ad(k,i)   + cp  *xheso_ad(k,i)
     xzo_ad(k,i)    = xzo_ad(k,i)    + grav   *xheso_ad(k,i)

     xqol_ad(k,i) = xqol_ad(k,i) + hvap*xheo_ad(k,i)
     xtol_ad(k,i) = xtol_ad(k,i) + cp  *xheo_ad(k,i)
     xzo_ad(k,i)  = xzo_ad(k,i)  + grav   *xheo_ad(k,i)
  enddo
!
! Adjoint of xheo and xheso calculation below k=kmax
  do i = 1,im
     do k = kmax(i)-1,1,-1
!
!       Redo nlm calculation
        po    = half * (p(k,i) + p(k+1,i))
        call fpvsx_ad(xtol(k,i),es0,adt,es0_ad,.false.)
        esl = 10.0_r_kind*es0
!
!       Adjoint of tlm code
        esl_ad = zero
        es0_ad = zero

        xqesol_ad(k,i) = xqesol_ad(k,i) + hvap *xheso_ad(k,i)
        xtol_ad(k,i)   = xtol_ad(k,i)   + cp   *xheso_ad(k,i)
        xzo_ad(k+1,i)  = xzo_ad(k+1,i)  + half*grav*xheso_ad(k,i)
        xzo_ad(k,i)    = xzo_ad(k,i)    + half*grav*xheso_ad(k,i)

        xqol_ad(k,i)   = xqol_ad(k,i)   + hvap *xheo_ad(k,i)
        xtol_ad(k,i)   = xtol_ad(k,i)   + cp   *xheo_ad(k,i)
        xzo_ad(k+1,i)  = xzo_ad(k+1,i)  + half*grav*xheo_ad(k,i)
        xzo_ad(k,i)    = xzo_ad(k,i)    + half*grav*xheo_ad(k,i)

        esl_ad = esl_ad - eps*esl/((po+epsm1*esl)**2) &
             *epsm1*xqesol_ad(k,i)
        esl_ad = esl_ad &
             + eps*xqesol_ad(k,i)/(po+epsm1*esl)

        es0_ad = es0_ad + 10*esl_ad

        adt = zero
        call fpvsx_ad(xtol(k,i),es0,adt,es0_ad,adjoint)
        xtol_ad(k,i) = xtol_ad(k,i) + adt
     end do
  end do
!
! Adjoint of xtol and xqol calculations
  do i = 1,im
     do k = kmax(i)-1,1,-1
!
!       Redo nlm calculations
        dz     = half * (xzo(k+1,i) - xzo(k,i))
        dp     = half * (p(k+1,i) - p(k,i))
        call fpvsx_ad(xto(k+1,i),es0,adt,es0_ad,.false.)
        es     = 10.0_r_kind*es0
        pprime = p(k+1,i) + epsm1 * es
        qs     = eps * es / pprime
        dqsdp  = - qs / pprime
        desdt  = ES * (FACTOR1 / xto(k+1,i) + FACTOR2 / (xto(k+1,i)**2))
        dqsdt  = qs * p(k+1,i) * desdt / (es * pprime)
        gamma  = el2orc * xqeso(k+1,i) / (xto(k+1,i)**2)
        dt     = (grav * dz + hvap*dqsdp*dp) / (cp*(one+gamma))
        dq     = dqsdt * dt + dqsdp * dp
!           
!       Initialize local ajm variables to zero
        dq_ad     = zero
        dt_ad     = zero
        dqsdp_ad  = zero
        dqsdt_ad  = zero
        gamma_ad  = zero
        dz_ad     = zero
        pprime_ad = zero
        es_ad     = zero
        desdt_ad  = zero
        qs_ad     = zero
        es0_ad    = zero
!
!       Adjoint of tlm code
        if (xqol0(k,i)>zero) then
           xqol0_ad(k,i) = xqol0_ad(k,i) + xqol_ad(k,i)
        else
           xqol0_ad(k,i) = zero
        endif
        xqo_ad(k+1,i) = xqo_ad(k+1,i) + xqol0_ad(k,i)
        dq_ad = dq_ad + xqol0_ad(k,i)

        xto_ad(k+1,i) = xto_ad(k+1,i) + xtol_ad(k,i)
        dt_ad = dt_ad + xtol_ad(k,i)

        dqsdp_ad = dqsdp_ad + dq_ad*dp
        dt_ad    = dt_ad + dqsdt*dq_ad
        dqsdt_ad = dqsdt_ad + dq_ad*dt

        gamma_ad = gamma_ad - &
             (grav*dz+hvap*dqsdp*dp)/((cp*(one+gamma))**2) * &
             cp*dt_ad
        dqsdp_ad = dqsdp_ad + hvap*dt_ad*dp/(cp*(one+gamma))
        dz_ad    = dz_ad + grav*dt_ad/(cp*(one+gamma))

        xto_ad(k+1,i) = xto_ad(k+1,i) - &
             two*el2orc*xqeso(k+1,i)/(xto(k+1,i)**3) * gamma_ad
        xqeso_ad(k+1,i) = xqeso_ad(k,i) + &
             el2orc*gamma_ad/(xto(k+1,i)**2)

        rterm = one/(es*pprime)
        pprime_ad = pprime_ad - (qs*p(k+1,i)*desdt)*rterm**2 * &
             es*dqsdt_ad
        es_ad     = es_ad - (qs*p(k+1,i)*desdt)*rterm**2 * &
             dqsdt_ad*pprime
        desdt_ad  = desdt_ad + p(k+1,i)*qs*dqsdt_ad*rterm
        qs_ad     = qs_ad + p(k+1,i)*dqsdt_ad*desdt*rterm

        es_ad = es_ad + &
             desdt_ad*(factor1/xto(k+1,i) + factor2/(xto(k+1,i)**2))
        xto_ad(k+1,i) = xto_ad(k+1,i) - desdt_ad * &
             (es/(xto(k+1,i)**2)) * (factor1+two*factor2/xto(k+1,i))

        pprime_ad = pprime_ad + (qs/(pprime**2))*dqsdp_ad
        qs_ad     = qs_ad - dqsdp_ad/pprime

        pprime_ad = pprime_ad - (eps*es/(pprime**2))*qs_ad
        es_ad     = es_ad + eps*qs_ad/pprime
            
        es_ad = es_ad + epsm1*pprime_ad

        es0_ad = es0_ad + 10*es_ad

        adt = zero
        call fpvsx_ad(xto(k+1,i),es0,adt,es0_ad,adjoint)
        xto_ad(k+1,i) = xto_ad(k+1,i) + adt

        xzo_ad(k,i)   = xzo_ad(k,i)   - half*dz_ad
        xzo_ad(k+1,i) = xzo_ad(k+1,i) + half*dz_ad
            
     end do
  end do

! Adjoint of integration of hydrostatic equation to get xzo
  do i = 1,im
     do k = kmax(i),2,-1
!          
!       Redo nlm calculations
        DLNSIG = LOG(SL(K,i) / SL(K-1,i))
        term1  = dlnsig * rd / grav
!
!       Adjoint of tlm code
        term2_ad = zero

        term2_ad = term2_ad - term1*xzo_ad(k,i)
        xzo_ad(k-1,i) = xzo_ad(k-1,i) + xzo_ad(k,i)

        xtvo_ad(k-1,i) = xtvo_ad(k-1,i) + half*term2_ad
        xtvo_ad(k,i)      = xtvo_ad(k,i)      + half*term2_ad
     end do
  end do
  do i = 1,im
     dlnsig = log(sl(1,i))
     xtvo_ad(1,i) = xtvo_ad(1,i) - dlnsig*rd/grav * xzo_ad(1,i)
  end do
!
! Adjoint of xqeso and xtvo calculations
  do i = 1,im
     do k = kmax(i),1,-1
        call fpvsx_ad(xto(k,i),es0,adt,es0_ad,.false.)
        es = 10.0_r_kind*es0

        es_ad  = zero
        es0_ad = zero

        xqo_ad(k,i) = xqo_ad(k,i) + delta*xto(k,i)*xtvo_ad(k,i)
        xto_ad(k,i) = xto_ad(k,i) + delta*xqo(k,i)*xtvo_ad(k,i)
        xto_ad(k,i) = xto_ad(k,i) + xtvo_ad(k,i)

        es_ad = es_ad &
             - eps*es/((p(k,i)+epsm1*es)**2) * &
             epsm1*xqeso_ad(k,i)
        es_ad = es_ad &
             + eps*xqeso_ad(k,i)/(p(k,i)+epsm1*es)

        es0_ad = es0_ad + 10*es_ad
        adt = zero
        call fpvsx_ad(xto(k,i),es0,adt,es0_ad,adjoint)
        xto_ad(k,i) = xto_ad(k,i) + adt
     end do
  end do

!
! Adjoint of xto,xqo calculation
  do i = 1,im
     do k = kmax(i),1,-1
        dellat_ad = zero

        if ( k <= ktcon(i) .and. ktcon(i)/=1 )  then
           to_ad(k,i) = to_ad(k,i) + xto_ad(k,i)
           dellat_ad  = dellat_ad + xto_ad(k,i)*mbdt

           dellaq_ad(k,i) = dellaq_ad(k,i) - hvap*dellat_ad/cp
           dellah_ad(k,i) = dellah_ad(k,i) + dellat_ad/cp
               
           if (xqo0(k,i)>zero) then
              xqo0_ad(k,i) = xqo0_ad(k,i) + xqo_ad(k,i)
           else
              xqo0_ad(k,i) = zero
           endif
           qo0_ad(k,i)    = qo0_ad(k,i) + xqo0_ad(k,i)
           dellaq_ad(k,i) = dellaq_ad(k,i) + xqo0_ad(k,i)*mbdt
        else
           to_ad(k,i) = to_ad(k,i) + xto_ad(k,i)
           qo0_ad(k,i)= qo0_ad(k,i)+ xqo_ad(k,i)
        endif
     end do
  end do
!
! Adjoint of dellah and dellaq calculations at cloud top
  do i = 1,im
     if (ktcon(i)/=1) then
        indx = ktcon(i)
!
!       Redo nlm calculations
        dp   = 100.0_r_kind * psfc(i)*del(indx,i)
        dv1  = heol(indx-1,i)
        dvq1 = qol(indx-1,i)
        dvu1 = uol(indx-1,i)
        dvv1 = vol(indx-1,i)
!
!       Adjoint of tlm code
        dvv1_ad = zero
        dvu1_ad = zero
        dvq1_ad = zero
        dv1_ad  = zero

        qlko_ktcon_ad(i) = qlko_ktcon_ad(i)       + &
             eta(indx-1,i)*dellal_ad(i) *grav/dp
        eta_ad(indx-1,i) = eta_ad(indx-1,i) + &
             dellal_ad(i)    *qlko_ktcon(i)*grav/dp


        dvv1_ad = dvv1_ad &
             - eta(indx-1,i) * dellav_ad(indx,i) * (grav/dp)
        vcko_ad(indx-1,i) = vcko_ad(indx-1,i) &
             + eta(indx-1,i) * dellav_ad(indx,i) * (grav/dp)
        eta_ad (indx-1,i) = eta_ad (indx-1,i) &
             + dellav_ad(indx,i)*(vcko(indx-1,i)-dvv1)*(grav/dp)
        vol_ad (indx-1,i) = vol_ad (indx-1,i) + dvv1_ad

        dvu1_ad = dvu1_ad &
             - eta(indx-1,i) * dellau_ad(indx,i) * (grav/dp)
        ucko_ad(indx-1,i) = ucko_ad(indx-1,i) &
             + eta(indx-1,i) * dellau_ad(indx,i) * (grav/dp)
        eta_ad (indx-1,i) = eta_ad (indx-1,i) &
             + dellau_ad(indx,i)*(ucko(indx-1,i)-dvu1)*(grav/dp)
        uol_ad (indx-1,i) = uol_ad (indx-1,i) + dvu1_ad

        dvq1_ad = dvq1_ad &
             - eta(indx-1,i) * dellaq_ad(indx,i) * (grav/dp) 
        qcko_ad(indx-1,i) = qcko_ad(indx-1,i) &
             + eta(indx-1,i) * dellaq_ad(indx,i) * (grav/dp)
        eta_ad (indx-1,i) = eta_ad (indx-1,i) &
             + dellaq_ad(indx,i)*(qcko(indx-1,i)-dvq1)*(grav/dp)
        qol_ad (indx-1,i) = qol_ad (indx-1,i) + dvq1_ad

        dv1_ad = dv1_ad &
             - eta(indx-1,i) * dellah_ad(indx,i) * (grav/dp)
        hcko1_ad(indx-1,i) = hcko1_ad(indx-1,i) &
             + eta(indx-1,i) * dellah_ad(indx,i) * (grav/dp)
        eta_ad  (indx-1,i) = eta_ad  (indx-1,i) &
             + dellah_ad(indx,i)*(hcko1(indx-1,i)-dv1)*(grav/dp)
        heol_ad (indx-1,i) = heol_ad (indx-1,i) + dv1_ad

     endif
  end do
!
! Adjoint of dellah and dellaq calculations within cloud
!
  do i = 1,im
     do k = kmax(i)-1,2,-1
        if ( k < ktcon(i) .and. ktcon(i)/=1 ) then
!
!          Redo nlm calculations
           aup = one
           if (k <= kb(i)) aup = zero
           adw = one
           if (k > jmin(i)) adw = zero
!
           dv1   = heol(k,i)
           dv2   = half * (heol(k,i) + heol(k+1,i))
           dv3   = heol(k-1,i)
           dv1q  = qol(k,i)
           dv2q  = half * (qol(k,i) + qol(k+1,i))
           dv3q  = qol(k-1,i)

           DV1u  = uOL(k,i)
           DV2u  = half * (uOL(k,i) + uOL(k+1,i))
           DV3u  = uOL(k-1,i)
           DV1v  = vOL(k,i)
           DV2v  = half * (vOL(k,i) + vOL(k+1,i))
           DV3v  = vOL(k-1,i)

           dp    = 100.0_r_kind * psfc(i) * del(k,i)
           detau = eta(k,i) - eta(k-1,i)
           detad = etad(k,i) - etad(k-1,i)              
!
           term1 = aup*eta(k,i) - adw*edto1(i)*etad(k,i)
           term2 = aup*detau
           term3 = aup*eta(k-1,i) - adw*edto1(i)*etad(k-1,i)
           term4 = adw*edto1(i)*detad
           termq = half*(qrcdo(k,i)+qrcdo(k-1,i))
!
!          Initialize loop local adjoint variables to zero
           termq_ad = zero
           term4_ad = zero
           term3_ad = zero
           term2_ad = zero
           term1_ad = zero
           dv3v_ad  = zero
           dv2v_ad  = zero
           dv1v_ad  = zero
           dv3u_ad  = zero
           dv2u_ad  = zero
           dv1u_ad  = zero
           dv3q_ad  = zero
           dv2q_ad  = zero
           dv1q_ad  = zero
           dv3_ad   = zero
           dv2_ad   = zero
           dv1_ad   = zero
           detad_ad = zero
           detau_ad = zero
!
!          Adjoint of tlm code.
           vcdo_ad(i) = vcdo_ad(i) + (grav/dp)*term4*dellav_ad(k,i)
           term4_ad   = term4_ad   + (grav/dp)*dellav_ad(k,i)*vcdo(i)
           dv2v_ad    = dv2v_ad    - (grav/dp)*term2*dellav_ad(k,i)
           term2_ad   = term2_ad   - (grav/dp)*dellav_ad(k,i)*dv2v
           dv3v_ad    = dv3v_ad    - (grav/dp)*term3*dellav_ad(k,i)
           term3_ad   = term3_ad   - (grav/dp)*dellav_ad(k,i)*dv3v
           dv1v_ad    = dv1v_ad    + (grav/dp)*term1*dellav_ad(k,i)
           term1_ad   = term1_ad   + (grav/dp)*dellav_ad(k,i)*dv1v

           ucdo_ad(i) = ucdo_ad(i) + (grav/dp)*term4*dellau_ad(k,i)
           term4_ad   = term4_ad   + (grav/dp)*dellau_ad(k,i)*ucdo(i)
           dv2u_ad    = dv2u_ad    - (grav/dp)*term2*dellau_ad(k,i)
           term2_ad   = term2_ad   - (grav/dp)*dellau_ad(k,i)*dv2u
           dv3u_ad    = dv3u_ad    - (grav/dp)*term3*dellau_ad(k,i)
           term3_ad   = term3_ad   - (grav/dp)*dellau_ad(k,i)*dv3u
           dv1u_ad    = dv1u_ad    + (grav/dp)*term1*dellau_ad(k,i)
           term1_ad   = term1_ad   + (grav/dp)*dellau_ad(k,i)*dv1u
               
           termq_ad   = termq_ad   + (grav/dp)*term4*dellaq_ad(k,i)
           term4_ad   = term4_ad   + (grav/dp)*termq*dellaq_ad(k,i)
           dv2q_ad    = dv2q_ad    - (grav/dp)*term2*dellaq_ad(k,i)
           term2_ad   = term2_ad   - (grav/dp)*dv2q *dellaq_ad(k,i)
           dv3q_ad    = dv3q_ad    - (grav/dp)*term3*dellaq_ad(k,i)
           term3_ad   = term3_ad   - (grav/dp)*dv3q *dellaq_ad(k,i)
           dv1q_ad    = dv1q_ad    + (grav/dp)*term1*dellaq_ad(k,i)
           term1_ad   = term1_ad   + (grav/dp)*dv1q *dellaq_ad(k,i)
!
           hcdo_ad(i) = hcdo_ad(i) + (grav/dp)*term4  *dellah_ad(k,i)
           term4_ad   = term4_ad   + (grav/dp)*hcdo(i)*dellah_ad(k,i)
           dv2_ad     = dv2_ad     - (grav/dp)*term2*dellah_ad(k,i)
           term2_ad   = term2_ad   - (grav/dp)*dv2  *dellah_ad(k,i)
           dv3_ad     = dv3_ad     - (grav/dp)*term3*dellah_ad(k,i)
           term3_ad   = term3_ad   - (grav/dp)*dv3  *dellah_ad(k,i)
           dv1_ad     = dv1_ad     + (grav/dp)*term1*dellah_ad(k,i)
           term1_ad   = term1_ad   + (grav/dp)*dv1  *dellah_ad(k,i)
!
           qrcdo_ad(k-1,i) = qrcdo_ad(k-1,i) + half*termq_ad
           qrcdo_ad(k,i)   = qrcdo_ad(k,i)      + half*termq_ad

           detad_ad    = detad_ad + adw*edto1(i)*term4_ad
           edto1_ad(i) = edto1_ad(i) + adw*term4_ad*detad

           etad_ad(k-1,i)  = etad_ad(k-1,i) - adw*edto1(i)*term3_ad
           edto1_ad(i)     = edto1_ad(i) - adw*term3_ad*etad(k-1,i)
           eta_ad(k-1,i)   = eta_ad(k-1,i) + aup*term3_ad
           
           detau_ad = detau_ad + aup*term2_ad

           etad_ad(k,i) = etad_ad(k,i) - adw*edto1(i)*term1_ad
           edto1_ad(i)  = edto1_ad(i)  - adw*term1_ad*etad(k,i)
           eta_ad(k,i)  = eta_ad(k,i) + aup*term1_ad

           etad_ad(k-1,i) = etad_ad(k-1,i) - detad_ad
           etad_ad(k,i)   = etad_ad(k,i)   + detad_ad

           eta_ad(k-1,i) = eta_ad(k-1,i) - detau_ad
           eta_ad(k,i)   = eta_ad(k,i)   + detau_ad

           uol_ad(k-1,i) = uol_ad(k-1,i) + dv3u_ad
           uol_ad(k+1,i) = uol_ad(k+1,i) + half*dv2u_ad
           uol_ad(k,i)   = uol_ad(k,i)   + half*dv2u_ad
           uol_ad(k,i)   = uol_ad(k,i)   + dv1u_ad

           vol_ad(k-1,i) = vol_ad(k-1,i) + dv3v_ad
           vol_ad(k+1,i) = vol_ad(k+1,i) + half*dv2v_ad
           vol_ad(k,i)   = vol_ad(k,i)   + half*dv2v_ad
           vol_ad(k,i)   = vol_ad(k,i)   + dv1v_ad

           qol_ad(k-1,i) = qol_ad(k-1,i) + dv3q_ad
           qol_ad(k+1,i) = qol_ad(k+1,i) + half*dv2q_ad
           qol_ad(k,i)   = qol_ad(k,i)   + half*dv2q_ad
           qol_ad(k,i)   = qol_ad(k,i)   + dv1q_ad

           heol_ad(k-1,i) = heol_ad(k-1,i) + dv3_ad
           heol_ad(k+1,i) = heol_ad(k+1,i) + half*dv2_ad
           heol_ad(k,i)   = heol_ad(k,i)   + half*dv2_ad
           heol_ad(k,i)   = heol_ad(k,i)   + dv1_ad
!
        endif
     end do
  end do
!     
! Adjoint of dellah and dellaq calculations at k=1
  do i = 1,im
     if (ktcon(i)/=1) then
        dp = 100.0_r_kind * psfc(i) * del(1,i)

        vol_ad(1,i) = vol_ad(1,i) &
             - edto1(i)*etad(1,i)*dellav_ad(1,i)*grav/dp
        vcdo_ad(i) = vcdo_ad(i) &
             + edto1(i)*etad(1,i)*dellav_ad(1,i)*grav/dp
        etad_ad(1,i) = etad_ad(1,i) &
             + edto1(i)*dellav_ad(1,i)*(vcdo(i)-vol(1,i))*grav/dp
        edto1_ad(i) = edto1_ad(i) &
             + dellav_ad(1,i)*etad(1,i)*(vcdo(i)-vol(1,i))*grav/dp

        uol_ad(1,i) = uol_ad(1,i) &
             - edto1(i)*etad(1,i)*dellau_ad(1,i)*grav/dp
        ucdo_ad(i) = ucdo_ad(i) &
             + edto1(i)*etad(1,i)*dellau_ad(1,i)*grav/dp
        etad_ad(1,i) = etad_ad(1,i) &
             + edto1(i)*dellau_ad(1,i)*(ucdo(i)-uol(1,i))*grav/dp
        edto1_ad(i) = edto1_ad(i) &
             + dellau_ad(1,i)*etad(1,i)*(ucdo(i)-uol(1,i))*grav/dp

        qol_ad(1,i) = qol_ad(1,i) &
             - edto1(i)*etad(1,i)*dellaq_ad(1,i)*grav/dp
        qcdo_ad(i) = qcdo_ad(i) &
             + edto1(i)*etad(1,i)*dellaq_ad(1,i)*grav/dp
        etad_ad(1,i) = etad_ad(1,i) &
             + edto1(i)*dellaq_ad(1,i)*(qcdo(i)-qol(1,i))*grav/dp
        edto1_ad(i) = edto1_ad(i) &
             + dellaq_ad(1,i)*etad(1,i)*(qcdo(i)-qol(1,i))*grav/dp
            
        heol_ad(1,i) = heol_ad(1,i) & 
             - edto1(i)*etad(1,i)*dellah_ad(1,i)*grav/dp
        hcdo_ad(i) = hcdo_ad(i) &
             + edto1(i)*etad(1,i)*dellah_ad(1,i)*grav/dp
        etad_ad(1,i) = etad_ad(1,i) &
             + edto1(i)*dellah_ad(1,i)*(hcdo(i)-heol(1,i))*grav/dp
        edto1_ad(i) = edto1_ad(i) &
             + dellah_ad(1,i)*etad(1,i)*(hcdo(i)-heol(1,i))*grav/dp
     endif
  end do
!
! Adjoint of cloud mass flux threshold function based
! on cloud work function after downdraft
  do i = 1,im
     term_ad = zero
     term = 1.e3_r_kind*aa1(i)
     if (abs(term) < 100._r_kind) then
        term_ad = dftanh(term)*tcwfdn_ad(i)
     else
        if (term < -100._r_kind) then
           term_ad = zero
        else
           term_ad = zero
        endif
     endif
     aa1_ad(i) = aa1_ad(i) + 1.e3_r_kind*term_ad
  end do
!
! Adjoint of downdraft cloudwork function calculation
  do i = 1,im
     do k = 1,kmax(i)-1
        if (dwnflg2(i) .and. k < jmin(i)) then
!
!          Redo nlm calculations
           GAMMA = EL2ORC * QESOL(k+1,i) / TOL(k+1,i)**2
           DHH   = HCDO(I)
           DT    = TOL(k+1,i)
           DG    = GAMMA
           DH    = HESOL(k+1,i)
           dz    = zo(k,i) - zo(k+1,i)
!
           term1 = edto1(i)
           term2 = dz
           term3 = grav/(cp*dt)
           term4 = (dhh-dh)
           term5 = one/(one+dg)
           term6 = one + delta*cp*dg*dt/hvap
!     
           dqs = qesol(k+1,i) - qol(k+1,i)
!
!          Adjoint of tlm code
           dqs_ad = zero
           dz_ad  = zero
           dt_ad  = zero
           dg_ad  = zero
           dh_ad  = zero
           dhh_ad = zero
           if (dqs > zero) then
              dqs_ad = dqs_ad + edto1(i)*dz*grav*delta * aa1_ad(i)
              dz_ad  = dz_ad + edto1(i)*grav*delta*dqs * aa1_ad(i)
              edto1_ad(i) = edto1_ad(i) + dz*grav*delta*dqs * aa1_ad(i)
           endif
                  
           qol_ad(k+1,i)   = qol_ad  (k+1,i) - dqs_ad
           qesol_ad(k+1,i) = qesol_ad(k+1,i) + dqs_ad
! 
           term1_ad = term2*term3*term4*term5*term6 * aa1_ad(i)
           term2_ad = term1*term3*term4*term5*term6 * aa1_ad(i)
           term3_ad = term1*term2*term4*term5*term6 * aa1_ad(i)
           term4_ad = term1*term2*term3*term5*term6 * aa1_ad(i)
           term5_ad = term1*term2*term3*term4*term6 * aa1_ad(i)
           term6_ad = term1*term2*term3*term4*term5 * aa1_ad(i)
! 
           dt_ad  = (delta*cp/hvap) * dg*term6_ad
           dg_ad  = (delta*cp/hvap) * term6_ad*dt
           dg_ad  = dg_ad - one/(one+dg)**2 * term5_ad
           dh_ad  = dh_ad - term4_ad
           dhh_ad = dhh_ad + term4_ad
           dt_ad  = dt_ad - one*grav/(cp*dt*dt) * term3_ad
           dz_ad  = dz_ad + term2_ad
           edto1_ad(i) = edto1_ad(i) + term1_ad
!
           zo_ad(k,i)     = zo_ad(k,i)   + dz_ad
           zo_ad(k+1,i)   = zo_ad(k+1,i) - dz_ad
           hesol_ad(k+1,i)= hesol_ad(k+1,i) + dh_ad
           gamma_ad       = dg_ad
           tol_ad(k+1,i)  = tol_ad(k+1,i) + dt_ad
           hcdo_ad(i)     = hcdo_ad(i) + dhh_ad
           tol_ad(k+1,i)  = tol_ad(k+1,i) &
                - two*el2orc*qesol(k+1,i)/tol(k+1,i)**3 * gamma_ad
           qesol_ad(k+1,i)= qesol_ad(k+1,i) + &
                el2orc*gamma_ad/tol(k+1,i)**2
        endif
     end do
  end do
!
! Adjoint of final edto1 calculation
  do i = 1,im
     edtmax = edtmaxl
     if (nint(slimsk(i)) == 0) edtmax = edtmaxs
     if (dwnflg2(i)) then
        if (pwevo(i) < zero) then
           if (edto10(i) > edtmax) then
!             do nothing since tlm is zero
              edto1_ad(i) = zero
           endif
           pwevo_ad(i) = pwevo_ad(i) &
                + edto(i)*pwavo(i)/pwevo(i)**2 * edto1_ad(i)
           pwavo_ad(i) = pwavo_ad(i) &
                - edto(i)*edto1_ad(i)/pwevo(i)
           edto_ad(i) = edto_ad(i) &
                - edto1_ad(i)*pwavo(i)/pwevo(i)
        else
!          do nothing since tlm is zero
           edto1_ad(i) = zero
        endif
     else
!       do nothing since tlm is zero             
        edto1_ad(i) = zero
     endif
  enddo
!
! Adjoint of qrcdo --> qcdo transfer
  do i = 1,im
     qrcdo_ad(1,i) = qrcdo_ad(1,i) + qcdo_ad(i)
  end do
!
! Adjoint of pwevo, pwdo, qrcdo calculation
  do i = 1,im
     do k = 1,kmax(i)-1
        if (k < jmin(i)) then
!
!          Recalculate part of foward model
           dq    = qesol(k,i)
           dt    = tol(k,i)
           gamma = el2orc * dq / dt**2
           dh    = hcdo(i) - hesol(k,i)
           term5 = (one/hvap) * gamma
           term6 = one/(one+gamma)
           detad = etad(k+1,i) - etad(k,i)               
!
!          Adjoint of forward tlm
           term1_ad = zero
           term2_ad = zero
           term3_ad = zero
           pwdo_ad(k,i) = pwdo_ad(k,i) + pwevo_ad(i)

           if (k < jmin(i)-1) then
              term1_ad = term1_ad + pwdo_ad(k,i)
              term2_ad = term2_ad - pwdo_ad(k,i)
              term3_ad = term3_ad - pwdo_ad(k,i)

              qrcdo_ad(k+1,i) = qrcdo_ad(k+1,i) + detad*half*term3_ad
              qrcdo_ad(k,i) = qrcdo_ad(k,i) + detad*half*term3_ad
              detad_ad = +half*(qrcdo(k,i)+qrcdo(k+1,i))*term3_ad
              qrcdo_ad(k,i) = qrcdo_ad(k,i) + etad(k,i)*term2_ad
              etad_ad(k,i) = etad_ad(k,i) + qrcdo(k,i)*term2_ad
              qcdo_ad(i) =qcdo_ad(i) +term1_ad*etad(k+1,i)
              etad_ad(k+1,i) =etad_ad(k+1,i) +term1_ad*qcdo(i)
           endif

           etad_ad(k,i)   = etad_ad(k,i)      - detad_ad
           etad_ad(k+1,i) = etad_ad(k+1,i) + detad_ad

           dh_ad    = term5 * term6 * qrcdo_ad(k,i)
           term6_ad = term5 * qrcdo_ad(k,i) * dh
           term5_ad = qrcdo_ad(k,i) * term6 * dh
           dq_ad    = qrcdo_ad(k,i)

           gamma_ad = -one/(one+gamma)**2 * term6_ad
           gamma_ad = gamma_ad + (one/hvap) * term5_ad
           term5_ad = zero
           term6_ad = zero

           hesol_ad(k,i) = hesol_ad(k,i) - dh_ad
           hcdo_ad(i) = hcdo_ad(i) + dh_ad
           dh_ad = zero

           dt_ad = -two*el2orc*dq/dt**3 * gamma_ad
           dq_ad = dq_ad + el2orc*gamma_ad/dt**2
           gamma_ad = zero

           tol_ad(k,i) = tol_ad(k,i) + dt_ad
           qesol_ad(k,i) = qesol_ad(k,i) + dq_ad
           dt_ad = zero
           dq_ad = zero

        endif
     end do
  end do

!
! Adjoint of hcdo and qrcdo transfer
  do i = 1,im
     jmn = jmin(i)
     heol_ad(jmn,i)  = heol_ad(jmn,i) + hcdo_ad(i)
     qesol_ad(jmn,i) = qesol_ad(jmn,i) + qrcdo_ad(jmn,i)
     uol_ad(jmn,i)   = uol_ad(jmn,i) + ucdo_ad(i)
     vol_ad(jmn,i)   = vol_ad(jmn,i) + vcdo_ad(i)
  end do
!
! Adjoint of etad (downdraft entrainment) calculation at k=1
  k = 1
  do i = 1,im
     dz_ad = zero
     if (kbdtr(i) > 1) then
!
!       Redo nlm
        dz = half * (zo(2,i)-zo(1,i))
!
!       Adjoint of tlm code.
        dz_ad = dz_ad + &
             etad(k+1,i)*exp(xlamd(i)*dz) * xlamd(i)*etad_ad(k,i)
        xlamd_ad(i) = xlamd_ad(i) + &
             etad(k+1,i)*exp(xlamd(i)*dz) * dz*etad_ad(k,i)
        etad_ad(k+1,i) = etad_ad(k+1,i) + &
             exp(xlamd(i)*dz) * etad_ad(k,i)

        zo_ad(1,i) = zo_ad(1,i) - half*dz_ad
        zo_ad(2,i) = zo_ad(2,i) + half*dz_ad
     endif
  end do
!
! Adjoint of etad (downdraft entrainment) calculation above k=1
  do i = 1,im
     do k = 2,kbmax(i)
        dz_ad = zero
        if (k < kbdtr(i)) then
!
!          Redo nlm
           dz = half * (zo(k+1,i)-zo(k-1,i))
!
!          Adjoint of tlm code.
           dz_ad = dz_ad + &
                etad(k+1,i)*exp(xlamd(i)*dz) * xlamd(i)*etad_ad(k,i)
           xlamd_ad(i) = xlamd_ad(i) + &
                etad(k+1,i)*exp(xlamd(i)*dz) * dz*etad_ad(k,i)
           etad_ad(k+1,i) = etad_ad(k+1,i) + &
                exp(xlamd(i)*dz) * etad_ad(k,i)
           
           zo_ad(k-1,i) = zo_ad(k-1,i) - half*dz_ad
           zo_ad(k+1,i) = zo_ad(k+1,i) + half*dz_ad
        endif
     end do
  end do
!
! Adjoint of xlamd calculation
  do i = 1,im
     dz_ad = zero
     beta = betas
     if (nint(slimsk(i)) == 1) beta = betal
     if (kbdtr(i) > 1) then
!
!       Redo nlm calculations
        dz = half * (zo(kbdtr(i),i) + zo(kbdtr(i)-1,i)) &
             - zo(1,i)
!
!       Adjoint code
        dz_ad = dz_ad - log(beta)/(dz*dz)*xlamd_ad(i)
        zo_ad(1,i) = zo_ad(1,i) - dz_ad
        zo_ad(kbdtr(i)-1,i) = zo_ad(kbdtr(i)-1,i) + half*dz_ad
        zo_ad(kbdtr(i),i)      = zo_ad(kbdtr(i),i)      + half*dz_ad
     endif
  end do
!
! Adjoint of upper and lower bounds on precipitation efficiency.
  do i = 1,im
     edt_ad(i) = edt_ad(i) + edtx_ad(i)
     edt_ad(i) = edt_ad(i) + edto_ad(i)
     if (edt0(i) > 0.9_r_kind) then
        edt_ad(i) = zero
     elseif (edt0(i) < zero) then
        edt_ad(i) = zero
     endif
  end do
!
! Adjoint for precipitation efficiency (edt) calculation
  do i = 1,im
     e1_ad   = zero
     term_ad = zero
     dz_ad   = zero
!
!    Redo nlm calculations
     dz = zo(ktcon(i),i)-zo(kb(i),i)
     if (abs(dz) > tiny_r_kind) then
        term = 1.e3_r_kind * vshear(i)/dz
!
!       Adjoint code
        e1_ad   = e1_ad - edt_ad(i)

        term_ad = term_ad + pcpeff3*three*term**2 * e1_ad
        term_ad = term_ad + pcpeff2*two*term * e1_ad
        term_ad = term_ad + pcpeff1 * e1_ad

        dz_ad = dz_ad - 1.e3_r_kind*vshear(i)/(dz*dz) * term_ad
        vshear_ad(i) = vshear_ad(i) + 1.e3_r_kind*term_ad/dz
     else
        term_ad = zero
        dz_ad = zero
     endif

     zo_ad(kb(i),i) = zo_ad(kb(i),i) - dz_ad
     zo_ad(ktcon(i),i) = zo_ad(ktcon(i),i) + dz_ad

  end do
!
! Adjoint of windshear calculation
  do i = 1,im
     do k = kmax(i),1,-1
        shear_ad = zero
        termv_ad = zero
        termu_ad = zero
        if (k >= kb(i) .and. k <= ktcon(i)) then
!
!          Adjoint code for vshear
           shear_ad = shear_ad + vshear_ad(i)
!
!          Redo nlm calculations
           termu = (uol(k+1,i) - uol(k,i))**2
           termv = (vol(k+1,i) - vol(k,i))**2
!
!          Adjoint code for uol and vol
           if (termu+termv > tiny_r_kind) then
              termv_ad = termv_ad + &
                   rcs(i)*(half/sqrt(termu+termv)) * shear_ad
              termu_ad = termu_ad + &
                   rcs(i)*(half/sqrt(termu+termv)) * shear_ad
           else
              shear_ad = zero
              termu_ad = zero
              termv_ad = zero
           endif

           vol_ad(k,i)      = vol_ad(k,i)      &
                - two*(vol(k+1,i)-vol(k,i))*termv_ad
           vol_ad(k+1,i) = vol_ad(k+1,i) &  
                + two*(vol(k+1,i)-vol(k,i))*termv_ad

           uol_ad(k,i)      = uol_ad(k,i)      &
                - two*(uol(k+1,i)-uol(k,i))*termu_ad
           uol_ad(k+1,i) = uol_ad(k+1,i) &  
                + two*(uol(k+1,i)-uol(k,i))*termu_ad

        endif
     end do
  end do
!
! Adjoint of cloud mass flux threshold function based
! on updraft cloud work function
  do i = 1,im
     term_ad = zero
     term = 1.e3_r_kind*aa10(i)
     if (abs(term) < 100._r_kind) then
        term_ad = dftanh(term)*tcwfup_ad(i)
     else
        if (term < -100.0_r_kind) then
           term_ad = zero
        else
           term_ad = zero
        endif
     endif
     aa1_ad(i) = aa1_ad(i) + 1.e3_r_kind*term_ad
  end do
!
!
! Adjoint of updraft cwf calculation
  do i=1,im
     do k=kmax(i),1,-1
        if (k > kbcon(i).and.k <= ktcon(i)) then
!
!          Recalculate part of forward model
           dz1   = zo(k,i) - zo(k-1,i)
           gamma = el2orc * qesol(k-1,i) / (tol(k-1,i)**2)
           rfact =  one + delta*cp*gamma*tol(k-1,i)/hvap
           term1 = dz1*grav/cp
           term2 = one/tol(k-1,i)
           term3 = dbyo(k-1,i)
           term4 = one/(one+gamma)
           term5 = rfact
           dqs = qesol(k-1,i) - qol(k-1,i)
!
!          Adjoint of forward model tlm
           dqs_ad = zero
           dz1_ad = zero
           if (dqs > zero) then
              dqs_ad   = dqs_ad + dz1*grav*delta * aa1_ad(i)
              dz1_ad   = dz1_ad + aa1_ad(i) * grav*delta*dqs
           endif

           qol_ad(k-1,i)   = qol_ad  (k-1,i) - dqs_ad
           qesol_ad(k-1,i) = qesol_ad(k-1,i) + dqs_ad

           term5_ad = term1*term2*term3*term4*aa1_ad(i)
           term4_ad = term1*term2*term3*term5*aa1_ad(i)
           term3_ad = term1*term2*term4*term5*aa1_ad(i)
           term2_ad = term1*term3*term4*term5*aa1_ad(i)
           term1_ad = term2*term3*term4*term5*aa1_ad(i)

           rfact_ad = term5_ad
           gamma_ad = -one/((one+gamma)**2) * term4_ad
           dbyo_ad(k-1,i) = dbyo_ad(k-1,i) + term3_ad
           tol_ad(k-1,i) = tol_ad(k-1,i) - &
                (one/(tol(k-1,i)**2)) * term2_ad
           dz1_ad = dz1_ad + term1_ad*grav/cp

           tol_ad(k-1,i) = tol_ad(k-1,i) + &
                delta*cp*gamma*rfact_ad/hvap
           gamma_ad = gamma_ad + &
                delta*cp*rfact_ad*tol(k-1,i)/hvap

           tol_ad(k-1,i) = tol_ad(k-1,i) - &
                two*el2orc*qesol(k-1,i)/(tol(k-1,i)**3)* &
                gamma_ad
           qesol_ad(k-1,i) = qesol_ad(k-1,i) + &
                el2orc*gamma_ad/(tol(k-1,i)**2)

           zo_ad(k,i)   = zo_ad(k,i)   + dz1_ad
           zo_ad(k-1,i) = zo_ad(k-1,i) - dz1_ad

        endif
     end do
  end do

! Adjoint of section is for cloud liquid water
  if (ncloud > 0) then
     do i=1,im
        k = ktcon(i)
        if (cnvflg(i)) then

!          Recompute nlm
           gamma = EL2ORC * QESOl(K,i) / (TOl(K,i)**2)
           QRCH = QESOl(K,i) &
                + GAMMA * DBYO(K,i) / (HVAP * (one + GAMMA))
           DQ = qcko00(K-1,i) - QRCH

!          Adjoint of tangent linear model
           dq_ad   = zero
           qrch_ad = zero
           gamma_ad= zero
           
           if (dq > zero) then
              qrch_ad = qrch_ad + qcko_ad(k-1,i)
              dq_ad   = dq_ad + qlko_ktcon_ad(i)
           endif

           qrch_ad = qrch_ad - dq_ad
           qcko_ad(k-1,i) = qcko_ad(k-1,i) + dq_ad

           gamma_ad = gamma_ad + &
                qrch_ad*dbyo(k,i)/(hvap*(one+gamma)*(one+gamma))
           dbyo_ad(k,i) = dbyo_ad(k,i) + &
                gamma*qrch_ad/(hvap*(one+gamma))
           qesol_ad(k,i) = qesol_ad(k,i) + qrch_ad

           tol_ad(k,i) = tol_ad(k,i) - gamma_ad * &
                two*el2orc*qesol(k,i)/(tol(k,i)**3)
           qesol_ad(k,i) = qesol_ad(k,i) + &
                el2orc*gamma_ad/(tol(k,i)**2)
        endif
     end do
  endif

!
! Adjoint of pwavo, pwo, qcko, and aa1 calculation
!
  do i = 1,im
     do k = kmax(i),1,-1
        if (k > kb(i).and.k < ktcon(i)) then
!
!          Recalculate part of forward model
           factor = eta(k-1,i)/eta(k,i)
           onemf  = one - factor
           temp   = qcko0(k,i)
           gamma  = el2orc*qesol(k,i)/(tol(k,i)**2) 
           qrch   = qesol(k,i) + &
                gamma*dbyo(k,i)/(hvap*(one+gamma))
           dq = dqk(k,i)
!
           if (dq > zero) then
!
!             Redo nlm
              dz   = half*(zo(k+1,i)-zo(k-1,i))
              dz1  = zo(k,i) - zo(k-1,i)
              etah = half*(eta(k,i)+eta(k-1,i))
              qlk  = dq/(eta(k,i) + etah*c0*dz)
!
!             Adjoint of tlm
              temp_ad = zero
              qc_ad = qcko_ad(k,i)
              qrch_ad = qc_ad
              qlk_ad = qc_ad

              pwo_ad(k,i) = pwo_ad(k,i) + pwavo_ad(i)
                   
              qlk_ad  = qlk_ad + etah*c0*dz*pwo_ad(k,i)
              dz_ad   = etah*c0*pwo_ad(k,i)*qlk
              etah_ad = pwo_ad(k,i)*c0*dz*qlk

              qlk_ad = qlk_ad  - dz1*grav*aa1_ad(i)
              dz1_ad = -one*aa1_ad(i)*grav*qlk
                   
              dz_ad  = dz_ad - &
                   (dq/((eta(k,i)+etah*c0*dz)**2))* &
                   etah*c0*qlk_ad
              etah_ad = etah_ad - &
                   (dq/((eta(k,i)+etah*c0*dz)**2))* &
                   qlk_ad*c0*dz
              eta_ad(k,i) = eta_ad(k,i) - &
                   (dq/((eta(k,i)+etah*c0*dz)**2))* &
                   qlk_ad
              dq_ad = qlk_ad/(eta(k,i)+etah*c0*dz)
              
              eta_ad(k-1,i) = eta_ad(k-1,i) + half*etah_ad
              eta_ad(k,i)   = eta_ad(k,i)      + half*etah_ad
                   
              zo_ad(k-1,i) = zo_ad(k-1,i) - dz1_ad
              zo_ad(k,i)   = zo_ad(k,i)      + dz1_ad
              zo_ad(k-1,i) = zo_ad(k-1,i) - half*dz_ad
              zo_ad(k+1,i) = zo_ad(k+1,i) + half*dz_ad

           else
              qc_ad   = zero
              qrch_ad = zero
              qlk_ad  = zero
              dz_ad   = zero
              etah_ad = zero
              dz1_ad  = zero
              dq_ad   = zero
              temp_ad = qcko_ad(k,i)
           endif
              
           qrch_ad = qrch_ad - eta(k,i)*dq_ad
           temp_ad = temp_ad + eta(k,i)*dq_ad
           eta_ad(k,i) = eta_ad(k,i) + dq_ad*(temp-qrch)
!
           gamma_ad = -gamma*dbyo(k,i)/((hvap*(one+gamma))**2)* &
                hvap*qrch_ad
           dbyo_ad(k,i) = dbyo_ad(k,i) + &
                gamma*qrch_ad/(hvap*(one+gamma))
           gamma_ad = gamma_ad + &
                qrch_ad*dbyo(k,i)/(hvap*(one+gamma))
           qesol_ad(k,i) = qesol_ad(k,i) + qrch_ad
           
           tol_ad(k,i) = tol_ad(k,i) - &
                el2orc*2*qesol(k,i)/(tol(k,i)**3)*gamma_ad
           qesol_ad(k,i) = qesol_ad(k,i) + &
                el2orc/(tol(k,i)**2)*gamma_ad
              
           qol_ad(k+1,i) = qol_ad(k+1,i) + onemf*half*temp_ad
           qol_ad(k,i)   = qol_ad(k,i)   + onemf*half*temp_ad
           onemf_ad = temp_ad*half*(qol(k,i)+qol(k+1,i))
           qcko_ad(k-1,i) = qcko_ad(k-1,i) + factor*temp_ad
           factor_ad = temp_ad*qcko(k-1,i)
              
           factor_ad = factor_ad - onemf_ad
                   
           eta_ad(k,i)      = eta_ad(k,i)      - &
                (eta(k-1,i)/(eta(k,i)**2))*factor_ad
           eta_ad(k-1,i) = eta_ad(k-1,i) + (one/eta(k,i))*factor_ad

        endif
     end do
  end do
!
! Adjoint of qol to qcko transfer
  do i = 1,im
     indx = kb(i)
     qol_ad(indx,i) = qol_ad(indx,i) + qcko_ad(indx,i)
  end do

!
! Adjoint of hcko,hckod transfer to hcko1 and dbyo calculation
  do i = 1,im
     do k = kmax(i)-1,2,-1
        if (k > kb(i)) then
           hesol_ad(k,i) = hesol_ad(k,i) - dbyo_ad(k,i)
           hcko1_ad(k,i) = hcko1_ad(k,i) + dbyo_ad(k,i)
        endif
        if (k > kb(i).and.k <= kbcon(i)) then
           hcko_ad(k,i) = hcko_ad(k,i) + hcko1_ad(k,i)
        elseif (k > kbcon(i).and.k <= ktcon(i)) then
           if (.not.dwnflg(i)) then
              hcko_ad(k,i) = hcko_ad(k,i) + hcko1_ad(k,i)
           else
              hckod_ad(k,i) = hckod_ad(k,i) + hcko1_ad(k,i)
           endif
        elseif (k > ktcon(i)) then
           hcko_ad(k,i) = hcko_ad(k,i) + hcko1_ad(k,i)
        endif
     end do
  end do

!
! Adjoint of modification of hcko1 by detrainment process
  do i = 1,im
     do k = kmax(i)-1,2,-1
        if (k > kbcon(i) .and. k <= ktcon(i)) then
!
!          Redo nlm calculations
           factor = eta(k-1,i)/eta(k,i)
           onemf  = one - factor
           fuv    = etau(k-1,i)/etau(k,i)
           onemfu = one - fuv
           heol2  = half*(heol(k,i) + heol(k+1,i))
           uol2   = half*(uol(k,i)+uol(k+1,i))
           vol2   = half*(vol(k,i)+vol(k+1,i))
!
!          Adjoint code
           fuv_ad    = zero
           onemfu_ad = zero
           onemf_ad  = zero
           factor_ad = zero
           heol2_ad  = zero
           vol2_ad   = zero
           uol2_ad   = zero

           vol2_ad   = vol2_ad + onemfu*vckod_ad(k,i)
           onemfu_ad = onemfu_ad + vckod_ad(k,i)*vol2
           vckod_ad(k-1,i) = vckod_ad(k-1,i) + fuv*vckod_ad(k,i)
           fuv_ad = fuv_ad + vckod_ad(k,i)*vckod(k-1,i)

           uol2_ad   = uol2_ad + onemfu*ukcod_ad(k,i)
           onemfu_ad = onemfu_ad + ukcod_ad(k,i)*uol2
           ukcod_ad(k-1,i) = ukcod_ad(k-1,i) + fuv*ukcod_ad(k,i)
           fuv_ad = fuv_ad + ukcod_ad(k,i)*uckod(k-1,i)
           
           heol2_ad = heol2_ad + onemf*hckod_ad(k,i)
           onemf_ad = onemf_ad + heol2*hckod_ad(k,i)
           hckod_ad(k-1,i) = hckod_ad(k-1,i) + factor*hckod_ad(k,i)
           factor_ad = factor_ad + hckod(k-1,i)*hckod_ad(k,i)
               
           heol_ad(k,i)   = heol_ad(k,i)   + half*heol2_ad
           heol_ad(k+1,i) = heol_ad(k+1,i) + half*heol2_ad

           fuv_ad = fuv_ad - onemfu_ad
           etau_ad(k-1,i) = etau_ad(k-1,i) + fuv_ad/etau(k,i)
           etau_ad(k,i) = etau_ad(k,i) -  fuv_ad * &
                etau(k-1,i)/(etau(k,i)**2)

           factor_ad = factor_ad - onemf_ad

           eta_ad(k,i) = eta_ad(k,i) - &
                factor_ad*(eta(k-1,i)/(eta(k,i)**2))
           eta_ad(k-1,i) = eta_ad(k-1,i) + factor_ad/eta(k,i)
        endif
     end do
  end do

!
! Adjoint of hcko to hckod transfer at kbcon
  do i = 1,im
     if (dwnflg(i)) then
        indx = kbcon(i)
        hcko_ad(indx,i) = hcko_ad(indx,i) + hckod_ad(indx,i)
        ucko_ad(indx,i) = ucko_ad(indx,i) + ukcod_ad(indx,i)
        vcko_ad(indx,i) = vcko_ad(indx,i) + vckod_ad(indx,i)
     endif
  end do
!
! Adjoint of eta calculation
  do i = 1,im
     if (dwnflg(i)) then
        do k = kmax(i)-1,2,-1
           if (k > jmin(i) .and. k <= kt2(i)) then

!             Redo nlm calculation                  
              dz = half*(zo(k+1,i)-zo(k-1,i))

!             Adjoint code
              dz_ad = zero

              dz_ad = dz_ad + etau(k-1,i) * &
                   (xlamdet(i)+xlambu)*etau_ad(k,i)

              xlamdet_ad(i) = xlamdet_ad(i) + & 
                   etau(k-1,i)*etau_ad(k,i)*dz

              etau_ad(k-1,i) = etau_ad(k-1,i) + & 
                   etau_ad(k,i)*(one+(xlamdet(i)+xlambu)*dz)

              dz_ad = dz_ad + eta(k-1,i)*xlamdet(i)*eta_ad(k,i)

              xlamdet_ad(i) = xlamdet_ad(i) + &
                   eta(k-1,i)*eta_ad(k,i)*dz

              eta_ad(k-1,i) = eta_ad(k-1,i) + &
                   eta_ad(k,i)*(one+xlamdet(i)*dz)

              zo_ad(k-1,i) = zo_ad(k-1,i) - half*dz_ad
              zo_ad(k+1,i) = zo_ad(k+1,i) + half*dz_ad
           endif
        end do
     else
        do k=kmax(i)-1,2,-1
           if (k > jmin(i) .and. k <= ktcon0(i)) then

!             Redo nlm calculation                  
              DZ = half * (ZO(k+1,i) - ZO(k-1,i))

!             Adjoint code
              dz_ad = zero

              dz_ad = dz_ad + etau(k-1,i)*xlambu*etau_ad(k,i)
              etau_ad(k-1,i) = etau_ad(k-1,i) + &
                   etau_ad(k,i)*(one+xlambu*dz)
           endif
        end do
     endif
  end do


!
!
! Adjoint of modification of cloud properties below cloud
! base via the entrainment process
  do i = 1,im
     tem1 = hcko(jmin(i),i) - hesol(kt2(i),i)
     tem2 = sumz(kt2(i),i) * hesol(kt2(i),i) - sumh(kt2(i),i)

     tem1_ad = zero
     tem2_ad = zero

     if (sumz(kt2(i),i)  >  0.000001_r_kind) then
        term = 2.3_r_kind/sumz(kt2(i),i)
        if (xlamdet1(i) < term) then
           xlamdet1_ad(i) = xlamdet1_ad(i) + xlamdet_ad(i)
        else
           term_ad = zero
           term_ad = term_ad + xlamdet_ad(i)
           sumz_ad(kt2(i),i) = sumz_ad(kt2(i),i) - term_ad * &
                2.3_r_kind/(sumz(kt2(i),i)**2)
        endif
     else
        xlamdet1_ad(i) = xlamdet1_ad(i) + xlamdet_ad(i)
     endif

     if (xlamdet0(i) >= zero) then
        xlamdet0_ad(i) = xlamdet0_ad(i) + xlamdet1_ad(i)
     else
        xlamdet0_ad(i) = zero
     endif

     if (abs(tem2)  >  0.000001_r_kind) then
        tem1_ad = zero
        tem2_ad = zero
        tem1_ad = tem1_ad + xlamdet0_ad(i)/tem2
        tem2_ad = tem2_ad - tem1/(tem2**2) * xlamdet0_ad(i)
     else
        tem1_ad = zero
        tem2_ad = zero
     endif

     sumz_ad(kt2(i),i) =sumz_ad(kt2(i),i) + tem2_ad*hesol(kt2(i),i)
     hesol_ad(kt2(i),i)=hesol_ad(kt2(i),i)+ tem2_ad*sumz(kt2(i),i)
     sumh_ad(kt2(i),i) =sumh_ad(kt2(i),i) - tem2_ad
         
     hesol_ad(kt2(i),i) = hesol_ad(kt2(i),i) - tem1_ad
     hcko_ad(jmin(i),i)  = hcko_ad(jmin(i),i) + tem1_ad

  end do

  do i=1,im
     do k=kmax(i)-1,2,-1
        if (k > jmin(i) .and. k <= ktcon0(i)) then

           heol_ad(k,i)      = heol_ad(k,i) + &
                half * (zo(k+1,i) - zo(k-1,i)) * sumh_ad(k,i)
           zo_ad(k-1,i)   = zo_ad(k-1,i)   - half*sumh_ad(k,i)*heol(k,i)
           zo_ad(k+1,i)   = zo_ad(k+1,i)   + half*sumh_ad(k,i)*heol(k,i)
           sumh_ad(k-1,i) = sumh_ad(k-1,i) + sumh_ad(k,i)
           
           zo_ad(k-1,i)   = zo_ad(k-1,i)   - half*sumz_ad(k,i)
           zo_ad(k+1,i)   = zo_ad(k+1,i)   + half*sumz_ad(k,i)
           sumz_ad(k-1,i) = sumz_ad(k-1,i) + sumz_ad(k,i)
        endif
     end do
  end do

! Adjoint of cloud property below cloud base is modified by the 
! entrainment process
  do i=1,im
     do k=kmax(i)-1,2,-1

!       Adjoint of last if-then block in tlm
        if (k > kbcon(i)) then
           hcko_ad(kbcon(i),i) = hcko_ad(kbcon(i),i) + hcko_ad(k,i)
           ucko_ad(kbcon(i),i) = ucko_ad(kbcon(i),i) + ucko_ad(k,i)
           vcko_ad(kbcon(i),i) = vcko_ad(kbcon(i),i) + vcko_ad(k,i)
        endif
!
!       Adjoint of first if-then block in tlm
        if (k > kb(i).and.k <= kbcon(i)) then
!
!          Redo nlm calculations
           factor = eta0(k-1,i)/eta0(k,i)
           onemf  = one - factor
!
!          Adjoint code
           factor_ad = zero
           onemf_ad = zero

           vol_ad(k+1,i) = vol_ad(k+1,i) + onemf*half*vcko_ad(k,i)
           vol_ad(k,i)   = vol_ad(k,i)      + onemf*half*vcko_ad(k,i)
           onemf_ad = onemf_ad + vcko_ad(k,i)*half* (vol(k,i)+vol(k+1,i))
           vcko_ad(k-1,i) = vcko_ad(k-1,i) + factor*vcko_ad(k,i)
           factor_ad = factor_ad + vcko_ad(k,i)*vcko0(k-1,i)

           uol_ad(k+1,i) = uol_ad(k+1,i) + onemf*half*ucko_ad(k,i)
           uol_ad(k,i)   = uol_ad(k,i)   + onemf*half*ucko_ad(k,i)
           onemf_ad = onemf_ad + ucko_ad(k,i)*half* (uol(k,i)+uol(k+1,i))
           ucko_ad(k-1,i) = ucko_ad(k-1,i) + factor*ucko_ad(k,i)
           factor_ad = factor_ad + ucko_ad(k,i)*ucko0(k-1,i)

           heol_ad(k+1,i) = heol_ad(k+1,i) + onemf*half*hcko_ad(k,i)
           heol_ad(k,i)   = heol_ad(k,i)      + onemf*half*hcko_ad(k,i)
           onemf_ad = onemf_ad + hcko_ad(k,i)*half*(heol(k,i)+heol(k+1,i))
           hcko_ad(k-1,i) = hcko_ad(k-1,i) + factor*hcko_ad(k,i)
           factor_ad = factor_ad + hcko_ad(k,i)*hcko0(k-1,i)
               

           factor_ad = factor_ad - onemf_ad
               
           eta_ad(k,i) = eta_ad(k,i) - factor_ad* &
                eta0(k-1,i)/(eta0(k,i)*eta0(k,i))
           eta_ad(k-1,i) = eta_ad(k-1,i) + factor_ad/eta0(k,i)
        endif
     end do
  end do
!
! Adjoint of hkbo --> hcko transfer
  do i = 1,im
     indx = kb(i)
     hkbo_ad(i) = hkbo_ad(i) + hcko_ad(indx,i)
     ukbo_ad(i) = ukbo_ad(i) + ucko_ad(indx,i)
     vkbo_ad(i) = vkbo_ad(i) + vcko_ad(indx,i)
  end do
!
! Adjoint of updraft mass flux calculation
  do i = 1,im
     dz_ad = zero
     if (kb(i)==1 .and. kbcon(i) > 1) then
        dz = half*(zo(2,i)-zo(1,i))
        eta_ad(1,i) = eta_ad(1,i) + etau_ad(1,i)
        dz_ad = dz_ad - &
             eta0(2,i)*exp(-xlamb(i)*dz)*xlamb(i)* &
             eta_ad(1,i)
        xlamb_ad(i) = xlamb_ad(i) - &
             eta0(2,i)*exp(-xlamb(i)*dz)*dz*eta_ad(1,i)
        eta_ad(2,i) = eta_ad(2,i) + eta_ad(1,i)* &
             exp(-xlamb(i)*dz)
        zo_ad(1,i) = zo_ad(1,i) - half*dz_ad
        zo_ad(2,i) = zo_ad(2,i) + half*dz_ad
     endif
  end do
  do i = 1,im
     do k = 2,kbmax(i)
        dz_ad = zero
        if (k < kbcon(i).and.k >= kb(i)) then
           dz = half * (zo(k+1,i) - zo(k-1,i))
           eta_ad(k,i) = eta_ad(k,i) + etau_ad(k,i)
           dz_ad = dz_ad - eta0(k+1,i)*exp(-xlamb(i)*dz) * &
                xlamb(i)*eta_ad(k,i)
           xlamb_ad(i) = xlamb_ad(i) - &
                eta0(k+1,i)*exp(-xlamb(i)*dz)*dz*eta_ad(k,i)
           eta_ad(k+1,i) = eta_ad(k+1,i) + eta_ad(k,i)* &
                exp(-xlamb(i)*dz)
           zo_ad(k-1,i) = zo_ad(k-1,i) - half*dz_ad
           zo_ad(k+1,i) = zo_ad(k+1,i) + half*dz_ad
        endif
     end do
  end do
!
! Adjoint of determination of entrainment rate between kb and kbcon
  do i = 1,im
!
!    Recalculate nlm
     alpha = alphas
     if (nint(slimsk(i)) == 1) alpha = alphal
     if (kb(i)==1) then
        dz = half * (zo(kbcon(i),i)+zo(kbcon(i)-1,i)) - zo(1,i)
     else
        dz = half*(zo(kbcon(i),i)+zo(kbcon(i)-1,i)) - &
             half*(zo(kb(i),i)  +zo(kb(i)-1,i))
     endif
!
!    Adjoint of tlm code.
     dz_ad = zero
     if (kbcon(i)/=kb(i)) then
        dz_ad = dz_ad + two*log(alpha)/(dz*dz) * xlamb_ad(i)
     else
        dz_ad = zero
     endif
     if (kb(i)==1) then
        zo_ad(1,i)          = zo_ad(1,i)             - dz_ad
        zo_ad(kbcon(i)-1,i) = zo_ad(kbcon(i)-1,i) + half*dz_ad
        zo_ad(kbcon(i),i)   = zo_ad(kbcon(i),i)      + half*dz_ad
     else
        zo_ad(kb(i)-1,i)    = zo_ad(kb(i)-1,i)    - half*dz_ad
        zo_ad(kb(i),i)      = zo_ad(kb(i),i)         - half*dz_ad
        zo_ad(kbcon(i),i)   = zo_ad(kbcon(i),i)      + half*dz_ad
        zo_ad(kbcon(i)-1,i) = zo_ad(kbcon(i)-1,i) + half*dz_ad
     endif
  end do
!
! Adjoint of dot to pdot transfer at kbcon
  do i = 1,im
     dot0_ad(kbcon(i),i) = dot0_ad(kbcon(i),i) + 10.0_r_kind*pdot_ad(i)
  end do
!
! Adjoint of heol --> hkbo transfer
  do i = 1,im
     indx = kb(i)
     heol_ad(indx,i) = heol_ad(indx,i) + hkbo_ad(i)
     uol_ad(indx,i)  = uol_ad(indx,i) + ukbo_ad(i)
     vol_ad(indx,i)  = vol_ad(indx,i) + vkbo_ad(i)
  end do
!
! Adjoint of heol and hesol calculation at kmax(i)
  do i = 1,im
     k=kmax(i)
     heo_ad(k,i)  = heo_ad(k,i) + heol_ad(k,i)
     heso_ad(k,i) = heso_ad(k,i) + hesol_ad(k,i)
  end do
!
! Adjoint of heol and hesol calculation below kmax      
  do i = 1,im

     do k = kmax(i)-1,1,-1
!
!       Make necessary nlm calculations
        PO  = half * (P(k,i) + P(k+1,i))
        call fpvsx_ad(tol(k,i),es0,adt,es0_ad,.false.)
        esl = 10.0_r_kind*es0

!
!       Zero local adjoint variables
        es0_ad = zero
        esl_ad = zero
!
!       Adjoint of tlm code
        vo_ad(k+1,i)  = vo_ad(k+1,i) + half*vol_ad(k,i)
        vo_ad(k,i)    = vo_ad(k,i)   + half*vol_ad(k,i)

        uo_ad(k+1,i)  = uo_ad(k+1,i) + half*uol_ad(k,i)
        uo_ad(k,i)    = uo_ad(k,i)   + half*uol_ad(k,i)

        qesol_ad(k,i) = qesol_ad(k,i)+ hvap*hesol_ad(k,i)
        tol_ad(k,i)   = tol_ad(k,i)  + cp*hesol_ad(k,i)
        zo_ad(k+1,i)  = zo_ad(k+1,i) + half*grav*hesol_ad(k,i)
        zo_ad(k,i)    = zo_ad(k,i)   + half*grav*hesol_ad(k,i)

        qol_ad(k,i)   = qol_ad(k,i)  + hvap*heol_ad(k,i)
        tol_ad(k,i)   = tol_ad(k,i)  + cp*heol_ad(k,i)
        zo_ad(k+1,i)  = zo_ad(k+1,i) + half*grav*heol_ad(k,i)
        zo_ad(k,i)    = zo_ad(k,i)   + half*grav*heol_ad(k,i)

        esl_ad = esl_ad - eps*esl/((po+epsm1*esl)**2) * &
             epsm1 * qesol_ad(k,i)
        esl_ad = esl_ad + eps*qesol_ad(k,i)/(po+epsm1*esl)

        es0_ad = es0_ad + 10*esl_ad

        adt = zero
        call fpvsx_ad(tol(k,i),es0,adt,es0_ad,adjoint)
        tol_ad(k,i) = tol_ad(k,i) + adt

     end do
  end do


!
! Adjoint of tol and qol calculation
  do i = 1,im
     do k = kmax(i)-1,1,-1
!
!       Make necessary nlm calculations
        DZ    = half * (ZO(k+1,i) - ZO(k,i))
        DP    = half * (P(k+1,i) - P(k,i))
!       ES    = 10.0_r_kind * FESB(TO(k+1,i))
        call fpvsx_ad(to(k+1,i),es0,adt,es0_ad,.false.)
        es = 10.0_r_kind*es0
        PPRIME= P(k+1,i) + EPSM1 * ES
        QS    = EPS * ES / PPRIME
        DQSDP = - QS / PPRIME
        DESDT = ES * (FACTOR1 / TO(k+1,i) + FACTOR2 / (TO(k+1,i)**2))
        DQSDT = QS * P(k+1,i) * DESDT / (ES * PPRIME)
        GAMMA = EL2ORC * QESO(k+1,i) / (TO(k+1,i)**2)
        DT    = (GRAV * DZ + HVAP * DQSDP * DP) / (CP * (one + GAMMA))
!
!       Initialize local ajm variables to zero
        es0_ad    = zero
        dq_ad     = zero
        dt_ad     = zero
        dqsdp_ad  = zero
        dqsdt_ad  = zero
        gamma_ad  = zero
        dz_ad     = zero
        pprime_ad = zero
        es_ad     = zero
        desdt_ad  = zero
        qs_ad     = zero
!
!       Adjoint of tlm code
        if (qol0(k,i)>zero) then
           qol0_ad(k,i) = qol0_ad(k,i) + qol_ad(k,i)
        else
           qol0_ad(k,i) = zero
        endif

        qo_ad(k+1,i) = qo_ad(k+1,i) + qol0_ad(k,i)
        dq_ad = dq_ad + qol0_ad(k,i)

        to_ad(k+1,i) = to_ad(k+1,i) + tol_ad(k,i)
        dt_ad = dt_ad + tol_ad(k,i)

        dqsdp_ad = dqsdp_ad + dq_ad*dp
        dt_ad    = dt_ad + dqsdt*dq_ad
        dqsdt_ad = dqsdt_ad + dq_ad*dt

        gamma_ad = gamma_ad - &
             (grav*dz+hvap*dqsdp*dp)/((cp*(one+gamma))**2) * &
             cp*dt_ad
        dqsdp_ad = dqsdp_ad + hvap*dt_ad*dp/(cp*(one+gamma))
        dz_ad    = dz_ad + grav*dt_ad/(cp*(one+gamma))

        to_ad(k+1,i) = to_ad(k+1,i) - &
             two*el2orc*qeso(k+1,i)/(to(k+1,i)**3) * gamma_ad
        qeso_ad(k+1,i) = qeso_ad(k,i) + & 
             el2orc*gamma_ad/(to(k+1,i)**2)

        rterm = one/(es*pprime)
        pprime_ad = pprime_ad - (qs*p(k+1,i)*desdt)*rterm**2 * es*dqsdt_ad
        es_ad     = es_ad - (qs*p(k+1,i)*desdt)*rterm**2 * dqsdt_ad*pprime
        desdt_ad  = desdt_ad + p(k+1,i)*qs*dqsdt_ad*rterm
        qs_ad     = qs_ad + p(k+1,i)*dqsdt_ad*desdt*rterm

        to_ad(k+1,i) = to_ad(k+1,i) + desdt_ad * es * &
             (-one*factor1/(to(k+1,i)**2) - two*factor2/(to(k+1,i)**3))
        es_ad = es_ad + desdt_ad * &
             (factor1/to(k+1,i) + factor2/(to(k+1,i)**2))

        pprime_ad = pprime_ad + (qs/(pprime**2))*dqsdp_ad
        qs_ad     = qs_ad - dqsdp_ad/pprime

        pprime_ad = pprime_ad - (eps*es/(pprime**2))*qs_ad
        es_ad     = es_ad + eps*qs_ad/pprime
 
        es_ad = es_ad + epsm1*pprime_ad

        es0_ad = es0_ad + 10.0_r_kind*es_ad
        adt=zero
        call fpvsx_ad(to(k+1,i),es0,adt,es0_ad,.false.)
        to_ad(k+1,i) = to_ad(k+1,i) + adt

        zo_ad(k,i)   = zo_ad(k,i)      - half*dz_ad
        zo_ad(k+1,i) = zo_ad(k+1,i) + half*dz_ad
 
     end do
  end do
!            
!
! Adjoint of moist static energy calculation
  do i = 1,im
     do k = 1,kmax(i)
        qeso_ad(k,i) = qeso_ad(k,i) + hvap*heso_ad(k,i)
        to_ad(k,i)   = to_ad(k,i) + cp*heso_ad(k,i)
        zo_ad(k,i)   = zo_ad(k,i) + grav*heso_ad(k,i)
        qo_ad(k,i)   = qo_ad(k,i) + hvap*heo_ad(k,i)
        to_ad(k,i)   = to_ad(k,i) + cp*heo_ad(k,i)
        zo_ad(k,i)   = zo_ad(k,i) + grav*heo_ad(k,i)
     end do
  end do
!
!
! Adjoint of zo calculation
  do i = 1,im
     do k = kmax(i),2,-1
        dlnsig = log(sl(k,i)/sl(k-1,i))
        term1  = dlnsig * rd / grav
        zo_ad(k-1,i) = zo_ad(k-1,i) + zo_ad(k,i)
        term2_ad = -one*term1*zo_ad(k,i)
        tvo_ad(k-1,i) = tvo_ad(k-1,i) + half*term2_ad
        tvo_ad(k,i)   = tvo_ad(k,i)      + half*term2_ad
     end do
  end do
!
  do i = 1,im
     dlnsig = log(sl(1,i))
     tvo_ad(1,i) = tvo_ad(1,i) - (dlnsig*rd/grav)*zo_ad(1,i)
  end do
!
!
! Adjoint of initial qeso and tvo calculation.
  do i = 1,im
     do k = 1,kmax(i)
        call fpvsx_ad(to(k,i),es0,adt,es0_ad,.false.)
        es = 10.0_r_kind*es0

        qo_ad(k,i) = qo_ad(k,i) + delta*to(k,i)*tvo_ad(k,i)
        to_ad(k,i) = to_ad(k,i) + delta*qo(k,i)*tvo_ad(k,i)
        to_ad(k,i) = to_ad(k,i) + tvo_ad(k,i)

        es_ad = zero
        es_ad = es_ad + qeso_ad(k,i) * &
             eps*p(k,i)/(p(k,i)+epsm1*es)**2

        es0_ad = zero
        es0_ad = 10.0_r_kind*es_ad

        adt=zero
        call fpvsx_ad(to(k,i),es0,adt,es0_ad,adjoint)
        to_ad(k,i) = to_ad(k,i) + adt

     end do
  end do
!
! Adjoint of initial 0 --> o transfer 


  do i = 1,im
     do k = 1,km
        q0_ad(k,i)   = q0_ad(k,i) + qo0_ad(k,i)
        cwm0_ad(k,i) = cwm0_ad(k,i) + cwmo_ad(k,i)
        v0_ad(k,i)   = v0_ad(k,i) + vo_ad(k,i)
        u0_ad(k,i)   = u0_ad(k,i) + uo_ad(k,i)

        if (q0(k,i)>zero) then
           q0_ad(k,i) = q0_ad(k,i) + qo_ad(k,i)
        else
           q0_ad(k,i) = zero
        endif
        t0_ad(k,i)   = t0_ad(k,i) + to_ad(k,i)
     end do
  end do

!
! End of routine
  return

! Define internal functions
contains

  real(r_kind) function ftanh(x)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ftanh
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-07  lueken - added subprogram doc block
!
!   input argument list:
!    x
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
    implicit none

    real(r_kind),intent(in   ) :: x

    ftanh = half*(one + tanh(x))
  end function ftanh

  real(r_kind) function dftanh(x)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    dftanh
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-07  lueken - added subprogram doc block
!
!   input argument list:
!    x
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
    implicit none

    real(r_kind),intent(in   ) :: x

    dftanh = half/(cosh(x)**2)
  end function dftanh

end subroutine nlmsas_ad_im_ix_

end module nlmsas_ad_mod
