
      module tuv_bndry
      use addfld_module,only: addfld
      implicit none
!
! Output: Lower boundary conditions for t,u,v:
!
      real,allocatable,dimension(:,:) ::  ! (lon0:lon1,lat0:lat1)
     |  tbnd,  ! tn lower boundary (not in use as of Mar, 2004)
     |  ubnd,  ! un lower boundary (used in duv)
     |  vbnd   ! vn lower boundary (used in duv)
      contains
!-----------------------------------------------------------------------
      subroutine tuvbnd(tn,tn_nm,un,un_nm,vn,vn_nm,w,z,barm,cp,
     |  hdt,hdu,hdv,gwt,gwu,gwv,cool_imp,cool_exp,
     |  lev0,lev1,lon0,lon1,lat0,lat1)
!
! Calculate lower boundaries of t,u,v.
!
      use params_module,only: dz,nlonp4,nlat
      use cons_module,only: shapiro,dtx2inv,expz,gask,racs,grav,re_inv,
     |  dtr,cor,tanphi=>tn,kut
      use init_module,only: glat
      use qrj_module,only: qtotal
      use mgw_module,only: rayk
      use input_module,only: ncep_ncfile
      use fields_module,only: tlbc,ulbc,vlbc,tlbc_nm,ulbc_nm,vlbc_nm
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in)::
     |  tn,       ! temperature                   (f(:,nj+nt  ))
     |  tn_nm,    ! temperature, time n-1         (f(:,nj+ntnm))
     |  un,       ! zonal wind                    (f(:,nj+nu))
     |  un_nm,    ! zonal wind, time n-1          (f(:,nj+nunm))
     |  vn,       ! meridional wind               (f(:,nj+nv))
     |  vn_nm,    ! meridional wind, time n-1     (f(:,nj+nvnm))
     |  w,        ! omega (for vertical velocity) (f(:,nj+nw))
     |  z,        ! geopotential height           (f(:,nj+nz))
     |  barm,     ! mean mol weight               (f(:,nj+nms))
     |  cp,       ! specific heat                 (f(:,ncp))
     |  hdt,      ! horizontal diffusion of T     (f(:,nqdh))
     |  hdu,      ! horizontal diffusion of U     (f(:,nflh))
     |  hdv,      ! horizontal diffusion of V     (f(:,nfph))
     |  gwt,      ! gravity wave heating (mgw.F)  (f(:,ngwvt))
     |  gwu,      ! gravity wave U tendency       (f(:,ngwvu))
     |  gwv,      ! gravity wave V tendency       (f(:,ngwvv))
     |  cool_exp, ! explicit cooling (radcool.F)  (f(:,nwte))
     |  cool_imp  ! implicit cooling (radcool.F)  (f(:,nwti))
!
! Local:
      integer :: k,i,j,lat,lonbeg,lonend,kutuv(nlat)
      integer ::i0,i1,nk,nkm1,ni ! for addfld
      real,dimension(lon0:lon1) :: 
     |  advect,    ! horizontal advection (t)
     |  advecu,    ! horizontal advection (u)
     |  advecv,    ! horizontal advection (v)
     |  cpmbar,    ! cp*barm
     |  shapiro2,  ! t shapiro stage 2
     |  ushapiro2, ! u shapiro stage 2
     |  vshapiro2, ! v shapiro stage 2
     |  qterm1,    ! 1./(2*dt)*t(n-1)-advect + q/cp
     |  qterm2,    ! 1./(2*dt)+cool_imp+w*r/(cp*mbar)
     |  heatlbc,   ! total heating at lower boundary
     |  a11,a12,a21,a22, ! final terms for ubnd,vbnd
     |  deta       ! a11*a22-a12*a21
      real :: shapiro1(lon0-2:lon1+2)   ! t shapiro stage 1
      real :: ushapiro1(lon0-2:lon1+2)  ! u shapiro stage 1
      real :: vshapiro1(lon0-2:lon1+2)  ! v shapiro stage 1
      real,dimension(1,lon0:lon1) :: 
     |  zl,zp             ! derivatives of z at bottom boundary
      real :: rlat,eqfric ! for rayleigh friction
      real,dimension(lev0:lev1,lon0:lon1) :: ! diag only (addfld calls)
     |  advect_diag,advecu_diag,advecv_diag,heatlbc_diag,qterm1_diag, 
     |  qterm2_diag,smooth1,smooth2,tbnd_diag,zl_diag,zp_diag,
     |  usmooth1,vsmooth1,usmooth2,vsmooth2,a11_diag,a12_diag,
     |  a21_diag,a22_diag,ubnd_diag,vbnd_diag
      real,parameter :: hlat = 20./57.295 ! for rayleigh friction
!
! For addfld calls:
      i0 = lon0 ; i1 = lon1 ; nk = lev1-lev0+1 ; nkm1 = nk-1
      ni = i1-i0+1
!
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = lon1-2
!
! Latitude scan:
      do lat=lat0,lat1
!
! Equatorial rayleigh friction (tunable) (s.a. hlat above):
        rlat = glat(lat)*dtr
!       eqfric = 5.E-7*exp(-(rlat/hlat)**2)
	eqfric = 1.e-12
!
! 5/4/06 btf: new call to advecl for lbc:
        call adveclbc(tlbc,ulbc,vlbc,advect,lev0,lev1,lon0,lon1,
     |    lat0,lat1,lat)

!       call xferdiag(advect,ni,advect_diag,nk,ni)
!       call addfld('TADVEC  ',' ',' ',advect_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Add in vertical advection
        do i=lon0,lon1
          advect(i) = advect(i)+2.*w(lev0,i,lat)*
     |      (tn(lev0,i,lat)-tlbc(i,lat))/dz
          cpmbar(i) = cp(lev0,i,lat)*barm(lev0,i,lat)
        enddo ! i=lon0,lon1
!
! Shapiro smoother, stage 1
        do i=lon0-2,lon1+2
          shapiro1(i) = tlbc_nm(i,lat)-shapiro*
     |         (tlbc_nm(i,lat+2)+tlbc_nm(i,lat-2)-
     |      4.*(tlbc_nm(i,lat+1)+tlbc_nm(i,lat-1))+
     |      6.*tlbc_nm(i,lat))
        enddo ! i=lon0-2,lon1+2
!
! Shapiro smoother, stage 2: 
        shapiro2 = 0. ! whole array
        do i=lonbeg,lonend
          shapiro2(i)= shapiro1(i)-shapiro*
     |         (shapiro1(i+2)+shapiro1(i-2)-
     |      4.*(shapiro1(i+1)+shapiro1(i-1))+6.*shapiro1(i))
        enddo ! i=lonbeg,lonend
!
! Save single level diags to secondary history:
!       smooth1 = 0. ; smooth2 = 0.
!       call xferdiag(shapiro1(lon0),ni,smooth1,nk,ni)
!       call xferdiag(shapiro2(lon0),ni,smooth2,nk,ni)
!       call addfld('SMOOTH1 ',' ',' ',smooth1,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('SMOOTH2 ',' ',' ',smooth2,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)

! tgcm24 adds periodic points to shapiro2 (s3) at this point.

!
! Save input fields to secondary history:
! (3d fields, but here only bottom 2 levels lev0 and lev0+1 are used)
!   gwt: gravity wave heating from gw_drag
!   hdt: horizontal diffusiong from hdif3
!   coolexp: explicit cooling from radcool
!
! For comparison with tgcm24, gravity wave tendencies must be turned
!   off in both models. timegcm1 has only new parameterization from
!   waccm/cam, and tgcm24 has only Hanli's old parameterization.
!   See calc_gw logical flag in mgw.F to turn off gwt,u,v.
!
!       call addfld('GWT     ',' ',' ',gwt(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('HDT     ',' ',' ',hdt(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('COOLEXP ',' ',' ',cool_exp(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)

        do i=lon0,lon1
          qterm1(i) = dtx2inv*shapiro2(i)-advect(i)
          heatlbc(i) = qtotal(lev0,i,lat)
!
! Add heating from fourth order diffusion and gravity waves:
          heatlbc(i) = heatlbc(i)+1.5*(hdt(lev0,i,lat)+gwt(lev0,i,lat))-
     |      .5*(hdt(lev0+1,i,lat)+gwt(lev0+1,i,lat))
!
! Add explicit newtonian cooling:
          heatlbc(i) = heatlbc(i)-(1.5*cool_exp(lev0,i,lat)/expz(lev0)-
     |      .5*cool_exp(lev0+1,i,lat)/expz(lev0+1))
!
! qterm1 = qterm1 + q/cp
          qterm1(i) = qterm1(i)+heatlbc(i)/cp(lev0,i,lat)
        enddo ! i=lon0,lon1
!
! Save single level heating and qterm1 to secondary history:
!       call xferdiag(heatlbc,ni,heatlbc_diag,nk,ni)
!       call xferdiag(qterm1 ,ni,qterm1_diag ,nk,ni)
!       call addfld('HEATLBC ',' ',' ',heatlbc_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('QTERM1  ',' ',' ',qterm1_diag ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! qterm2 = W*R/(CP*MBAR) (adiabatic heating)
        do i=lon0,lon1
          qterm2(i) = w(lev0,i,lat)*gask/cpmbar(i)
!
! qterm2 = qterm2 + 1./(2.*dt)+cool_imp
          qterm2(i) = dtx2inv+1.5*cool_imp(lev0,i,lat)-
     |      0.5*cool_imp(lev0+1,i,lat)+qterm2(i)
!
! Final tn bottom boundary:
          if (len_trim(ncep_ncfile) <= 0) 
     |      tbnd(i,lat) = qterm1(i)/qterm2(i)
        enddo ! i=lon0,lon1
!
! Save single level qterm2 and tbnd to secondary history:
!       call xferdiag(qterm2        ,ni,qterm2_diag,nk,ni)
!       call xferdiag(tbnd(lon0,lat),ni,tbnd_diag  ,nk,ni)
!       call addfld('QTERM2  ',' ',' ',qterm2_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('TBND    ',' ',' ',tbnd_diag  ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! Derive boundary conditions for U,V.
!
! 5/4/06 btf: New advecl calls for lbc:
        call adveclbc(ulbc,ulbc,vlbc,advecu,lev0,lev1,lon0,lon1,
     |    lat0,lat1,lat)
        call adveclbc(vlbc,ulbc,vlbc,advecv,lev0,lev1,lon0,lon1,
     |    lat0,lat1,lat)
!
! Save to secondary history:
!       call xferdiag(advecu,ni,advecu_diag,nk,ni)
!       call xferdiag(advecv,ni,advecv_diag,nk,ni)
!       call addfld('UADVEC  ',' ',' ',advecu_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('VADVEC  ',' ',' ',advecv_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Longitudinal and latitudinal derivatives of geopotential bottom boundary:
        call dldp(z(lev0:lev0,:,lat-2:lat+2),zl,zp,lev0,lev0,lon0,lon1,
     |    lat)
!       write(6,"(/,'tuvbnd after dldp: lat=',i3,' zl=',/,(6e12.4))")
!    |    lat,zl
!       write(6,"(' zp=',/,(6e12.4))") zp
!
! Save to secondary history:
!       call xferdiag(zl,ni,zl_diag,nk,ni)
!       call xferdiag(zp,ni,zp_diag,nk,ni)
!       call addfld('ZL_DLDP ',' ',' ',zl_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('ZP_DLDP ',' ',' ',zp_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! zl = G/(RAD*COS(PHI))*DZ/DLAMDA
! zp = G/RA*DZ/DPHI
! (formerly sub glpl in tgcm24)
!
        do i=lonbeg,lonend
          zl(1,i) = zl(1,i)*grav*racs(lat)
          zp(1,i) = zp(1,i)*grav*re_inv
        enddo ! i=lonbeg,lonend
!
! 5/19/04 btf: May need to do periodic points on zl,zp here. 
! Try setting to zero for now to avoid FP exception:
!
        if (lon0==1) then
          zl(1,1:2) = 0.
          zp(1,1:2) = 0.
        endif
        if (lon1==nlonp4) then
          zl(1,nlonp4-1:nlonp4) = 0.
          zp(1,nlonp4-1:nlonp4) = 0.
        endif
!       write(6,"(/,'tuvbnd after glpl calc: lat=',i3,' racs(lat)=',
!    |    e12.4,' zl=',/,(6e12.4))") lat,racs(lat),zl
!       write(6,"(' zp=',/,(6e12.4))") zp
!
! Save to secondary history:
!       call xferdiag(zl,ni,zl_diag,nk,ni)
!       call xferdiag(zp,ni,zp_diag,nk,ni)
!       call addfld('ZL_LBC  ',' ',' ',zl_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('ZP_LBC  ',' ',' ',zp_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Add in vertical advection:
        do i=lon0,lon1
          advecu(i) = advecu(i)+2.*w(lev0,i,lat)*
     |      (un(lev0,i,lat)-ulbc(i,lat))/dz
          advecv(i) = advecv(i)+2.*w(lev0,i,lat)*
     |      (vn(lev0,i,lat)-vlbc(i,lat))/dz
!
! Add horizontal derivatives of geopotential:
          advecu(i) = advecu(i)+zl(1,i)
          advecv(i) = advecv(i)+zp(1,i)
!
! Fourth order horizontal diffusion:
!   hdu,v: zonal,meridional diffusion from hdif3 (hdif.F)
!   gwu,v: gravity wave tendencies from mgw (mgw.F).
!          (gw may be turned off, check calc_gw)
!
!         advecu(i) = advecu(i)-
!    |      1.5*(hdu(lev0  ,i,lat)+gwu(lev0,  i,lat))+
!    |      0.5*(hdu(lev0+1,i,lat)+gwu(lev0+1,i,lat))
!         advecv(i) = advecv(i)-
!    |      1.5*(hdv(lev0  ,i,lat)+gwv(lev0,  i,lat))+
!    |      0.5*(hdv(lev0+1,i,lat)+gwv(lev0+1,i,lat))
!
          advecu(i) = advecu(i)-
     |      1.5*(hdu(lev0  ,i,lat))+
     |      0.5*(hdu(lev0+1,i,lat))
          advecv(i) = advecv(i)-
     |      1.5*(hdv(lev0  ,i,lat))+
     |      0.5*(hdv(lev0+1,i,lat))
        enddo ! i=lon0,lon1
!
! Save to secondary history:
!       call xferdiag(advecu,ni,advecu_diag,nk,ni)
!       call xferdiag(advecv,ni,advecv_diag,nk,ni)
!       call addfld('UADVEC1 ',' ',' ',advecu_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('VADVEC1 ',' ',' ',advecv_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Shapiro smoother for u,v, stage 1
        do i=lon0-2,lon1+2
          ushapiro1(i) = ulbc_nm(i,lat)-shapiro*
     |         (ulbc_nm(i,lat+2)+ulbc_nm(i,lat-2)-
     |      4.*(ulbc_nm(i,lat+1)+ulbc_nm(i,lat-1))+
     |      6.* ulbc_nm(i,lat))
          vshapiro1(i) = vlbc_nm(i,lat)-shapiro*
     |         (vlbc_nm(i,lat+2)+vlbc_nm(i,lat-2)-
     |      4.*(vlbc_nm(i,lat+1)+vlbc_nm(i,lat-1))+
     |      6.* vlbc_nm(i,lat))
        enddo ! i=lon0-2,lon1+2

!       usmooth1 = 0. ; vsmooth1 = 0.
!       call xferdiag(ushapiro1(lon0),ni,usmooth1,nk,ni)
!       call xferdiag(vshapiro1(lon0),ni,vsmooth1,nk,ni)
!       call addfld('USMOOTH1',' ',' ',usmooth1,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('VSMOOTH1',' ',' ',vsmooth1,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Shapiro smoother for u,v, stage 2:
        ushapiro2 = 0. ! whole array
        vshapiro2 = 0. ! whole array
        do i=lonbeg,lonend
          ushapiro2(i)= ushapiro1(i)-shapiro*
     |         (ushapiro1(i+2)+ushapiro1(i-2)-
     |      4.*(ushapiro1(i+1)+ushapiro1(i-1))+6.*ushapiro1(i))
          vshapiro2(i)= vshapiro1(i)-shapiro*
     |         (vshapiro1(i+2)+vshapiro1(i-2)-
     |      4.*(vshapiro1(i+1)+vshapiro1(i-1))+6.*vshapiro1(i))
        enddo ! i=lonbeg,lonend

!       usmooth2 = 0. ; vsmooth2 = 0.
!       call xferdiag(ushapiro2(lon0),ni,usmooth2,nk,ni)
!       call xferdiag(vshapiro2(lon0),ni,vsmooth2,nk,ni)
!       call addfld('USMOOTH2',' ',' ',usmooth2,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('VSMOOTH2',' ',' ',vsmooth2,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)

! tgcm24 adds periodic points to ushapiro2,vspairo2 at this point.

!
! advecu = U(N-1)/(2.*DT)-advecu
! advecv = V(N-1)/(2.*DT)-advecv
!
        do i=lon0,lon1
          advecu(i) = ushapiro2(i)*dtx2inv-advecu(i)             ! s1
          advecv(i) = vshapiro2(i)*dtx2inv-advecv(i)             ! s2
        enddo ! i=lon0,lon1

!       call xferdiag(advecu,ni,advecu_diag,nk,ni)
!       call xferdiag(advecv,ni,advecv_diag,nk,ni)
!       call addfld('UADVEC2 ',' ',' ',advecu_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('VADVEC2 ',' ',' ',advecv_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Final terms for ubnd,vbnd
        do i=lon0,lon1
          a11(i) = dtx2inv+(rayk(lev0)+eqfric)**1.5/             ! s3
     |                     (rayk(lev0+1)+eqfric)**0.5 
!         a12(i) = -(cor(lat)+un(lev1,i,lat)*tanphi(lat)*re_inv) ! s4
          a12(i) = -(cor(lat)+ulbc(i,lat)*tanphi(lat)*re_inv)    ! s4
          a21(i) = -a12(i)                                       ! s5
          a22(i) = dtx2inv+(rayk(lev0)+eqfric)**1.5/             ! s6
     |                     (rayk(lev0+1)+eqfric)**0.5 
!
! deta = DET(A) = A11*A22-A12*A21
          deta(i) = a11(i)*a22(i)-a12(i)*a21(i)                  ! s7
!
! Solve for ubnd,vbnd:
          ubnd(i,lat)= ( a22(i)*advecu(i)-a12(i)*advecv(i))/deta(i)
          vbnd(i,lat)= (-a21(i)*advecu(i)+a11(i)*advecv(i))/deta(i)
        enddo ! i=lon0,lon1

!       call xferdiag(a11,ni,a11_diag,nk,ni)
!       call xferdiag(a12,ni,a12_diag,nk,ni)
!       call xferdiag(a21,ni,a21_diag,nk,ni)
!       call xferdiag(a22,ni,a22_diag,nk,ni)
!       call addfld('A11     ',' ',' ',a11_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('A12     ',' ',' ',a12_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('A21     ',' ',' ',a21_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('A22     ',' ',' ',a22_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Save u,v lower boundaries to secondary history:
!       call xferdiag(ubnd(lon0,lat),ni,ubnd_diag,nk,ni)
!       call xferdiag(vbnd(lon0,lat),ni,vbnd_diag,nk,ni)
!       call addfld('UBND    ',' ',' ',ubnd_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('VBND    ',' ',' ',vbnd_diag,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! End latitude scan:
      enddo ! lat=lat0,lat1

!
! Save 2d pre-filtered u,v lower boundaries to secondary history:
!     call addfld('UBND',' ',' ',ubnd,
!    |  'lon',lon0,lon1,'lat',lat0,lat1,0)
!     call addfld('VBND',' ',' ',vbnd,
!    |  'lon',lon0,lon1,'lat',lat0,lat1,0)
! 
! 8/14/06 btf: filter (fourier smooth) u,v boundaries:
! 10/18/06 btf: add filter to timegcm1.3, as in tiegcmdr10
!
! kut(nlat) is use-associated from cons.F. If different
! values are desired, change local kutuv here:
!
      do j=1,nlat
        kutuv(j) = kut(j)
      enddo
 
      call filter_uvbnd(ubnd,lev0,lev1,lon0,lon1,lat0,lat1,kutuv)
      call filter_uvbnd(vbnd,lev0,lev1,lon0,lon1,lat0,lat1,kutuv)
!
! Save 2d filtered u,v lower boundaries to secondary history:
!     call addfld('UBND_FIL',' ',' ',ubnd,
!    |  'lon',lon0,lon1,'lat',lat0,lat1,0)
!     call addfld('VBND_FIL',' ',' ',vbnd,
!    |  'lon',lon0,lon1,'lat',lat0,lat1,0)

      end subroutine tuvbnd
!-----------------------------------------------------------------------
      subroutine dldp(f,xl,xp,lev0,lev1,lon0,lon1,lat)
!
! Derivatives of f are returned in xl,xp. Bottom boundary only.
!
      use cons_module,only: dlamda_2div3, dlamda_1div12,
     |                      dphi_2div3,   dphi_1div12
      use params_module,only: nlonp4
!
! Args:
      implicit none
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev0,lon0-2:lon1+2,lat-2:lat+2),
     |  intent(in) :: f
      real,dimension(lev0:lev0,lon0:lon1),intent(out) :: xl,xp
!
! Local:
      integer :: i,lonbeg,lonend
!
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = lon1-2
!
      do i=lonbeg,lonend
!
! Longitudinal derivative:
        xl(lev0,i) = dlamda_2div3 *(f(lev0,i+1,lat)-f(lev0,i-1,lat))-
     |            dlamda_1div12*(f(lev0,i+2,lat)-f(lev0,i-2,lat)) 
!
! Latitudinal derivative:
        xp(lev0,i) = dphi_2div3 *(f(lev0,i,lat+1)-f(lev0,i,lat-1))-
     |               dphi_1div12*(f(lev0,i,lat+2)-f(lev0,i,lat-2))

!       write(6,"('dldp: lat=',i3,' i=',i3,' xl=',/,(6e12.4))")
!    |    lat,i,xl(:,i)

      enddo ! i=lonbeg,lonend

      end subroutine dldp
!-----------------------------------------------------------------------
      subroutine xferdiag(fin,idin, fout,idout1,idout2)
!
! Transfer 1d array fin to 2d fout. First dimension of fout (usually zp)
!   will be redundant.
! This is used to transfer (lon) arrays to diagnostic (lev,lon) arrays 
!   for addfld.
!
      integer,intent(in) :: idin,idout1,idout2
      real,intent(in) :: fin(idin)
      real,intent(out) :: fout(idout1,idout2)
!
      integer :: k
      do k=1,idout1
        fout(k,:) = fin(:)
      enddo
      end subroutine xferdiag
!-----------------------------------------------------------------------
      subroutine alloc_tuvbnd(lon0,lon1,lat0,lat1)
!
! Args:
      integer,intent(in) :: lon0,lon1,lat0,lat1
!
! Local:
      integer :: istat
!
! Allocate lower boundaries for t,u,v: 
!
      allocate(tbnd(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"(/,'>>> alloc_tuvbnd: error allocating',
     |  ' tbnd: stat=',i3)") istat
!
      allocate(ubnd(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"(/,'>>> alloc_tuvbnd: error allocating',
     |  ' ubnd: stat=',i3)") istat
!
      allocate(vbnd(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"(/,'>>> alloc_tuvbnd: error allocating',
     |  ' vbnd: stat=',i3)") istat

      end subroutine alloc_tuvbnd
!-----------------------------------------------------------------------
      subroutine filter_uvbnd(fout,lev0,lev1,lon0,lon1,lat0,lat1,kut)
!
! Filter updated u,v boundaries:
! Maintain k-dimension in fik and fkij, even tho fout does not have
!   a k-dim, for the calls to gather, scatter, and filter. 
!
      use params_module,only: nlat,nlonp4,nlon,spval
      use filter_module,only: filter
      use addfld_module,only: addfld ! add fields to secondary histories
#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(lon0:lon1,lat0:lat1)
!
! 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_uvbnd' ; 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) ! k-dim in fkij is redundant 
        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(1,i,j) ! k-levels are redundant in fkij
        enddo
      enddo
!
#ifdef VT
!     code = 131 ; state = 'filter_uvbnd' ; activity='Filtering'
      call vtend(131,ier)
#endif
      end subroutine filter_uvbnd
!-----------------------------------------------------------------------
      end module tuv_bndry
