!
      subroutine duv(tn,tn_upd,tn_nm,un,vn,un_nm,vn_nm,w_upd,o2,o1,n2,
     |  barm,z,hdu,hdv,ui,vi,lxx,lyy,lxy,lyx,km,un_upd,unm_upd,vn_upd,
     |  vnm_upd,lev0,lev1,lon0,lon1,lat0,lat1)
!
! This software is part of the NCAR TIE-GCM.  Use is governed by the 
! Open Source Academic Research License Agreement contained in the file 
! tiegcmlicense.txt.
!
! Advance neutral velocities for current time step:
!
      use params_module,only: nlonp4,dz,nlat,spval,nlonp2,nlon
      use init_module,only: iter,igetgswm,iday
      use lbc,only: u_lbc,v_lbc
      use cons_module,only: dtx2inv,expz,p0,gask,grav,
     |  xmue,cor,tanphi=>tn,re,dtsmooth_div2,dtsmooth,kut,
     |  rmassinv_n2,rmassinv_o2,rmassinv_o1
      use fields_module,only: tlbc,ulbc,vlbc,tlbc_nm,ulbc_nm,vlbc_nm
      use addfld_module,only: addfld
#ifdef MPI
      use mpi_module,only: mp_periodic_f3d,mp_periodic_f2d,
     |  mp_bndlons_f2d,mp_bndlats_f2d
#endif
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
!
! Inputs at full subdomain:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in)::
     |  tn,      ! neutral temperature (deg K)
     |  tn_upd,  ! updated neutral temperature (from dt.F)
     |  tn_nm,   ! neutral temperature from time n-1
     |  un,      ! neutral zonal velocity
     |  vn,      ! neutral meridional velocity
     |  un_nm,   ! neutral zonal velocity at time n-1
     |  vn_nm,   ! neutral meridional velocity at time n-1
     |  w_upd,   ! updated vertical velocity (swdot.F)
     |  o2,      ! o2 (for mbar only)
     |  o1,      ! o1 (for mbar only)
     |  n2,      ! n2 (for mbar only)
     |  barm,    ! mean molecular weight
     |  z,       ! geopotential height
     |  hdu,     ! horizontal diffusion of U (hdif3 in hdif.F)
     |  hdv,     ! horizontal diffusion of V (hdif3 in hdif.F)
     |  ui,      ! zonal ion velocity input
     |  vi,      ! meridional ion velocity input
     |  lxx,     ! xx ion drag coefficient
     |  lyy,     ! yy ion drag coefficient
     |  lxy,     ! xy ion drag coefficient
     |  lyx,     ! yx ion drag coefficient
     |  km       ! molecular viscosity (cpktkm.F)
!
! Outputs at full subdomain:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out)::
     |  un_upd,  ! updated zonal velocity
     |  unm_upd, ! updated zonal velocity at time n-1
     |  vn_upd,  ! updated meridional velocity
     |  vnm_upd  ! updated meridional velocity at time n-1
!
! Local:
      integer :: k,kk,i,j,lat,lonbeg,lonend
      integer :: nk,nkm1,nlevs
      real :: 
     |  unlbc(lon0:lon1), ! un lower boundary (t1)
     |  vnlbc(lon0:lon1)  ! vn lower boundary (t2)
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  g,           ! g**2*(kM+XMUE)*MBAR*/(P0*R*T*DS**2)
     |  dwdz,        ! exp(-s)*w(k+1/2)/(2.*Ds)
     |  tni          ! tn at interfaces
      real,dimension(lev0:lev1,lon0:lon1,lat0:lat1) ::
     |  advec_un,    ! zonal advection (output of sub advec)
     |  advec_vn,    ! meridional advection (output of sub advec)
     |  zl,zp,       ! output from glp (and dldp)
     |  unm_smooth,  ! un at time n-1, smoothed
     |  vnm_smooth   ! vn at time n-1, smoothed
      real :: flbc(lon0-2:lon1+2,lat0-2:lat1+2,4) ! ulbc,vlbc,ulbc_nm,vlbc_n
      real :: rtxmue(lev0:lev1) ! sqrt(xmue)
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  eddyvisc,    ! eddy viscosity (?)
     |  mbar         ! mean molecular weight (slightly different than barm)
!
! For diagnostic, debug:
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  unlbc_diag,   ! un lbc redundant in vertical, for diagnostic
     |  vnlbc_diag,   ! vn lbc redundant in vertical, for diagnostic
     |  ss_un,ss_vn,  ! also for diagnostics
     |  pp11,pp12,pp21,pp22, qq11,qq12,qq21,qq22, rr11,rr12,rr21,rr22
!
! pp, qq, rr, beta AND gamma are (2 x 2) matrices defined for 
!   k = 3/2, K-1/2, 1 and i = 3, IMAX+2, 1
      real,dimension(2,2,lev0:lev1,lon0:lon1) :: 
     |  pp,qq,rr,beta,gamma 
!
! ss, xx and yy are two component vectors similarly defined.
      real,dimension(2,lev0:lev1,lon0:lon1) :: 
     |  ss,xx,yy
!
! For addfld:
      nk = lev1-lev0+1
      nkm1 = nk-1
      nlevs = nk
!
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = nlonp4-2
!
! First latitude scan:
      do lat=lat0,lat1
!
! Horizontal advection:
!     subroutine advec(f,hadvec,lev0,lev1,lon0,lon1,lat)
!     real,dimension(lev0:lev1,lon0-2:lon1+2,lat-2:lat+2),intent(in) :: 
!    |  f ! input field with ghost cells for finite differencing
!     real,dimension(lev0:lev1,lon0:lon1),intent(out) :: hadvec
!      
        call advec(un(:,:,lat-2:lat+2),advec_un(:,:,lat), 
     |    lev0,lev1,lon0,lon1,lat)
        call advec(vn(:,:,lat-2:lat+2),advec_vn(:,:,lat), 
     |    lev0,lev1,lon0,lon1,lat)

!       call addfld('ADVEC_U0','ADVEC_U0',' ',advec_un(:,:,lat),
!    |    'ilev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('ADVEC_V0','ADVEC_V0',' ',advec_vn(:,:,lat),
!    |    'ilev',lev0,lev1,'lon',lon0,lon1,lat)
!
! End first latitude scan:
      enddo ! lat=lat0,lat1
!
! Horizontal pressure forcing output of glp is in zl,zp. Pass full 
! subdomains for derivatives in lat and lon (see sub dldp, called by glp):
!
      call glp(tn,tn_upd,tn_nm,barm,z,zl,zp,
     |  lev0,lev1,lon0,lon1,lat0,lat1)
!
! Second latitude scan:
      do lat=lat0,lat1

       call addfld('ZL','ZL',' ',zl(:,:,lat),
     |    'ilev',lev0,lev1,'lon',lon0,lon1,lat)
       call addfld('ZP','ZP',' ',zp(:,:,lat),
     |    'ilev',lev0,lev1,'lon',lon0,lon1,lat)

        do i=lon0,lon1
          do k=lev0,lev1-1
            advec_un(k,i,lat) = advec_un(k,i,lat)+
     |        zl(k,i,lat)-hdu(k,i,lat)
            advec_vn(k,i,lat) = advec_vn(k,i,lat)+
     |        zp(k,i,lat)-hdv(k,i,lat)
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1

!       call addfld('ADVEC_U1','ADVEC_U1',' ',advec_un(:,:,lat),
!    |     'ilev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('ADVEC_V1','ADVEC_V1',' ',advec_vn(:,:,lat),
!    |     'ilev',lev0,lev1,'lon',lon0,lon1,lat)
!
! End second lat scan:
      enddo ! lat=lat0,lat1
!
! Smooth un,vn at time n-1 with 2-part shapiro smoother.
! Input un_nm,vn_nm at full subdomain.
! Output in unm_smooth,vnm_smooth at subdomain.
!
      call smooth(un_nm,unm_smooth,lev0,lev1,lon0,lon1,lat0,lat1,0)
      call smooth(vn_nm,vnm_smooth,lev0,lev1,lon0,lon1,lat0,lat1,0)
!
! Start third latitude scan:
      do lat=lat0,lat1

!       call addfld('UNSMOOTH','UNSMOOTH',' ',unm_smooth(:,:,lat),
!    |     'ilev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('VNSMOOTH','VNSMOOTH',' ',vnm_smooth(:,:,lat),
!    |     'ilev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('UN_NM_DUV','UN_NM_DUV',' ',un_nm(:,lon0:lon1,lat),
!    |     'ilev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('VN_NM_DUV','VN_NM_DUV',' ',vn_nm(:,lon0:lon1,lat),
!    |     'ilev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Lower boundaries u_lbc,v_lbc were calculated by sub tuvz_lbc (lbc.F)
        unlbc(:) = u_lbc(lon0:lon1,lat)
        vnlbc(:) = v_lbc(lon0:lon1,lat)
        do k=lev0,lev1
          unlbc_diag(k,:) = unlbc(:)
          vnlbc_diag(k,:) = vnlbc(:)
        enddo
!
!       call addfld('UNLBC_DUV','UNLBC_DUV',' ',unlbc_diag,
!    |    'lon',lon0,lon1,'lat',lat0,lat1,0)
!       call addfld('VNLBC_DUV','VNLBC_DUV',' ',vnlbc_diag,
!    |    'lon',lon0,lon1,'lat',lat0,lat1,0)

        do i=lonbeg,lonend
          do k=lev0,lev1-1
            advec_un(k,i,lat) = advec_un(k,i,lat)-
     |        dtx2inv*unm_smooth(k,i,lat)
            advec_vn(k,i,lat) = advec_vn(k,i,lat)-
     |        dtx2inv*vnm_smooth(k,i,lat)
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1

!       call addfld('ADVEC_U2','ADVEC_U2',' ',advec_un(:,:,lat),
!    |    'ilev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('ADVEC_V2','ADVEC_V2',' ',advec_vn(:,:,lat),
!    |    'ilev',lev0,lev1,'lon',lon0,lon1,lat)

        do i=lon0,lon1
          do k=lev0,lev1-1
            ss(1,k,i) = expz(k)*(.25*((lxx(k,i,lat)+lxx(k+1,i,lat))*
     |        (ui(k,i,lat)+ui(k+1,i,lat))+(lxy(k,i,lat)+lxy(k+1,i,lat))*
     |        (vi(k,i,lat)+vi(k+1,i,lat)))-advec_un(k,i,lat))
            ss(2,k,i) = expz(k)*(.25*((lyy(k,i,lat)+lyy(k+1,i,lat))*
     |        (vi(k,i,lat)+vi(k+1,i,lat))-(lyx(k,i,lat)+lyx(k+1,i,lat))*
     |        (ui(k,i,lat)+ui(k+1,i,lat)))-advec_vn(k,i,lat))
            ss_un(k,i) = ss(1,k,i)
            ss_vn(k,i) = ss(2,k,i)
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1

!       call addfld('SS_UN',' ',' ',ss_un,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('SS_VN',' ',' ',ss_vn,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)

!
! xmue = eddy viscosity (cons.F). rtxmue = sqrt(xmue)
        do k=lev0+1,lev1-1          
          rtxmue(k) = sqrt(xmue(k-1,iday)*xmue(k,iday))
        enddo ! k=lev0+1,lev1-1          
!
! tni = tn at interfaces:
        do i=lon0,lon1
          do k=lev0+1,lev1-1          
            tni(k,i) = .5*(tn(k-1,i,lat)+tn(k,i,lat))
          enddo ! k=lev0+1,lev1-1          
        enddo ! i=lon0,lon1
!
! Boundaries:
        rtxmue(lev0) = sqrt(xmue(lev0,iday)**3/
     |                      xmue(lev0+1,iday))
        rtxmue(lev1) = sqrt(xmue(lev1-1,iday)**3/
     |                      xmue(lev1-2,iday))
!
        do i=lon0,lon1
          do k=lev0,lev1-1
            mbar(k,i) = 1./(o2(k,i,lat)*rmassinv_o2 + 
     |        o1(k,i,lat)*rmassinv_o1 + n2(k,i,lat)*rmassinv_n2)
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1
!
        do i=lon0,lon1
          do k=lev0,lev1-1          
            eddyvisc(k,i) = p0*expz(k)*gask*tn(k,i,lat)*xmue(k,iday)/
     |                      (mbar(k,i)*grav**2)
          enddo ! k=lev0+1,lev1-1          
!
! Upper boundary:
          eddyvisc(lev1,i) = sqrt(eddyvisc(lev1-1,i)**3/             ! s3 top
     |                       eddyvisc(lev1-2,i))
!            
! Loop from top down:
          k = lev1+1
          do kk=lev0+1,lev1
            k = k-1 
            eddyvisc(k,i) = sqrt(eddyvisc(k-1,i)*eddyvisc(k,i))      ! s3
          enddo ! kk=lev0+1,lev1-1
!
! Bottom boundary in tlbc (from t_lbc of lbc.F):
          eddyvisc(lev0,i) = eddyvisc(lev0,i)**2/eddyvisc(lev0+1,i)  ! s3
        enddo ! i=lon0,lon1
!
        do i=lon0,lon1
          tni(lev0,i) = tlbc_nm(i,lat)    ! Lower boundary is in tlbc_nm
          tni(lev1,i) = tn(lev1-1,i,lat)

! G = g*KM/(P0*H*Ds**2) = g**2*(kM+XMUE)*MBAR*/(P0*R*T*DS**2)
          do k=lev0,lev1
            g(k,i) = grav**2*(km(k,i,lat)+eddyvisc(k,i))*barm(k,i,lat)/
     |        (p0*gask*tni(k,i)*dz**2)
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1

!       call addfld('TNI'  ,' ',' ',tni,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('DUV_G',' ',' ',g  ,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)

!
! dwdz = exp(-s)*w(k+1/2)/(2.*Ds)
        do i=lon0,lon1
          do k=lev0,lev1-1
            dwdz(k,i) = 0.5*expz(k)*(w_upd(k,i,lat)+w_upd(k+1,i,lat))/
     |        (2.*dz)
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1

!       call addfld('DWDZ',' ',' ',dwdz,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)

!
! P = -(G(k-1/2,n) + exp(-s)*w(k,n)/   |1., 0.|
!                            (2.*Ds)) *|0., 1.|
!                                      |_    _|
!
! R = -(G(k+1/2,n) - exp(-s)*w(k,n)/   |1., 0.|
!                            (2.*Ds)) *|0., 1.|
!                                      |_    _|
!  for k = 3/2, K-1/2, 1
!
        do i=lon0,lon1
          do k=lev0,lev1-1
            pp(1,1,k,i) = -(g(k  ,i)+dwdz(k,i))
            rr(1,1,k,i) = -(g(k+1,i)-dwdz(k,i))
            pp(2,2,k,i) = pp(1,1,k,i)
            rr(2,2,k,i) = rr(1,1,k,i)
            pp(1,2,k,i) = 0.
            rr(1,2,k,i) = 0.
            pp(2,1,k,i) = 0.
            rr(2,1,k,i) = 0.
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1
!
! Now calculate Q from (5) for k = 3/2, K-1/2, 1
!  Q(1,1,k) = G(k-1/2) + G(k+1/2) + exp(-s)*(1./(2.*Dt) + lamxx(k))
!    and
!  Q(2,2,k) = G(k-1/2) + G(k+1/2) + exp(-s)*(1./(2.*Dt) + lamyy(k))
!
        do i=lon0,lon1
          do k=lev0,lev1-1
            qq(1,1,k,i) = g(k,i) + g(k+1,i)
            qq(2,2,k,i) = qq(1,1,k,i)
            qq(1,1,k,i) = qq(1,1,k,i) + 
     |        expz(k)*(dtx2inv+.5*(lxx(k,i,lat)+lxx(k+1,i,lat)))
            qq(2,2,k,i) = qq(2,2,k,i) +
     |        expz(k)*(dtx2inv+.5*(lyy(k,i,lat)+lyy(k+1,i,lat)))
!
!     Q(1,2) = -exp(-s) * (f+u/r*tan(theta)-lamxy)
!      and 
!     Q(2,1) = exp(-s) * (f+u/r*tan(theta)-lamyx)
!
            qq(1,2,k,i) = cor(lat)+un(k,i,lat)/re*tanphi(lat)
            qq(2,1,k,i) = qq(1,2,k,i)
            qq(1,2,k,i) = -expz(k)*(qq(1,2,k,i)-
     |        .5*(lxy(k,i,lat)+lxy(k+1,i,lat)))
            qq(2,1,k,i) =  expz(k)*(qq(2,1,k,i)-
     |        .5*(lyx(k,i,lat)+lyx(k+1,i,lat)))
          enddo ! k=lev0,lev1-1
!
! Boundaries at k=3/2 and k=k-1/2
! Lower boundary:
!     Q(3/2) = Q(3/2) - P(3/2)
!     S(3/2) = S(3/2) -2.*P(3/2)*Vb
!     P(3/2) = 0.
!
          qq(1,1,lev0,i) = qq(1,1,lev0,i)-pp(1,1,lev0,i)
          qq(1,2,lev0,i) = qq(1,2,lev0,i)-pp(1,2,lev0,i)
          qq(2,1,lev0,i) = qq(2,1,lev0,i)-pp(1,2,lev0,i)
          qq(2,2,lev0,i) = qq(2,2,lev0,i)-pp(2,2,lev0,i)
          ss(1,lev0,i) = ss(1,lev0,i)-2.*(pp(1,1,lev0,i)*unlbc(i)+
     |                                    pp(1,2,lev0,i)*vnlbc(i))
          ss(2,lev0,i) = ss(2,lev0,i)-2.*(pp(2,1,lev0,i)*unlbc(i)+
     |                                    pp(2,2,lev0,i)*vnlbc(i))
          ss_un(lev0,i) = ss(1,lev0,i)
          ss_vn(lev0,i) = ss(2,lev0,i)
          pp(1,1,lev0,i) = 0.
          pp(1,2,lev0,i) = 0.
          pp(2,1,lev0,i) = 0.
          pp(2,2,lev0,i) = 0.
!
! Upper boundary:
!     Q(K-1/2) = Q(K-1/2) + R(K-1/2)
!     R(K-1/2) = 0.
!
          qq(1,1,lev1-1,i) = qq(1,1,lev1-1,i)+rr(1,1,lev1-1,i)
          qq(1,2,lev1-1,i) = qq(1,2,lev1-1,i)+rr(1,2,lev1-1,i)
          qq(2,1,lev1-1,i) = qq(2,1,lev1-1,i)+rr(2,1,lev1-1,i)
          qq(2,2,lev1-1,i) = qq(2,2,lev1-1,i)+rr(2,2,lev1-1,i)
          rr(1,1,lev1-1,i) = 0.
          rr(1,2,lev1-1,i) = 0.
          rr(2,1,lev1-1,i) = 0.
          rr(2,2,lev1-1,i) = 0.
        enddo ! i=lon0,lon1

!       pp11(:,:) = pp(1,1,:,:) ; pp12(:,:) = pp(1,2,:,:)
!       pp21(:,:) = pp(2,1,:,:) ; pp22(:,:) = pp(2,2,:,:)
!       qq11(:,:) = qq(1,1,:,:) ; qq12(:,:) = qq(1,2,:,:)
!       qq21(:,:) = qq(2,1,:,:) ; qq22(:,:) = qq(2,2,:,:)
!       rr11(:,:) = rr(1,1,:,:) ; rr12(:,:) = rr(1,2,:,:)
!       rr21(:,:) = rr(2,1,:,:) ; rr22(:,:) = rr(2,2,:,:)
!       call addfld('PP11',' ',' ',pp11,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('PP12',' ',' ',pp12,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('PP21',' ',' ',pp21,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('PP22',' ',' ',pp22,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('QQ11',' ',' ',qq11,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('QQ12',' ',' ',qq12,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('QQ21',' ',' ',qq21,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('QQ22',' ',' ',qq22,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('RR11',' ',' ',rr11,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('RR12',' ',' ',rr12,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('RR21',' ',' ',rr21,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('RR22',' ',' ',rr22,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('SS_UN',' ',' ',ss_un,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('SS_VN',' ',' ',ss_vn,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Invoke tridiagonal (2x2) block matrix solver for U and V: 
! pp,qq,rr,ss are input, beta,gamma,yy,xx are output.
! (final solution is in xx)
!
        call blktri(pp,qq,rr,ss, beta,gamma,yy,xx,
     |    lev0,lev1,lon0,lon1,lat)
!
! Save updated U and V:
        do i=lonbeg,lonend 
          do k=lev0,lev1-1
            un_upd(k,i,lat) = xx(1,k,i)
            vn_upd(k,i,lat) = xx(2,k,i)
          enddo ! k=lev0,lev1-1
!
! Put spval in top nlevp1 level:
          un_upd(lev1,i,lat) = spval
          vn_upd(lev1,i,lat) = spval
!
! Lower boundary is in ulbc,vlbc (fields.F):
          ulbc_nm(i,lat) = ulbc(i,lat)  ! LB for un_nm(itc) = un(itp) 
          vlbc_nm(i,lat) = vlbc(i,lat)  ! LB for vn_nm(itc) = vn(itp) 
          ulbc(i,lat) = unlbc(i)        ! LB for un(itc)
          vlbc(i,lat) = vnlbc(i)        ! LB for vn(itc) 
        enddo ! i=lonbeg,lonend 

!       call addfld('UNLBC',' ',' ',unlbc_diag,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('VNLBC',' ',' ',vnlbc_diag,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)

!       call addfld('UN_SOLV',' ',' ',un_upd(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('VN_SOLV',' ',' ',vn_upd(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! End third latitude scan:
      enddo ! lat=lat0,lat1

!     call addfld('UNLBC_IJ','UNLBC_IJ','[cm/s]',unlbc_ij,
!    |  'lon',lon0,lon1,'lat',lat0,lat1,0)
!     call addfld('VNLBC_IJ','VNLBC_IJ','[cm/s]',vnlbc_ij,
!    |  'lon',lon0,lon1,'lat',lat0,lat1,0)
!
! Fourier smoothing of U and V:
      call filter_uv(un_upd(lev0:lev1-1,:,:),lev0,lev1-1,lon0,lon1,
     |  lat0,lat1,kut,'UN')
      call filter_uv(vn_upd(lev0:lev1-1,:,:),lev0,lev1-1,lon0,lon1,
     |  lat0,lat1,kut,'VN')
!
! Fourth latitude scan:
      do lat=lat0,lat1
!       call addfld('UN_FILT',' ',' ',un_upd(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('VN_FILT',' ',' ',vn_upd(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Smooth updated un,vn:
        do i=lon0,lon1 
          do k=lev0,lev1-1
            unm_upd(k,i,lat) = dtsmooth_div2*(un_nm(k,i,lat)+
     |        un_upd(k,i,lat)) + dtsmooth*un(k,i,lat)
            vnm_upd(k,i,lat) = dtsmooth_div2*(vn_nm(k,i,lat)+
     |        vn_upd(k,i,lat)) + dtsmooth*vn(k,i,lat)
          enddo ! k=lev0,lev1-1
!
! put spval in top nlevp1 level:
          unm_upd(lev1,i,lat) = spval
          vnm_upd(lev1,i,lat) = spval
        enddo ! i=lon0,lon1 
!       call addfld('UN_NMOUT',' ',' ',unm_upd(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('VN_NMOUT',' ',' ',vnm_upd(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! End fourth latitude scan:
      enddo ! lat=lat0,lat1
#ifdef MPI
!
! Boundary (halo) latitudes:
      flbc(:,:,1) = ulbc(:,:)
      flbc(:,:,2) = vlbc(:,:)
      flbc(:,:,3) = ulbc_nm(:,:)
      flbc(:,:,4) = vlbc_nm(:,:)
      call mp_bndlats_f2d(flbc,lon0,lon1,lat0,lat1,4)
      ulbc(:,:)    = flbc(:,:,1)
      vlbc(:,:)    = flbc(:,:,2)
      ulbc_nm(:,:) = flbc(:,:,3)
      vlbc_nm(:,:) = flbc(:,:,4)
!
! Boundary (halo) longitudes:
      flbc(:,:,1) = ulbc(:,:)
      flbc(:,:,2) = vlbc(:,:)
      flbc(:,:,3) = ulbc_nm(:,:)
      flbc(:,:,4) = vlbc_nm(:,:)
      call mp_bndlons_f2d(flbc,lon0,lon1,lat0,lat1,4)
      ulbc(:,:)    = flbc(:,:,1)
      vlbc(:,:)    = flbc(:,:,2)
      ulbc_nm(:,:) = flbc(:,:,3)
      vlbc_nm(:,:) = flbc(:,:,4)
!
! Periodic points for lbc:
      flbc(:,:,1) = ulbc(:,:)
      flbc(:,:,2) = vlbc(:,:)
      flbc(:,:,3) = ulbc_nm(:,:)
      flbc(:,:,4) = vlbc_nm(:,:)
      call mp_periodic_f2d(flbc(lon0:lon1,lat0:lat1,1:4),
     |       lon0,lon1,lat0,lat1,4)
      ulbc(:,:)    = flbc(:,:,1)
      vlbc(:,:)    = flbc(:,:,2)
      ulbc_nm(:,:) = flbc(:,:,3)
      vlbc_nm(:,:) = flbc(:,:,4)
#else
      do j=lat0,lat1
        do i=1,2
          ulbc(i,j)        = ulbc(nlon+i,j)
          ulbc(nlonp2+i,j) = ulbc(i+2,j) 
          vlbc(i,j)        = vlbc(nlon+i,j)
          vlbc(nlonp2+i,j) = vlbc(i+2,j)  
          ulbc_nm(i,j)        = ulbc_nm(nlon+i,j)
          ulbc_nm(nlonp2+i,j) = ulbc_nm(i+2,j) 
          vlbc_nm(i,j)        = vlbc_nm(nlon+i,j)
          vlbc_nm(nlonp2+i,j) = vlbc_nm(i+2,j)  
        enddo
      enddo    
#endif
!
!     do lat=lat0,lat1
!       call addfld('UN_FINAL',' ',' ',un_upd(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('VN_FINAL',' ',' ',vn_upd(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     enddo ! lat=lat0,lat1
      end subroutine duv
!-----------------------------------------------------------------------
      subroutine glp(tn,tn_upd,tn_nm,barm,z,zl,zp,
     |  lev0,lev1,lon0,lon1,lat0,lat1)
!
! Horizontal pressure forcing for un,vn: 
!
      use params_module,only: nlonp4,nlat
      use cons_module,only: dz,dzgrav,grav,racs,re_inv
      use init_module,only: iter
      use addfld_module,only: addfld
#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
!
! Inputs at full subdomain:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in)::
     |  tn,      ! neutral temperature (deg K)
     |  tn_upd,  ! updated neutral temperature (from dt.F)
     |  tn_nm,   ! neutral temperature from time n-1
     |  barm,    ! mean molecular weight
     |  z        ! geopotential
!
! Outputs at subdomain:
      real,dimension(lev0:lev1,lon0:lon1,lat0:lat1),intent(out) :: 
     |  zl,zp
!
! Local:
      integer :: k,i,nk,nkm1,nlevs,lat,lonbeg,lonend
      complex :: expt1,expt2,expta
      real,parameter :: wt=0.225
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2) ::
     |  ztmp
      real,dimension(lev0:lev1,lon0-2:lon1+2) ::
     |  tbar,barmi,dztbar
      real :: ftmp(lev0:lev1,lon0:lon1,lat0:lat1,2)
!
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = nlonp4-2

      nk = lev1-lev0+1 ; nkm1=nk-1 ; nlevs = nk ! for addfld
!
! Latitude scan:
!
      do lat=lat0,lat1
        do i=lon0,lon1
          do k=lev0,lev1-1
            tbar(k,i) = tn(k,i,lat)*(-2.)+tn_nm(k,i,lat)     ! s6
            tbar(k,i) = (tbar(k,i)+tn_upd(k,i,lat))*wt
            tbar(k,i) = tbar(k,i)+tn(k,i,lat)
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1

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

        do i=lon0,lon1
          do k=lev0,lev1-1
            barmi(k,i) = .5*(barm(k,i,lat)+barm(k+1,i,lat))    ! s7
            dztbar(k,i) = (dz/dzgrav) * (tbar(k,i)/barmi(k,i)) ! s7
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1

!       call addfld('DZTBAR',' ',' ',dztbar(:,lon0:lon1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! In earlier model versions, geopotential bottom boundary was recalculated
! here. In this code, it is just copied from z, where it is already available:
        do i=lon0,lon1
          ztmp(lev0,i,lat) = z(lev0,i,lat)
        enddo ! i=lon0,lon1
!
! This will differ from z(..[itp,itc]) because tn_upd is in dztbar.
        do i=lon0,lon1
          do k=lev0,lev1-1
            ztmp(k+1,i,lat) = dztbar(k,i)+ztmp(k,i,lat)
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1
      enddo
!
#ifdef MPI
!
! Exchange boundary longitudes in ztmp for 4th order differencing in dldp. 
! Not necessary if ztmp has been defined at full subdomain, including lon 
! halos, but as of 4/02, this was not working (s.a. comments above)
!
      call mp_periodic_f3d(ztmp(:,lon0:lon1,lat0:lat1),
     |  lev0,lev1,lon0,lon1,lat0,lat1,1)
      call mp_bndlons_f3d(ztmp,nlevs,lon0,lon1,lat0,lat1,1,0)
#else
!
! Set periodic points for serial run:
      ztmp(:,1:2,:) = ztmp(:,lonend-2:lonend-1,:)
      ztmp(:,lonend+1:lonend+2,:) = ztmp(:,3:4,:)
#endif
!     do lat=lat0,lat1
!       call addfld('ZTMP',' ',' ',ztmp(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     enddo ! lat=lat0,lat1
!
! Calculate latitudinal and longitudinal derivatives in zl,zp:
! (dldp input ztmp is defined at full subdomain)
!
      call dldp(ztmp,z,zl,zp,lev0,lev1,lon0,lon1,lat0,lat1)
!
! Complete zl,zp at subdomain:
!
      do lat=lat0,lat1
        do i=lonbeg,lonend
          do k=lev0,lev1-1
            zl(k,i,lat) = (zl(k,i,lat)+zl(k+1,i,lat))*
     |        (.5*grav*racs(lat))
            zp(k,i,lat) = (zp(k,i,lat)+zp(k+1,i,lat))*
     |        (.5*grav*re_inv)
          enddo ! k=lev0,lev1-1
        enddo ! i=lonbeg,lonend
!
! Set periodic points zero to avoid NaNS fpe in duv.
! (earlier versions actually set periodic points here)
!        if (lon0==1) then
!          zl(:,lon0:lon0+1,lat) = 0.
!          zp(:,lon0:lon0+1,lat) = 0.
!        endif
!        if (lon1==nlonp4) then
!          zl(:,lon1-1:lon1,lat) = 0.
!          zp(:,lon1-1:lon1,lat) = 0.
!        endif
      enddo ! lat=lat0,lat1
!
! This part is modified to restore the periodic boundary conditions
! for zl and zp. zl and zp are horizontal pressure gradients. NaNs occurred before
! because values at periodic points (lon=1,2 and 75,76) of tn_upd in dt.F were not given. 
! They are used in subroutine dldp to preform fouth order finite difference. 
! Now the halo longtudinal periodic boundary for tn_upd is set in dt.F, 
! so the forth order differencing in dldp can be made betwen lon=3 to 74 
! (for 5 degree resolution), and so the values of periodic 
! points of zl and zp are now no longer NaNs.
!                                                   Wenbin  04/24/08
#ifdef MPI
      ftmp(:,:,:,1) = zl(:,lon0:lon1,lat0:lat1)
      ftmp(:,:,:,2) = zp(:,lon0:lon1,lat0:lat1)

      call mp_periodic_f3d(ftmp,lev0,lev1,lon0,lon1,lat0,lat1,2)

      zl(:,lon0:lon1,lat0:lat1) = ftmp(:,:,:,1)
      zp(:,lon0:lon1,lat0:lat1) = ftmp(:,:,:,2)
#else
!
! Set periodic points for serial run:
      zl(:,1:2,:) = zl(:,lonend-2:lonend-1,:)
      zl(:,lonend+1:lonend+2,:) = zp(:,3:4,:)
      zp(:,1:2,:) = zl(:,lonend-2:lonend-1,:)
      zp(:,lonend+1:lonend+2,:) = zp(:,3:4,:)
#endif
      end subroutine glp
!-----------------------------------------------------------------------
      subroutine dldp(zin,z,zl,zp,lev0,lev1,lon0,lon1,lat0,lat1)
      use params_module,only: nlonp4
      use cons_module,only: dlamda_1div12,dlamda_2div3,
     |  dphi_2div3,dphi_1div12
      use addfld_module,only: addfld
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
!
! Inputs at full subdomain:
! (Note new Z with new temperatures from glp is available only at the
!  current latitude, not adjacent latitudes. Therefore, the new Z 
!  from glp is used in longitudinal derivatives, whereas the current
!  Z (w/o updated temperature) is used in the latitudinal derivatives.
!  See ztmp in glp.)
!
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in)::
     |  zin,   ! new     Z (with updated temperature) from glp for zl
     |  z      ! current Z (with current temperature) from glp for zp
!
! Output at subdomain:
      real,dimension(lev0:lev1,lon0:lon1,lat0:lat1),intent(out) :: 
     |  zl,zp  ! output at subdomain
!
! Local:
      integer :: lonbeg,lonend,k,i,lat,nk
!
      nk = lev1-lev0+1
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = nlonp4-2
!
! Longitudinal derivatives:
      do lat=lat0,lat1

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

        do i=lonbeg,lonend
          do k=lev0,lev1
            zl(k,i,lat) = 
     |        dlamda_2div3 *(zin(k,i+1,lat)-zin(k,i-1,lat))-
     |        dlamda_1div12*(zin(k,i+2,lat)-zin(k,i-2,lat))
          enddo ! k=lev0,lev1

!         write(6,"('dldp zl: lat=',i2,' i=',i2,' lonbeg,end=',2i3)")
!    |      lat,i,lonbeg,lonend
!         write(6,"('zin(:,i+1,lat)=',/,(6e12.4))") zin(:,i+1,lat)
!         write(6,"('zin(:,i-1,lat)=',/,(6e12.4))") zin(:,i-1,lat)
!         write(6,"('zin(:,i+2,lat)=',/,(6e12.4))") zin(:,i+2,lat)
!         write(6,"('zin(:,i-2,lat)=',/,(6e12.4))") zin(:,i-2,lat)
!         write(6,"('zl(:,i,lat)=',/,(6e12.4))") zl(:,i,lat)

        enddo ! i=lonbeg,lonend

!
! Latitudinal derivatives:
        do i=lonbeg,lonend
          do k=lev0,lev1
            zp(k,i,lat) = 
     |        dphi_2div3 *(z(k,i,lat+1)-z(k,i,lat-1))-
     |        dphi_1div12*(z(k,i,lat+2)-z(k,i,lat-2))
          enddo ! k=lev0,lev1

!         write(6,"('dldp zp: lat=',i2,' i=',i2,' zin=',/,(6e12.4))") 
!    |      lat,i,z(:,i,lat)
!         write(6,"('dldp zp: lat=',i2,' i=',i2,' zin=',/,(6e12.4))") 
!    |      lat,i,zin(:,i,lat)
!         write(6,"('dldp zp: lat=',i2,' i=',i2)") lat,i
!         write(6,"('zin(:,i,lat+1)=',/,(6e12.4))") zin(:,i,lat+1)
!         write(6,"('zin(:,i,lat-1)=',/,(6e12.4))") zin(:,i,lat-1)
!         write(6,"('zin(:,i,lat+2)=',/,(6e12.4))") zin(:,i,lat+2)
!         write(6,"('zin(:,i,lat-2)=',/,(6e12.4))") zin(:,i,lat-2)
!         write(6,"('zp(:,i,lat)=',/,(6e12.4))") zp(:,i,lat)

        enddo ! i=lonbeg,lonend

!       call addfld('Z_JM2',' ',' ',z(:,lon0:lon1,lat-2),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('Z_JM1',' ',' ',z(:,lon0:lon1,lat-1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('Z_JP1',' ',' ',z(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('Z_JP2',' ',' ',z(:,lon0:lon1,lat+2),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('DLDP_ZL',' ',' ',zl(:,:,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       call addfld('DLDP_ZP',' ',' ',zp(:,:,lat),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)

      enddo ! lat=lat0,lat1
      end subroutine dldp
!-----------------------------------------------------------------------
      subroutine blktri(a,b,c,f,beta,gamma,y,x,lev0,lev1,lon0,lon1,lat)
!
! This procedure solves (I2-I1+1) tridiagonal block matrix
! systems in which all blocks are 2 x 2 matrices.
!
      use params_module,only: nlonp4 
      implicit none
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(2,2,lev0:lev1,lon0:lon1),intent(in) ::
     |  a,b,c
      real,dimension(2,lev0:lev1,lon0:lon1),intent(in) ::
     |  f
!
! Output args:
      real,dimension(2,2,lev0:lev1,lon0:lon1),intent(out) ::
     |  beta,gamma
      real,dimension(2,lev0:lev1,lon0:lon1),intent(out) ::
     |  x,y
!
! Local:
      integer :: k,i,lonbeg,lonend
!
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = nlonp4-2
!
! Lower boundary:
      do i=lonbeg,lonend
!
! Y(1,lev0,i) = determinant(B(lev0))
        y(1,lev0,i) = b(1,1,lev0,i)*b(2,2,lev0,i)-
     |                b(1,2,lev0,i)*b(2,1,lev0,i)
!
! BETA(K1) = B(K1)**(-1)
        beta(1,1,lev0,i) =  b(2,2,lev0,i)/y(1,lev0,i)
        beta(1,2,lev0,i) = -b(1,2,lev0,i)/y(1,lev0,i)
        beta(2,1,lev0,i) = -b(2,1,lev0,i)/y(1,lev0,i)
        beta(2,2,lev0,i) =  b(1,1,lev0,i)/y(1,lev0,i)
!
! Y(K1) = BETA(K1)*F(K1)
        y(1,lev0,i) = beta(1,1,lev0,i)*f(1,lev0,i)+
     |                beta(1,2,lev0,i)*f(2,lev0,i)
        y(2,lev0,i) = beta(2,1,lev0,i)*f(1,lev0,i)+
     |                beta(2,2,lev0,i)*f(2,lev0,i)
!
! Now deal with levels (K1+1),K2,1
        do k=lev0+1,lev1-1
!
! GAMMA(K-1) = BETA(K-1)*C(K-1)
          gamma(1,1,k-1,i) = beta(1,1,k-1,i)*c(1,1,k-1,i)+
     |                       beta(1,2,k-1,i)*c(2,1,k-1,i)          
          gamma(1,2,k-1,i) = beta(1,1,k-1,i)*c(1,2,k-1,i)+
     |                       beta(1,2,k-1,i)*c(2,2,k-1,i)          
          gamma(2,1,k-1,i) = beta(2,1,k-1,i)*c(1,1,k-1,i)+
     |                       beta(2,2,k-1,i)*c(2,1,k-1,i)          
          gamma(2,2,k-1,i) = beta(2,1,k-1,i)*c(1,2,k-1,i)+
     |                       beta(2,2,k-1,i)*c(2,2,k-1,i)          
!
! GAMMA(K) = B(K) - A(K)*GAMMA(K-1)
          gamma(1,1,k,i) = b(1,1,k,i)-a(1,1,k,i)*gamma(1,1,k-1,i)-
     |                                a(1,2,k,i)*gamma(2,1,k-1,i)
          gamma(1,2,k,i) = b(1,2,k,i)-a(1,1,k,i)*gamma(1,2,k-1,i)-
     |                                a(1,2,k,i)*gamma(2,2,k-1,i)
          gamma(2,1,k,i) = b(2,1,k,i)-a(2,1,k,i)*gamma(1,1,k-1,i)-
     |                                a(2,2,k,i)*gamma(2,1,k-1,i)
          gamma(2,2,k,i) = b(2,2,k,i)-a(2,1,k,i)*gamma(1,2,k-1,i)-
     |                                a(2,2,k,i)*gamma(2,2,k-1,i)
!
! Y(1,k,i) = determinant(GAMMA(K))
          y(1,k,i) = gamma(1,1,k,i)*gamma(2,2,k,i)-
     |               gamma(1,2,k,i)*gamma(2,1,k,i)
!
! BETA(K) = GAMMA(K)**(-1)
          beta(1,1,k,i) =  gamma(2,2,k,i)/y(1,k,i)
          beta(1,2,k,i) = -gamma(1,2,k,i)/y(1,k,i)
          beta(2,1,k,i) = -gamma(2,1,k,i)/y(1,k,i)
          beta(2,2,k,i) =  gamma(1,1,k,i)/y(1,k,i)
!
!  X(K) = F(K) - A(K)*Y(K-1)
          x(1,k,i) = f(1,k,i)-a(1,1,k,i)*y(1,k-1,i)-
     |                        a(1,2,k,i)*y(2,k-1,i)
          x(2,k,i) = f(2,k,i)-a(2,1,k,i)*y(1,k-1,i)-
     |                        a(2,2,k,i)*y(2,k-1,i)
!
! Y(K) = BETA(K)*X(K)
          y(1,k,i) = beta(1,1,k,i)*x(1,k,i)+beta(1,2,k,i)*x(2,k,i)
          y(2,k,i) = beta(2,1,k,i)*x(1,k,i)+beta(2,2,k,i)*x(2,k,i)
        enddo ! k=lev0+1,lev1-1
      enddo ! i=lonbeg,lonend
!
      do i=lonbeg,lonend
!
!  X(K2) = Y(K2)
        x(1,lev1-1,i) = y(1,lev1-1,i)
        x(2,lev1-1,i) = y(2,lev1-1,i)
!
! Backward sweep to determine final solution, X(K) for k=k2,k1,-1
        do k=lev1-2,lev0,-1
          x(1,k,i) = y(1,k,i)-gamma(1,1,k,i)*x(1,k+1,i)-
     |                        gamma(1,2,k,i)*x(2,k+1,i)
          x(2,k,i) = y(2,k,i)-gamma(2,1,k,i)*x(1,k+1,i)-
     |                        gamma(2,2,k,i)*x(2,k+1,i)
        enddo ! k=lev1-2,lev1-1,-1
      enddo ! i=lon0,lon1
      end subroutine blktri
!-----------------------------------------------------------------------
      subroutine filter_uv(fout,lev0,lev1,lon0,lon1,lat0,lat1,kut,name)
!
! Filter updated W omega:
!
      use params_module,only: nlat,nlonp4,nlon
      use filter_module,only: filter
      use addfld_module,only: addfld
#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)
      character(len=*),intent(in) :: name
!
! VT vampir tracing:
!
#ifdef VT
#include <VT.inc>
#endif
!
! Local:
      integer :: i,j,nlevs
      real :: fik(nlonp4,lev0:lev1),fkij(lev0:lev1,nlonp4,lat0:lat1)
      real :: fmin,fmax
!
#ifdef VT
!     code = 131 ; state = 'filter_uv' ; activity='Filtering'
      call vtbegin(131,ier)
#endif
!
      nlevs = lev1-lev0+1
!
! Save subdomain of current task in fkij:
      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,name)
#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,name)

!     write(6,"('filter uv after mp_scatterlons: lon0,1=',2i4,' lat0,1='
!    |  ,2i4,' fkij(:,lon0:lon1,:) min,max=',2e12.4)") lon0,lon1,
!    |  lat0,lat1,minval(fkij(:,lon0:lon1,:)),
!    |            maxval(fkij(:,lon0:lon1,:))
#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_uv' ; activity='Filtering'
      call vtend(131,ier)
#endif
      end subroutine filter_uv