subroutine calctends_tl(fields,fields_dt,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    calctends_tl       tlm of calctends
!   prgmmr: kleist           org: np20                date: 2005-09-29
!
! abstract: TLM of routine that compute tendencies for pressure, Tv, u, v 
!
! program history log:
!   2005-09-29  kleist
!   2005-10-17  kleist - changes to improve computational efficiency
!   2005-11-21  kleist - add tracer tendencies, use new module
!   2006-04-12  treadon - replace sigi with bk5
!   2006-04-21  kleist - add divergence tendency parts
!   2006-07-31  kleist - changes to use ps instead of ln(ps)
!   2006-09-21  kleist - add rescaling to divergence tendency formulation
!   2006-10-04  rancic - correct bug in tracer advection terms
!   2007-04-16  kleist - move constraint specific items elsewhere
!   2007-05-08  kleist - add bits for fully generalized vertical coordinate
!   2007-06-21  rancic - add pbl 
!   2007-07-26  cucurull - add 3d pressure pri in argument list ;
!                          move getprs_tl outside calctends_tl;
!                          call getprs_horiz_tl;
!                          remove ps from argument list
!   2007-08-08  derber - optimize
!   2008-06-05  safford - rm unused uses
!   2009-08-20  parrish - replace curvfct with curvx, curvy.  this allows tendency computation to
!                          work for any general orthogonal coordinate.
!   2010-11-03  derber - moved threading calculations to gridmod and modified
!   2011-05-01  todling - cwmr no longer in guess-grids; use metguess bundle now
!   2011-12-02  zhu     - add safe-guard for the case when there is no entry in the metguess table
!   2012-03-11  tong    - added condition to calcuate cw_t
!   2013-10-19  todling - revamp interface (pass all in bundles); derivatives and
!                         guess fields also in bundles
!   2013-10-28  todling - rename p3d to prse
!   2014-06-05  eliu    - move location to get cw index(icw) from sv table; add condition
!                         to get pointer for cw 
!
! usage:
!   input argument list:
!     fields    - bundle holding relevant fields
!     mype      - task id
!
!   output argument list:
!     fields_dt - bundle holding related time tendencies
!
!   notes:
!     TLM check performed & verified on 2005-09-29 by d. kleist
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds,only: r_kind,i_kind
  use gridmod, only: lat2,lon2,nsig,istart,nlat,idvc5,bk5,&
      wrf_nmm_regional,nems_nmmb_regional,eta2_ll,regional,nthreads,jtstart,jtstop
  use constants, only: zero,half,two,rd,rcp
  use derivsmod, only: cwgues
  use tendsmod, only: what9,prsth9,r_prsum9,prdif9,r_prdif9,pr_xsum9,pr_xdif9,&
      pr_ysum9,pr_ydif9,curvx,curvy,coriolis
  use guess_grids, only: ntguessig,ges_teta,ges_prsi
  use gsi_metguess_mod, only: gsi_metguess_get,gsi_metguess_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundledup
  use gsi_bundlemod, only: gsi_bundleinquire
  use gsi_bundlemod, only: gsi_bundledestroy
  use mpeu_util, only: die, getindex
  use derivsmod, only: gsi_xderivative_bundle
  use derivsmod, only: gsi_yderivative_bundle
  implicit none

! Declare passed variables
  type(gsi_bundle) :: fields
  type(gsi_bundle) :: fields_dt
  integer(i_kind),intent(in) :: mype

! Declare local variables
  character(len=*),parameter::myname='calctends_tl'
  real(r_kind),dimension(:,:,:),pointer :: u_t,v_t,t_t,q_t,oz_t,cw_t
  real(r_kind),dimension(:,:,:),pointer :: p_t
  real(r_kind),dimension(:,:,:),pointer :: pri
  real(r_kind),dimension(:,:,:),pointer :: u,v,t,q,oz,cw
  real(r_kind),dimension(:,:,:),pointer :: u_x,u_y,v_x,v_y,t_x,t_y
  real(r_kind),dimension(:,:,:),pointer :: q_x,q_y,oz_x,oz_y,cw_x,cw_y
  real(r_kind),dimension(:,:)  ,pointer :: ps_x,ps_y
  real(r_kind),dimension(:,:,:),pointer :: pri_x,pri_y
  real(r_kind),dimension(lat2,lon2,nsig+1):: prsth,what
  real(r_kind),dimension(lat2,lon2,nsig):: prsum,prdif,pr_xsum,pr_xdif,&
       pr_ysum,pr_ydif
  real(r_kind),dimension(lat2,lon2):: sumkm1,sumvkm1,sum2km1,sum2vkm1
  real(r_kind),dimension(lat2,lon2,nsig):: t_thor9

  real(r_kind) tmp,tmp2,tmp3,sumk,sumvk,sum2k,sum2vk,uduvdv
  integer(i_kind) i,j,k,ix,it,kk,istatus,n_actual_clouds,ier,icw
  logical ihave_xtra_derivatives
  character(len=10),allocatable,dimension(:)::fvars2d,fvars3d

  type(gsi_bundle) :: xderivative
  type(gsi_bundle) :: yderivative

  real(r_kind),pointer,dimension(:,:,:) :: ges_u=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_v=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_tv=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_q =>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_oz=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_cwmr=>NULL()

  real(r_kind),pointer,dimension(:,:,:) :: ges_u_lon=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_v_lon=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_tv_lon=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_q_lon =>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_oz_lon=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_cw_lon=>NULL()

  real(r_kind),pointer,dimension(:,:,:) :: ges_u_lat=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_v_lat=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_tv_lat=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_q_lat =>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_oz_lat=>NULL()
  real(r_kind),pointer,dimension(:,:,:) :: ges_cw_lat=>NULL()

! linearized about guess solution, so set it flag accordingly
  it=ntguessig

  if(fields%n2d>0) allocate(fvars2d(fields%n2d))
  if(fields%n3d>0) allocate(fvars3d(fields%n3d))
  call gsi_bundleinquire (fields,'shortnames::2d',fvars2d,istatus)
  call gsi_bundleinquire (fields,'shortnames::3d',fvars3d,istatus)
  icw=getindex(fvars3d,'cw')

  ier=0
  call gsi_bundlegetpointer(fields,'u',   u,   istatus);ier=istatus+ier
  call gsi_bundlegetpointer(fields,'v',   v,   istatus);ier=istatus+ier
  call gsi_bundlegetpointer(fields,'tv',  t,   istatus);ier=istatus+ier
  call gsi_bundlegetpointer(fields,'q',   q,   istatus);ier=istatus+ier
  call gsi_bundlegetpointer(fields,'oz' , oz,  istatus);ier=istatus+ier
  if (icw>0) &  
  call gsi_bundlegetpointer(fields,'cw' , cw,  istatus);ier=istatus+ier 
  call gsi_bundlegetpointer(fields,'prse',pri, istatus);ier=istatus+ier
  if(ier/=0) then
     write(6,*) myname,': pointers not found on input, ier=', ier
     call stop2(999)
  endif

  ier=0
  call gsi_bundlegetpointer(fields_dt,'u',   u_t, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(fields_dt,'v',   v_t, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(fields_dt,'tv',  t_t, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(fields_dt,'q',   q_t, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(fields_dt,'oz' , oz_t,istatus);ier=istatus+ier
  if (icw>0) &
  call gsi_bundlegetpointer(fields_dt,'cw' , cw_t,istatus);ier=istatus+ier       
  call gsi_bundlegetpointer(fields_dt,'prse',p_t ,istatus);ier=istatus+ier
  if(ier/=0) then
     write(6,*) myname,': pointers not found on tendency, ier=', ier
     call stop2(999)
  endif

  ier=0
  call gsi_bundlegetpointer(gsi_metguess_bundle(it),'u' ,ges_u, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_metguess_bundle(it),'v' ,ges_v, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_metguess_bundle(it),'tv',ges_tv,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_metguess_bundle(it),'q', ges_q ,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_metguess_bundle(it),'oz',ges_oz,istatus);ier=istatus+ier
  if(ier/=0) then
     write(6,*) myname, ': pointers not found in met-guess, ier=', ier
     call stop2(999)
  endif

  ier=0
  call gsi_bundlegetpointer(gsi_xderivative_bundle(it),'u' ,ges_u_lon, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_xderivative_bundle(it),'v' ,ges_v_lon, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_xderivative_bundle(it),'tv',ges_tv_lon,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_xderivative_bundle(it),'q', ges_q_lon ,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_xderivative_bundle(it),'oz',ges_oz_lon,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_xderivative_bundle(it),'cw',ges_cw_lon,istatus);ier=istatus+ier
  if(ier/=0) then
     write(6,*) myname, ': pointers not found in lon-derivatives, ier=', ier
     call stop2(999)
  endif

  ier=0
  call gsi_bundlegetpointer(gsi_yderivative_bundle(it),'u' ,ges_u_lat, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_yderivative_bundle(it),'v' ,ges_v_lat, istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_yderivative_bundle(it),'tv',ges_tv_lat,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_yderivative_bundle(it),'q', ges_q_lat ,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_yderivative_bundle(it),'oz',ges_oz_lat,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(gsi_yderivative_bundle(it),'cw',ges_cw_lat,istatus);ier=istatus+ier
  if(ier/=0) then
     write(6,*) myname, ': pointers not found in lat-derivatives, ier=', ier
     call stop2(999)
  endif


! Get pointer to could water mixing ratio
  call gsi_metguess_get('clouds::3d',n_actual_clouds,ier)
  if (n_actual_clouds>0) then
     call gsi_bundlegetpointer (gsi_metguess_bundle(it),'cw',ges_cwmr,istatus)
     if (istatus/=0) then 
        if (regional) then 
           ges_cwmr => cwgues     ! temporily, revise after moist physics is ready
        else
           call die('setuppcp','cannot get pointer to cwmr, istatus =',istatus)
        end if
     end if
  else
     ges_cwmr => cwgues
  end if

! preliminaries:
  ihave_xtra_derivatives=.false. 
  ier=0
  call gsi_bundledup ( fields, xderivative, 'lon-derivatives', istatus )
  ier=ier+istatus
  call gsi_bundledup ( fields, yderivative, 'lat-derivatives', istatus )
  ier=ier+istatus
  if(ier==0) then
    ihave_xtra_derivatives=.true. 
  else
    call die(myname,': trouble creating temporary space for derivatives',ier)
  endif

  call get_derivatives( fields, xderivative, yderivative )

  call gsi_bundlegetpointer(xderivative,'ps'  , ps_x,istatus)
  call gsi_bundlegetpointer(xderivative,'u'   ,  u_x,istatus)
  call gsi_bundlegetpointer(xderivative,'v'   ,  v_x,istatus)
  call gsi_bundlegetpointer(xderivative,'tv'  ,  t_x,istatus)
  call gsi_bundlegetpointer(xderivative,'q'   ,  q_x,istatus)
  call gsi_bundlegetpointer(xderivative,'oz'  , oz_x,istatus)
  call gsi_bundlegetpointer(xderivative,'cw'  , cw_x,istatus)
  call gsi_bundlegetpointer(xderivative,'prse',pri_x,istatus)

  call gsi_bundlegetpointer(yderivative,'ps'  , ps_y,istatus)
  call gsi_bundlegetpointer(yderivative,'u'   ,  u_y,istatus)
  call gsi_bundlegetpointer(yderivative,'v'   ,  v_y,istatus)
  call gsi_bundlegetpointer(yderivative,'tv'  ,  t_y,istatus)
  call gsi_bundlegetpointer(yderivative,'q'   ,  q_y,istatus)
  call gsi_bundlegetpointer(yderivative,'oz'  , oz_y,istatus)
  call gsi_bundlegetpointer(yderivative,'cw'  , cw_y,istatus)
  call gsi_bundlegetpointer(yderivative,'prse',pri_y,istatus)

  call getprs_horiz_tl(ps_x,ps_y,pri,pri_x,pri_y)

!$omp parallel do private(i,j,k,kk,tmp,tmp2,uduvdv, &
!$omp                  tmp3,sumk,sumvk,sum2k,sum2vk,ix)
  do kk=1,nthreads

    do k=1,nsig
      do j=jtstart(kk),jtstop(kk) 
        do i=1,lat2
           prsum  (i,j,k)=pri  (i,j,k)+pri  (i,j,k+1)
           prdif  (i,j,k)=pri  (i,j,k)-pri  (i,j,k+1)
           pr_xsum(i,j,k)=pri_x(i,j,k)+pri_x(i,j,k+1)
           pr_xdif(i,j,k)=pri_x(i,j,k)-pri_x(i,j,k+1)
           pr_ysum(i,j,k)=pri_y(i,j,k)+pri_y(i,j,k+1)
           pr_ydif(i,j,k)=pri_y(i,j,k)-pri_y(i,j,k+1)
        end do
      end do
    end do

!   Compute horizontal part of tendency for 3d pressure

    do j=jtstart(kk),jtstop(kk) 
      do i=1,lat2
        prsth(i,j,nsig+1)=zero
      end do
    end do
    do k=nsig,1,-1
      do j=jtstart(kk),jtstop(kk) 
        do i=1,lat2
           prsth(i,j,k)=prsth(i,j,k+1) - ( u(i,j,k)*pr_xdif9(i,j,k) + &
              ges_u(i,j,k)*pr_xdif(i,j,k) + v(i,j,k)*pr_ydif9(i,j,k) + &
              ges_v(i,j,k)*pr_ydif(i,j,k) + &
              (u_x(i,j,k) + v_y(i,j,k))*(prdif9(i,j,k)) + &
              (ges_u_lon(i,j,k) + ges_v_lat(i,j,k))*(prdif(i,j,k)) )
        end do
      end do
    end do

!   Get horizontal part of temperature tendency for vertical velocity term

    do k=1,nsig
      do j=jtstart(kk),jtstop(kk) 
        do i=1,lat2
           tmp=rd*ges_tv(i,j,k)*r_prsum9(i,j,k)
           t_t(i,j,k)=-u(i,j,k)*ges_tv_lon(i,j,k) - ges_u(i,j,k)*t_x(i,j,k) - &
              v(i,j,k)*ges_tv_lat(i,j,k) - ges_v(i,j,k)*t_y(i,j,k) +          &
              tmp*rcp*( ges_u(i,j,k)*pr_xsum(i,j,k) + &
              u(i,j,k)*pr_xsum9(i,j,k) + &
              ges_v(i,j,k)*pr_ysum(i,j,k) + &
              v(i,j,k)*pr_ysum9(i,j,k) + &
              prsth(i,j,k) + prsth(i,j,k+1)) + &
              rcp*( ges_u(i,j,k)*pr_xsum9(i,j,k) + &
              ges_v(i,j,k)*pr_ysum9(i,j,k) + &
              prsth9(i,j,k)+prsth9(i,j,k+1) ) * &
              ( rd*t(i,j,k)*r_prsum9(i,j,k) - tmp*prsum(i,j,k)*r_prsum9(i,j,k) )
 
        end do
      end do
    end do

!   calculate vertical velocity term:  z(dp/dz) (zero at top/bottom interfaces)
!   if running global, and there is a c(k) coefficient, we call the vvel subroutine

    if ( (.not.regional) .AND. (idvc5==3)) then

!     Get horizontal part of temperature tendency for vertical velocity term

      do k=1,nsig
        do j=jtstart(kk),jtstop(kk) 
           do i=1,lat2
              tmp=rd*ges_tv(i,j,k)*r_prsum9(i,j,k)
  
              t_thor9(i,j,k)=-ges_u(i,j,k)*ges_tv_lon(i,j,k) - &
                   ges_v(i,j,k)*ges_tv_lat(i,j,k)
              t_thor9(i,j,k)=t_thor9(i,j,k) -tmp*rcp * ( ges_u(i,j,k)*pr_xsum9(i,j,k) + &
                 ges_v(i,j,k)*pr_ysum9(i,j,k) + &
                 prsth9(i,j,k) + prsth9(i,j,k+1) )
           end do
        end do
      end do
      call getvvel_tl(t,t_t,t_thor9,prsth,prdif,what,jtstart(kk),jtstop(kk))
    else
      if(wrf_nmm_regional.or.nems_nmmb_regional) then
        do k=2,nsig
           do j=jtstart(kk),jtstop(kk) 
              do i=1,lat2
                 what(i,j,k)=prsth(i,j,k)-eta2_ll(k)*prsth(i,j,1)
              end do
           end do
        end do
      else
        do k=2,nsig
           do j=jtstart(kk),jtstop(kk) 
              do i=1,lat2
                 what(i,j,k)=prsth(i,j,k)-bk5    (k)*prsth(i,j,1)
              end do
           end do
        end do
      end if
    end if

!   top/bottom boundary condition:

    do j=jtstart(kk),jtstop(kk) 
      do i=1,lat2
        what(i,j,1)=zero
        what(i,j,nsig+1)=zero
      enddo
    enddo


!   load actual dp/dt

    do k=1,nsig+1
      do j=jtstart(kk),jtstop(kk) 
        do i=1,lat2
           p_t(i,j,k)=prsth(i,j,k)-what(i,j,k)
        end do
      end do
    end do

!   before big k loop, zero out the km1 summation arrays

    do j=jtstart(kk),jtstop(kk) 
      do i=1,lat2
        sumkm1  (i,j)=zero
        sum2km1 (i,j)=zero
        sumvkm1 (i,j)=zero
        sum2vkm1(i,j)=zero
      end do
    end do

!   Compute terms for tendencies of wind components & Temperature

    do k=1,nsig
      do j=jtstart(kk),jtstop(kk) 
        do i=1,lat2
           uduvdv=two*(ges_u(i,j,k)*u(i,j,k) + ges_v(i,j,k)*v(i,j,k))
           u_t(i,j,k)=-u(i,j,k)*ges_u_lon(i,j,k) - ges_u(i,j,k)*u_x(i,j,k) - &
              v(i,j,k)*ges_u_lat(i,j,k) - ges_v(i,j,k)*u_y(i,j,k) + &
              coriolis(i,j)*v(i,j,k) + curvx(i,j)*uduvdv
           v_t(i,j,k)=-u(i,j,k)*ges_v_lon(i,j,k) - ges_u(i,j,k)*v_x(i,j,k) - &
              v(i,j,k)*ges_v_lat(i,j,k) - ges_v(i,j,k)*v_y(i,j,k) - &
              coriolis(i,j)*u(i,j,k) + curvy(i,j)*uduvdv

           tmp=rd*ges_tv(i,j,k)*r_prsum9(i,j,k)
           tmp2=rd*t(i,j,k)*r_prsum9(i,j,k)

           u_t(i,j,k) = u_t(i,j,k)-tmp*( pr_xsum(i,j,k) - &
              (prsum(i,j,k)*pr_xsum9(i,j,k)*r_prsum9(i,j,k)) ) - &
              tmp2*pr_xsum9(i,j,k)
 
           v_t(i,j,k) = v_t(i,j,k)-tmp*( pr_ysum(i,j,k) - &
              (prsum(i,j,k)*pr_ysum9(i,j,k)*r_prsum9(i,j,k)) ) - &
              tmp2*pr_ysum9(i,j,k)

!     vertical flux terms

           if (k > 1) then
              tmp=half*what(i,j,k)*r_prdif9(i,j,k)
              tmp2=half*what9(i,j,k)*r_prdif9(i,j,k)

              u_t(i,j,k) = u_t(i,j,k) - tmp*(ges_u (i,j,k-1)-ges_u (i,j,k)) - &
                 tmp2*( (u(i,j,k-1)-u(i,j,k)) - (ges_u (i,j,k-1)-ges_u (i,j,k))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))          
              v_t(i,j,k) = v_t(i,j,k) - tmp*(ges_v (i,j,k-1)-ges_v (i,j,k)) - &
                 tmp2*( (v(i,j,k-1)-v(i,j,k)) - (ges_v (i,j,k-1)-ges_v (i,j,k))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
              t_t(i,j,k) = t_t(i,j,k) - tmp*(ges_tv(i,j,k-1)-ges_tv(i,j,k)) - &
                 tmp2*( (t(i,j,k-1)-t(i,j,k)) - (ges_tv(i,j,k-1)-ges_tv(i,j,k))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
           end if
           if (k < nsig) then
              tmp=half*what(i,j,k+1)*r_prdif9(i,j,k)
              tmp2=half*what9(i,j,k+1)*r_prdif9(i,j,k)
              u_t(i,j,k) = u_t(i,j,k) - tmp*(ges_u (i,j,k)-ges_u (i,j,k+1)) - &
                 tmp2*( (u(i,j,k)-u(i,j,k+1)) - (ges_u (i,j,k)-ges_u (i,j,k+1))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
              v_t(i,j,k) = v_t(i,j,k) - tmp*(ges_v (i,j,k)-ges_v (i,j,k+1)) - &
                 tmp2*( (v(i,j,k)-v(i,j,k+1)) - (ges_v (i,j,k)-ges_v (i,j,k+1))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
              t_t(i,j,k) = t_t(i,j,k) - tmp*(ges_tv(i,j,k)-ges_tv(i,j,k+1)) - &
                 tmp2*( (t(i,j,k)-t(i,j,k+1)) - (ges_tv(i,j,k)-ges_tv(i,j,k+1))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
           end if
         end do   !end do i
      end do     !end do j

!     first sum to level k-1     

      do j=jtstart(kk),jtstop(kk) 
        do i=1,lat2
           tmp=rd*t(i,j,k)*r_prsum9(i,j,k)
           tmp2=rd*ges_tv(i,j,k)*r_prsum9(i,j,k)
           tmp3=prdif9(i,j,k)*r_prsum9(i,j,k)

           sumk = sumkm1(i,j) + ( tmp - tmp2*prsum(i,j,k)*r_prsum9(i,j,k) ) * &
              ( pr_xdif9(i,j,k) - ( pr_xsum9(i,j,k)*tmp3 ) )
           sumk = sumk + tmp2*( pr_xdif(i,j,k) - &
              (tmp3 *pr_xsum(i,j,k) + pr_xsum9(i,j,k)*( ( &
              prdif(i,j,k) - tmp3* prsum(i,j,k) )*r_prsum9(i,j,k) ) ) )

           sumvk = sumvkm1(i,j) + ( tmp - tmp2*prsum(i,j,k)*r_prsum9(i,j,k) ) * &
              ( pr_ydif9(i,j,k) - ( pr_ysum9(i,j,k)*tmp3 ) )
           sumvk = sumvk + tmp2*( pr_ydif(i,j,k) - &
              (tmp3*pr_ysum(i,j,k) + pr_ysum9(i,j,k)*( ( &
              prdif(i,j,k) - tmp3* prsum(i,j,k) )*r_prsum9(i,j,k) ) ) )

           sum2k  = sum2km1 (i,j) + t_x(i,j,k)*tmp3 + &
              ges_tv_lon(i,j,k)*( (prdif(i,j,k) - &
              tmp3*prsum(i,j,k))*r_prsum9(i,j,k))

           sum2vk = sum2vkm1(i,j) + t_y(i,j,k)*tmp3 + &
              ges_tv_lat(i,j,k)*( (prdif(i,j,k) - &
              tmp3*prsum(i,j,k))*r_prsum9(i,j,k))

           u_t(i,j,k) = u_t(i,j,k) - sumkm1 (i,j) - rd*sum2km1 (i,j) - &
              sumk  - rd*sum2k
           v_t(i,j,k) = v_t(i,j,k) - sumvkm1(i,j) - rd*sum2vkm1(i,j) - &
              sumvk - rd*sum2vk
 
! load up the km1 arrays for next k loop

           sumkm1  (i,j)=sumk
           sumvkm1 (i,j)=sumvk
           sum2km1 (i,j)=sum2k
           sum2vkm1(i,j)=sum2vk

        end do !end do i
      end do   !end do j
    end do  !end do k

    call turbl_tl(ges_prsi(1,1,1,it),ges_tv,ges_teta(1,1,1,it),&
               u,v,pri,t,u_t,v_t,t_t,jtstart(kk),jtstop(kk))

    if(.not.wrf_nmm_regional.and..not.nems_nmmb_regional)then
      do k=1,nsig

!       Zero out time derivatives at poles
        do j=jtstart(kk),jtstop(kk) 
           do i=1,lat2
              ix=istart(mype+1)+i-2
              if (ix == 1 .or. ix == nlat) then
                 u_t(i,j,k)=zero
                 v_t(i,j,k)=zero
              end if
           end do
        end do

      end do  !end do k 
    end if
    do k=1,nsig
      do j=jtstart(kk),jtstop(kk) 
        do i=1,lat2
!          tracer advection terms
           q_t (i,j,k) = -u(i,j,k)*ges_q_lon    (i,j,k) - ges_u(i,j,k)*q_x (i,j,k) - &
              v(i,j,k)*ges_q_lat    (i,j,k) - ges_v(i,j,k)*q_y (i,j,k)
           oz_t(i,j,k) = -u(i,j,k)*ges_oz_lon   (i,j,k) - ges_u(i,j,k)*oz_x(i,j,k) - &
              v(i,j,k)*ges_oz_lat   (i,j,k) - ges_v(i,j,k)*oz_y(i,j,k)
           if (icw>0) then
              cw_t(i,j,k) = -u(i,j,k)*ges_cw_lon(i,j,k) - ges_u(i,j,k)*cw_x(i,j,k) - &
                 v(i,j,k)*ges_cw_lat(i,j,k) - ges_v(i,j,k)*cw_y(i,j,k)
           end if
           if(k > 1)then
              tmp=half*what(i,j,k)*r_prdif9(i,j,k)
              tmp2=half*what9(i,j,k)*r_prdif9(i,j,k)
              q_t (i,j,k) = q_t (i,j,k) - tmp*(ges_q   (i,j,k-1)-ges_q   (i,j,k)) - &
                 tmp2*( (q (i,j,k-1)-q (i,j,k)) - (ges_q   (i,j,k-1)-ges_q   (i,j,k))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
              oz_t(i,j,k) = oz_t(i,j,k) - tmp*(ges_oz  (i,j,k-1)-ges_oz  (i,j,k)) - &
                 tmp2*( (oz(i,j,k-1)-oz(i,j,k)) - (ges_oz  (i,j,k-1)-ges_oz  (i,j,k))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
              if (icw>0) then 
                 cw_t(i,j,k) = cw_t(i,j,k) - tmp*(ges_cwmr(i,j,k-1)-ges_cwmr(i,j,k)) - &
                    tmp2*( (cw(i,j,k-1)-cw(i,j,k)) - (ges_cwmr(i,j,k-1)-ges_cwmr(i,j,k))* &
                    prdif(i,j,k)*r_prdif9(i,j,k))
              end if
           end if
           if(k < nsig)then
              tmp=half*what(i,j,k+1)*r_prdif9(i,j,k)
              tmp2=half*what9(i,j,k+1)*r_prdif9(i,j,k)
              q_t (i,j,k) = q_t (i,j,k) - tmp*(ges_q   (i,j,k)-ges_q   (i,j,k+1)) - &
                 tmp2*( (q (i,j,k)-q (i,j,k+1)) - (ges_q   (i,j,k)-ges_q   (i,j,k+1))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
              oz_t(i,j,k) = oz_t(i,j,k) - tmp*(ges_oz  (i,j,k)-ges_oz  (i,j,k+1)) - &
                 tmp2*( (oz(i,j,k)-oz(i,j,k+1)) - (ges_oz  (i,j,k)-ges_oz  (i,j,k+1))* &
                 prdif(i,j,k)*r_prdif9(i,j,k))
              if (icw>0) then
                 cw_t(i,j,k) = cw_t(i,j,k) - tmp*(ges_cwmr(i,j,k)-ges_cwmr(i,j,k+1)) - &
                    tmp2*( (cw(i,j,k)-cw(i,j,k+1)) - (ges_cwmr(i,j,k)-ges_cwmr(i,j,k+1))* &
                    prdif(i,j,k)*r_prdif9(i,j,k))
              end if
           end if
        end do
      end do
    end do

  end do
  if (ihave_xtra_derivatives) then
     deallocate(fvars2d,fvars3d)
     call gsi_bundledestroy(yderivative,ier)
     call gsi_bundledestroy(xderivative,ier)
  endif
! end of threading loop
  return
end subroutine calctends_tl