subroutine turbl_ad(pges,tges,oges,u,v,prs,t,termu,termv,termt,jstart,jstop)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    turbl_ad
!
!   prgrmmr:
!
! abstract:      Adjoint of turbl_tl
!
! program history log:
!   2008-04-02  safford -- add subprogram doc block, rm unused vars and uses
!
!   input argument list:
!     pges     -
!     tges     -
!     oges     -
!     prs      -
!     t        -
!     u        -
!     v        -
!     termu    -
!     termv    -
!     termt    -
!
!   output argument list:
!     prs      -
!     t        -
!     u        -
!     v        -
!     termu    -
!     termv    -
!     termt    -
!
! attributes:
!   language:  f90
!   machine:
!
!$$$

  use constants, only: rd_over_cp,two,rd_over_g,half,zero,one,three,grav
  use kinds, only: r_kind,i_kind 
  use gridmod, only: lat2,lon2,nsig,nsig_hlf
  use turblmod, only: use_pbl
  use turblmod, only: dudz,dvdz,dodz,ri,rf,kar0my20,zi,km,kh,sm,sh
  use turblmod, only: lmix,dudtm,dvdtm,dtdtm,rdzi,rdzl
  use turblmod, only: a0my20,c0my20,d0my20, &
                      f7my20,f8my20,karmy20
  use turblmod, only: eps_m
  use turblmod, only: fsm_my20,fsh_my20
  use turblmod, only: ri_int
  implicit none


! Declare passed variables
  real(r_kind),dimension(lat2,lon2,nsig+1),intent(in   ) :: pges
  real(r_kind),dimension(lat2,lon2,nsig)  ,intent(in   ) :: tges,oges
  real(r_kind),dimension(lat2,lon2,nsig+1),intent(inout) :: prs
  real(r_kind),dimension(lat2,lon2,nsig)  ,intent(inout) :: t,u,v
  real(r_kind),dimension(lat2,lon2,nsig)  ,intent(inout) :: termu,termv,termt
  integer(i_kind)                         ,intent(in   ) :: jstart,jstop

! Declare local variables
  real(r_kind),dimension(nsig_hlf):: t_bck,o_bck
  real(r_kind),dimension(nsig_hlf):: dudz_bck,dvdz_bck,dodz_bck,ri_bck,rf_bck
  real(r_kind),dimension(nsig_hlf):: rdudz_bck,rdvdz_bck,sm_bck,sh_bck
  real(r_kind),dimension(nsig_hlf):: rdzl_bck,rdzi_bck,t_tl,u_tl,v_tl,o_tl,dzi_tl,zl_tl
  real(r_kind),dimension(nsig_hlf):: rssq,rofbck,rshbck
  real(r_kind),dimension(nsig_hlf+1):: km_bck,kh_bck,p_bck,zi_bck
  real(r_kind),dimension(nsig_hlf+1):: km_tl,kh_tl,p_tl,zi_tl
  real(r_kind),dimension(2:nsig_hlf):: dzl_tl,dodz_tl,dudz_tl,dvdz_tl,ri_tl
  real(r_kind),dimension(2:nsig_hlf):: rf_tl,sh_tl,sm_tl,lmix_tl
  integer(i_kind),dimension(nsig_hlf):: lssq
  real(r_kind):: px,rpx,a1,a2,ax,bx,ssq,rtbck
  real(r_kind):: alph,beta,gamm,gam1,gam2,alph1,alph2,delt,bet
  real(r_kind):: rdzibk,rdzlbk,kmbk,khbk,rpbck,hrdzbk,ardzbk
  real(r_kind):: kmaz_bck,khaz_bck,kmaz_tl,khaz_tl
  integer(i_kind) i,j,k

  if(.not. use_pbl)return

  do i=1,lat2
     do j=jstart,jstop

! background fields
        do k=1,nsig_hlf
           t_bck(k)=tges(i,j,k)
           p_bck(k)=pges(i,j,k)
           o_bck(k)=oges(i,j,k)
           rdzi_bck(k)=rdzi(i,j,k)
           rdzl_bck(k)=rdzl(i,j,k)
           zi_bck(k)=zi(i,j,k)
        end do
        p_bck(nsig_hlf+1) =pges(i,j,nsig_hlf+1)
        zi_bck(nsig_hlf+1)=zi  (i,j,nsig_hlf+1)

        do k=1,nsig_hlf
           dodz_bck(k)=dodz(i,j,k)
           dudz_bck(k)=dudz(i,j,k)
           dvdz_bck(k)=dvdz(i,j,k)
           ri_bck(k)=ri(i,j,k)
           rf_bck(k)=rf(i,j,k)
           sm_bck(k)=sm(i,j,k)
           sh_bck(k)=sh(i,j,k)
           km_bck(k)=km(i,j,k)
           kh_bck(k)=kh(i,j,k)
        end do
        km_bck(nsig_hlf+1)=zero
        kh_bck(nsig_hlf+1)=zero

        do k=2,nsig_hlf
           ssq=dudz_bck(k)**2+dvdz_bck(k)**2
           lssq(k)=1
           if(ssq < eps_m) then
              lssq(k)=0
              ssq = eps_m
           end if
           rssq(k)=one/ssq
           rdudz_bck(k)=dudz_bck(k)*rssq(k)
           rdvdz_bck(k)=dvdz_bck(k)*rssq(k)
           rofbck(k)=one/(one-rf_bck(k))
           rshbck(k)=one/sh_bck(k)
        end do 

! initialize perturbations

        do k=nsig_hlf,2,-1
           dzl_tl(k)=zero
           lmix_tl(k)=zero
           sm_tl(k)=zero; sh_tl(k)=zero
           ri_tl(k)=zero; rf_tl(k)=zero
           dudz_tl(k)=zero; dvdz_tl(k)=zero; dodz_tl(k)=zero
        enddo
        do k=nsig_hlf,1,-1
           dzi_tl(k)=zero
           zi_tl(k)=zero;   zl_tl(k)=zero
           km_tl(k)=zero; kh_tl(k)=zero
           t_tl(k)=zero
           o_tl(k)=zero; u_tl(k)=zero; v_tl(k)=zero; p_tl(k)=zero
        end do
        p_tl(nsig_hlf+1)=zero; zi_tl(nsig_hlf+1)=zero
        km_tl(nsig_hlf+1)=zero; kh_tl(nsig_hlf+1)=zero

        
! adjoint of update of perturbation tendencies

        do k=nsig_hlf,1,-1
           kmaz_tl=zero
           khaz_tl=zero
           rdzibk=rdzi_bck(k)
           ax=t_bck(k)/o_bck(k)
           hrdzbk=half*rdzibk
           ardzbk=ax*hrdzbk
           kmaz_bck=(km_bck(k)+km_bck(k+1))*hrdzbk
           khaz_bck=(kh_bck(k)+kh_bck(k+1))*ardzbk
           if(k>1) then
              kmaz_tl=-dudz_bck(k)*termu(i,j,k)-dvdz_bck(k)*termv(i,j,k)
              khaz_tl=-dodz_bck(k)*termt(i,j,k)
              dudz_tl(k)=-kmaz_bck*termu(i,j,k)
              dvdz_tl(k)=-kmaz_bck*termv(i,j,k)
              dodz_tl(k)=-khaz_bck*termt(i,j,k)
           end if
           if(k<nsig_hlf) then
              kmaz_tl= dudz_bck(k+1)*termu(i,j,k)+dvdz_bck(k+1)*termv(i,j,k)+kmaz_tl
              khaz_tl= dodz_bck(k+1)*termt(i,j,k)+khaz_tl
              dudz_tl(k+1)= kmaz_bck*termu(i,j,k)+dudz_tl(k+1)
              dvdz_tl(k+1)= kmaz_bck*termv(i,j,k)+dvdz_tl(k+1)
              dodz_tl(k+1)= khaz_bck*termt(i,j,k)+dodz_tl(k+1)
           end if
           dzi_tl(k)= -rdzibk*(dudtm(i,j,k)*termu(i,j,k)+&
                               dvdtm(i,j,k)*termv(i,j,k)+&
                               dtdtm(i,j,k)*termt(i,j,k) )
           if(k>1) then
              km_tl(k)  =hrdzbk*kmaz_tl
              kh_tl(k)  =ardzbk*khaz_tl
           end if
           if(k<nsig_hlf) then
              km_tl(k+1)=hrdzbk*kmaz_tl+km_tl(k+1)
              kh_tl(k+1)=ardzbk*khaz_tl+kh_tl(k+1)
           end if
        end do


! adjoint of perturbation of km and kh

        do k=nsig_hlf,2,-1
           kmbk=km_bck(k)
           khbk=kh_bck(k)
           if(ri_int(i,j,k)==1) then
              km_tl(k)=zero
              kh_tl(k)=zero
           end if
           bx=half/sm_bck(k)
           bet=-half*rofbck(k)
           ax=two/lmix(i,j,k)
           delt= kmbk*ax
           gam1= kmbk*rdudz_bck(k)
           gam2= kmbk*rdvdz_bck(k)
           beta= kmbk*bet
           alph= kmbk*three*bx
           sm_tl(k)  =km_tl(k)*alph  
           rf_tl(k)  =km_tl(k)*beta  
           dudz_tl(k)=km_tl(k)*gam1+dudz_tl(k)
           dvdz_tl(k)=km_tl(k)*gam2+dvdz_tl(k)
           lmix_tl(k)=km_tl(k)*delt          
 
           delt= khbk*ax
           gam1= khbk*rdudz_bck(k)
           gam2= khbk*rdvdz_bck(k)
           beta= khbk*bet
           alph2=khbk*rshbck(k)
           alph= khbk *bx
           sh_tl(k)  =kh_tl(k)*alph2  
           sm_tl(k)  =kh_tl(k)*alph+sm_tl(k)
           rf_tl(k)  =kh_tl(k)*beta+rf_tl(k)
           dudz_tl(k)=kh_tl(k)*gam1+dudz_tl(k)
           dvdz_tl(k)=kh_tl(k)*gam2+dvdz_tl(k)
           lmix_tl(k)=kh_tl(k)*delt+lmix_tl(k)      
        end do


! adjoint of perturbation of flux Richardson number, sh, sm and lmix

        do k=nsig_hlf,2,-1
           zi_tl(k)=karmy20/(one+kar0my20(i,j)*(zi_bck(k)-zi_bck(1)))**2 *lmix_tl(k)  
           if(ri_int(i,j,k)==1) then
              sh_tl(k)=zero
              rf_tl(k)=zero
              sm_tl(k)=zero
           end if 
           sh_tl(k)=sh_tl(k)+sm_bck(k)*rshbck(k)           *sm_tl(k) 
           rf_tl(k)=rf_tl(k)+fsm_my20*sh_bck(k)   &     
                            /(f7my20-f8my20*rf_bck(k))**2  *sm_tl(k) &
                            +fsh_my20*rofbck(k)**2         *sh_tl(k)
           ri_tl(k)=        a0my20*(one-(two*ri_bck(k)-c0my20)/ &
                   (two*sqrt(ri_bck(k)*(ri_bck(k)-c0my20)+d0my20))) *rf_tl(k)
        end do


! adjoint of perturbation of Richardson number

        do k=nsig_hlf,2,-1
           if(ri_int(i,j,k)==1) ri_tl(k) = zero
           rtbck=one/(t_bck(k)+t_bck(k-1))
           alph=-ri_bck(k)*two*lssq(k)
           alph1=alph*rdudz_bck(k)
           alph2=alph*rdvdz_bck(k)
           beta= grav*rssq(k)*rtbck*two
           gamm=-ri_bck(k)*rtbck
           dudz_tl(k)=ri_tl(k)*alph1 + dudz_tl(k)
           dvdz_tl(k)=ri_tl(k)*alph2 + dvdz_tl(k)
           dodz_tl(k)=ri_tl(k)*beta  + dodz_tl(k)
           t_tl(k)   =ri_tl(k)*gamm  + t_tl(k)
           t_tl(k-1) =ri_tl(k)*gamm 
        end do


! adjoint of perturbation of vertical derivatives


        do k=nsig_hlf,2,-1
           rdzlbk=rdzl_bck(k)
           o_tl(k  )=o_tl(k  )+rdzlbk*dodz_tl(k)
           o_tl(k-1)=         -rdzlbk*dodz_tl(k)
           v_tl(k  )=v_tl(k  )+rdzlbk*dvdz_tl(k)
           v_tl(k-1)=         -rdzlbk*dvdz_tl(k)
           u_tl(k  )=u_tl(k  )+rdzlbk*dudz_tl(k)
           u_tl(k-1)=         -rdzlbk*dudz_tl(k)
           dzl_tl(k)=         -rdzlbk*dvdz_bck(k)*dvdz_tl(k) &
                              -rdzlbk*dodz_bck(k)*dodz_tl(k) &
                              -rdzlbk*dudz_bck(k)*dudz_tl(k)
        end do


! adjoint of perturbation of heights

        do k=nsig_hlf,2,-1
           zl_tl(k  )=+dzl_tl(k)+zl_tl(k  )
           zl_tl(k-1)=-dzl_tl(k) !+zl_tl(k-1)
        end do

        do k=nsig_hlf,1,-1
           zi_tl(k+1)=zl_tl(k)*half  + zi_tl(k+1)
           zi_tl(k     )=zl_tl(k)*half  + zi_tl(k  )
        end do

        do k=nsig_hlf,1,-1
           zi_tl(k) =zi_tl(k+1)+zi_tl(k)
           dzi_tl(k)=zi_tl(k+1)+dzi_tl(k)
        end do
        zi_tl(1)=zero    


! adjoint of perturbation fields

        do k=nsig_hlf,1,-1
           rpbck=one/(p_bck(k)+p_bck(k+1))
           rpx = two*rpbck
           a1 = rd_over_g*rpx
           a2 = a1*t_bck(k)*rpx
           gamm=-a1*(p_bck(k+1)-p_bck(k))
           beta=-a2*p_bck(k)
           alph= a2*p_bck(k+1)
           px=rd_over_cp*rpbck
           alph1=o_bck(k)/t_bck(k)
           alph2=-o_bck(k)*px
           t_tl(k  )= gamm*dzi_tl(k)+o_tl(k)*alph1 + t_tl(k  )
           p_tl(k+1)= beta*dzi_tl(k)+o_tl(k)*alph2 + p_tl(k+1)
           p_tl(k     )= alph*dzi_tl(k)+o_tl(k)*alph2 
        end do


        do k=nsig_hlf,1,-1
           t(i,j,k)=t_tl(k)+t(i,j,k)
           u(i,j,k)=u_tl(k)+u(i,j,k)
           v(i,j,k)=v_tl(k)+v(i,j,k)
           prs(i,j,k)=p_tl(k)+prs(i,j,k)
        end do
        prs(i,j,nsig_hlf+1)=p_tl(nsig_hlf+1)+prs(i,j,nsig_hlf+1)

     end do 
  end do 

  return
end subroutine turbl_ad