 
      module dt_module
      use addfld_module,only: addfld
      implicit none
      contains
!-----------------------------------------------------------------------
      subroutine dt(tn,tn_nm,un,vn,o2,o1,ox,n2,barm,cp,kt,km,hdt,qji_tn,
     |  gwt,difkt,cool_imp,cool_exp,w_upd,tn_upd,tn_nm_upd,
     |  lev0,lev1,lon0,lon1,lat0,lat1)
!
! Advance neutral temperature at current latitude: 
! updates tn(itc) and tn_nm(itc)
! updates tlbc and tlbc_nm, corresponding to tn(itc) and tn(itp)= tn_nm(itc)
! LB is locally calculated and stored in tnlbc
!
      use params_module,only: nlonp4,dz,nlat,spval
      use input_module,only: step,ncep_ncfile,ncep_reanalysis,
     |  ecmwf_ncfile
      use init_module,only: iter,igetgswm,istep
      use cons_module,only: freq_semidi,tbound,shapiro,dtx2inv,expz,
     |  rmassinv_o2,rmassinv_o1,rmassinv_n2,tsurplus,p0,boltz,avo,grav,
     |  gask,expzmid,expzmid_inv,dtsmooth,dtsmooth_div2,freq_3m3,kut
      use qrj_module,only: qtotal ! qtotal(nlevp1,lon0:lon1,lat0:lat1)
      use chemrates_module,only: rkm12
      use zatmos_module,only: zatmos_tn ! zatmos temperature (nlat)
      use ncep_module,only: tncep       ! ncep temperature (nlonp4,nlat)
      use ncep_rean_module,only: t_ncep
      use ecmwf_module,only: t_ecmwf
      use fields_module,only: tlbc,ulbc,vlbc,tlbc_nm
      use lbc,only: t_lbc
      use diags_module,only: 
     |  mkdiag_DEN,           ! total density
     |  mkdiag_QJOULE,        ! joule heating
     |  mkdiag_QHORDIF,       ! heating from horizontal diffusion
     |  mkdiag_QGRAVW,        ! gravity wave heating
     |  mkdiag_QADIAB,        ! adiabatic heating
     |  mkdiag_QHORADV,       ! horizontal advection heating
     |  mkdiag_QVERADV,       ! vertical advection heating
     |  mkdiag_QDISSIP,       ! dissipative heating
     |  mkdiag_QTOTAL         ! total heating
#ifdef MPI
      use mpi_module,only: mp_bndlons_f3d, mp_periodic_f3d
#endif
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
!
! Full subdomains:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in)::
     |  tn,    ! neutral temperature (deg K)
     |  tn_nm, ! neutral temperature, time n-1 
     |  un,    ! neutral zonal velocity (cm/sec)
     |  vn,    ! neutral zonal velocity (cm/sec)
     |  o2,    ! molecular oxygen (mmr)
     |  o1,    ! atomic oxygen (mmr)
     |  ox,    ! OX family (mmr)
     |  n2,    ! molecular nitrogen (mmr)
     |  barm,  ! mean molecular weight
     |  cp,    ! specific heat (ergs/deg/gm)           (sub cpktkm)
     |  kt,    ! molecular diffusion (ergs/cm/deg/sec) (sub cpktkm)
     |  km,    ! molecular viscosity (gm/cm/sec)       (sub cpktkm)
     |  hdt,   ! horizontal diffusion of tn (from sub hdif3, hdif.F)
     |  qji_tn,! joule heating for tn (from sub qjoule_tn, qjoule.F)
     |  gwt,   ! gravity wave heating (mgw.F)
     |  difkt, ! eddy thermal diffusion
     |  cool_imp, ! implicit cooling (newton.F)
     |  cool_exp, ! explicit cooling (newton.F)
     |  w_upd  ! updated vertical velocity (swdot.F)
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  tn_upd,   ! updated tn (output)
     |  tn_nm_upd ! updated tn at time n-1 (output)
!
! Local:
      integer :: k,i,j,lonbeg,lonend,lat,kutt(nlat)
      integer :: i0,i1,nk,nkm1,nlevs
!     complex :: expt,expt2,expt3m3
!     real :: rstep
      real :: tnlbc(lon0:lon1,lat0:lat1) ! lower boundary condition
      real,parameter :: joulefac = 1.0   ! Joule heating multiplication factor
                                         ! (was 2.0 in timegcm1-2dev7, 1.5 in tiegcm1.9)
!
! Local at 2d:
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  cptn,       ! cp*exp(-s)*(V.del(T(n))-1/(2*DT))*T(n-1) (k+1/2)
     |  mbar,       ! mean molecular weight
     |  qm,         ! heating due to molecular diffusion
     |  total_heat, ! total heating
     |  dudz,       ! du/dz(k)
     |  dvdz,       ! du/dz(k)
     |  g,          ! g*KT/(p0*H*Ds**2)
     |  f,          ! g*eps/(p0*2*H*Ds)
     |  h,          ! scale height R*T/(M*g) (cm)
     |  rho,        ! density
     |  tni,        ! tn at interfaces
     |  p, q, r,    ! coefficients for tridiagonal solver
     |  rhs,        ! right hand side of trsolv
     |  qpart,      ! part of q coeff
     |  tnlbc_diag, ! tnlbc redundant in vertical
     |  radcool     ! radiative cooling diagnostic
!
! Local at 3d (tnsmooth needs lat dimension only for sub smooth):
      real,dimension(lev0:lev1,lon0:lon1,lat0:lat1) ::
     |  tnsmooth,   ! zonal and meridional smoothing of tn_nm
     |  advec_tn    ! horizontal advection (output of sub advec)
!
! 4/17/06 btf:
! Local diagnostics for thermodynamic balance (turn on/off with thermo_diag)
!
      logical,parameter :: thermo_diag=.true.
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  verdif,     ! vertical temperature diffusion
     |  grwvht,     ! gravity wave vertical heat transfer
     |  dtds  ,     ! vertical temperature gradient
     |  sht   ,     ! scale height
     |  hflxm ,     ! heat flux from molecular diffusion
     |  coef  ,     ! coefficient
     |  dhflxm,     ! gradient of horizontal molecular heat flux
     |  radnet,     ! net heating and cooling
     |  xmbr  ,     ! mean mass
     |  xkmolt,     !
     |  totfab,     !
     |  hflxe ,     !
     |  dhflxe      !
!
! These diags need lat dimension because they are calculated in the 2nd
! latitude scan, but are needed for diags totfab and radnet at the end of
! the 3rd latitude scan.
!
      real,dimension(lev0:lev1,lon0:lon1,lat0:lat1) ::
     |  qsolar,     ! total solar heating diagnostic
     |  qjoule,     ! joule heating diagnostic
     |  hordif,     ! horizontal temperature diffusion
     |  horadv,     ! horizontal temperature advection
     |  veradv,     ! vertical temperature advection
     |  htflxm,     !
     |  adiaht,     ! adiabatic heating
     |  htflxe      !
!
      nk = lev1-lev0+1
      nkm1 = nk-1
      nlevs = nk
      i0 = lon0 ; i1 = lon1
      lonbeg = lon0-2
      if (lon0==1) lonbeg = 1
      lonend = lon1+2
      if (lon1==nlonp4) lonend = nlonp4
!
!     rstep = float(step)
!     expt = cexp(ci*freq_semidi*rstep*iter)
!     expt2 = cexp(ci*.5*freq_semidi*rstep*iter)
!     expt3m3 = cexp(ci*freq_3m3*rstep*iter)
!
! First latitude scan for dt:
      do lat=lat0,lat1
!
! TN lower boundary is in t_lbc (lbc.F):
        do i=lon0,lon1
          tnlbc(i,lat) = t_lbc(i,lat)
        enddo
!
! Horizontal advection (pass k vs i slices at full task subdomain 
! longitudes, and the 5 latitudes centered over the current latitude).
!
      call advec(tn(:,:,lat-2:lat+2),advec_tn(:,:,lat),
     |  lev0,lev1,lon0,lon1,lat)

!     call addfld('HADVECTN',' ',' ',advec_tn(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)

        if (thermo_diag) horadv(:,:,lat) = advec_tn(:,:,lat)
!
! Vertical advection. Sub advecv adds vertical advection to advec_tn.
      call advecv(tn(:,:,lat),tnlbc(:,lat),advec_tn(:,:,lat),
     |  lev0,lev1,lon0,lon1,lat)

      if (thermo_diag) 
     |  veradv(:,:,lat) = advec_tn(:,:,lat) - horadv(:,:,lat)

!     call addfld('ADVEC_TN',' ',' ',advec_tn(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)

        if (thermo_diag) then
          horadv(:,:,lat) = horadv(:,:,lat)*86400.
          veradv(:,:,lat) = veradv(:,:,lat)*86400.

          call mkdiag_QHORADV('QHORADV',horadv(:,:,lat),
     |      lev0,lev1,lon0,lon1,lat)
          call mkdiag_QVERADV('QVERADV',veradv(:,:,lat),
     |      lev0,lev1,lon0,lon1,lat)

        endif
!
! End first latitude scan:
      enddo ! lat=lat0,lat1
!
! Shapiro smoother for tn at time n-1:
      call smooth(tn_nm,tnsmooth,lev0,lev1,lon0,lon1,lat0,lat1,0)
!
! Begin second latitude scan:
      do lat=lat0,lat1
!
! Set cptn and mbar (k+1/2):
! (Earlier versions apparently assumed zero periodic points for
!  tnsmooth, since they were not set in smoothing. See sub smooth,
!  where the periodic points are set to zero to avoid NaNS fpe
!  in the following loop)
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          cptn(k,i) = .5*(cp(k,i,lat)+cp(k+1,i,lat))*expz(k)*
     |      (advec_tn(k,i,lat)-dtx2inv*tnsmooth(k,i,lat))
          mbar(k,i) = 1./(o2(k,i,lat)*rmassinv_o2 + 
     |      ox(k,i,lat)*rmassinv_o1 + n2(k,i,lat)*rmassinv_n2)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('CP'   ,' ',' ',cp(:,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('TNSMOOTH'   ,' ',' ',tnsmooth(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('CPTN0',' ',' ',cptn,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('MBAR' ,' ',' ',mbar,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Total heat sources are in total_heat (s5).
      do i=lon0,lon1
        do k=lev0,lev1-1
!
! Solar heating from qrj:
          total_heat(k,i) = .5*(qtotal(k,i,lat)+qtotal(k+1,i,lat))
          qsolar(k,i,lat) = total_heat(k,i)*86400./cp(k,i,lat)
!
! Add heating from 4th order horizontal diffusion (hdt from sub hdif3):
          total_heat(k,i) = total_heat(k,i)+hdt(k,i,lat)
!
! Add gravity wave heating (gw may be turned off, check calc_gw):
          total_heat(k,i) = total_heat(k,i)+gwt(k,i,lat)
        enddo ! k=lev0,lev1-1
        qsolar(lev1,i,lat) = spval
      enddo ! i=lon0,lon1
!
! Save heating from horizontal diffusion and gravity wave heating (deg/day):
!
      hordif(:,:,lat)=hdt(:,lon0:lon1,lat)*86400./cp(:,lon0:lon1,lat)
      grwvht(:,:) = gwt(:,lon0:lon1,lat)*86400./cp(:,lon0:lon1,lat)

      call mkdiag_QHORDIF('QHORDIF',hordif(:,:,lat),lev0,lev1,
     |  lon0,lon1,lat)
      call mkdiag_QGRAVW('QGRAVW',grwvht,lev0,lev1,lon0,lon1,lat)
!
! Add joule heating (see sub qjoule_tn, qjoule.F):
      do i=lon0,lon1
        do k=lev0,lev1-1
          total_heat(k,i) = total_heat(k,i)+qji_tn(k,i,lat)*joulefac
        enddo ! k=lev0,lev1-1
        total_heat(lev1,i) = spval
      enddo ! i=lon0,lon1
!
! Save joule heating rate in deg/day (see mkdiag_xxx calls below):
      do i=lon0,lon1
        do k=lev0,lev1-1
          qjoule(k,i,lat) = qji_tn(k,i,lat)*86400./cp(k,i,lat)
        enddo ! k=lev0,lev1-1
        qjoule(lev1,i,lat) = spval
      enddo ! i=lon0,lon1
!
      if (thermo_diag) then
!
! Init to missing value to cover undefined boundaries and near-boundary
! levels at top and bottom:
        dtds = spval
        sht = spval
        hflxm = spval
        coef = spval
        xmbr = spval
        xkmolt = spval
        adiaht(:,:,lat) = spval
        hflxe = spval
        dhflxm = spval
        htflxm(:,:,lat) = spval
        dhflxe = spval
        htflxe(:,:,lat) = spval
        do i=lon0,lon1
          do k=lev0+1,lev1-1
            dtds(k,i) = tn(k+1,i,lat)-tn(k-1,i,lat)
            sht(k,i) = gask*tn(k,i,lat)/(mbar(k,i)*grav)
            hflxm(k,i) = kt(k,i,lat)*dtds(k,i)/sht(k,i)
            coef(k,i) = grav/(p0*cp(k,i,lat)*expz(k))
            xmbr(k,i) = mbar(k,i)
            xkmolt(k,i) = km(k,i,lat)
            adiaht(k,i,lat) = -gask*tn(k,i,lat)/(cp(k,i,lat)*mbar(k,i))*
     |        w_upd(k,i,lat)*86400.
            hflxe(k,i) = difkt(k,i,lat)*gask*tn(k,i,lat)*p0*expz(k)*
     |        cp(k,i,lat)*dtds(k,i)/(sht(k,i)*mbar(k,i)*grav*grav)+
     |        (difkt(k,i,lat)*gask*gask*tn(k,i,lat)*p0*expz(k)*
     |        tn(k,i,lat))/(sht(k,i)*mbar(k,i)**2*grav**2)
          enddo ! k=lev0,lev1-1
          do k=lev0+2,lev1-2
            dhflxm(k,i) = hflxm(k+1,i)-hflxm(k-1,i)
            htflxm(k,i,lat) = coef(k,i)*dhflxm(k,i)*86400.
            dhflxe(k,i) = hflxe(k+1,i)-hflxe(k-1,i)
            htflxe(k,i,lat) = coef(k,i)*dhflxe(k,i)*86400.
          enddo
          hflxm(lev1-1,i)      = spval
          dtds(lev1-1,i)       = spval
          dhflxm(lev1-2,i)     = spval
          htflxm(lev1-2,i,lat) = spval
        enddo ! i=lon0,lon1

!       call addfld('DTDS'  ,' ',' ',dtds  ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('SHT'   ,' ',' ',sht   ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('HFLXM' ,' ',' ',hflxm ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('COEF'  ,' ',' ',coef  ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('DHFLXM',' ',' ',dhflxm,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('DHFLXE',' ',' ',dhflxe,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('HTFLXM',' ',' ',htflxm(:,:,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('HTFLXE',' ',' ',htflxe(:,:,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('XMBR'  ,' ',' ',xmbr,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('XKMOLT',' ',' ',xkmolt,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('ADIAHT',' ',' ',adiaht(:,:,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)

        call mkdiag_QADIAB('QADIAB',adiaht(:,:,lat),
     |    lev0,lev1,lon0,lon1,lat)

      endif ! thermo_diag
!
! Add heating due to molecular diffusion:
! du/dz and dv/dz (s10, s11):
      do i=lon0,lon1
        do k=lev0+1,lev1-2
          dudz(k,i) = (un(k+1,i,lat)-un(k-1,i,lat))/(2.*dz)  ! s10
          dvdz(k,i) = (vn(k+1,i,lat)-vn(k-1,i,lat))/(2.*dz)  ! s11
        enddo ! k=lev0+1,lev1-2
!
! Lower boundary is in ulbc,vlbc:
        dudz(1,i) = (un(1,i,lat)+1./3.*un(2,i,lat)-4./3.*
     |    ulbc(i,lat))/dz
        dvdz(1,i) = (vn(1,i,lat)+1./3.*vn(2,i,lat)-4./3.*
     |    vlbc(i,lat))/dz
!
! Upper boundary:
        dudz(lev1-1,i) = dudz(lev1-2,i)/3.
        dvdz(lev1-1,i) = dvdz(lev1-2,i)/3.
!
! qm = heating due to molecular diffusion:
! (km = molecular viscosity from sub cpktkm)
        do k=lev0,lev1-1
          qm(k,i) = grav**2*mbar(k,i)*.5*(km(k,i,lat)+km(k+1,i,lat))/
     |      (p0*gask*expz(k)*tn(k,i,lat))*(dudz(k,i)**2+dvdz(k,i)**2)
!
! Add qm to total heating:
          total_heat(k,i) = total_heat(k,i)+qm(k,i)
!
! Complete cptn:
! -cp*exp(-s)*(T(k,n-1)/(2*Dt) - V.del(T(k,n)) +Q/cp)
!
          cptn(k,i) = cptn(k,i)-expz(k)*total_heat(k,i) ! s1
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          total_heat(k,i) = total_heat(k,i)*86400./cp(k,i,lat)
        enddo
      enddo
!
! Save additional heating diagnostics in deg K/day.
! (note that QJOULE (erg/g/s) may have been saved from qjoule.F)
!
      call mkdiag_QJOULE('QJOULE',qjoule(:,:,lat),
     |  lev0,lev1,lon0,lon1,lat)
      call mkdiag_QDISSIP('QDISSIP',qm,lev0,lev1,lon0,lon1,lat)
      call mkdiag_QTOTAL('QTOTAL',total_heat,
     |  lev0,lev1,lon0,lon1,lat)
!
! H = R*T/(M*g)                            (s4)
! rho = p0*exp(-s)*M/(R*T)                 (s5)
! tni = T                                  (s6)
!
! Levels 2 through lev1-1:
      do i=lon0,lon1
        do k=lev0+1,lev1-1
          tni(k,i) = .5*(tn(k-1,i,lat)+tn(k,i,lat))
          h(k,i) = gask*tni(k,i)/barm(k,i,lat)
          rho(k,i) = p0*expzmid_inv*expz(k)/h(k,i)
          h(k,i) = h(k,i)/grav
        enddo ! k=lev0+1,lev1-1
!
! Boundaries:
!       tni(lev0,i) = tn(lev1,i,lat)  ! bottom boundary is in tn(lev1,..)
        tni(lev0,i) = tlbc(i,lat)     ! Lower boundary is in tlbc = tn(itp)
        tni(lev1,i) = tn(lev1-1,i,lat)
        h(lev0,i) = gask*tni(lev0,i)/barm(lev0,i,lat)
        h(lev1,i) = gask*tni(lev1,i)/barm(lev1,i,lat)
        rho(lev0,i) = p0*expzmid_inv*expz(lev0)/h(lev0,i)
        rho(lev1,i) = p0*expzmid_inv*expz(lev1-1)/h(lev1,i)
        h(lev0,i) = h(lev0,i)/grav
        h(lev1,i) = h(lev1,i)/grav
!
! G = g*(kT + H**2*rho*cp*kE)/(p0*H*Ds**2) (s2)
! F = g*(kE*H**3*rho*g/T)/(p0*2*H*Ds)      (s3)
!
        do k=lev0,lev1-1 
          g(k,i) = grav*(kt(k,i,lat)+h(k,i)**2*rho(k,i)*cp(k,i,lat)*
     |      difkt(k,i,lat))/(p0*h(k,i)*dz**2)
          f(k,i)=grav**2*difkt(k,i,lat)*h(k,i)**2*rho(k,i)/
     |      (tni(k,i)*p0*2.*dz)
        enddo ! k=lev0,lev1-1 
      enddo ! i=lon0,lon1
!
!     call addfld('DIFKT' ,' ',' ',difkt(:,i0:i1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('TNI'   ,' ',' ',tni,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('H_COEF',' ',' ',h  ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('DEN'   ,' ',' ',rho,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('G_COEF',' ',' ',g  ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('F_COEF',' ',' ',f  ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Save diagnostic density:
      call mkdiag_DEN('DEN',rho,lev0,lev1,lon0,lon1,lat)
!
! Coefficients for trsolv:
! Levels 3/2 through K-3/2
      do i=lon0,lon1
        do k=lev0,lev1-2
          p(k,i) = g(k,i)-f(k,i)
          q(k,i) = -g(k,i)-g(k+1,i) - f(k,i)+f(k+1,i)
          r(k,i) = g(k+1,i) + f(k+1,i)
          rhs(k,i) = cptn(k,i)
        enddo ! k=lev0,lev1-2
! Level k-1/2
        p(lev1-1,i) =  g(lev1-1,i)-f(lev1-1,i)
        q(lev1-1,i) = -g(lev1-1,i)-f(lev1-1,i)
        r(lev1-1,i) = 0.
        rhs(lev1-1,i) = cptn(lev1-1,i)
      enddo ! i=lon0,lon1

!     call addfld('P_COEF0' ,' ',' ',p,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('Q_COEF0' ,' ',' ',q,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('R_COEF0' ,' ',' ',r,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('RHS0'    ,' ',' ',rhs,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('COOL_IMP',' ',' ',cool_imp(:,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('COOL_EXP',' ',' ',cool_exp(:,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! qpart = cp*(1/(2*Dt)+ai/cp+w*R/(cp*M))
      do i=lon0,lon1
        do k=lev0,lev1-1
          qpart(k,i) =                                             ! s13
     |      .5*(cp(k,i,lat)+cp(k+1,i,lat))*(dtx2inv+cool_imp(k,i,lat))+
     |      .5*(w_upd(k,i,lat)+w_upd(k+1,i,lat))*gask/mbar(k,i) 
          rhs(k,i) = rhs(k,i)+cool_exp(k,i,lat)                    ! s12
          q(k,i) = q(k,i)-expz(k)*qpart(k,i)                       ! s10
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('QPART'  ,' ',' ',qpart,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('Q_COEF1',' ',' ',q    ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('RHS1'   ,' ',' ',rhs  ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     do k=lev0,lev1
!       tnlbc_diag(k,:) = tnlbc(:,lat)
!     enddo
!     call addfld('TNLBC',' ',' ',tnlbc_diag,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Lower boundary:
      do i=lon0,lon1
        q(lev0,i) = q(lev0,i)-p(lev0,i)
!
! Diffs in rhs lbc ??:
        rhs(lev0,i) = rhs(lev0,i)-2.*p(lev0,i)*tnlbc(i,lat)
        p(lev0,i) = 0.
      enddo ! i=lon0,lon1

!     call addfld('P_COEF2',' ',' ',p,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('Q_COEF2',' ',' ',q,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('R_COEF2',' ',' ',r,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('RHS2'   ,' ',' ',rhs,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)

!
! Solve tridiagonal system for new tn:
!     subroutine trsolv(a,b,c,f,x,lev0,lev1,k1,k2,lon0,lon1,lonmax,lat,
!    |  idebug)
!
      call trsolv(p,q,r,rhs,tn_upd(:,lon0:lon1,lat),lev0,lev1,
     |  lev0,lev1-1,lon0,lon1,nlonp4,lat,0)

!     call addfld('TN_SOLV',' ',' ',tn_upd(:,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! End second latitude scan:
      enddo ! lat=lat0,lat1
!
! Set kut for wave filtering according to dlat (2.5 or 5.0):
!     call set_wave_filter(36,kut_5,nlat,kutt)
!
! Filter updated tn:
      call filter_tn(tn_upd,lev0,lev1,lon0,lon1,lat0,lat1,kut)
!
! Third latitude scan:
      do lat=lat0,lat1
!       call addfld('TN_FILT',' ',' ',tn_upd(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Smooth updated tn:
        do i=lon0,lon1 
          do k=lev0,lev1-1
            tn_nm_upd(k,i,lat) = dtsmooth_div2*(tn_nm(k,i,lat)+
     |        tn_upd(k,i,lat)) + dtsmooth*tn(k,i,lat)
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1 
!       call addfld('TN_NMOUT',' ',' ',tn_nm_upd(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! 10/10/06 btf: upper boundary top slot nlevp1 <- lev1-1
        tn_upd(lev1,lon0:lon1,lat)    = tn_upd(lev1-1,lon0:lon1,lat)
        tn_nm_upd(lev1,lon0:lon1,lat) = tn_nm_upd(lev1-1,lon0:lon1,lat)
!
! save LB for tlb c= tn(itp) in tlbc_nm
! Save lower boundary for tn(itc) in tlbc for history:
! Lower boundary is in tlbc (fields.F):
        tlbc_nm(lon0:lon1,lat) = tlbc(lon0:lon1,lat)
        tlbc(lon0:lon1,lat)    = tnlbc(:,lat)
!
! Convective adjustment (timegcm only):
        call conv_adjust(tn_upd(:,:,lat),cp(:,:,lat),barm(:,:,lat),
     |    lev0,lev1,lon0,lon1,lat)
!
! Tn must be at least 100 K:
        do i=lonbeg,lonend
          do k=lev0,lev1
            if (tn_upd(k,i,lat) < 100.) tn_upd(k,i,lat) = 100.
          enddo
        enddo
!
#ifdef MPI
!
! 5/4/06 btf: These calls opened up as per Wenbin's suggestion.
!
! Define halo longitudes in tn_upd, for use by duv:
        call mp_bndlons_f3d(tn_upd,nlevs,lon0,lon1,lat0,lat1,1,0)
!
! Periodic points:
        call mp_periodic_f3d(tn_upd(:,lon0:lon1,lat0:lat1),
     |    lev0,lev1,lon0,lon1,lat0,lat1,1)
#endif
!
! Save final temperature:
!       call addfld('TN_FINAL',' ',' ',tn_upd(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Save total radiative cooling (deg/day):
        do i=lon0,lon1
          do k=lev0,lev1-1
            radcool(k,i) = cool_exp(k,i,lat) / expz(k) +
     |        cool_imp(k,i,lat) * cp(k,i,lat) * tn_upd(k,i,lat)
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1
        do i=lon0,lon1
          do k=lev0,lev1-1
            radcool(k,i) = radcool(k,i)*86400./cp(k,i,lat) ! deg/day
          enddo ! k=lev0,lev1-1
          radcool(lev1,i) = spval
        enddo ! i=lon0,lon1

        if (thermo_diag) then
          do i=lon0,lon1
            do k=lev0,lev1-1
              radnet(k,i) = qsolar(k,i,lat)+qjoule(k,i,lat)-radcool(k,i)
              totfab(k,i) = qsolar(k,i,lat)-radcool(k,i)+htflxm(k,i,lat)
     |          +adiaht(k,i,lat)-horadv(k,i,lat)-veradv(k,i,lat)+
     |           hordif(k,i,lat)+htflxe(k,i,lat)+qjoule(k,i,lat)
            enddo ! k=lev0,lev1-1
            totfab(lev0:lev0+1,i) = spval
            totfab(lev1-2:lev1,i) = spval
          enddo ! i=lon0,lon1
          radnet(lev1,:) = spval
	  call addfld('RADNET',' ',' ',radnet,
     |      'lev',lev0,lev1,'lon',lon0,lon1,lat)
	  call addfld('TOTFAB',' ',' ',totfab,
     |      'lev',lev0,lev1,'lon',lon0,lon1,lat)
        endif
!
!       call addfld('RADCOOL',' ','deg/day',radcool,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('COOL_IMP',' ',' ',cool_imp(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('COOL_EXP',' ',' ',cool_exp(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! End third lat scan:
      enddo ! lat=lat0,lat1
      end subroutine dt
!-----------------------------------------------------------------------
      subroutine filter_tn(fout,lev0,lev1,lon0,lon1,lat0,lat1,kut)
!
! Filter updated TN:
!
      use params_module,only: nlat,nlonp4,nlon
      use filter_module,only: filter
#ifdef MPI
      use mpi_module,only: mp_gatherlons_f3d,mp_scatterlons_f3d,mytidi
      implicit none
#else
      implicit none
      integer :: mytidi=0
#endif
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1,kut(nlat)
      real,intent(inout) :: fout(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2)
!
! VT vampir tracing:
!
#ifdef VT
#include <VT.inc>
#endif
!
! Local:
      integer :: i,j,nlevs,nlons,nlats
      real :: fik(nlonp4,lev0:lev1),fkij(lev0:lev1,nlonp4,lat0:lat1)
      real :: fmin,fmax
!
#ifdef VT
!     code = 131 ; state = 'filter_tn' ; activity='Filtering'
      call vtbegin(131,ier)
#endif
!
      nlevs = lev1-lev0+1
      nlons = lon1-lon0+1
      nlats = lat1-lat0+1
!
! Define lons in w_ki from current task:
      fkij = 0.
      do j=lat0,lat1
        do i=lon0,lon1
          fkij(:,i,j) = fout(:,i,j)
        enddo
      enddo ! j=lat0,lat1
!
#ifdef MPI
!
! Gather longitudes into tasks in first longitude column of task table
!   (leftmost of each j-row) for global fft. (i.e., tasks with mytidi==0 
!   gather lons from other tasks in that row). This includes all latitudes.
!
      call mp_gatherlons_f3d(fkij,lev0,lev1,lon0,lon1,lat0,lat1,1)
#endif
!
! Only leftmost tasks at each j-row of tasks does the global filtering:
      if (mytidi==0) then
!
! Define 2d array with all longitudes for filter at each latitude:
        latscan: do j=lat0,lat1
          if (kut(j) >= nlon/2) cycle latscan
          do i=1,nlonp4
            fik(i,:) = fkij(:,i,j)
          enddo ! i=1,nlonp4
!
! Remove wave numbers > kut(lat):
          call filter(fik,lev0,lev1,kut(j),j)
!
! Return filtered array to fkij:
          do i=1,nlonp4
            fkij(:,i,j) = fik(i,:)
          enddo ! i=1,nlonp4
        enddo latscan ! j=lat0,lat1
      endif ! mytidi==0
#ifdef MPI
!
! Now leftmost task at each j-row must redistribute filtered data
! back to other tasks in the j-row (mytidi>0,mytidj) (includes latitude):
!
      call mp_scatterlons_f3d(fkij,lev0,lev1,lon0,lon1,lat0,lat1,1)
#endif
!
! Return filtered array to fout at current task longitudes and latitudes:
      do j=lat0,lat1
        do i=lon0,lon1
          fout(:,i,j) = fkij(:,i,j)
        enddo
      enddo
!
#ifdef VT
!     code = 131 ; state = 'filter_tn' ; activity='Filtering'
      call vtend(131,ier)
#endif
      end subroutine filter_tn
!-----------------------------------------------------------------------
      subroutine conv_adjust(tnupd,cp,barm,
     |    lev0,lev1,lon0,lon1,lat)
      use params_module,only: nlonp4
      use cons_module,only: expz,p0,gask
!
! Convective adjustment to tn (time-gcm only).
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(inout) ::
     |  tnupd
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  cp,   ! specific heat
     |  barm  ! mean molecular weight
!
! Local:
      integer :: i,k,lonbeg,lonend,nk,nkm1
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  pi,
     |  q,
     |  beta,
     |  tnadj,
     |  theta
      integer :: doadjust(lon0:lon1)
!
      nk = lev1-lev0+1
      nkm1 = nk-1
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = nlonp4-2
!
!   pi(1) = (p(1))**(g/cp(1))
!   pi(k) = pi(k-1)*(p(k)/p(k-1))**(g/cp(k-1))
!
      pi = 0. ! init
      do i=lonbeg,lonend
        pi(lev0,i) = 1.
!
! Apparent bug in tgcm24: NMSK is init to NJ+NMS+1 (k==2), but is
!   not incremented in the DO 2 loop (analogous to this k-loop).
! As of 3/17/04, the loop below is set to duplicate tgcm24 results.
!   To correct the bug, use the statement below which uses barm(k,i) 
!   rather than the statement that uses barm(lev0+1,i).
! However, I am not sure why barm is included at all in calculation 
!   of pi (see expression (6) for II (here named pi) in Akmaev (1991), 
!   cited in sub adjust below)
! 3/25/04: Using corrected tgcm24 bug, i.e., increment barm
!          (tgcm24 corrected in ~foster/timegcm/tgcm24/modsrc.valid/tstadj.F) 
!
        do k=lev0+1,lev1-1
          pi(k,i) = pi(k-1,i)*                                   ! s1
!    |      (expz(k)/expz(k-1))**(gask/(cp(k+1,i)*barm(lev0+1,i))) ! tgcm24 bug
     |      (expz(k)/expz(k-1))**(gask/(cp(k,i)*barm(k,i)))        ! bug fixed
        enddo ! k=lev0+1,lev1-1
      enddo ! i=lonbeg,lonend
!     call addfld('TADJ_CP',' ',' ',cp(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('TADJ_MB',' ',' ',barm(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('TADJ_PI',' ',' ',pi,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! q(1) = pi(1)*(p(1)-p(2))
! q(k) = pi(k)*(p(k-1)-p(k+1))
!
      q = 0. ! init
      do i=lonbeg,lonend
        q(lev0,i)   = pi(lev0,i)*p0*(expz(lev0)-expz(lev0+1))      ! s2 lbc
        q(lev1-1,i) = pi(lev1-1,i)*p0*(expz(lev1-2)-expz(lev1-1))  ! s2 ubc
        do k=lev0+1,lev1-2
          q(k,i) = pi(k,i)*p0*(expz(k-1)-expz(k+1))              ! s2
        enddo ! k=lev0-1,lev1-1
      enddo ! i=lonbeg,lonend
!     call addfld('TADJ_Q',' ',' ',q,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! beta(k) = 1./pi(k)
! theta = T/pi = potential temperature
!
      beta = 0. ; tnadj = 0. ; theta = 0.
      do i=lonbeg,lonend
        do k=lev0,lev1-1
          beta(k,i) = 1./pi(k,i)                                 ! s1
          tnadj(k,i) = tnupd(k,i)                                ! s3
          theta(k,i) = beta(k,i)*tnadj(k,i)                      ! s4
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lonend
!     call addfld('TADJBETA',' ',' ',beta     ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('TNPREADJ',' ',' ',tnadj    ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('ADJTHETA',' ',' ',theta    ,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Identify which columns require convective adjustment:
!
      doadjust = 0  ! init whole array (t1 in tgcm24)
      do i=lonbeg,lonend
        do k=lev0,lev1-2
          if (theta(k+1,i) < theta(k,i)) doadjust(i) = 1 ! adjust this column
        enddo ! k=lev0,lev1-2
      enddo ! i=lon0,lon1
!
! Perform convective ajustment where needed:
!
      do i=lonbeg,lonend
        if (doadjust(i)==1) then
!         write(6,"('conv_adjust calling adjust: i=',i3,' lat=',i3,
!    |      ' beta=',/,(6e12.4))")
!    |      i,lat,beta(:,i)
! 5/17/04 btf: Changed call to adjust to avoid div by 0 with beta(lev1)==0:
!         call adjust(tnadj(:,i),q(:,i),beta(:,i),lev0,lev1)
          call adjust(tnadj(:,i),q(:,i),beta(:,i),lev0,lev1-1) 
!
! Apply newly adjusted temperature:
          do k=lev0,lev1-1
            tnupd(k,i) = tnadj(k,i)
          enddo
        endif ! doadjust(i)==1
      enddo ! i=lon0,lon1
      end subroutine conv_adjust
!-----------------------------------------------------------------------
      subroutine adjust(temp,q,beta,lev0,lev1)
!
! Perform convective adjustment on temperature.
! See Akmaev, R.A., Monthly Weather Review Note, October 1991, p.2499-2504
!
! 3/17/04 btf: 
!   Attempted to disentangle spaghetti code in old source (tgcm24/adjust.F), 
!   but gave up. It is reproduced below with minor changes.
!
! Args:
      integer,intent(in) :: lev0,lev1
      real,intent(inout) :: temp(lev0:lev1)
      real,intent(in) :: q(lev0:lev1),beta(lev0:lev1)
!
! Local:
      integer :: k,nn,l,kmax
      real :: ss,tt,ttheta
      real,dimension(lev1-lev0+1) :: s, t, theta
      integer :: n(lev1-lev0+1)
!
      kmax = lev1-lev0+1

    1 K = lev0
      N(1) = lev0
      L = lev0+1
      THETA(lev0) = BETA(lev0)*TEMP(lev0)
    2 NN = 1
      TTHETA = BETA(L)*TEMP(L)
    3 IF(THETA(K).LE.TTHETA)GO TO 6
      IF(NN.GT.1)GO TO 4
      SS = Q(L)
      TT = SS*TTHETA
    4 IF(N(K).GT.1)GO TO 5
      S(K) = Q(L-NN)
      T(K) = S(K)*THETA(K)
    5 NN = NN+N(K)
      S(K) = SS+S(K)
      SS = S(K)
      T(K) = TT+T(K)
      TT = T(K)
      TTHETA = TT/SS
      IF(K.EQ.1)GO TO 7
      K = K-1
      GO TO 3
    6 K = K+1
    7 IF(L.EQ.KMAX)GO TO 8
      L = L+1
      N(K) = NN
      THETA(K) = TTHETA
      GO TO 2
    8 IF(NN.EQ.1)GO TO 11
    9 TEMP(L) = TTHETA/BETA(L)
      IF(NN.EQ.1)GO TO 11
   10 L = L-1
      NN = NN-1
      GO TO 9
   11 IF(K.EQ.1)RETURN
      K = K-1
      L = L-1
      NN = N(K)
      TTHETA = THETA(K)
      GO TO 8

      end subroutine adjust
!-----------------------------------------------------------------------
      end module dt_module
