#include "dims.h"
!
      subroutine hdif3(lat,ixt)
      use cons_module,only: imax,imaxp2,len2,kmax,imaxp4,long
      implicit none
!
! Calculate horizontal diffusion terms (phys fields NFLH->NPDHNAT)
! in f-array, from kld terms in fg-array (kld terms were saved 
! globally in fg (also global nrh is in fnrh)) by hdif1 and hdif2).
!
#include "params.h"
#include "fcom.h"
#include "fgcom.h"
#include "vscr.h"
#include "index.h"
!
! Args:
      integer,intent(in) :: lat,ixt
!
! Local:
      integer :: i,k,n
      real :: fmin,fmax
!
! S2=-(T+T0)/M=-1./S4
      do k=1,kmax
        do i=1,imaxp4
          s2(i,k) = -1. / fnrh(i,k,lat-1)
        enddo
      enddo
!
! FLH (use kldu):
      call lsqdsq(lat,fg(2,nkldu+1,lat+1,ixt),fg(2,nkldu+1,lat,ixt),
     |                fg(2,nkldu+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,nflh) = s1(i+2,1)*s2(i+2,1)
      enddo
!     call addfsech('HDU',' ',' ',f(1,nflh),zimxp,zkmxp,zkmx,lat)
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,nflh+k) = f(i+imax,nflh+k)
          f(i+imaxp2,nflh+k) = f(i+2,nflh+k)
        enddo
      enddo
      f(:,nflh+kmax) = 0.
!     call addfsech('HDU',' ',' ',f(1,nflh),zimxp,zkmxp,zkmx,lat)
!
! FPH (use kldv):
      call lsqdsq(lat,fg(2,nkldv+1,lat+1,ixt),fg(2,nkldv+1,lat,ixt),
     |                fg(2,nkldv+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,nfph) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,nfph+k) = f(i+imax,nfph+k)
          f(i+imaxp2,nfph+k) = f(i+2,nfph+k)
        enddo
      enddo
      f(:,nfph+kmax) = 0.
!     call addfsech('HDV',' ',' ',f(1,nfph),zimxp,zkmxp,zkmx,lat)
!
! QDH (use kldt):
!
      call lsqdsq(lat,fg(2,nkldt+1,lat+1,ixt),fg(2,nkldt+1,lat,ixt),
     |                fg(2,nkldt+1,lat-1,ixt),s1(3,1),long)
      do i=1,len2
        s3(i,1) = (f(i,ncp)+f(i,ncp+1))*0.5
      enddo
      do i=1,long
        f(i+2,nqdh) = s1(i+2,1)*s2(i+2,1)*s3(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,nqdh+k) = f(i+imax,nqdh+k)
          f(i+imaxp2,nqdh+k) = f(i+2,nqdh+k)
        enddo
      enddo
      f(:,nqdh+kmax) = 0.
!     call addfsech('HDT',' ',' ',f(1,nqdh),zimxp,zkmxp,zkmx,lat)
!
! PSDH (use kldps):
      call lsqdsq(lat,fg(2,nkldps+1,lat+1,ixt),fg(2,nkldps+1,lat,ixt),
     |                fg(2,nkldps+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npsdh) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npsdh+k) = f(i+imax,npsdh+k)
          f(i+imaxp2,npsdh+k) = f(i+2,npsdh+k)
        enddo
      enddo
      f(:,npsdh+kmax) = 0.
!     call addfsech('HDO2',' ',' ',f(1,npsdh),zimxp,zkmxp,zkmx,lat)
!
! PSDH2 (use kldp2):
      call lsqdsq(lat,fg(2,nkldp2+1,lat+1,ixt),fg(2,nkldp2+1,lat,ixt),
     |                fg(2,nkldp2+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npsdh2) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npsdh2+k) = f(i+imax,npsdh2+k)
          f(i+imaxp2,npsdh2+k) = f(i+2,npsdh2+k)
        enddo
      enddo
      f(:,npsdh2+kmax) = 0.
!     call addfsech('HDO1',' ',' ',f(1,npsdh2),zimxp,zkmxp,zkmx,lat)
      return
      end