subroutine precpd_ad( im, ix, km, dt, del, sl, ps, rhc, q_in, &
     cwm_in, t_in, q_out, cwm_out, t_out, rn_out, q_in_ad, cwm_in_ad, &
     t_in_ad, q_out_ad, cwm_out_ad, t_out_ad, rn_out_ad, &
     adjoint)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    precpd_ad    compute precipitation processes from suspended cloud water/ice
!     prgmmr:    treadon     org: np23                date: 2003-12-18
!
! abstract:  This routine parameterizes the effect of precipitation processes from
!            suspended cloud water and ice on temperature, moisture, and cloud 
!            condensate.  For more details regarding the algorithm, please refer
!            to
!               Zhao and Carr (1997), Monthly Weather Review (August)
!               Sundqvist et al., (1989) Monthly Weather review. (August)
!
!
! program history log:
!   1995-01-01  zhao - original routine
!   1998-10-10  moorthi,pan - modified and rewritten for GFS application
!   2003-12-18  treadon - add adjoint code
!   2004-06-14  treadon - reformat documenation
!   2004-08-04  treadon - add only on use declarations; add intent in/out
!   2006-04-12  treadon - change del and sl from 1d to 2d arrays
!   2008-04-24  safford - rm unused vars
!
!   input argument list:
!       im        - inner dimension over which calculation is made
!       ix        - maximum inner dimension
!       km        - number of vertical levels
!       dt        - time step in seconds
!       del       - pressure layer thickness (bottom to top)
!       sl        - sigma level
!       ps        - surface pressure (centibars)
!       rhc       - critical relative humidity threshold
!       q_in      - specific humidity
!       cwm_in    - condensate mixing ratio
!       t_in      - temperature
!       q_out_ad   - specific humidity perturbation
!       cwm_out_ad - cloud condensate mixing perturbation
!       t_out_ad   - temperature perturbation
!       rn_out_ad  - precipitation perturbation
!       adjoint   - logical flag (.false.=forward model only, .true.=forward and ajoint)
!
!   output argument list:
!       q_out     - q following grid scale precipitation
!       cwm_out   - cloud condensate mixing ratio following precipitation
!       t_out     - temperature following precipitation
!       rn_out    - precipitation over one time step, dt
!       q_in_ad    - change in rain rate with respect to moisture
!       cwm_in_ad  - change in rain rate with respect to cloud condensate mixing ratio
!       t_in_ad    - change in rain rate with respect to temperature
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
!
!***************************************************************
!==============================================
! all entries are defined explicitly
!==============================================
  use kinds, only: r_kind,i_kind
  use constants, only: cmr,cws,half,epsm1,h300,rrow,two,&
       hfus,hvap,hsub,rcp,h1000,ke2,zero,one,ttp,eps,epsq,&
       grav,climit
  implicit none
!==============================================
! define arguments
!==============================================
  logical        ,intent(in   ):: adjoint
  integer(i_kind),intent(in   ):: ix
  integer(i_kind),intent(in   ):: km
  integer(i_kind),intent(in   ):: im
  real(r_kind)   ,intent(inout):: cwm_in_ad(km,ix)
  real(r_kind)   ,intent(inout):: cwm_out_ad(km,ix)
  real(r_kind)   ,intent(inout):: q_in_ad(km,ix)
  real(r_kind)   ,intent(inout):: q_out_ad(km,ix)
  real(r_kind)   ,intent(inout):: rn_out_ad(im)
  real(r_kind)   ,intent(inout):: t_in_ad(km,ix)
  real(r_kind)   ,intent(inout):: t_out_ad(km,ix)
  real(r_kind)   ,intent(in   ):: cwm_in(km,ix)
  real(r_kind)   ,intent(  out):: cwm_out(km,ix)
  real(r_kind)   ,intent(in   ):: del(km,ix)
  real(r_kind)   ,intent(in   ):: dt
  real(r_kind)   ,intent(in   ):: ps(im)
  real(r_kind)   ,intent(in   ):: q_in(km,ix)
  real(r_kind)   ,intent(  out):: q_out(km,ix)
  real(r_kind)   ,intent(in   ):: rhc(km,ix)
  real(r_kind)   ,intent(  out):: rn_out(im)
  real(r_kind)   ,intent(in   ):: sl(km,ix)
  real(r_kind)   ,intent(in   ):: t_in(km,ix)
  real(r_kind)   ,intent(  out):: t_out(km,ix)
!==============================================
! define local variables
!==============================================
  real(r_kind) aa2
  real(r_kind) amaxcmr_ad
  real(r_kind) amaxcms_ad
  real(r_kind) amaxps_ad
  real(r_kind) amaxrq_ad
  real(r_kind) ccr_ad
  real(r_kind) cwmin_ad
  real(r_kind) cwmk_ad
  real(r_kind) cwmks_ad
  real(r_kind) erk_ad
  real(r_kind) err_ad
  real(r_kind) ers_ad
  real(r_kind) es_ad
  real(r_kind) expf_ad
  real(r_kind) fi_ad
  real(r_kind) ppr0_ad
  real(r_kind) ppr1_ad
  real(r_kind) ppr2_ad
  real(r_kind) pps0_ad
  real(r_kind) pps1_ad
  real(r_kind) pps2_ad
  real(r_kind) praut_ad
  real(r_kind) praut0_ad
  real(r_kind) precrk0_ad
  real(r_kind) precrk1_ad
  real(r_kind) precrl1k_ad
  real(r_kind) precrl_0_ad
  real(r_kind) precrl_1_ad
  real(r_kind) precrl_2_ad
  real(r_kind) precrl_3_ad
  real(r_kind) precsk0_ad
  real(r_kind) precsk1_ad
  real(r_kind) precsl1k_ad
  real(r_kind) precsl_0_ad
  real(r_kind) precsl_1_ad
  real(r_kind) precsl_2_ad
  real(r_kind) precsl_3_ad
  real(r_kind) psaci_ad
  real(r_kind) psaut_ad
  real(r_kind) psm1_ad
  real(r_kind) psm2_ad
  real(r_kind) psm3_ad
  real(r_kind) q2_ad
  real(r_kind) qin_ad
  real(r_kind) qk_ad
  real(r_kind) qs_ad
  real(r_kind) qs0_ad
  real(r_kind) rprs_ad
  real(r_kind) rq_ad
  real(r_kind) rqkll_ad
  real(r_kind) t2_ad
  real(r_kind) tem1_ad
  real(r_kind) tem2_ad
  real(r_kind) tem3_ad
  real(r_kind) tem4_ad
  real(r_kind) tem5_ad
  real(r_kind) term_ad
  real(r_kind) term1_ad
  real(r_kind) term2r_ad
  real(r_kind) term2s_ad
  real(r_kind) term3r_ad
  real(r_kind) term3s_ad
  real(r_kind) term4_ad
  real(r_kind) term5_ad
  real(r_kind) term6_ad
  real(r_kind) tin_ad
  real(r_kind) tmt0_ad
  real(r_kind) tmt0k_ad
  real(r_kind) ww0_ad
  real(r_kind) ww1_ad
  real(r_kind) ww2_ad
  real(r_kind) wwn_ad
  real(r_kind) amaxcmr
  real(r_kind) amaxcms
  real(r_kind) amaxps
  real(r_kind) amaxrq
  real(r_kind) c00
  real(r_kind) ccr
  logical comput,comput0(ix)
  real(r_kind) conde
  real(r_kind) condt
  real(r_kind) const
  real(r_kind) cr
  real(r_kind) crs1
  real(r_kind) crs2
  real(r_kind) csm1
  real(r_kind) cwmin
  real(r_kind) cwmk
  real(r_kind) cwmks
  real(r_kind) dtcp
  real(r_kind) dum1
  real(r_kind) dum2
  real(r_kind) erk
  real(r_kind) err
  real(r_kind) ers
  real(r_kind) es
  real(r_kind) expf
  real(r_kind) fi
  integer(i_kind) i
  integer(i_kind) iwl
  integer(i_kind) iwl1,iwl1k(km+1,ix)
  integer(i_kind) k
  real(r_kind) ppr0
  real(r_kind) ppr1
  real(r_kind) ppr2
  real(r_kind) pps0
  real(r_kind) pps1
  real(r_kind) pps2
  real(r_kind) praut
  real(r_kind) praut0
  real(r_kind) precrk0
  real(r_kind) precrk1
  real(r_kind) precrl1k,precrl1ki(km+1,ix)
  real(r_kind) precrl_0
  real(r_kind) precrl_1
  real(r_kind) precrl_2
  real(r_kind) precrl_3
  real(r_kind) precsk0
  real(r_kind) precsk1
  real(r_kind) precsl1k,precsl1ki(km+1,ix)
  real(r_kind) precsl_0
  real(r_kind) precsl_1
  real(r_kind) precsl_2
  real(r_kind) precsl_3
  real(r_kind) pres
  real(r_kind) psaci
  real(r_kind) psaut
  real(r_kind) psm1
  real(r_kind) psm2
  real(r_kind) psm3
  real(r_kind) q2
  real(r_kind) qin
  real(r_kind) qk
  real(r_kind) qs
  real(r_kind) qs0
  real(r_kind) rconde
  real(r_kind) rdt
  real(r_kind) rhci(km,ix)
  real(r_kind) rke2
  real(r_kind) rprs
  real(r_kind) rq
  real(r_kind) rqkll
  real(r_kind) t2
  real(r_kind) tem
  real(r_kind) tem1
  real(r_kind) tem2
  real(r_kind) tem3
  real(r_kind) tem4
  real(r_kind) tem5
  real(r_kind) term
  real(r_kind) term1
  real(r_kind) term2r
  real(r_kind) term2s
  real(r_kind) term3r
  real(r_kind) term3s
  real(r_kind) term4
  real(r_kind) term5
  real(r_kind) term6
  real(r_kind) tin
  real(r_kind) tmt0
  real(r_kind) tmt0k
  real(r_kind) wmin
  real(r_kind) ww0
  real(r_kind) ww1
  real(r_kind) ww2
  real(r_kind) wwn
  real(r_kind) zaodt
  real(r_kind) onem10,onem20

!----------------------------------------------
! RESET LOCAL ADJOINT VARIABLES
!----------------------------------------------
  onem10 = 1.0e-10_r_kind
  onem20 = 1.0e-20_r_kind
  amaxcmr_ad = zero
  amaxcms_ad = zero
  amaxps_ad = zero
  amaxrq_ad = zero
  ccr_ad = zero
  cwmin_ad = zero
  cwmk_ad = zero
  cwmks_ad = zero
  erk_ad = zero
  err_ad = zero
  ers_ad = zero
  es_ad = zero
  expf_ad = zero
  fi_ad = zero
  ppr0_ad = zero
  ppr1_ad = zero
  ppr2_ad = zero
  pps0_ad = zero
  pps1_ad = zero
  pps2_ad = zero
  praut_ad = zero
  praut0_ad = zero
  precrk0_ad = zero
  precrk1_ad = zero
  precrl1k_ad = zero
  precrl_0_ad = zero
  precrl_1_ad = zero
  precrl_2_ad = zero
  precrl_3_ad = zero
  precsk0_ad = zero
  precsk1_ad = zero
  precsl1k_ad = zero
  precsl_0_ad = zero
  precsl_1_ad = zero
  precsl_2_ad = zero
  precsl_3_ad = zero
  psaci_ad = zero
  psaut_ad = zero
  psm1_ad = zero
  psm2_ad = zero
  psm3_ad = zero
  q2_ad = zero
  qin_ad = zero
  qk_ad = zero
  qs_ad = zero
  qs0_ad = zero
  rprs_ad = zero
  rq_ad = zero
  rqkll_ad = zero
  t2_ad = zero
  tem1_ad = zero
  tem2_ad = zero
  tem3_ad = zero
  tem4_ad = zero
  tem5_ad = zero
  term_ad = zero
  term1_ad = zero
  term2r_ad = zero
  term2s_ad = zero
  term3r_ad = zero
  term3s_ad = zero
  term4_ad = zero
  term5_ad = zero
  term6_ad = zero
  tin_ad = zero
  tmt0_ad = zero
  tmt0k_ad = zero
  ww0_ad = zero
  ww1_ad = zero
  ww2_ad = zero
  wwn_ad = zero

!----------------------------------------------
! ROUTINE BODY
!----------------------------------------------
!----------------------------------------------
! FUNCTION AND TAPE COMPUTATIONS
!----------------------------------------------
  rdt = one/dt
  zaodt = 800._r_kind*rdt
  csm1 = 5.e-8_r_kind*zaodt
  crs1 = 5.e-6_r_kind*zaodt
  crs2 = 6.666e-10_r_kind*zaodt
  cr = 0.0005_r_kind*zaodt
  aa2 = 0.00125_r_kind*zaodt
  rke2 = ke2*sqrt(rdt)
  dtcp = dt*rcp
  c00 = 0.0001_r_kind*dt
  do i=1,im
     do k = 1, km
        rhci(k,i) = one/rhc(k,i)
     end do
  end do
  do i = 1, im

!new
     comput0(i) =  .false.
     do k = 1, km
        tem = 0.00001_r_kind*sl(k,i)*ps(i)*0.01_r_kind
        if (cwm_in(k,i) > tem) then
           comput0(i) =  .true.
        endif
     end do
!new

     iwl1      = 0
     precrl1k  = zero
     precsl1k  = zero
     rn_out(i) = zero
     const     = ps(i)*(h1000*dt/grav)
!new
     if (comput0(i)) then
!new
        do k = km, 1, -1
           wmin     = 0.00001_r_kind*sl(k,i)*ps(i)*0.01_r_kind
           precrl_0 = precrl1k
           precsl_0 = precsl1k
           
!          Save forward model values for use in adjoint model
           precrl1ki(k,i) = precrl_0
           precsl1ki(k,i) = precsl_0
           iwl1k(k,i)     = iwl1

           iwl   = 0
           tin   = t_in(k,i)
           qin   = q_in(k,i)
           cwmin = cwm_in(k,i)
           pres  = ps(i)*sl(k,i)
           if (precrl_0 > zero) then
              precrk0 = precrl_0
           else
              precrk0 = zero
           endif
           if (precsl_0 > zero) then
              precsk0 = precsl_0
           else
              precsk0 = zero
           endif
           if (cwmin > climit) then
              wwn = cwmin
           else
              wwn = climit
           endif
           if (wwn > climit .or. precrk0+precsk0 > zero) then
              comput =  .true. 
           else
              comput =  .false. 
           endif
           if (comput) then
              conde = const*del(k,i)
              condt = conde*rdt
              rconde = one/conde
              if (qin > epsq) then
                 qk = qin
              else
                 qk = epsq
              endif
              tmt0 = tin-ttp
              call fpvsx_ad(tin,es,dum1,dum2,.false.)
              qs0  = eps*es/(pres+epsm1*es)
              if (qs0 > epsq) then
                 qs = qs0
              else
                 qs = epsq
              endif
              if (tmt0 < (-15._r_kind)) then
                 fi = qk-rhc(k,i)*qs
                 if (fi > zero .or. wwn > climit) then
                    iwl = 1
                 else
                    iwl = 0
                 endif
              else if (tmt0 >= zero) then
                 iwl = 0
              else
                 iwl = 0
                 if (iwl1 == 1 .and. wwn > climit) then
                    iwl = 1
                 endif
              endif
              if (qs <= onem10) then
                 rq = zero
              else
                 rq = qk/qs
              endif
              if (rq < rhc(k,i)) then
                 ccr = zero
              else if (rq >= one) then
                 ccr = one
              else
                 if (rq < one) then
                    rqkll = rq
                 else
                    rqkll = one
                 endif
                 ccr = one-sqrt((one-rqkll)/(one-rhc(k,i)))
              endif
              if (ccr > zero) then
                 ww0 = cwmin
                 if (ww0 > zero) then
                    cwmk = ww0
                 else
                    cwmk = zero
                 endif
                 if (iwl == 1) then
                    term1 = cwmk-wmin
                    if (term1 > zero) then
                       amaxcms = term1
                    else
                       amaxcms = zero
                    endif
                    expf   = dt*exp(0.025_r_kind*tmt0)
                    term2s = 0.0004_r_kind*expf*amaxcms
                    if (term2s < cwmk) then
                       psaut = term2s
                    else
                       psaut = cwmk
                    endif
                    ww1 = ww0-psaut
                    if (ww1 > zero) then
                       cwmks = ww1
                    else
                       cwmks = zero
                    endif
                    term3s = aa2*expf*precsl_0*cwmks
                    if (term3s < cwmks) then
                       psaci = term3s
                    else
                       psaci = cwmks
                    endif
                    ww2      = ww1-psaci
                    precsl_1 = precsl_0+(ww0-ww2)*condt
                    precrl_1 = precrl_0
                 else
                    amaxcmr = cwmk
                    tem1    = precsl_0+precrl_0
                    term2r  = 268._r_kind-tin
                    if (term2r > zero) then
                       term3r = term2r
                    else
                       term3r = zero
                    endif
                    if (term3r < 20._r_kind) then
                       term4 = term3r
                    else
                       term4 = 20._r_kind
                    endif
                    tem2 = term4
                    tem3 = (one+h300*sqrt(tem1*rdt))*(one+half*sqrt(tem2))
                    if (ccr > 0.01_r_kind) then
                       term5 = ccr
                    else
                       term5 = 0.01_r_kind
                    endif
                    tem4  = amaxcmr*cmr*tem3/term5
                    term6 = tem4*tem4
                    if (term6 < 50._r_kind) then
                       tem5 = term6
                    else
                       tem5 = 50._r_kind
                    endif
                    praut = c00*tem3*amaxcmr*(one-exp(-tem5))
                    if (praut < cwmk) then
                       praut0 = praut
                    else
                       praut0 = cwmk
                    endif
                    ww2      = ww0-praut0
                    precrl_1 = precrl_0+(ww0-ww2)*condt
                    precsl_1 = precsl_0
                 endif
              else
                 ww2      = cwmin
                 precrl_1 = precrl_0
                 precsl_1 = precsl_0
              endif
              if (tmt0 > (-30._r_kind)) then
                 tmt0k = tmt0
              else
                 tmt0k = -30._r_kind
              endif
              if (precrl_1 > zero) then
                 precrk1 = precrl_1
              else
                 precrk1 = zero
              endif
              if (precsl_1 > zero) then
                 precsk1 = precsl_1
              else
                 precsk1 = zero
              endif
              term = rhc(k,i)-rq
              if (term > zero) then
                 amaxrq = term*conde
              else
                 amaxrq = zero
              endif
              ppr0 = rke2*amaxrq*sqrt(precrk1)
              if (tmt0 >= zero) then
                 pps0 = zero
              else
                 pps0 = (crs1+crs2*tmt0k)*amaxrq*precsk1*rhci(k,i)
              endif
              erk = precrk1+precsk1
              if (rq >= onem10) then
                 erk = amaxrq*qk*rdt/rq
              endif
              if (ppr0+pps0 > abs(erk)) then
                 rprs = erk/(precrk1+precsk1)
                 ppr1 = precrk1*rprs
                 pps1 = precsk1*rprs
              else
                 ppr1 = ppr0
                 pps1 = pps0
              endif
              if (ppr1 < precrk1) then
                 ppr2 = ppr1
              else
                 ppr2 = precrk1
              endif
              if (pps1 < precsk1) then
                 pps2 = pps1
              else
                 pps2 = precsk1
              endif
              err      = ppr2*rconde
              ers      = pps2*rconde
              precrl_2 = precrl_1-ppr2
              precsl_2 = precsl_1-pps2
              if (tmt0 > zero) then
                 if (precsl_2 > zero) then
                    amaxps = precsl_2
                 else
                    amaxps = zero
                 endif
                 psm1 = csm1*tmt0*tmt0*amaxps
                 if (ww2 > zero) then
                    psm2 = cws*cr*ww2*amaxps
                 else
                    psm2 = zero
                 endif
                 ppr0 = (psm1+psm2)*conde
                 if (ppr0 > amaxps) then
                    ppr1 = amaxps
                    psm3 = amaxps*rconde
                 else
                    ppr1 = ppr0
                    psm3 = psm1
                 endif
                 precrl_3 = precrl_2+ppr1
                 precsl_3 = precsl_2-ppr1
              else
                 psm3     = zero
                 precrl_3 = precrl_2
                 precsl_3 = precsl_2
              endif
              t2       = tin-dtcp*(hvap*err+hsub*ers+hfus*psm3)
              q2       = qin+dt*(err+ers)
           else
              t2       = tin
              q2       = qin
              precrl_3 = precrl_0
              precsl_3 = precsl_0
              ww2      = cwmin
           endif
           iwl1 = iwl
           if (precrl_3 > zero) then
              precrl1k = precrl_3
           else
              precrl1k = zero
           endif
           if (precsl_3 > zero) then
              precsl1k = precsl_3
           else
              precsl1k = zero
           endif
           if (ww2 < zero) then
              q_out(k,i)   = q2+ww2
              t_out(k,i)   = t2-hvap*rcp*ww2
              cwm_out(k,i) = zero
           else
              q_out(k,i)   = q2
              t_out(k,i)   = t2
              cwm_out(k,i) = ww2
           endif
        end do
        rn_out(i) = (precrl1k+precsl1k)*rrow

!new - no compute case
     else
        do k = km, 1, -1
           if (cwm_in(k,i) < zero) then
              q_out(k,i)   = q_in(k,i)+cwm_in(k,i)
              t_out(k,i)   = t_in(k,i)-hvap*rcp*cwm_in(k,i)
              cwm_out(k,i) = zero
           else
              q_out(k,i)   = q_in(k,i)
              t_out(k,i)   = t_in(k,i)
              cwm_out(k,i) = cwm_in(k,i)
           endif
        end do
        rn_out(i) = zero
     endif
     
!     End of loop over profiles         
  end do
  
  if (.not.adjoint) return

!----------------------------------------------
! ADJOINT COMPUTATIONS
!----------------------------------------------
  do i = im, 1, -1

     if (comput0(i)) then
        
        const        = ps(i)*(h1000*dt/grav)
        precrl1k_ad  = precrl1k_ad+rn_out_ad(i)*rrow
        precsl1k_ad  = precsl1k_ad+rn_out_ad(i)*rrow
        rn_out_ad(i) = zero
        do k = 1, km

!          Load saved forward model values
           precrl1k = precrl1ki(k,i)
           precsl1k = precsl1ki(k,i)
           iwl1     = iwl1k(k,i)
           
!          Recompute forward model values
           wmin     = 0.00001_r_kind*sl(k,i)*ps(i)*0.01_r_kind
           precrl_0 = precrl1k
           precsl_0 = precsl1k
           tin      = t_in(k,i)
           qin      = q_in(k,i)
           cwmin    = cwm_in(k,i)
           pres     = ps(i)*sl(k,i)
           if (precrl_0 > zero) then
              precrk0 = precrl_0
           else
              precrk0 = zero
           endif
           if (precsl_0 > zero) then
              precsk0 = precsl_0
           else
              precsk0 = zero
           endif
           if (cwmin > climit) then
              wwn = cwmin
           else
              wwn = climit
           endif
           if (wwn > climit .or. precrk0+precsk0 > zero) then
              comput =  .true. 
           else
              comput =  .false. 
           endif
           if (comput) then
              conde = const*del(k,i)
              condt = conde*rdt
              if (qin > epsq) then
                 qk = qin
              else
                 qk = epsq
              endif
              tmt0 = tin-ttp
              call fpvsx_ad(tin,es,dum1,dum2,.false.)
              qs0  = eps*es/(pres+epsm1*es)
              if (qs0 > epsq) then
                 qs = qs0
              else
                 qs = epsq
              endif
              if (tmt0 < (-15._r_kind)) then
                 fi = qk-rhc(k,i)*qs
                 if (fi > zero .or. wwn > climit) then
                    iwl = 1
                 else
                    iwl = 0
                 endif
              else if (tmt0 >= zero) then
                 iwl = 0
              else
                 iwl = 0
                 if (iwl1 == 1 .and. wwn > climit) then
                    iwl = 1
                 endif
              endif
              if (qs <= onem10) then
                 rq = zero
              else
                 rq = qk/qs
              endif
              if (rq < rhc(k,i)) then
                 ccr = zero
              else if (rq >= one) then
                 ccr = one
              else
                 if (rq < one) then
                    rqkll = rq
                 else
                    rqkll = one
                 endif
                 ccr = one-sqrt((one-rqkll)/(one-rhc(k,i)))
              endif
              if (ccr > zero) then
                 ww0 = cwmin
                 if (ww0 > zero) then
                    cwmk = ww0
                 else
                    cwmk = zero
                 endif
                 if (iwl == 1) then
                    term1 = cwmk-wmin
                    if (term1 > zero) then
                       amaxcms = term1
                    else
                       amaxcms = zero
                    endif
                    expf   = dt*exp(0.025_r_kind*tmt0)
                    term2s = 0.0004_r_kind*expf*amaxcms
                    if (term2s < cwmk) then
                       psaut = term2s
                    else
                       psaut = cwmk
                    endif
                    ww1 = ww0-psaut
                    if (ww1 > zero) then
                       cwmks = ww1
                    else
                       cwmks = zero
                    endif
                    term3s = aa2*expf*precsl_0*cwmks
                    if (term3s < cwmks) then
                       psaci = term3s
                    else
                       psaci = cwmks
                    endif
                    ww2      = ww1-psaci
                    precsl_1 = precsl_0+(ww0-ww2)*condt
                    precrl_1 = precrl_0
                 else
                    amaxcmr  = cwmk
                    tem1     = precsl_0+precrl_0
                    term2r   = 268._r_kind-tin
                    if (term2r > zero) then
                       term3r = term2r
                    else
                       term3r = zero
                    endif
                    if (term3r < 20._r_kind) then
                       term4 = term3r
                    else
                       term4 = 20._r_kind
                    endif
                    tem2 = term4
                    tem3 = (one+h300*sqrt(tem1*rdt))*(one+half*sqrt(tem2))
                    if (ccr > 0.01_r_kind) then
                       term5 = ccr
                    else
                       term5 = 0.01_r_kind
                    endif
                    tem4  = amaxcmr*cmr*tem3/term5
                    term6 = tem4*tem4
                    if (term6 < 50._r_kind) then
                       tem5 = term6
                    else
                       tem5 = 50._r_kind
                    endif
                    praut = c00*tem3*amaxcmr*(one-exp(-tem5))
                    if (praut < cwmk) then
                       praut0 = praut
                    else
                       praut0 = cwmk
                    endif
                    ww2      = ww0-praut0
                    precrl_1 = precrl_0+(ww0-ww2)*condt
                    precsl_1 = precsl_0
                 endif
              else
                 ww2      = cwmin
                 precrl_1 = precrl_0
                 precsl_1 = precsl_0
              endif
              if (tmt0 > (-30._r_kind)) then
                 tmt0k = tmt0
              else
                 tmt0k = -30._r_kind
              endif
              if (precrl_1 > zero) then
                 precrk1 = precrl_1
              else
                 precrk1 = zero
              endif
              if (precsl_1 > zero) then
                 precsk1 = precsl_1
              else
                 precsk1 = zero
              endif
              term = rhc(k,i)-rq
              if (term > zero) then
                 amaxrq = term*conde
              else
                 amaxrq = zero
              endif
!            ppr0 = rke2*amaxrq*sqrt(precrk1)
              if (precrk1>zero) then
                 ppr0 = rke2*amaxrq*sqrt(precrk1)
              else
                 ppr0 = zero
              endif
              
              if (tmt0 >= zero) then
                 pps0 = zero
              else
                 pps0 = (crs1+crs2*tmt0k)*amaxrq*precsk1*rhci(k,i)
              endif
              erk = precrk1+precsk1
              if (rq >= onem10) then
                 erk = amaxrq*qk*rdt/rq
              endif
              if (ppr0+pps0 > abs(erk)) then
                 rprs = erk/(precrk1+precsk1)
                 ppr1 = precrk1*rprs
                 pps1 = precsk1*rprs
              else
                 ppr1 = ppr0
                 pps1 = pps0
              endif
              if (ppr1 < precrk1) then
                 ppr2 = ppr1
              else
                 ppr2 = precrk1
              endif
              if (pps1 < precsk1) then
                 pps2 = pps1
              else
                 pps2 = precsk1
              endif
              precrl_2 = precrl_1-ppr2
              precsl_2 = precsl_1-pps2
              if (tmt0 > zero) then
                 if (precsl_2 > zero) then
                    amaxps = precsl_2
                 else
                    amaxps = zero
                 endif
                 psm1 = csm1*tmt0*tmt0*amaxps
                 if (ww2 > zero) then
                    psm2 = cws*cr*ww2*amaxps
                 else
                    psm2 = zero
                 endif
                 ppr0 = (psm1+psm2)*conde
                 if (ppr0 > amaxps) then
                    ppr1 = amaxps
                 else
                    ppr1 = ppr0
                 endif
                 precrl_3 = precrl_2+ppr1
                 precsl_3 = precsl_2-ppr1
              else
                 precrl_3 = precrl_2
                 precsl_3 = precsl_2
              endif
           else
              precrl_3 = precrl_0
              precsl_3 = precsl_0
              ww2      = cwmin
           endif
           
!          Adjoint model
           if (ww2 < zero) then
              cwm_out_ad(k,i) = zero
              t2_ad           = t2_ad+t_out_ad(k,i)
              ww2_ad          = ww2_ad-t_out_ad(k,i)*hvap*rcp
              t_out_ad(k,i)   = zero
              q2_ad           = q2_ad+q_out_ad(k,i)
              ww2_ad          = ww2_ad+q_out_ad(k,i)
              q_out_ad(k,i)   = zero
           else
              ww2_ad          = ww2_ad+cwm_out_ad(k,i)
              cwm_out_ad(k,i) = zero
              t2_ad           = t2_ad+t_out_ad(k,i)
              t_out_ad(k,i)   = zero
              q2_ad           = q2_ad+q_out_ad(k,i)
              q_out_ad(k,i)   = zero
           endif
           if (precsl_3 > zero) then
              precsl_3_ad = precsl_3_ad+precsl1k_ad
              precsl1k_ad = zero
           else
              precsl1k_ad = zero
           endif
           if (precrl_3 > zero) then
              precrl_3_ad = precrl_3_ad+precrl1k_ad
              precrl1k_ad = zero
           else
              precrl1k_ad = zero
           endif
           if (comput) then
              rconde  = one/conde
              err_ad  = err_ad+q2_ad*dt
              ers_ad  = ers_ad+q2_ad*dt
              qin_ad  = qin_ad+q2_ad
              q2_ad   = zero
              err_ad  = err_ad-t2_ad*dtcp*hvap
              ers_ad  = ers_ad-t2_ad*dtcp*hsub
              psm3_ad = psm3_ad-t2_ad*dtcp*hfus
              tin_ad  = tin_ad+t2_ad
              t2_ad   = zero
              if (tmt0 > zero) then
                 ppr1_ad     = ppr1_ad-precsl_3_ad
                 precsl_2_ad = precsl_2_ad+precsl_3_ad
                 precsl_3_ad = zero
                 ppr1_ad     = ppr1_ad+precrl_3_ad
                 precrl_2_ad = precrl_2_ad+precrl_3_ad
                 precrl_3_ad = zero
                 if (ppr0 > amaxps) then
                    amaxps_ad = amaxps_ad+psm3_ad*rconde
                    psm3_ad   = zero
                    amaxps_ad = amaxps_ad+ppr1_ad
                    ppr1_ad   = zero
                 else
                    psm1_ad   = psm1_ad+psm3_ad
                    psm3_ad   = zero
                    ppr0_ad   = ppr0_ad+ppr1_ad
                    ppr1_ad   = zero
                 endif
                 psm1_ad = psm1_ad+ppr0_ad*conde
                 psm2_ad = psm2_ad+ppr0_ad*conde
                 ppr0_ad = zero
                 if (ww2 > zero) then
                    amaxps_ad = amaxps_ad+psm2_ad*cws*cr*ww2
                    ww2_ad    = ww2_ad+psm2_ad*cws*cr*amaxps
                    psm2_ad   = zero
                 else
                    psm2_ad   = zero
                 endif
                 amaxps_ad = amaxps_ad+psm1_ad*csm1*tmt0*tmt0
                 tmt0_ad   = tmt0_ad+2*psm1_ad*csm1*tmt0*amaxps
                 psm1_ad   = zero
                 if (precsl_2 > zero) then
                    precsl_2_ad = precsl_2_ad+amaxps_ad
                    amaxps_ad   = zero
                 else
                    amaxps_ad   = zero
                 endif
              else
                 precsl_2_ad = precsl_2_ad+precsl_3_ad
                 precsl_3_ad = zero
                 precrl_2_ad = precrl_2_ad+precrl_3_ad
                 precrl_3_ad = zero
                 psm3_ad     = zero
              endif
              pps2_ad     = pps2_ad-precsl_2_ad
              precsl_1_ad = precsl_1_ad+precsl_2_ad
              precsl_2_ad = zero
              ppr2_ad     = ppr2_ad-precrl_2_ad
              precrl_1_ad = precrl_1_ad+precrl_2_ad
              precrl_2_ad = zero
              pps2_ad     = pps2_ad+ers_ad*rconde
              ers_ad      = zero
              ppr2_ad     = ppr2_ad+err_ad*rconde
              err_ad      = zero
              if (pps1 < precsk1) then
                 pps1_ad    = pps1_ad+pps2_ad
                 pps2_ad    = zero
              else
                 precsk1_ad = precsk1_ad+pps2_ad
                 pps2_ad    = zero
              endif
              ppr0 = rke2*amaxrq*sqrt(precrk1)
              if (ppr0+pps0 > abs(erk)) then
                 ppr1 = precrk1*rprs
              else
                 ppr1 = ppr0
              endif
              if (ppr1 < precrk1) then
                 ppr1_ad    = ppr1_ad+ppr2_ad
                 ppr2_ad    = zero
              else
                 precrk1_ad = precrk1_ad+ppr2_ad
                 ppr2_ad    = zero
              endif
              ppr0 = rke2*amaxrq*sqrt(precrk1)
              if (ppr0+pps0 > abs(erk)) then
                 precsk1_ad = precsk1_ad+pps1_ad*rprs
                 rprs_ad    = rprs_ad+pps1_ad*precsk1
                 pps1_ad    = zero
                 precrk1_ad = precrk1_ad+ppr1_ad*rprs
                 rprs_ad    = rprs_ad+ppr1_ad*precrk1
                 ppr1_ad    = zero
                 erk_ad     = erk_ad+rprs_ad/(precrk1+precsk1)
                 precrk1_ad = precrk1_ad-rprs_ad*(erk/((precrk1+precsk1)* &
                      (precrk1+precsk1)))
                 precsk1_ad = precsk1_ad-rprs_ad*(erk/((precrk1+precsk1)* &
                      (precrk1+precsk1)))
                 rprs_ad    = zero
              else
                 pps0_ad    = pps0_ad+pps1_ad
                 pps1_ad    = zero
                 ppr0_ad    = ppr0_ad+ppr1_ad
                 ppr1_ad    = zero
              endif
              if (rq >= onem10) then
                 amaxrq_ad = amaxrq_ad+erk_ad*(qk*rdt/rq)
                 qk_ad     = qk_ad+erk_ad*(amaxrq*rdt/rq)
                 rq_ad     = rq_ad-erk_ad*(amaxrq*qk*rdt/(rq*rq))
                 erk_ad    = zero
              endif
              precrk1_ad = precrk1_ad+erk_ad
              precsk1_ad = precsk1_ad+erk_ad
              erk_ad     = zero
              if (tmt0 >= zero) then
                 pps0_ad    = zero
              else
                 amaxrq_ad  = amaxrq_ad+pps0_ad*(crs1+crs2*tmt0k)*precsk1* &
                      rhci(k,i)
                 precsk1_ad = precsk1_ad+pps0_ad*(crs1+crs2*tmt0k)*amaxrq* &
                      rhci(k,i)
                 tmt0k_ad   = tmt0k_ad+pps0_ad*crs2*amaxrq*precsk1*rhci(k,i)
                 pps0_ad    = zero
              endif
!            amaxrq_ad  = amaxrq_ad+ppr0_ad*rke2*sqrt(precrk1)
!            precrk1_ad = precrk1_ad+ppr0_ad*rke2*amaxrq*one/(two* &
!                 sqrt(precrk1))
!            ppr0_ad    = zero
              if (precrk1>zero) then
                 amaxrq_ad  = amaxrq_ad+ppr0_ad*rke2*sqrt(precrk1)
                 precrk1_ad = precrk1_ad+ppr0_ad*rke2*amaxrq*one/(two* &
                      sqrt(precrk1))
                 ppr0_ad    = zero
              else
                 ppr0_ad    = zero
              endif
              
              if (term > zero) then
                 term_ad   = term_ad+amaxrq_ad*conde
                 amaxrq_ad = zero
              else
                 amaxrq_ad = zero
              endif
              rq_ad   = rq_ad-term_ad
              term_ad = zero
              if (precsl_1 > zero) then
                 precsl_1_ad = precsl_1_ad+precsk1_ad
                 precsk1_ad  = zero
              else
                 precsk1_ad  = zero
              endif
              if (precrl_1 > zero) then
                 precrl_1_ad = precrl_1_ad+precrk1_ad
                 precrk1_ad  = zero
              else
                 precrk1_ad  = zero
              endif
              if (tmt0 > (-30._r_kind)) then
                 tmt0_ad  = tmt0_ad+tmt0k_ad
                 tmt0k_ad = zero
              else
                 tmt0k_ad = zero
              endif
              if (ccr > zero) then
                 if (iwl == 1) then
                    precrl_0_ad = precrl_0_ad+precrl_1_ad
                    precrl_1_ad = zero
                    precsl_0_ad = precsl_0_ad+precsl_1_ad
                    ww0_ad      = ww0_ad+precsl_1_ad*condt
                    ww2_ad      = ww2_ad-precsl_1_ad*condt
                    precsl_1_ad = zero
                    psaci_ad    = psaci_ad-ww2_ad
                    ww1_ad      = ww1_ad+ww2_ad
                    ww2_ad      = zero
                    if (term3s < cwmks) then
                       term3s_ad = term3s_ad+psaci_ad
                       psaci_ad  = zero
                    else
                       cwmks_ad  = cwmks_ad+psaci_ad
                       psaci_ad  = zero
                    endif
                    cwmks_ad    = cwmks_ad+term3s_ad*aa2*expf*precsl_0
                    expf_ad     = expf_ad+term3s_ad*aa2*precsl_0*cwmks
                    precsl_0_ad = precsl_0_ad+term3s_ad*aa2*expf*cwmks
                    term3s_ad   = zero
                    if (ww1 > zero) then
                       ww1_ad   = ww1_ad+cwmks_ad
                       cwmks_ad = zero
                    else
                       cwmks_ad = zero
                    endif
                    psaut_ad = psaut_ad-ww1_ad
                    ww0_ad   = ww0_ad+ww1_ad
                    ww1_ad   = zero
                    if (term2s < cwmk) then
                       term2s_ad = term2s_ad+psaut_ad
                       psaut_ad  = zero
                    else
                       cwmk_ad   = cwmk_ad+psaut_ad
                       psaut_ad  = zero
                    endif
                    amaxcms_ad = amaxcms_ad+0.0004_r_kind*term2s_ad*expf
                    expf_ad    = expf_ad+0.0004_r_kind*term2s_ad*amaxcms
                    term2s_ad  = zero
                    tmt0_ad    = tmt0_ad+0.025_r_kind*expf_ad*dt*exp(0.025_r_kind*tmt0)
                    expf_ad    = zero
                    if (term1 > zero) then
                       term1_ad   = term1_ad+amaxcms_ad
                       amaxcms_ad = zero
                    else
                       amaxcms_ad = zero
                    endif
                    cwmk_ad     = cwmk_ad+term1_ad
                    term1_ad    = zero
                 else
                    precsl_0_ad = precsl_0_ad+precsl_1_ad
                    precsl_1_ad = zero
                    precrl_0_ad = precrl_0_ad+precrl_1_ad
                    ww0_ad      = ww0_ad+precrl_1_ad*condt
                    ww2_ad      = ww2_ad-precrl_1_ad*condt
                    precrl_1_ad = zero
                    praut0_ad   = praut0_ad-ww2_ad
                    ww0_ad      = ww0_ad+ww2_ad
                    ww2_ad      = zero
                    if (praut < cwmk) then
                       praut_ad  = praut_ad+praut0_ad
                       praut0_ad = zero
                    else
                       cwmk_ad   = cwmk_ad+praut0_ad
                       praut0_ad = zero
                    endif
                    amaxcmr_ad = amaxcmr_ad+praut_ad*c00*tem3*(one-exp(-tem5))
                    tem3_ad    = tem3_ad+praut_ad*c00*amaxcmr*(one-exp(-tem5))
                    tem5_ad    = tem5_ad+praut_ad*c00*tem3*amaxcmr*exp(-tem5)
                    praut_ad   = zero
                    if (term6 < 50._r_kind) then
                       term6_ad = term6_ad+tem5_ad
                       tem5_ad  = zero
                    else
                       tem5_ad  = zero
                    endif
                    tem4_ad    = tem4_ad+2*term6_ad*tem4
                    term6_ad   = zero
                    amaxcmr_ad = amaxcmr_ad+tem4_ad*(cmr*tem3/term5)
                    tem3_ad    = tem3_ad+tem4_ad*(amaxcmr*cmr/term5)
                    term5_ad   = term5_ad-tem4_ad*(amaxcmr*cmr*tem3/(term5*term5))
                    tem4_ad    = zero
                    if (ccr > 0.01_r_kind) then
                       ccr_ad   = ccr_ad+term5_ad
                       term5_ad = zero
                    else
                       term5_ad = zero
                    endif

!old - code below can result in divide by zero
!                tem1_ad = tem1_ad+tem3_ad*h300*one/(two*sqrt(tem1*rdt))* &
!                     rdt*(one+half*sqrt(tem2))
!                tem2_ad = tem2_ad+tem3_ad*(one+h300*sqrt(tem1*rdt))*half* &
!                     (one/(two*sqrt(tem2)))
                    if (abs(tem1)>onem20) then 
                       tem1_ad = tem1_ad+tem3_ad*h300*one/ &
                            (two*sqrt(tem1*rdt))*rdt*(one+half*sqrt(tem2))
                    else
                       tem1_ad = tem1_ad+tem3_ad*h300*one/ &
                            (two*sqrt(onem20*rdt))*rdt*(one+half*sqrt(tem2))
                    endif
                    if (abs(tem2)>onem20) then
                       tem2_ad = tem2_ad+tem3_ad*(one+h300*sqrt(tem1*rdt))* &
                            half*(one/(two*sqrt(tem2)))
                    else
                       tem2_ad = tem2_ad+tem3_ad*(one+h300*sqrt(tem1*rdt))* &
                            half*(one/(two*sqrt(onem20)))
                    endif
!new
                    
                    tem3_ad  = zero
                    term4_ad = term4_ad+tem2_ad
                    tem2_ad  = zero
                    if (term3r < 20._r_kind) then
                       term3r_ad = term3r_ad+term4_ad
                       term4_ad  = zero
                    else
                       term4_ad  = zero
                    endif
                    if (term2r > zero) then
                       term2r_ad = term2r_ad+term3r_ad
                       term3r_ad = zero
                    else
                       term3r_ad = zero
                    endif
                    tin_ad      = tin_ad-term2r_ad
                    term2r_ad   = zero
                    precrl_0_ad = precrl_0_ad+tem1_ad
                    precsl_0_ad = precsl_0_ad+tem1_ad
                    tem1_ad     = zero
                    cwmk_ad     = cwmk_ad+amaxcmr_ad
                    amaxcmr_ad  = zero
                 endif
                 if (ww0 > zero) then
                    ww0_ad  = ww0_ad+cwmk_ad
                    cwmk_ad = zero
                 else
                    cwmk_ad = zero
                 endif
                 cwmin_ad    = cwmin_ad+ww0_ad
                 ww0_ad      = zero
              else
                 precsl_0_ad = precsl_0_ad+precsl_1_ad
                 precsl_1_ad = zero
                 precrl_0_ad = precrl_0_ad+precrl_1_ad
                 precrl_1_ad = zero
                 cwmin_ad    = cwmin_ad+ww2_ad
                 ww2_ad      = zero
              endif
              if (rq < rhc(k,i)) then
                 ccr_ad = zero
              else if (rq >= one) then
                 ccr_ad = zero
              else
                 rqkll_ad = rqkll_ad+ccr_ad*(one/(two*sqrt((one-rqkll)/ &
                      (one-rhc(k,i))))/(one-rhc(k,i)))
                 ccr_ad   = zero
                 if (rq < one) then
                    rq_ad    = rq_ad+rqkll_ad
                    rqkll_ad = zero
                 else
                    rqkll_ad = zero
                 endif
              endif
              if (qs <= onem10) then
                 rq_ad = zero
              else
                 qk_ad = qk_ad+rq_ad/qs
                 qs_ad = qs_ad-rq_ad*(qk/(qs*qs))
                 rq_ad = zero
              endif
              if (tmt0 < (-15._r_kind)) then
                 qk_ad = qk_ad+fi_ad
                 qs_ad = qs_ad-fi_ad*rhc(k,i)
                 fi_ad = zero
              endif
              if (qs0 > epsq) then
                 qs0_ad = qs0_ad+qs_ad
                 qs_ad  = zero
              else
                 qs_ad  = zero
              endif
              es_ad   = es_ad+qs0_ad*(eps/(pres+epsm1*es)-eps*es*epsm1/ &
                   ((pres+epsm1*es)*(pres+epsm1*es)))
              qs0_ad  = zero
              call fpvsx_ad( tin,es,tin_ad,es_ad,adjoint )
              es_ad   = zero
              tin_ad  = tin_ad+tmt0_ad
              tmt0_ad = zero
              if (qin > epsq) then
                 qin_ad = qin_ad+qk_ad
                 qk_ad  = zero
              else
                 qk_ad  = zero
              endif
           else
              cwmin_ad    = cwmin_ad+ww2_ad
              ww2_ad      = zero
              precsl_0_ad = precsl_0_ad+precsl_3_ad
              precsl_3_ad = zero
              precrl_0_ad = precrl_0_ad+precrl_3_ad
              precrl_3_ad = zero
              qin_ad      = qin_ad+q2_ad
              q2_ad       = zero
              tin_ad      = tin_ad+t2_ad
              t2_ad       = zero
           endif
           if (cwmin > climit) then
              cwmin_ad = cwmin_ad+wwn_ad
              wwn_ad   = zero
           else
              wwn_ad   = zero
           endif
           if (precsl_0 > zero) then
              precsl_0_ad = precsl_0_ad+precsk0_ad
              precsk0_ad  = zero
           else
              precsk0_ad  = zero
           endif
           if (precrl_0 > zero) then
              precrl_0_ad = precrl_0_ad+precrk0_ad
              precrk0_ad  = zero
           else
              precrk0_ad  = zero
           endif
           cwm_in_ad(k,i) = cwm_in_ad(k,i)+cwmin_ad
           cwmin_ad       = zero
           q_in_ad(k,i)   = q_in_ad(k,i)+qin_ad
           qin_ad         = zero
           t_in_ad(k,i)   = t_in_ad(k,i)+tin_ad
           tin_ad         = zero
           precsl1k_ad    = precsl1k_ad+precsl_0_ad
           precsl_0_ad    = zero
           precrl1k_ad    = precrl1k_ad+precrl_0_ad
           precrl_0_ad    = zero
        end do
        
     else
        rn_out_ad(i) = zero
        do k = km, 1, -1
           if (cwm_in(k,i) < zero) then
              cwm_out_ad(k,i) = zero
              cwm_in_ad(k,i)  = cwm_in_ad(k,i)-t_out_ad(k,i)*hvap*rcp
              t_in_ad(k,i)    = t_in_ad(k,i)+t_out_ad(k,i)
              t_out_ad(k,i)   = zero
              cwm_in_ad(k,i)  = cwm_in_ad(k,i)+q_out_ad(k,i)
              q_in_ad(k,i)    = q_in_ad(k,i)+q_out_ad(k,i)
              q_out_ad(k,i)   = zero
           else
              cwm_in_ad(k,i)  = cwm_in_ad(k,i)+cwm_out_ad(k,i)
              cwm_out_ad(k,i) = zero
              t_in_ad(k,i)    = t_in_ad(k,i)+t_out_ad(k,i)
              t_out_ad(k,i)   = zero
              q_in_ad(k,i)    = q_in_ad(k,i)+q_out_ad(k,i)
              q_out_ad(k,i)   = zero
           endif
        end do
     endif
     rn_out_ad(i) = zero
     precsl1k_ad  = zero
     precrl1k_ad  = zero
  end do
  
  return
end subroutine precpd_ad


