subroutine calctends_model_tl(u,v,t,q,oz,cw,pri, &
   phi_x,phi_y,u_x,u_y,v_x,v_y,t_x,t_y,pri_x,pri_y,q_x,q_y,oz_x,oz_y,cw_x,cw_y, &
   mype,u_t,v_t,t_t,q_t,oz_t,cw_t,ps_t)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    calctends_model_tl       tlm of calctends_model
!   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
!   2010-02-24  rancic - adjust for use in 4dvar perturbation model 
!
! usage:
!   input argument list:
!     u        - zonal wind on subdomain
!     v        - meridional wind on subdomain
!     t        - virtual temperature on subdomain
!     q        - q on subdomain
!     oz       - ozone on subdomain
!     cw       - cloud water mixing ratio on subdomain 
!     pri      - 3d pressure
!     u_x      - zonal derivative of u
!     u_y      - meridional derivative of u
!     v_x      - zonal derivative of v
!     v_y      - meridional derivative of v
!     t_x      - zonal derivative of t
!     t_y      - meridional derivative of t
!!!     ps_x     - zonal derivative of ps
!!!     ps_y     - meridional derivative of ps
!     pri_x    - zonal gradient of pri
!     pri_y    - meridional gradient of por
!     q_x      - zonal derivative of q
!     q_y      - meridional derivative of q
!     oz_x     - zonal derivative of ozone
!     oz_y     - meridional derivative of ozone
!     cw_x     - zonal derivative of cloud water
!     cw_y     - meridional derivative of cloud water
!     mype     - task id
!     tracer   - logical flag if true tracer time derivatives calculated
!
!   output argument list:
!     u_t      - time tendency of u
!     v_t      - time tendency of v
!     t_t      - time tendency of t
!!!     p_t      - time tendency of 3d pressure
!     p_t      - time tendency of surface pressure
!     q_t      - time tendency of q
!     oz_t     - time tendency of ozone
!     cw_t     - time tendency of cloud water
!
!   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,nlon,nlat,rlats,idvc5,bk5,&
     wrf_nmm_regional,eta2_ll,regional,region_lat,region_lon,jstart,corlats
  use constants, only: ione,zero,half,one,two,rd,rcp,grav,rearth,pi,omega
  use tendsmod, only: coriolis,ctph0,stph0,tlm0
  use tends4pertmod, only: curvfct,time_step
  use nonlinmod, only: bck_u,bck_v,bck_tv,bck_q,bck_oz,bck_cw, &
     bck_u_lon,bck_u_lat,bck_v_lon,bck_v_lat,bck_tvlat,bck_tvlon, &
     bck_qlon,bck_qlat,bck_ozlon,bck_ozlat,bck_cwlon,bck_cwlat, &
     what_bck,prsth_bck,prsum_bck,r_prsum_bck,prdif_bck,r_prdif_bck,pr_xsum_bck,&
     pr_ysum_bck,rdtop_bck
  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2,nsig)     ,intent(in   ) :: u,v,t,u_x,u_y,v_x,v_y,&
     t_x,t_y,q,oz,cw,q_x,q_y,oz_x,oz_y,cw_x,cw_y,phi_x,phi_y
  real(r_kind),dimension(lat2,lon2,nsig+ione),intent(in   ) :: pri,pri_x,pri_y
  integer(i_kind)                            ,intent(in   ) :: mype
  real(r_kind),dimension(lat2,lon2,nsig)     ,intent(  out) :: u_t,v_t,t_t,q_t,oz_t,cw_t
  real(r_kind),dimension(lat2,lon2)          ,intent(  out) :: ps_t

! Declare local variables
  real(r_kind),dimension(lat2,lon2,nsig+ione):: prsth,what
  real(r_kind),dimension(lat2,lon2,nsig):: prsum,prdif,pr_xsum,pr_ysum
  real(r_kind),dimension(lat2,lon2,nsig):: t_thor9
!m-----------------------------------------------------------------------------B
  real(r_kind),dimension(lat2,lon2,nsig):: prdifu,prdifv
  real(r_kind),dimension(lat2,lon2,nsig):: pgf_x,pgf_y,pgf_xx,pgf_yy
  real(r_kind),dimension(lat2,lon2,nsig):: pgf1_x,pgf1_y,pgf1_xx,pgf1_yy
  real(r_kind),dimension(lat2,lon2,nsig):: rdtop,div,prdifu_x,prdifv_y2
  real(r_kind),dimension(lat2,lon2):: rdrdx2,rdrdy2
  real(r_kind) psum,psum9,psump1,psum9p1
  real(r_kind) rr,dlam,dphi,wpdar,pdif
  real(r_kind) relm,crlm,aph,sph,cph,cc,tph
!m-----------------------------------------------------------------------------E
  real(r_kind) tmp,tmp2,tmp9,tmp9u,tmp9v,tmp9t,tmp9q,tmp9oz,tmp9cw
  integer(i_kind) i,j,k,ix,jx
  integer(i_kind) :: jjstart,jjstop
  integer(i_kind) :: nth,tid,omp_get_num_threads,omp_get_thread_num

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  real(r_kind), parameter:: dts=600.,eps_damp=0.2,fac_nk_damp=2.6
  integer(i_kind) k_top,nk_damp,k_damp
  real(r_kind) dampwt,rdampwt,pihalf,arg,rnk_damp
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

! linearized about guess solution, so set it flag accordingly
  jjstart=1
  jjstop=lon2

! preliminaries:

  do k=1,nsig
    do j=jjstart,jjstop
      do i=1,lat2
        prsum(i,j,k)=pri(i,j,k)+pri(i,j,k+ione)
        prdif(i,j,k)=pri(i,j,k)-pri(i,j,k+ione)
        pr_xsum(i,j,k)=pri_x(i,j,k)+pri_x(i,j,k+ione)
        pr_ysum(i,j,k)=pri_y(i,j,k)+pri_y(i,j,k+ione)
      end do
    end do
  end do

! constants
  if(wrf_nmm_regional) then
    do j=1,lon2
      jx=j+jstart(mype+ione)-2_i_kind
      jx=min(max(1,jx),nlon)
      do i=1,lat2
        ix=istart(mype+ione)+i-2_i_kind
        ix=min(max(ix,1),nlat)
        coriolis(i,j)=two*omega*sin(region_lat(ix,jx))
        relm=region_lon(ix,jx)-tlm0
        crlm=cos(relm)
        aph=region_lat(ix,jx)
        sph=sin(aph)
        cph=cos(aph)
        cc=cph*crlm
        tph=asin(ctph0*sph-stph0*cc)
        curvfct(i,j)=tan(tph)/rearth
      end do
    end do
  else
    do j=1,lon2
      do i=1,lat2
        ix=istart(mype+ione)+i-2_i_kind
        ix=min(max(ix,2),nlat-1)
        coriolis(i,j)=corlats(ix)
        curvfct(i,j)=tan(rlats(ix))/rearth
      end do
    end do
  end if

! Preliminaries for 0)


  rr=rd/rearth**2
  dlam=two*pi/nlon
  dphi=pi/nlat      
!!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
!!  wpdar=  grav*time_step  * 0.09 !m>>>> test
!!  wpdar=  grav*time_step  * 0.0 !m>>>> test
!!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM

    do j=1,lon2
      do i=1,lat2
        ix=istart(mype+ione)+i-2_i_kind
        ix=min(max(ix,2),nlat-ione)
        rdrdx2(i,j)=rr/(cos(rlats(ix))*dlam)**2
        rdrdy2(i,j)=rr/dphi**2
      end do
    end do


! 0) TLM of A-grid modification for divergence 

  do k=1,nsig
    do j=1,lon2
      do i=1,lat2
        rdtop(i,j,k)=-rd*t(i,j,k)*r_prsum_bck(i,j,k) &
                     -rdtop_bck(i,j,k)*r_prsum_bck(i,j,k)*prsum(i,j,k)
        pgf_x(i,j,k)=-rdtop_bck(i,j,k)*pr_xsum(i,j,k) &
                     -rdtop(i,j,k)*pr_xsum_bck(i,j,k)+phi_x(i,j,k)
        pgf_y(i,j,k)=-rdtop_bck(i,j,k)*pr_ysum(i,j,k) &
                     -rdtop(i,j,k)*pr_ysum_bck(i,j,k)+phi_y(i,j,k)
      end do
    end do
  end do


!!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
!!    pgf1_x(:,:,:)=zero
!!    pgf1_y(:,:,:)=zero

!!  call mp_compact_dlon2(pgf_x,pgf_xx,.false.,nsig,mype)
!!  call mp_compact_dlat2(pgf_y,pgf_yy,.false.,nsig,mype)


!!  do k=1,nsig
!!    do j=1,lon2-ione
!!      do i=1,lat2
!!          psum9=prsum_bck(i,j,k)
!!          psum9p1=prsum_bck(i,j+ione,k)
!!          psum=prsum(i,j,k)
!!          psump1=prsum(i,j+ione,k)
!!          tmp9=one/(psum9+psum9p1)
!!        pgf1_x(i,j,k)=rdrdx2(i,j)*tmp9*( &
!!          (psum9p1-psum9)*(t(i,j+ione,k)+t(i,j,k))+&
!!          two*tmp9*(bck_tv(i,j+ione,k)+bck_tv(i,j,k))* &
!!                    (psum9*psump1-psum9p1*psum) )
!!      end do
!!    end do
!!  end do


!!  do k=1,nsig
!!    do j=1,lon2
!!      do i=1,lat2-ione
!!          psum9=prsum_bck(i,j,k)
!!          psum9p1=prsum_bck(i+ione,j,k)
!!          psum=prsum(i,j,k)
!!          psump1=prsum(i+ione,j,k)
!!          tmp9=one/(psum9+psum9p1)
!!        pgf1_y(i,j,k)=rdrdy2(i,j)*tmp9*( &
!!          (psum9p1-psum9)*(t(i+ione,j,k)+t(i,j,k))+&
!!          two*tmp9*(bck_tv(i+ione,j,k)+bck_tv(i,j,k))* &
!!                    (psum9*psump1-psum9p1*psum) )
!!      end do
!!    end do
!!  end do


!!  do k=1,nsig
!!    do j=2,lon2-ione    
!!      do i=1,lat2
!!        pgf1_xx(i,j,k)=pgf1_x(i,j,k)-pgf1_x(i,j-ione,k)
!!      end do
!!    end do
!!  end do
!!  do k=1,nsig
!!    do j=1,lon2
!!      do i=2,lat2-ione   
!!        pgf1_yy(i,j,k)=pgf1_y(i,j,k)-pgf1_y(i-ione,j,k)
!!      end do
!!    end do
!!  end do
!!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM


  do k=1,nsig
    do j=1,lon2
      do i=1,lat2
        prdifu(i,j,k)=prdif_bck(i,j,k)*u(i,j,k)+prdif(i,j,k)*bck_u(i,j,k)
        prdifv(i,j,k)=prdif_bck(i,j,k)*v(i,j,k)+prdif(i,j,k)*bck_v(i,j,k)
      end do
    end do
  end do


  call mp_compact_dlon2(prdifu,prdifu_x,.false.,nsig,mype)
  call mp_compact_dlat2(prdifv,prdifv_y2,.true.,nsig,mype)


  div(:,:,:)=zero 
  do k=1,nsig
    do j=1,lon2
      do i=1,lat2
!!    do j=2,lon2-ione
!!      do i=2,lat2-ione
!!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
!!        div(i,j,k)=wpdar*( pgf1_xx(i,j,k) + pgf1_yy(i,j,k)- &
!!                           pgf_xx(i,j,k) -  pgf_yy(i,j,k) ) 
!!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
        div(i,j,k)=div(i,j,k)+prdifu_x(i,j,k)+prdifv_y2(i,j,k)
      end do
    end do
  end do


! 1) Compute horizontal part of tendency for 3d pressure
   do j=jjstart,jjstop
     do i=1,lat2
       prsth(i,j,nsig+ione)=zero
     end do
   end do
   do k=nsig,1,-ione
     do j=jjstart,jjstop
       do i=1,lat2
         prsth(i,j,k)=prsth(i,j,k+ione) - div(i,j,k)
       end do
     end do
   end do


! 1.1) Get horizontal part of temperature tendency for vertical velocity term
  do k=1,nsig
    do j=jjstart,jjstop
      do i=1,lat2
        tmp9=-rdtop_bck(i,j,k)
        tmp =-rdtop(i,j,k)
        t_t(i,j,k)=-u(i,j,k)*bck_tvlon(i,j,k)-bck_u(i,j,k)*t_x(i,j,k) - &
                    v(i,j,k)*bck_tvlat(i,j,k)-bck_v(i,j,k)*t_y(i,j,k) + &
          tmp9*rcp*(bck_u(i,j,k)*pr_xsum(i,j,k)+u(i,j,k)*pr_xsum_bck(i,j,k) + &
                    bck_v(i,j,k)*pr_ysum(i,j,k)+v(i,j,k)*pr_ysum_bck(i,j,k) + &
                    prsth(i,j,k) + prsth(i,j,k+ione) ) + &
          tmp *rcp*(bck_u(i,j,k)*pr_xsum_bck(i,j,k)+bck_v(i,j,k)*pr_ysum_bck(i,j,k)+&
                    prsth_bck(i,j,k)+prsth_bck(i,j,k+ione) ) 
      end do
    end do
  end do


! 2) 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 'getvvel'
! subroutine

  if ( (.not.regional) .AND. (idvc5.eq.3)) then
!   1.1) Get horizontal part of temperature tendency for vertical velocity term
    do k=1,nsig
      do j=jjstart,jjstop
        do i=1,lat2
          tmp=-rdtop_bck(i,j,k)
          t_thor9(i,j,k)=-bck_u(i,j,k)*bck_tvlon(i,j,k) - &
                          bck_v(i,j,k)*bck_tvlat(i,j,k)
          t_thor9(i,j,k)=t_thor9(i,j,k) - &
	     tmp*rcp * ( bck_u(i,j,k)*pr_xsum_bck(i,j,k) + &
                         bck_v(i,j,k)*pr_ysum_bck(i,j,k) + &
			 prsth_bck(i,j,k) + prsth_bck(i,j,k+ione) )
        end do
      end do
    end do
    call getvvel_tl(t,t_t,t_thor9,prsth,prdif,what,jjstart,jjstop)
  else
    do k=2,nsig
      do j=jjstart,jjstop
        do i=1,lat2
          if(wrf_nmm_regional) then
            what(i,j,k)=prsth(i,j,k)-eta2_ll(k)*prsth(i,j,1)
          else
            what(i,j,k)=prsth(i,j,k)-bk5(k)*prsth(i,j,1)
          end if
        end do
      end do
    end do
  end if

! top/bottom boundary condition:
    do j=jjstart,jjstop
      do i=1,lat2
        what(i,j,1)=zero
        what(i,j,nsig+ione)=zero
      enddo
    enddo

   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Shiqiu Peng added implicit Rayleigh damping of veritical velocity

!  pihalf=pi*half

!   k_top=nsig
!   nk_damp=k_top/fac_nk_damp+ione
!   k_damp=k_top-nk_damp
!   rnk_damp=one/nk_damp

!    do k=2,nsig
!      do j=1,lon2
!        do i=1,lat2
!            if(k.ge.k_damp) then
!              arg=pihalf*(k-k_damp)*rnk_damp
!              dampwt=eps_damp*dts**sin(arg)**2
!              rdampwt=one/(one+dampwt)
!            else 
!              rdampwt=one
!            end if
!          what(i,j,k)=what(i,j,k)*rdampwt
!        end do
!      end do
!    end do

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

! 3) load actual dp/dt
!!  do k=1,nsig+ione
    do j=jjstart,jjstop
      do i=1,lat2
!!        p_t(i,j,k)=prsth(i,j,k)-what(i,j,k)
        ps_t(i,j)=prsth(i,j,1)-what(i,j,1)
      end do
    end do
!!  end do


! 4) Compute terms for tendencies of wind components & Temperature
  do k=1,nsig
    do j=jjstart,jjstop
      do i=1,lat2
! horizontal part of momentum equations
        u_t(i,j,k)= &
           +coriolis(i,j)*v(i,j,k)  &
           -pgf_x(i,j,k) &
           +curvfct(i,j)*(bck_u(i,j,k)*v(i,j,k) + u(i,j,k)*bck_v(i,j,k)) &
           -u(i,j,k)*bck_u_lon(i,j,k)-bck_u(i,j,k)*u_x(i,j,k)  &
           -v(i,j,k)*bck_u_lat(i,j,k)-bck_v(i,j,k)*u_y(i,j,k) 
        v_t(i,j,k)= &
           -coriolis(i,j)*u(i,j,k)  &
           -pgf_y(i,j,k) &
           -two*curvfct(i,j)*bck_u(i,j,k)*u(i,j,k) &
           -u(i,j,k)*bck_v_lon(i,j,k)-bck_u(i,j,k)*v_x(i,j,k) &
           -v(i,j,k)*bck_v_lat(i,j,k)-bck_v(i,j,k)*v_y(i,j,k) 
      end do   !end do i
    end do     !end do j
  end do       !end do k


  do k=1,nsig
    do j=jjstart,jjstop
      do i=1,lat2
! horizontal advection of "tracer" quantities
        q_t(i,j,k)=-u(i,j,k)*bck_qlon(i,j,k)-bck_u(i,j,k)*q_x(i,j,k) &
                   -v(i,j,k)*bck_qlat(i,j,k)-bck_v(i,j,k)*q_y(i,j,k)
        oz_t(i,j,k)=-u(i,j,k)*bck_ozlon(i,j,k)-bck_u(i,j,k)*oz_x(i,j,k) &
                    -v(i,j,k)*bck_ozlat(i,j,k)-bck_v(i,j,k)*oz_y(i,j,k)
        cw_t(i,j,k)=-u(i,j,k)*bck_cwlon(i,j,k)-bck_u(i,j,k)*cw_x(i,j,k) &
                    -v(i,j,k)*bck_cwlat(i,j,k)-bck_v(i,j,k)*cw_y(i,j,k)
      end do   !end do i
    end do     !end do j
  end do       !end do k

! vertical flux terms
  do k=1,nsig
    do j=jjstart,jjstop
      do i=1,lat2
          pdif=prdif(i,j,k)*r_prdif_bck(i,j,k)
        if (k > 1) then
          tmp=half*what(i,j,k)*r_prdif_bck(i,j,k)
          tmp2=half*what_bck(i,j,k)*r_prdif_bck(i,j,k)
          tmp9u=bck_u(i,j,k-ione)-bck_u(i,j,k)
          tmp9v=bck_v(i,j,k-ione)-bck_v(i,j,k)
          tmp9t=bck_tv(i,j,k-ione)-bck_tv(i,j,k)
          tmp9q=bck_q(i,j,k-ione)-bck_q(i,j,k) 
          tmp9oz=bck_oz(i,j,k-ione)-bck_oz(i,j,k) 
          tmp9cw=bck_cw(i,j,k-ione)-bck_cw(i,j,k) 
        u_t(i,j,k) = u_t(i,j,k) - tmp*tmp9u - &
                     tmp2*( (u(i,j,k-ione)-u(i,j,k))-tmp9u*pdif )   
        v_t(i,j,k) = v_t(i,j,k) - tmp*tmp9v - &
                     tmp2*( (v(i,j,k-ione)-v(i,j,k))-tmp9v*pdif )
        t_t(i,j,k) = t_t(i,j,k)-tmp*tmp9t - &
                     tmp2*( (t(i,j,k-ione)-t(i,j,k))-tmp9t*pdif )
        q_t(i,j,k)= q_t(i,j,k) - tmp*tmp9q - &
                     tmp2*( (q(i,j,k-ione)-q(i,j,k))-tmp9q*pdif )
        oz_t(i,j,k)= oz_t(i,j,k)-tmp*tmp9oz - &
                     tmp2*( (oz(i,j,k-ione)-oz(i,j,k))-tmp9oz*pdif )
        cw_t(i,j,k)= cw_t(i,j,k)-tmp*tmp9cw - &
                     tmp2*( (cw(i,j,k-ione)-cw(i,j,k))-tmp9cw*pdif )
        end if
        if (k < nsig) then
          tmp=half*what(i,j,k+ione)*r_prdif_bck(i,j,k)
          tmp2=half*what_bck(i,j,k+ione)*r_prdif_bck(i,j,k)
          tmp9u=bck_u(i,j,k)-bck_u(i,j,k+ione)
          tmp9v=bck_v(i,j,k)-bck_v(i,j,k+ione)
          tmp9t=bck_tv(i,j,k)-bck_tv(i,j,k+ione)
          tmp9q=bck_q(i,j,k)-bck_q(i,j,k+ione) 
          tmp9oz=bck_oz(i,j,k)-bck_oz(i,j,k+ione) 
          tmp9cw=bck_cw(i,j,k)-bck_cw(i,j,k+ione) 
        u_t(i,j,k)= u_t(i,j,k) - tmp*tmp9u - &
                     tmp2*( (u(i,j,k)-u(i,j,k+ione))-tmp9u*pdif )
        v_t(i,j,k)= v_t(i,j,k) - tmp*tmp9v - &
                     tmp2*( (v(i,j,k)-v(i,j,k+ione))-tmp9v*pdif )
        t_t(i,j,k)= t_t(i,j,k)- tmp*tmp9t - &
                     tmp2*( (t(i,j,k)-t(i,j,k+ione))-tmp9t*pdif )
        q_t(i,j,k)= q_t(i,j,k) - tmp*tmp9q - &
                     tmp2*( (q(i,j,k)-q(i,j,k+ione)) - tmp9q*pdif )
        oz_t(i,j,k)= oz_t(i,j,k) - tmp*tmp9oz - &
                     tmp2*( (oz(i,j,k)-oz(i,j,k+ione))-tmp9oz*pdif )
        cw_t(i,j,k)= cw_t(i,j,k)-tmp*tmp9cw - &
                     tmp2*( (cw(i,j,k)-cw(i,j,k+ione))-tmp9cw*pdif )
        end if
      end do   !end do i
    end do     !end do j
  end do       !end do k

  call hdiff(u_x,u_y,v_x,v_y,t_x,t_y,u_t,v_t,t_t,mype)
  call sfcdrag_tl(u,v,t,pri,u_t,v_t)
 
  if(.not.wrf_nmm_regional)then
    do k=1,nsig
      do j=1,lon2
        do i=1,lat2
          ix=istart(mype+ione)+i-2_i_kind
          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

  return
end subroutine calctends_model_tl