subroutine calctends_model(u,v,t,q,oz,cw,pri,phi,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,z,mype,u_t,v_t,t_t,q_t,oz_t,cw_t,ps_t)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    calctends       calculate u,v,t,p tendencies
!   prgmmr: kleist           org: np20                date: 2005-09-29
!
! abstract: compute tendencies for pressure, wind, and virtual 
!           temperature
!
! 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
!   2006-07-31  kleist - changes to use ps instead of ln(ps)
!   2007-04-16  kleist - remove divergence tendency bits to outside
!   2007-05-08  kleist - add bits for fully generalized coordinate
!   2007-06-21  rancic - add pbl 
!   2007-07-02  derber - move calculation of z_x, z_y into routine
!   2007-07-26 cucurull - add pri in argument list, call getprs_horiz;
!                         move getprs outside calctends;
!                         remove ps from argument list
!   2007-08-08  derber - optimize, remove calculation of t_over* and dp_over* unless needed.
!   2008-06-05  safford - rm unused vars and 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
!     z        - sfc terrain height
!     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
!     pri_x    - zonal derivative of 3d pressure
!     pri_y    - meridional derivative of 3d pressure
!     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
!
!   output argument list:
!     u_t      - time tendency of u
!     v_t      - time tendency of v
!     t_t      - time tendency of t
!xxxx p_t      - time tendency of 3d prs  xxxxxxxxxxxxxxxxxx
!     ps_t     - replaced back by Rancic in order to increase efficiecy
!     q_t      - time tendency of q
!     oz_t     - time tendency of ozone
!     cw_t     - time tendency of cloud water
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds,only: r_kind,i_kind
  use gridmod, only: lat2,lon2,nsig,istart,rlats,nlat,idvc5,bk5,&
     jstart,region_lat,region_lon,eta2_ll,wrf_nmm_regional,nlon,regional,&
     corlats
  use constants, only: ione,zero,half,one,two,rearth,rd,rcp,omega,grav,cp,pi
  use tendsmod, only: what9,prsth9,r_prsum9,r_prdif9,prdif9,pr_xsum9,pr_xdif9,pr_ysum9,&
     pr_ydif9,coriolis,ctph0,stph0,tlm0
  use tends4pertmod, only: prsum9,curvfct,rdtop9
  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,phi_x,phi_y
  real(r_kind),dimension(lat2,lon2,nsig+ione)          ,intent(in   ):: pri_x,pri_y
  real(r_kind),dimension(lat2,lon2)          ,intent(in   ) :: z
  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
  real(r_kind),dimension(lat2,lon2,nsig+ione),intent(in   ) :: pri

! Declare local variables
  real(r_kind),dimension(lat2,lon2,nsig):: div
  real(r_kind),dimension(lat2,lon2,nsig):: prdif9u,prdif9v
  real(r_kind),dimension(lat2,lon2,nsig):: prdif9u_x,prdif9v_y2
  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):: rdrdx2,rdrdy2
  real(r_kind),dimension(nsig):: t_over_p,dp_over_p
  real(r_kind) tmp,tmp2,tmp9,count,count0,kap
  integer(i_kind) i,j,k,ix,jx
  real(r_kind) relm,crlm,aph,sph,cph,cc,tph
  real(r_kind) rr,dlam,dphi,wpdar
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  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
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
  real(r_kind),dimension(lat2):: dslam,rlat_deg
  real(r_kind),dimension(lon2):: rlon_deg
  integer(i_kind) unit_number
  character(len=4) ch_mype
  character(len=9) fname
!TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT

! NOTES:
!  - equations taken from NCEP Office Note 445 (Juang 2005)
!  - this is the nonlinear routine, which currently in the GSI is only
!    called based on the current guess solution.  As such, basic state
!    variables that are needed for the TLM are loaded here (in the *9
!    arrays)


  what9=zero

! 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

  rr=one/rearth**2
  dlam=two*pi/nlon
  dphi=pi/nlat      
  wpdar=  zero

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!TEST    TEST     TEST    TEST    TEST    TEST     TEST    TEST   TEST    TEST
!.............................................................................
!   write(ch_mype,'(i4.4)') mype
!   fname='dlat.'//ch_mype
!   unit_number=401+mype
!
!   do i=1,lat2
!     ix=istart(mype+1)+i-2
!     ix=min(max(ix,2),nlat-1)
!       dslam(i)=rearth*cos(rlats(ix))*dlam
!       rlat_deg(i)=rlats(ix)*180./pi
!   end do      
!   open(unit_number,file=fname,form='formatted')
!   do i=1,lat2
!     write(unit_number,1000) mype,dslam(i),rlat_deg(i),i
!   end do
! 1000 format(i6,2x,f15.4,2x,f10.4,2x,i5)
!   close(unit_number)
!   stop
!.............................................................................
!TEST    TEST     TEST    TEST    TEST    TEST     TEST    TEST   TEST    TEST
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

  do k=1,nsig
    do j=1,lon2
      do i=1,lat2
        pr_xsum9(i,j,k)=pri_x(i,j,k)+pri_x(i,j,k+ione)
        pr_xdif9(i,j,k)=pri_x(i,j,k)-pri_x(i,j,k+ione)
        pr_ysum9(i,j,k)=pri_y(i,j,k)+pri_y(i,j,k+ione)
        pr_ydif9(i,j,k)=pri_y(i,j,k)-pri_y(i,j,k+ione)
        prdif9u(i,j,k)=prdif9(i,j,k)*u(i,j,k)
        prdif9v(i,j,k)=prdif9(i,j,k)*v(i,j,k)
      end do
    end do
  end do

! 0) Define A-grid modification for divergence 
  do k=1,nsig
    do j=1,lon2
      do i=1,lat2
        rdtop9(i,j,k)=-rd*t(i,j,k)*r_prsum9(i,j,k)
        pgf_x(i,j,k)=-pr_xsum9(i,j,k)*rdtop9(i,j,k)+phi_x(i,j,k)
        pgf_y(i,j,k)=-pr_ysum9(i,j,k)*rdtop9(i,j,k)+phi_y(i,j,k)
      end do
    end do
  end do
   
!!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMB

!!  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-1
!!      do i=1,lat2
!!        pgf1_x(i,j,k)= rdrdx2(i,j)*( rd*(t(i,j+ione,k)+t(i,j,k)) * &
!!                         (prsum9(i,j+ione,k)-prsum9(i,j,k))/   &
!!                         (prsum9(i,j+ione,k)+prsum9(i,j,k)) + &
!!                           (phi(i,j+ione,k)-phi(i,j,k)) )
!!      end do
!!    end do
!!  end do

!!  do k=1,nsig
!!    do j=1,lon2
!!      do i=1,lat2-1
!!        pgf1_y(i,j,k)= rdrdy2(i,j)*( rd*(t(i+ione,j,k)+t(i,j,k)) * &
!!                         (prsum9(i+ione,j,k)-prsum9(i,j,k))/   &
!!                         (prsum9(i+ione,j,k)+prsum9(i,j,k)) + &
!!                            (phi(i+ione,j,k)-phi(i,j,k)) )
!!      end do
!!    end do
!!  end do


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

  call mp_compact_dlon2(prdif9u,prdif9u_x,.false.,nsig,mype)
  call mp_compact_dlat2(prdif9v,prdif9v_y2,.true.,nsig,mype)

  div(:,:,:)=zero 
  do k=1,nsig
    do j=1,lon2
      do i=1,lat2
!!    do j=2,lon2-1
!!      do i=2,lat2-1
!!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMB
!!        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) ) 
!!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMME
        div(i,j,k)=div(i,j,k)+prdif9u_x(i,j,k)+prdif9v_y2(i,j,k)
      end do
    end do
  end do
   

! 1) Compute horizontal part of tendency for 3d pressure (so dps/dt is the same
!    as prsth9(i,j,1) . . . also note that at the top, dp/dt=0
!    or: prsth9(i,j,nsig+ione)=0

  do j=1,lon2
    do i=1,lat2
      prsth9(i,j,nsig+ione)=zero
    end do
  end do
  do k=nsig,1,-1
    do j=1,lon2
      do i=1,lat2
         prsth9(i,j,k)=prsth9(i,j,k+ione) - div(i,j,k)
      end do
    end do
  end do   

! 1.1) Compute horizontal part of tendency for T (needed for vertical velocity in hybrid 
!      theta coordinates)

  do k=1,nsig
    do j=1,lon2
      do i=1,lat2
         tmp9=-rdtop9(i,j,k)
         t_t(i,j,k)=-u(i,j,k)*t_x(i,j,k) - v(i,j,k)*t_y(i,j,k)
         t_t(i,j,k)=t_t(i,j,k) + tmp9*rcp * ( u(i,j,k)*pr_xsum9(i,j,k) + &
             v(i,j,k)*pr_ysum9(i,j,k) + prsth9(i,j,k) + prsth9(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 vvel subroutine

  if ( (.not.regional) .AND. (idvc5.eq.3)) then
    call getvvel(t,t_t,prsth9,prdif9,what9,1,lon2)
  else
    do k=2,nsig
      do j=1,lon2
        do i=1,lat2
          if(wrf_nmm_regional) then
            what9(i,j,k)=prsth9(i,j,k)-eta2_ll(k)*prsth9(i,j,1)
          else
            what9(i,j,k)=prsth9(i,j,k)-bk5(k)*prsth9(i,j,1)
          end if
        end do
      end do
    end do
  end if ! end if on 

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

!   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
!          what9(i,j,k)=what9(i,j,k)*rdampwt
!        end do
!      end do
!    end do

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

! 3) load actual dp/dt here now, as prsth9 is reused in 
!    what9(i,k,1) & what9(i,j,nsig+ione) = zero
!    p_t(i,j,1) is the same as the surface pressure tendency

!!mr  do k=1,nsig+ione
    do j=1,lon2
      do i=1,lat2
!!mr      p_t(i,j,k)=prsth9(i,j,k)-what9(i,j,k)
        ps_t(i,j)=prsth9(i,j,1)-what9(i,j,1)
      end do
    end do
!!mr  end do
 
  do k=1,nsig
    do j=1,lon2
      do i=1,lat2

! horizontal part of momnetum equations
        u_t(i,j,k)= &
           -u(i,j,k)*u_x(i,j,k) - v(i,j,k)*u_y(i,j,k) &
           + coriolis(i,j)*v(i,j,k)  &
           + curvfct(i,j)*(u(i,j,k)*v(i,j,k))  &
           - pgf_x(i,j,k)

        v_t(i,j,k)= &
          -u(i,j,k)*v_x(i,j,k) - v(i,j,k)*v_y(i,j,k)   &
         - coriolis(i,j)*u(i,j,k) &
         - curvfct(i,j)*(u(i,j,k)*u(i,j,k))   &
         - pgf_y(i,j,k)

! horizontal advection of "tracer" quantities
        q_t(i,j,k) = -u(i,j,k)*q_x(i,j,k) - v(i,j,k)*q_y(i,j,k)
        oz_t(i,j,k) = -u(i,j,k)*oz_x(i,j,k) - v(i,j,k)*oz_y(i,j,k)
        cw_t(i,j,k) = -u(i,j,k)*cw_x(i,j,k) - v(i,j,k)*cw_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=1,lon2
      do i=1,lat2

! vertical flux terms
        if (k.gt.1) then
          tmp = half*what9(i,j,k)*r_prdif9(i,j,k)
          u_t(i,j,k) = u_t(i,j,k) - tmp*(u(i,j,k-1)-u(i,j,k))
          v_t(i,j,k) = v_t(i,j,k) - tmp*(v(i,j,k-1)-v(i,j,k))
          t_t(i,j,k) = t_t(i,j,k) - tmp*(t(i,j,k-1)-t(i,j,k))
          q_t(i,j,k) = q_t(i,j,k) - tmp*(q(i,j,k-1)-q(i,j,k))
          oz_t(i,j,k) = oz_t(i,j,k) - tmp*(oz(i,j,k-1)-oz(i,j,k))
          cw_t(i,j,k) = cw_t(i,j,k) - tmp*(cw(i,j,k-1)-cw(i,j,k))
        end if
        if (k.lt.nsig) then
          tmp = half*what9(i,j,k+ione)*r_prdif9(i,j,k)
          u_t(i,j,k) = u_t(i,j,k) - tmp*(u(i,j,k)-u(i,j,k+ione))
          v_t(i,j,k) = v_t(i,j,k) - tmp*(v(i,j,k)-v(i,j,k+ione))
          t_t(i,j,k) = t_t(i,j,k) - tmp*(t(i,j,k)-t(i,j,k+ione))
          q_t(i,j,k) = q_t(i,j,k) - tmp*(q(i,j,k)-q(i,j,k+ione))
          oz_t(i,j,k) = oz_t(i,j,k) - tmp*(oz(i,j,k)-oz(i,j,k+ione))
          cw_t(i,j,k) = cw_t(i,j,k) - tmp*(cw(i,j,k)-cw(i,j,k+ione))
        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(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