subroutine gfidi_hyb_resonan
     & (lon_dim_a,lon_dim_h,lons_lat,lat,
     &  dg,tg,ug,vg,ug2,vg2,ug_m,vg_m,ww,ww2,ww_m,rgs,dqdphi,dqdlam,qg,
     &  rcl,etamid,etaint,spdmax,deltim,
     &  dtdf,dtdl,
     &  td3_h1,td3_h2,rgt_h,ud3_h1,ud3_h2,vd3_h1,vd3_h2,
     &  sd_m,z_m,
     &  unlm,vnlm,tnlm,unla,vnla,tnla,
     &  dpdt_a,dpdt_d3_h1,dpdt_d3_h2,
     &  grad_gzlam,grad_gzphi,gz_grid,go_forward,batah,gg_tracers,tb_k,
     &  kdt,r_rt,omega_v)
!sk10302012 modifications for cont_eq_opt1
      use gfs_dyn_MACHINE, only : kind_evod
!mjr  use gfs_dyn_MACHINE, only : kind_phys
      use gfs_dyn_resol_def, only : latg2,levm1,levp1,levs,ntrac
      use gfs_dyn_coordinate_def, only : ak5,bk5,ck,dbk,sv_ecm,t_ecm,
     &                                   y_ecm
      use namelist_dynamics_def, only : ref_temp,herm_x,herm_y,herm_z,
     &                            lin_xyz,lin_xy,time_extrap_etadot,
     &                            settls_dep3ds,settls_dep3dg,
     &                            cont_eq_opt1
      use gfs_dyn_layout1, only : me
      use gfs_dyn_physcons, only : rerth => con_rerth,
     &                               rd => con_rd,
     &                               cp => con_cp,
     &                            omega => con_omega,
     &                             cvap => con_cvap,
     &                             grav => con_g,
     &                               fv => con_fvirt
      implicit none
      integer lon_dim_a,lon_dim_h,lons_lat
      integer kdt
      real coriol,rcl,rk,sinra,deltim,sinlat,speed2,omega_v
      real  ugsave,vgsave,wwsave, tref , rdtref, r_rt 
      real etamid(levs),etaint(levp1)
!     real etamid(levs),etaint(levp1),wgtw(7)
      real grad_gzlam(lon_dim_a),grad_gzphi(lon_dim_a)
      real    gz_grid(lon_dim_a)
      real spdmax(levs),
     1    dg(lon_dim_a,levs),   tg(lon_dim_a,levs),
     1  dtdf(lon_dim_a,levs), dtdl(lon_dim_a,levs),
     1  sd_m(lon_dim_a     ),  z_m(lon_dim_a,levs),
     1dpdt_a(lon_dim_a     ),
     1  unlm(lon_dim_a,levs),
     1  vnlm(lon_dim_a,levs), tnlm(lon_dim_a,levs),
     1  unla(lon_dim_a,levs),
     1  vnla(lon_dim_a,levs), tnla(lon_dim_a,levs),
     1   rgs(lon_dim_a,levs,ntrac),
     1  dqdphi(lon_dim_a), dqdlam(lon_dim_a), qg(lon_dim_a)
      real  ww(lon_dim_h,levs),ww2(lon_dim_h,levs),
     1      wfilt(lon_dim_h,levp1),
     1     ww_m(lon_dim_h,levs),
     1      zz(lon_dim_h,levp1),  
     1      ug(lon_dim_h,levs),        vg(lon_dim_h,levs),
     1     ug2(lon_dim_h,levs),       vg2(lon_dim_h,levs),
     1    ug_m(lon_dim_h,levs),      vg_m(lon_dim_h,levs),
     1   ud3_h1(lon_dim_h,levs),     vd3_h1(lon_dim_h,levs),
     1   ud3_h2(lon_dim_h,levs),     vd3_h2(lon_dim_h,levs),
     1     u_n(lon_dim_h,levs),       u_l(lon_dim_h,levs),
     1     v_n(lon_dim_h,levs),       v_l(lon_dim_h,levs),
     1     t_n(lon_dim_h,levs),       t_l(lon_dim_h,levs),
     1   dpdt_d3_h1(lon_dim_h,levs),   dpdt_d3_h2(lon_dim_h,levs),
     1   td3_h1(lon_dim_h,levs),td3_h2(lon_dim_h,levs),
     1   tb(lon_dim_a,levs),
     1   rgt_h(lon_dim_h,levs,ntrac)
      real pk5(lons_lat,levp1), dpk(lons_lat,levs),sv(levs)
      real
     &     dot(lons_lat,levp1),   
     &    advz(lons_lat,levs),      cg(lons_lat,levs),
     &      cb(lons_lat,levs),      db(lons_lat,levs),
     &    zlam(lons_lat,levs),    zphi(lons_lat,levs),
     &   worka(lons_lat,levs),   workb(lons_lat,levs),
     &   workc(lons_lat,levs),
     &    phiu(lons_lat,levs),    phiv(lons_lat,levs),
     &    uprs(lons_lat,levs),    vprs(lons_lat,levs),
     &    cofa(lons_lat,levs),    cofb(lons_lat,levs),
     &    alfa(lons_lat,levs),    rlnp(lons_lat,levs),
     &    px1u(lons_lat,levs),    px1v(lons_lat,levs),
     &    px2u(lons_lat,levs),    px2v(lons_lat,levs),
     &     px2(lons_lat,levs),
     &    px3u(lons_lat,levs),    px3v(lons_lat,levs),
     &    px4u(lons_lat,levs),    px4v(lons_lat,levs),
     &    px5u(lons_lat,levs),    px5v(lons_lat,levs),
     &    uphi(lons_lat,levs),    vphi(lons_lat,levs),
     &    sadx(lons_lat,levs),    sadk(lons_lat,levs),
     &    expq(lons_lat),         sadt(lons_lat,levs),
     &    dqdt(lons_lat),
     1    sd_n(lons_lat)     ,   z_n(lons_lat  ,levs),
     & z_n_sum(lons_lat),
     &    rdel(lons_lat,levs),   rdel2(lons_lat,levs),
     &  sumdel(lons_lat),          del(lons_lat,levs),
     &    rk1,rkr,                 dif(lons_lat)
      real clog2
      real rmin,rmax,delta,delta1,batah,batah2
      logical gg_tracers,test_mode,go_forward
      real rcoslat,coslat,tbar
      real tb_k(levs)           
      integer i,j,k,n,ifirst,lat,k_row,k_col,kk,kk1
      real, parameter :: cons0=0.d0, cons0p5=0.5d0, cons1=1.d0,
     &                   cons2=2.d0
!     save clog2,ifirst,delta,delta1
!     data ifirst /1/
!!      data wgtw /0.015625,-0.09375,0.234375,0.6875,
!!     .           0.234375,-0.09375,0.015625/

!!    data wgtw /-.03125,0.,.28125,.5,.28125,0.,-.03125/

!Moor   data wgtw /0.015625, 0.09375,0.234375,0.3125,
!    .             0.234375, 0.09375,0.015625/

      rk= rd /cp
!--------------------------------------------------
      call top2bot_array(lon_dim_a,lons_lat,dg)
      call top2bot_array(lon_dim_a,lons_lat,tg)
      call top2bot_array(lon_dim_h,lons_lat,ug)
      call top2bot_array(lon_dim_h,lons_lat,vg)
      call top2bot_array(lon_dim_a,lons_lat,dtdf)
      call top2bot_array(lon_dim_a,lons_lat,dtdl)

      call top2bot_ntrac(lon_dim_a,lons_lat,rgs)
!
      if(gg_tracers .and. kdt == 1)then          ! do only once at start
        do k=1,levs
          do j=1,lons_lat
            rgt_h(j,k,1) = max(0.0, rgs(j,k,1))
          enddo
        enddo
        if (ntrac > 1) then
          do n=2,ntrac
            do k=1,levs
              do j=1,lons_lat
                rgt_h(j,k,n) = rgs(j,k,n)
              enddo
            enddo
          enddo
        endif
      endif                                      ! tracers at kdt.eq.1

      if (.not. gg_tracers) then
        do n=1,ntrac
          do k=1,levs
            do j=1,lons_lat
              rgt_h(j,k,n) = rgs(j,k,n)
            enddo
          enddo
        enddo
      endif

      sinra  = sqrt(cons1-cons1/rcl)
      coriol = cons2*omega*sinra
      if (lat > latg2) coriol = -coriol
!     coriol=0.

      rcoslat = sqrt(rcl)      ! =1/cos(lat)
      coslat  = 1.0 / rcoslat  ! cos(lat)
      clog2   = log(cons2)     ! constant
      delta   = cvap/cp  ! check these cpv cpd (at const p for vapor and dry
      delta1  = delta-cons1
      rk1     = rk + 1.e0
      rkr     = 1.0/rk


!     ifirst=0

      do j=1,lons_lat
        expq(j)=exp(qg(j))
      enddo
      do k=1,levp1
        do j=1,lons_lat
          pk5(j,k) = ak5(k) + bk5(k)*expq(j)
        enddo
      enddo
      do k=1,levs
        do j=1,lons_lat
          dpk(j,k)   = pk5(j,k+1) - pk5(j,k)
          rdel(j,k)  = cons1/dpk(j,k)
          rdel2(j,k) = 0.5*rdel(j,k)*gz_grid(j)
        enddo
      enddo
      k = 1
      do j=1,lons_lat
        alfa(j,k)=clog2
        rlnp(j,1)=cons0
      enddo
      do  k=2,levs
        do j=1,lons_lat
          rlnp(j,k) = log( pk5(j,k+1)/pk5(j,k) )
          alfa(j,k) = cons1-( pk5(j,k)/dpk(j,k) )*rlnp(j,k)
        enddo
      enddo
      spdmax = 0.

      do  k=1,levs
        do j=1,lons_lat
          ug(j,k)   = ug(j,k)*rcoslat ! ug at input is scaled U=u*cos(theta)
          vg(j,k)   = vg(j,k)*rcoslat ! vg at input is scaled V=v*cos(theta)

          dtdl(j,k) = dtdl(j,k)*coslat !dtdl in=(dt/dlamda  )/cos**2(theta)
          dtdf(j,k) = dtdf(j,k)*coslat !dtdf in=(dt/d(theta))/cos**2(theta)

          speed2    = ug(j,k)*ug(j,k)+vg(j,k)*vg(j,k)

          if (speed2  > spdmax(k))  spdmax(k) = speed2 
        enddo
      enddo

      do j=1,lons_lat
        dqdlam(j) = dqdlam(j)*rcoslat                 
        dqdphi(j) = dqdphi(j)*rcoslat                 
      enddo
         
      do k=1,levs
        do j=1,lons_lat
          cg(j,k) = ug(j,k)*dqdlam(j) + vg(j,k)*dqdphi(j)
        enddo
      enddo

      k = 1
      do j=1,lons_lat
        db(j,k) = dg(j,k)*dpk(j,k)
        cb(j,k) = cg(j,k)*dbk(k)
      enddo

      do k=1,levm1
        do j=1,lons_lat
          db(j,k+1) = db(j,k) + dg(j,k+1)*dpk(j,k+1)
          cb(j,k+1) = cb(j,k) + cg(j,k+1)*dbk(k+1)
        enddo
      enddo

!sela dqdt=partial(ln(p_s)/partial(t)
      do j=1,lons_lat
        dqdt(j)      = -db(j,levs)/expq(j)-cb(j,levs) !old eulerian dqdt
        dot(j,    1) = cons0
        dot(j,levp1) = cons0
      enddo
!sela dot=( eta_dot*d(p)/d(eta) )(k+1/2) as in Eulerian 
      do k=1,levs-1
        do j=1,lons_lat
          dot(j,k+1) = -expq(j)*(bk5(k+1)*dqdt(j)+cb(j,k)) -db(j,k)
        enddo
      enddo
!sela ww=eta_dot
      do k=2,levs
        do j=1,lons_lat
          zz(j,k) = dot(j,k)
     &     *( etaint(k+1)-etaint(k-1) )/( pk5(j,k+1)-pk5(j,k-1) )
        enddo
      enddo
      do j=1,lons_lat
        zz(j,    1) = 0.0
        zz(j,levp1) = 0.0
      enddo

!     do j=1,lons_lat
!$$$       wfilt(j,    1)=0.0
!$$$       wfilt(j,levp1)=0.0
!     enddo
!$$$!    filter ertical velocities
!$$$!     k=2
!$$$        do i=1,lons_lat
!$$$         wfilt(i,2)=-zz(i,3)*wgtw(1)
!$$$     .              -zz(i,2)*wgtw(2)
!$$$! zz(i,1)=0         +zz(i,1)*wgtw(3)
!$$$     .              +zz(i,2)*wgtw(4)
!$$$     .              +zz(i,3)*wgtw(5)
!$$$     .              +zz(i,4)*wgtw(6)
!$$$     .              +zz(i,5)*wgtw(7)
!$$$        enddo
!$$$!     k=3
!$$$        do i=1,lons_lat
!$$$         wfilt(i,3)=-zz(i,2)*wgtw(1)
!$$$! zz(i,1)=0         -zz(i,1)*wgtw(2)
!$$$     .              +zz(i,2)*wgtw(3)
!$$$     .              +zz(i,3)*wgtw(4)
!$$$     .              +zz(i,4)*wgtw(5)
!$$$     .              +zz(i,5)*wgtw(6)
!$$$     .              +zz(i,6)*wgtw(7)
!$$$        enddo
!$$$
!$$$       do k=4,levs-2
!$$$        do i=1,lons_lat
!$$$         wfilt(i,k)=zz(i,k-3)*wgtw(1)
!$$$     .             +zz(i,k-2)*wgtw(2)
!$$$     .             +zz(i,k-1)*wgtw(3)
!$$$     .             +zz(i,k  )*wgtw(4)
!$$$     .             +zz(i,k+1)*wgtw(5)
!$$$     .             +zz(i,k+2)*wgtw(6)
!$$$     .             +zz(i,k+3)*wgtw(7)
!$$$        enddo
!$$$       enddo
!$$$!     k=levs-1
!$$$        do i=1,lons_lat
!$$$         wfilt(i,levs-1)= zz(i,levs-4)*wgtw(1)
!$$$     .                   +zz(i,levs-3)*wgtw(2)
!$$$     .                   +zz(i,levs-2)*wgtw(3)
!$$$     .                   +zz(i,levs-1)*wgtw(4)
!$$$     .                   +zz(i,levs  )*wgtw(5)
!$$$! zz(levs+1)=0.          +zz(i,levs+1)*wgtw(6)
!$$$     .                   -zz(i,levs  )*wgtw(7)
!$$$        enddo
!$$$!     k=levs   
!$$$        do i=1,lons_lat
!$$$         wfilt(i,levs  )= zz(i,levs-3)*wgtw(1)
!$$$     .                   +zz(i,levs-2)*wgtw(2)
!$$$     .                   +zz(i,levs-1)*wgtw(3)
!$$$     .                   +zz(i,levs  )*wgtw(4)
!$$$! zz(levs+1)=0           +zz(i,levs+1)*wgtw(5)
!$$$     .                   -zz(i,levs  )*wgtw(6)
!$$$     .                   -zz(i,levs-1)*wgtw(7)
!$$$        enddo
!sela do k=1,levs
!sela    do j=1,lons_lat
!sela       ww(j,k)=0.5*(wfilt(j,k)+wfilt(j,k+1) )
!sela    enddo
!$$$! fin   filter ertical velocities

      tref   = ref_temp
      rdtref = rd*tref
!
      do k=1,levs
        do j=1,lons_lat
          ww(j,k) = 0.5*(zz(j,k)+zz(j,k+1) )

          tb(j,k) = tb_k(k)*gz_grid(j)                          
!
          sadx(j,k) = ug(j,k)*grad_gzlam(j) + vg(j,k)*grad_gzphi(j)     
          advz(j,k) = sadx(j,k)*r_rt  
        enddo
      enddo
!
!sela check sign of rdel2
!
      do j=1,lons_lat
        sadk(j,1) = rdel2(j,1) * dot(j,2) * (tb_k(2)-tb_k(1))
      enddo
      k = levs
      do j=1,lons_lat
        sadk(j,k) = rdel2(j,k) * dot(j,k) * (tb_k(k)-tb_k(k-1))
      enddo

      do k=2,levm1
        do j=1,lons_lat
          sadk(j,k) = rdel2(j,k) * (dot(j,k+1) * (tb_k(k+1)-tb_k(k))
     &                           +  dot(j,k)   * (tb_k(k)-tb_k(k-1)))
        enddo
      enddo
      do k=1,levs
        do j=1,lons_lat
          sadt(j,k) = tb_k(k)*sadx(j,k) + sadk(j,k)
        enddo
      enddo
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
!!!for testing input fields=output or any locally computed fields=output
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      test_mode = .false.
      if (test_mode) then !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
        print*,' enter test mode in gfidi_tracers with lin temp.'
        do k=1,levs
          do j=1,lons_lat
           ud3_h1(j,k)=ug(j,k)*rcl !watch scaling for test
           vd3_h1(j,k)=vg(j,k)*rcl !watch scaling for test
            dpdt_d3_h1(j,k) = qg(j)
            tg(j,k)      =  sadt(j,k) !for testing with ncar graphics thru td3_h
          td3_h1(j,k)=tg(j,k)
          enddo
        enddo
        do j=1,lons_lat
          dpdt_a(j) = qg(j)
        enddo
        return
      endif
!      if (test_mode) return !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      batah2 = batah*0.5
      do k=1,levs
        sv(k) = batah2*sv_ecm(k)
      enddo

      z_n_sum = 0.
      do k=1,levs
        do j=1,lons_lat
          z_n(j,k)   = dbk(k)*( dqdt(j) + cg(j,k) + advz(j,k) )
          z_n_sum(j) = z_n_sum(j) + z_n(j,k)
        enddo
      enddo
      do j=1,lons_lat
        sd_n(j) = sv(1)*dg(j,1)
      enddo
      do k=2,levs
        do j=1,lons_lat
          sd_n(j) = sd_n(j) + sv(k)*dg(j,k)
        enddo
      enddo
!sk splitting of dpdt_d3 into dpdt_d3_h1 and dpdt_d3_h2
!sk to facilitate lin_xy = t or cont_eq_opt1 = t
      if (go_forward) then
        do j=1,lons_lat
          dpdt_a(j) = deltim*( 0.5*z_n_sum(j)+sd_n(j) ) ! at arrival points
          sd_m(j)   = sd_n(j)
        enddo
        if (cont_eq_opt1) then
         do k=1,levs
           do j=1,lons_lat
             dpdt_d3_h1(j,k) = deltim*( 0.5*z_n(j,k) )
             dpdt_d3_h2(j,k) = ( qg(j) + gz_grid(j)*r_rt )*dbk(k)
             z_m(j,k)     = z_n(j,k)
           enddo
         enddo
        elseif (lin_xy) then
         do k=1,levs
           do j=1,lons_lat
              dpdt_d3_h1(j,k) = qg(j)*dbk(k)
              dpdt_d3_h2(j,k) = deltim*( 0.5*z_n(j,k) )
     &                       + ( qg(j) + gz_grid(j)*r_rt )*dbk(k)
     &                       - dpdt_d3_h1(j,k)
              z_m(j,k)     = z_n(j,k)
           enddo
         enddo
        else  !default
        do k=1,levs
          do j=1,lons_lat
            dpdt_d3_h1(j,k) = deltim*( 0.5*z_n(j,k) )
     &                   + ( qg(j) + gz_grid(j)*r_rt )*dbk(k)
            dpdt_d3_h2(j,k) = 0.0
            z_m(j,k)     = z_n(j,k)
          enddo
        enddo
        endif

      else  ! go_forward

        if (cont_eq_opt1) then
         do k=1,levs
           do j=1,lons_lat
             dpdt_d3_h1(j,k) = deltim*( z_n(j,k)-0.5*z_m(j,k) )
             dpdt_d3_h2(j,k) = ( qg(j) + gz_grid(j)*r_rt )*dbk(k)
     &                       + dbk(k)*deltim*( sd_n(j) - sd_m(j) )
             z_m(j,k)     = z_n(j,k)
           enddo
         enddo
        elseif (lin_xy) then
         do k=1,levs
           do j=1,lons_lat
             dpdt_d3_h1(j,k) = qg(j)*dbk(k)
             dpdt_d3_h2(j,k) = deltim*( z_n(j,k)-0.5*z_m(j,k) )
     &                    + ( qg(j) + gz_grid(j)*r_rt )*dbk(k)
     &                    + dbk(k)*deltim*( sd_n(j) - sd_m(j) )
     &                    - dpdt_d3_h1(j,k)
             z_m(j,k)     = z_n(j,k)
           enddo
         enddo
        else !default
        do k=1,levs
          do j=1,lons_lat
            dpdt_d3_h1(j,k) = deltim*( z_n(j,k)-0.5*z_m(j,k) )
     &                   + ( qg(j) + gz_grid(j)*r_rt )*dbk(k)
     &                   + dbk(k)*deltim*( sd_n(j) - sd_m(j) )
            dpdt_d3_h2(j,k) = 0.0
            z_m(j,k)     = z_n(j,k)
          enddo
        enddo
        endif
        do j=1,lons_lat
          dpdt_a(j) = deltim*( 0.5*z_n_sum(j)+sd_n(j) ) ! at arrival points
          sd_m(j)   = sd_n(j)
        enddo
      endif !go_forward
!
      do j=1,lons_lat
        cofb(j,1)    = -rdel(j,1) * (alfa(j,1)*dbk(1))
        px2(j,levs)  = cons0
        px3u(j,levs) = cons0
        px3v(j,levs) = cons0
      enddo
!
      do k=2,levs
        kk  = levp1 - k + 1
        kk1 = kk + 1
        do j=1,lons_lat
          cofb(j,k)    = -rdel(j,k)*(bk5(k)*rlnp(j,k)+alfa(j,k)*dbk(k))

          px2(j,kk-1)  = px2(j,kk)
     &       - rd*(bk5(kk1)/pk5(j,kk1)-bk5(kk)/pk5(j,kk)) * tg(j,kk)
!
          px3u(j,kk-1) = px3u(j,kk) - rd*rlnp(j,kk)*dtdl(j,kk)
          px3v(j,kk-1) = px3v(j,kk) - rd*rlnp(j,kk)*dtdf(j,kk)
        enddo
      enddo
      do k=1,levs
        do j=1,lons_lat
          u_n(j,k) =  vg(j,k)*coriol ! ug,vg are not scaled in this version
          v_n(j,k) = -ug(j,k)*coriol ! ug,vg are not scaled in this version
!sela     u_n(j,k)=0.
!sela     v_n(j,k)=0.
!
          uprs(j,k) = cofb(j,k)*rd*tg(j,k)*expq(j)*dqdlam(j)
          vprs(j,k) = cofb(j,k)*rd*tg(j,k)*expq(j)*dqdphi(j)
!
          cofa(j,k) = -rdel(j,k)*(
     &                 bk5(k+1)*pk5(j,k)/pk5(j,k+1) - bk5(k)
     &              +  rlnp(j,k) * (bk5(k)-pk5(j,k)*dbk(k)*rdel(j,k)) )
!
          px1u(j,k) = -grad_gzlam(j)  ! grad_gz is expected to be true grad
          px1v(j,k) = -grad_gzphi(j)  ! grad_gz is expected to be true grad
!
          px2u(j,k) = px2(j,k)*expq(j)*dqdlam(j)
          px2v(j,k) = px2(j,k)*expq(j)*dqdphi(j)
        enddo
      enddo

!     do k=1,levs
!     do j=1,lons_lat
!      px3u(j,k)=px3u(j,k)/rcl
!      px3v(j,k)=px3v(j,k)/rcl
!     enddo
!     enddo

      do k=1,levs
        do j=1,lons_lat
          px4u(j,k) = -rd*alfa(j,k)*dtdl(j,k) ! in scaled version /rcl
          px4v(j,k) = -rd*alfa(j,k)*dtdf(j,k) ! in scaled version /rcl
!
          px5u(j,k) = -cofa(j,k)*rd*tg(j,k)*expq(j)*dqdlam(j)
          px5v(j,k) = -cofa(j,k)*rd*tg(j,k)*expq(j)*dqdphi(j)
!
          uphi(j,k) = px1u(j,k)+px2u(j,k)+px3u(j,k)+px4u(j,k)+px5u(j,k)
          vphi(j,k) = px1v(j,k)+px2v(j,k)+px3v(j,k)+px4v(j,k)+px5v(j,k)
!
          u_n(j,k) = u_n(j,k) + uphi(j,k) + uprs(j,k)
          v_n(j,k) = v_n(j,k) + vphi(j,k) + vprs(j,k)
!
          worka(j,k) = rk*tg(j,k) * (cons1+fv*rgt_h(j,k,1))
     &               / (cons1+delta1*rgt_h(j,k,1)) * rdel(j,k)
        enddo
      enddo

      do j=1,lons_lat
        workb(j,1)=
     &  alfa(j,1)*( dg(j,1)*dpk(j,1)+expq(j)*cb(j,1)*dbk(1) )
!
        workc(j,1) = expq(j)*cg(j,1)*dbk(1)
      enddo

      do k=2,levs
        do j=1,lons_lat
          workb(j,k) = rlnp(j,k) * (db(j,k-1)+expq(j)*cb(j,k-1))
     &    + alfa(j,k) * (dg(j,k)*dpk(j,k)+expq(j)*cg(j,k)*dbk(k))
!
          workc(j,k) = expq(j) * cg(j,k)
     &               * ( dbk(k) + ck(k)*rlnp(j,k)*rdel(j,k) )
        enddo
      enddo

      do k=1,levs
        do j=1,lons_lat
          t_n(j,k) = worka(j,k) * ( -workb(j,k) + workc(j,k))
        enddo
      enddo

!$$$  do k=1,levs
!$$$  do j=1,lons_lat
!$$$   u_l(j,k)= 0.
!$$$   v_l(j,k)= 0.
!$$$   t_l(j,k)= 0.
!$$$  enddo
!$$$  enddo

!$$$  do k_row=1,levs
!$$$   do k_col=1,levs
!$$$    do j=1,lons_lat
!$$$     u_l(j,k_row)=u_l(j,k_row)+y_ecm(k_row,k_col)*dtdl(j,k_col)
!$$$     v_l(j,k_row)=v_l(j,k_row)+y_ecm(k_row,k_col)*dtdf(j,k_col)
!$$$     t_l(j,k_row)=t_l(j,k_row)+t_ecm(k_row,k_col) * dg(j,k_col)
!$$$    enddo
!$$$   enddo
!$$$  enddo

      if ( kind_evod .eq. 8 ) then !------------------------------------
         call dgemm ('n', 't',
     &               lons_lat, levs, levs, cons1,
     &               dtdl(1,1), lon_dim_a,
     &               y_ecm(1,1), levs, cons0,
     &               u_l(1,1), lon_dim_h)
         call dgemm ('n', 't',
     &               lons_lat, levs, levs, cons1,
     &               dtdf(1,1), lon_dim_a,
     &               y_ecm(1,1), levs, cons0,
     &               v_l(1,1), lon_dim_h)
         call dgemm ('n', 't',
     &               lons_lat, levs, levs, cons1,
     &               dg(1,1), lon_dim_a,
     &               t_ecm(1,1), levs, cons0,
     &               t_l(1,1), lon_dim_h)
      else !------------------------------------------------------------
         call sgemm ('n', 't',
     &               lons_lat, levs, levs, cons1,
     &               dtdl(1,1), lon_dim_a,
     &               y_ecm(1,1), levs, cons0,
     &               u_l(1,1), lon_dim_h)
         call sgemm ('n', 't',
     &               lons_lat, levs, levs, cons1,
     &               dtdf(1,1), lon_dim_a,
     &               y_ecm(1,1), levs, cons0,
     &               v_l(1,1), lon_dim_h)
         call sgemm ('n', 't',
     &               lons_lat, levs, levs, cons1,
     &               dg(1,1), lon_dim_a,
     &               t_ecm(1,1), levs, cons0,
     &               t_l(1,1), lon_dim_h)
      endif !-----------------------------------------------------------

      do k=1,levs
        do j=1,lons_lat
          u_l(j,k) = rdtref*dqdlam(j) + u_l(j,k) ! in scaled version /rcl
          v_l(j,k) = rdtref*dqdphi(j) + v_l(j,k) ! in scaled version /rcl
        enddo
      enddo

      if (go_forward) then
        do k=1,levs
          do j=1,lons_lat
            unla(j,k) = u_n(j,k) + batah*u_l(j,k)
            vnla(j,k) = v_n(j,k) + batah*v_l(j,k)
            tnla(j,k) = t_n(j,k) + batah*t_l(j,k) - sadt(j,k)
!sela   ud3_h(j,k)=ug(j,k)+omega_v+
        ud3_h1(j,k) = ug(j,k)
        ud3_h2(j,k) = deltim*(0.5*unla(j,k)-batah2*u_l(j,k))       

        vd3_h1(j,k) = vg(j,k)
        vd3_h2(j,k) = deltim*(0.5*vnla(j,k)-batah2*v_l(j,k))

        td3_h1(j,k) = tg(j,k)
        td3_h2(j,k) = deltim*(0.5*tnla(j,k)-batah2*t_l(j,k)) -tb(j,k)

            unlm(j,k)  = unla(j,k)
            vnlm(j,k)  = vnla(j,k)
            tnlm(j,k)  = tnla(j,k)

            ug_m(j,k)  = ug(j,k) !for midpoint value in dep. point calc
            vg_m(j,k)  = vg(j,k) !for midpoint value in dep. point calc
            ww_m(j,k)  = ww(j,k)
          enddo
        enddo

      else  ! go_forward

        do k=1,levs
          do j=1,lons_lat
            unla(j,k) = u_n(j,k) + batah*u_l(j,k)
            vnla(j,k) = v_n(j,k) + batah*v_l(j,k)
            tnla(j,k) = t_n(j,k) + batah*t_l(j,k)-sadt(j,k)
!sela   ud3_h(j,k)=ug(j,k)+omega_v+
        ud3_h1(j,k) = ug(j,k)
        ud3_h2(j,k) = deltim*(unla(j,k)-0.5*unlm(j,k)-batah2*u_l(j,k))

        vd3_h1(j,k) = vg(j,k)
        vd3_h2(j,k) = deltim*(vnla(j,k)-0.5*vnlm(j,k)-batah2*v_l(j,k))

        td3_h1(j,k) = tg(j,k) 
        td3_h2(j,k) = deltim*
     &               (tnla(j,k)-0.5*tnlm(j,k)-batah2*t_l(j,k)) - tb(j,k)

            unlm(j,k)  = unla(j,k)
            vnlm(j,k)  = vnla(j,k)
            tnlm(j,k)  = tnla(j,k)
          enddo
        enddo

      endif ! go_forward


!my if any hermite interpolation then linxyz and cubwgtlin should not be used in slgscan 
!my separation of linear and nonlinear terms is not needed
!my also if lin_xyz is false do not separate
      if (herm_x .or. herm_y .or. herm_z .or. .not. lin_xyz) then
         do k=1,levs
           do j=1,lons_lat
                ud3_h1(j,k) = ud3_h1(j,k) + ud3_h2(j,k)
                vd3_h1(j,k) = vd3_h1(j,k) + vd3_h2(j,k)
                td3_h1(j,k) = td3_h1(j,k) + td3_h2(j,k)
             dpdt_d3_h1(j,k) = dpdt_d3_h1(j,k) + dpdt_d3_h2(j,k)
           enddo
         enddo
      endif
      if (.not.lin_xy .and. .not.cont_eq_opt1) then
         do k=1,levs
           do j=1,lons_lat
             dpdt_d3_h1(j,k) = dpdt_d3_h1(j,k) + dpdt_d3_h2(j,k)
           enddo
         enddo
      endif

!my12272012
      do k=1,levs
      do j=1,lons_lat
        ugsave    = ug(j,k)
        ug(j,k)   = 1.5*ug(j,k) - 0.5*ug_m(j,k)
        ug_m(j,k) = ugsave
        vgsave    = vg(j,k)
        vg(j,k)   = 1.5*vg(j,k) - 0.5*vg_m(j,k)
        vg_m(j,k) = vgsave
      enddo
      enddo
      if (time_extrap_etadot) then
        do k=1,levs
        do j=1,lons_lat
          wwsave    = ww(j,k)
          ww(j,k)   = 1.5*ww(j,k) - 0.5*ww_m(j,k)
          ww_m(j,k) = wwsave
        enddo
        enddo
      endif
!sk10052012
      if (settls_dep3ds.or.settls_dep3dg) then
        do k=1,levs
        do j=1,lons_lat
          ugsave    = ug(j,k)
          ug2(j,k)   = 2.0*ug(j,k) - ug_m(j,k)
          ug_m(j,k) = ugsave
          vgsave    = vg(j,k)
          vg2(j,k)   = 2.0*vg(j,k) - vg_m(j,k)
          vg_m(j,k) = vgsave
          wwsave    = ww(j,k)
          ww2(j,k)   = 2.0*ww(j,k) - ww_m(j,k)
          ww_m(j,k) = wwsave
        enddo
        enddo
      endif

      call top2bot_array(1,1,spdmax)


!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
!!!for testing input fields=output or any locally computed fields=output
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      test_mode =.false.
      if (test_mode) then !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
       print*,' enter test mode in gfidi_tracers with lin temp.'

!$$$         do j=1,lons_lat
!$$$           dpdt_a(j)=qg(j)
!$$$         enddo
        do k=1,levs
          do j=1,lons_lat

            unlm(j,k) = unla(j,k) * coslat 
            vnlm(j,k) = vnla(j,k) * coslat
            tnlm(j,k) = tnla(j,k) ! for testing with ncar graphics thru td3_h
           enddo
        enddo
      endif
       if(test_mode) return !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

      return
      end

      subroutine top2bot_array(lon_dim,lons_lat,array)
!mjr  use gfs_dyn_MACHINE, only : kind_phys
      use gfs_dyn_resol_def, only : levp1,levs
!
      implicit none
      integer lon_dim,lons_lat,i,k,kk
      real array(lon_dim,levs),tempo(lons_lat)
!
      do k=1,levs/2
        kk = levp1-k
        do i=1,lons_lat
          tempo(i)    = array(i,kk)
          array(i,kk) = array(i,k)
          array(i,k)  = tempo(i)
        enddo
      enddo
      return
      end
      subroutine top2bot_ntrac(lon_dim,lons_lat,array)
!mjr  use gfs_dyn_MACHINE, only : kind_phys
      use gfs_dyn_resol_def, only : levp1,levs,ntrac
!
      implicit none
      integer lon_dim,lons_lat,i,k,n,kk
      real array(lon_dim,levs,ntrac),tempo(lons_lat)
      do n=1,ntrac
        do k=1,levs/2
          kk = levp1-k
          do i=1,lons_lat
            tempo(i)      = array(i,kk,n)
            array(i,kk,n) = array(i,k,n)
            array(i,k,n)  = tempo(i)
          enddo
        enddo
      enddo
      return
      end