subroutine horizontal_deln( i_diff_order, nx,ny,nz, 
     &                                gamma1, gammaw, gammat, rns, 
     &                                fu, fv, fw, ft, u, v, w, t ) 
!
! Calculates del^(i_diff_order) of input fields and adds to input tendencies
! (after mult. by appropriate scalar factors).
!
! Operates on a single level.
!
! Assumes dx = dy.
!
! Fields (u,v,w,t) should be lagged from current time level.
!
! INPUTS: 
!	i_diff_order = order of diffusion (2 or 4)	
!	gamma1       = nondiml. diffusion coeff. for u,v
!	gammaw, gammat
!	             = same for w,t
!       rns	     = 1 / ( number of small steps )
! 	fu, fv, fw, ft
!	             = partial tendencies for u,v,w,t
! 	u, v, w, t   = arrays for fields
! OUTPUTS:
! 	fu, fv, fw, ft
!	             = as above, but diffusive tendency added
!
      implicit none
      integer i_diff_order, nx, ny, nz, i,j
      real gamma1, gammaw, gammat, rns,
     & fu(0:nx,ny), fv(nx,0:ny), fw(nx,ny), ft(nx,ny),
     & u (0:nx,ny), v (nx,0:ny), w (nx,ny), t (nx,ny),
     & tmp(-1:nx+1, -1:ny+1)

      tmp = 0

      if ( i_diff_order.eq.2 ) then
	!-- del^2 diffusion

        !-- u ... fill tmp array with periodic rows, cols (don't
        !         need periodic corners); then apply del^2
        tmp(1:nx,1:ny) = u(1:nx,1:ny)
        tmp(0,1:ny) = u(nx-1,1:ny); tmp(1:nx,0) = u(1:nx,ny-1); 

        do i = 1,nx-1; do j = 1,ny-1
          fu(i,j) =  fu(i,j) + rns * gamma1 * (
     &  tmp(i+1,j) + tmp(i,j+1) + tmp(i-1,j) + tmp(i,j-1)
     & -4.*tmp(i,j)                            )
        enddo; enddo

        fu(nx,:) = fu(1,:); fu(:,ny) = fu(:,1); fu(nx,ny) = fu(1,1);
        fu(0,:) = fu(nx-1,:)

        !-- v 
        tmp(1:nx,1:ny) = v(1:nx,1:ny)
        tmp(0,1:ny) = v(nx-1,1:ny); tmp(1:nx,0) = v(1:nx,ny-1); 

        do i = 1,nx-1; do j = 1,ny-1
          fv(i,j) = fv(i,j) + rns * gamma1 * (
     &  tmp(i+1,j) + tmp(i,j+1) + tmp(i-1,j) + tmp(i,j-1)
     & -4.*tmp(i,j)                          )
        enddo; enddo

        fv(nx,:) = fv(1,:); fv(:,ny) = fv(:,1); fv(nx,ny) = fv(1,1);
        fv(:,0) = fv(:,ny-1)

        !-- w
        tmp(1:nx,1:ny) = w(1:nx,1:ny)
        tmp(0,1:ny) = w(nx-1,1:ny); tmp(1:nx,0) = w(1:nx,ny-1); 

        do i = 1,nx-1; do j = 1,ny-1
          fw(i,j) = fw(i,j) + rns * gammaw * ( 
     &  tmp(i+1,j) + tmp(i,j+1) + tmp(i-1,j) + tmp(i,j-1)
     & -4.*tmp(i,j)                          )
        enddo; enddo

        fw(nx,:) = fw(1,:); fw(:,ny) = fw(:,1); fw(nx,ny) = fw(1,1);

        !-- t
        tmp(1:nx,1:ny) = t(1:nx,1:ny)
        tmp(0,1:ny) = t(nx-1,1:ny); tmp(1:nx,0) = t(1:nx,ny-1); 

        do i = 1,nx-1; do j = 1,ny-1
          ft(i,j) = ft(i,j) + gammat * ( 
     &  tmp(i+1,j) + tmp(i,j+1) + tmp(i-1,j) + tmp(i,j-1)
     & -4.*tmp(i,j)                    )
        enddo; enddo

        ft(nx,:) = ft(1,:); ft(:,ny) = ft(:,1); ft(nx,ny) = ft(1,1);

      else if ( i_diff_order.eq.4 ) then
	!-- del^4 diffusion

        !-- u ... fill tmp array with periodic rows, cols (don't need 
        !         periodic corners except i = j = 0); then apply del^4
        tmp(1:nx,1:ny) = u(1:nx,1:ny)
        tmp(-1,1:ny) = u(nx-2,1:ny); tmp(0,1:ny) = u(nx-1,1:ny); 
        tmp(nx+1,1:ny) = u(2,1:ny);
        tmp(1:nx,-1) = u(1:nx,ny-2); tmp(1:nx,0) = u(1:nx,ny-1); 
        tmp(1:nx,ny+1) = u(1:nx,2);
        tmp(0,0) = u(nx-1,ny-1)

        do i = 1,nx-1; do j = 1,ny-1
          fu(i,j) = fu(i,j) + rns * gamma1 * (
     &  tmp(i+2,j) + tmp(i,j+2) + tmp(i-2,j) + tmp(i,j-2)
     & +2. * (  tmp(i+1,j+1) + tmp(i-1,j+1) 
     &        + tmp(i-1,j-1) + tmp(i+1,j-1) )
     & -8. * (  tmp(i+1,j) + tmp(i,j+1) 
     &        + tmp(i-1,j) + tmp(i,j-1) )
     & +20.*tmp(i,j)                         )
        enddo; enddo

        fu(nx,:) = fu(1,:); fu(:,ny) = fu(:,1); fu(nx,ny) = fu(1,1);
        fu(0,:) = fu(nx-1,:)

        !-- v 
        tmp(1:nx,1:ny) = v(1:nx,1:ny)
        tmp(-1,1:ny) = v(nx-2,1:ny); tmp(0,1:ny) = v(nx-1,1:ny); 
        tmp(nx+1,1:ny) = v(2,1:ny);
        tmp(1:nx,-1) = v(1:nx,ny-2); tmp(1:nx,0) = v(1:nx,ny-1); 
        tmp(1:nx,ny+1) = v(1:nx,2);
        tmp(0,0) = v(nx-1,ny-1)

        do i = 1,nx-1; do j = 1,ny-1
          fv(i,j) = fv(i,j) + rns * gamma1 * (
     &  tmp(i+2,j) + tmp(i,j+2) + tmp(i-2,j) + tmp(i,j-2)
     & +2. * (  tmp(i+1,j+1) + tmp(i-1,j+1) 
     &        + tmp(i-1,j-1) + tmp(i+1,j-1) )
     & -8. * (  tmp(i+1,j) + tmp(i,j+1) 
     &        + tmp(i-1,j) + tmp(i,j-1) )
     & +20.*tmp(i,j)                         )
        enddo; enddo

        fv(nx,:) = fv(1,:); fv(:,ny) = fv(:,1); fv(nx,ny) = fv(1,1);
        fv(:,0) = fv(:,ny-1)

        !-- w 
        tmp(1:nx,1:ny) = w(1:nx,1:ny)
        tmp(-1,1:ny) = w(nx-2,1:ny); tmp(0,1:ny) = w(nx-1,1:ny); 
        tmp(nx+1,1:ny) = w(2,1:ny);
        tmp(1:nx,-1) = w(1:nx,ny-2); tmp(1:nx,0) = w(1:nx,ny-1); 
        tmp(1:nx,ny+1) = w(1:nx,2);
        tmp(0,0) = w(nx-1,ny-1)

        do i = 1,nx-1; do j = 1,ny-1
          fw(i,j) = fw(i,j) + rns * gammaw * (
     &  tmp(i+2,j) + tmp(i,j+2) + tmp(i-2,j) + tmp(i,j-2)
     & +2. * (  tmp(i+1,j+1) + tmp(i-1,j+1) 
     &        + tmp(i-1,j-1) + tmp(i+1,j-1) )
     & -8. * (  tmp(i+1,j) + tmp(i,j+1) 
     &        + tmp(i-1,j) + tmp(i,j-1) )
     & +20.*tmp(i,j)                         )
        enddo; enddo

        fw(nx,:) = fw(1,:); fw(:,ny) = fw(:,1); fw(nx,ny) = fw(1,1);

        !-- t 
        tmp(1:nx,1:ny) = t(1:nx,1:ny)
        tmp(-1,1:ny) = t(nx-2,1:ny); tmp(0,1:ny) = t(nx-1,1:ny); 
        tmp(nx+1,1:ny) = t(2,1:ny);
        tmp(1:nx,-1) = t(1:nx,ny-2); tmp(1:nx,0) = t(1:nx,ny-1); 
        tmp(1:nx,ny+1) = t(1:nx,2);
        tmp(0,0) = t(nx-1,ny-1)

        do i = 1,nx-1; do j = 1,ny-1
          ft(i,j) = ft(i,j) + gammat * (
     &  tmp(i+2,j) + tmp(i,j+2) + tmp(i-2,j) + tmp(i,j-2)
     & +2. * (  tmp(i+1,j+1) + tmp(i-1,j+1) 
     &        + tmp(i-1,j-1) + tmp(i+1,j-1) )
     & -8. * (  tmp(i+1,j) + tmp(i,j+1) 
     &        + tmp(i-1,j) + tmp(i,j-1) )
     & +20.*tmp(i,j)                   )
        enddo; enddo

        ft(nx,:) = ft(1,:); ft(:,ny) = ft(:,1); ft(nx,ny) = ft(1,1);

      else if ( i_diff_order.eq.5 ) then
	!-- del^4 in x, y separately 

        !-- u ... fill tmp array with periodic rows, cols (don't need 
        !         periodic corners except i = j = 0); then apply del^4
        tmp(1:nx,1:ny) = u(1:nx,1:ny)
        tmp(-1,1:ny) = u(nx-2,1:ny); tmp(0,1:ny) = u(nx-1,1:ny); 
        tmp(nx+1,1:ny) = u(2,1:ny);
        tmp(1:nx,-1) = u(1:nx,ny-2); tmp(1:nx,0) = u(1:nx,ny-1); 
        tmp(1:nx,ny+1) = u(1:nx,2);
        tmp(0,0) = u(nx-1,ny-1)

        do i = 1,nx-1; do j = 1,ny-1
          fu(i,j) = fu(i,j) + rns * gamma1 * (
     &  tmp(i+2,j) + tmp(i,j+2) + tmp(i-2,j) + tmp(i,j-2)
     & -4. * (  tmp(i+1,j) + tmp(i,j+1) 
     &        + tmp(i-1,j) + tmp(i,j-1) )
     & +12.*tmp(i,j)                         )
        enddo; enddo

        fu(nx,:) = fu(1,:); fu(:,ny) = fu(:,1); fu(nx,ny) = fu(1,1);
        fu(0,:) = fu(nx-1,:)

        !-- v 
        tmp(1:nx,1:ny) = v(1:nx,1:ny)
        tmp(-1,1:ny) = v(nx-2,1:ny); tmp(0,1:ny) = v(nx-1,1:ny); 
        tmp(nx+1,1:ny) = v(2,1:ny);
        tmp(1:nx,-1) = v(1:nx,ny-2); tmp(1:nx,0) = v(1:nx,ny-1); 
        tmp(1:nx,ny+1) = v(1:nx,2);
        tmp(0,0) = v(nx-1,ny-1)

        do i = 1,nx-1; do j = 1,ny-1
          fv(i,j) = fv(i,j) + rns * gamma1 * (
     &  tmp(i+2,j) + tmp(i,j+2) + tmp(i-2,j) + tmp(i,j-2)
     & -4. * (  tmp(i+1,j) + tmp(i,j+1) 
     &        + tmp(i-1,j) + tmp(i,j-1) )
     & +12.*tmp(i,j)                         )
        enddo; enddo

        fv(nx,:) = fv(1,:); fv(:,ny) = fv(:,1); fv(nx,ny) = fv(1,1);
        fv(:,0) = fv(:,ny-1)

        !-- w 
        tmp(1:nx,1:ny) = w(1:nx,1:ny)
        tmp(-1,1:ny) = w(nx-2,1:ny); tmp(0,1:ny) = w(nx-1,1:ny); 
        tmp(nx+1,1:ny) = w(2,1:ny);
        tmp(1:nx,-1) = w(1:nx,ny-2); tmp(1:nx,0) = w(1:nx,ny-1); 
        tmp(1:nx,ny+1) = w(1:nx,2);
        tmp(0,0) = w(nx-1,ny-1)

        do i = 1,nx-1; do j = 1,ny-1
          fw(i,j) = fw(i,j) + rns * gammaw * (
     &  tmp(i+2,j) + tmp(i,j+2) + tmp(i-2,j) + tmp(i,j-2)
     & -4. * (  tmp(i+1,j) + tmp(i,j+1) 
     &        + tmp(i-1,j) + tmp(i,j-1) )
     & +12.*tmp(i,j)                         )
        enddo; enddo

        fw(nx,:) = fw(1,:); fw(:,ny) = fw(:,1); fw(nx,ny) = fw(1,1);

        !-- t 
        tmp(1:nx,1:ny) = t(1:nx,1:ny)
        tmp(-1,1:ny) = t(nx-2,1:ny); tmp(0,1:ny) = t(nx-1,1:ny); 
        tmp(nx+1,1:ny) = t(2,1:ny);
        tmp(1:nx,-1) = t(1:nx,ny-2); tmp(1:nx,0) = t(1:nx,ny-1); 
        tmp(1:nx,ny+1) = t(1:nx,2);
        tmp(0,0) = t(nx-1,ny-1)

        do i = 1,nx-1; do j = 1,ny-1
          ft(i,j) = ft(i,j) + gammat * (
     &  tmp(i+2,j) + tmp(i,j+2) + tmp(i-2,j) + tmp(i,j-2)
     & -4. * (  tmp(i+1,j) + tmp(i,j+1) 
     &        + tmp(i-1,j) + tmp(i,j-1) )
     & +12.*tmp(i,j)                   )
        enddo; enddo

        ft(nx,:) = ft(1,:); ft(:,ny) = ft(:,1); ft(nx,ny) = ft(1,1);

      else 
        write(6,*) 'CHECK I_DIFF_ORDER.  SHOULD BE EITHER 2, 4 OR 5 !!!!'
        stop

      endif

      return
      end