#include "dims.h"
!
      subroutine hdif3(lat,ixt)
      use cons_module,only: len2,kmax,imax,imaxp2,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
      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.
!
! 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.
!
! 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.
!
! 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.
!
! 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.
!
! The following 11 minor species are for time-gcm only (not tiegcm):
!
! N4S (use kldn4s):
      call lsqdsq(lat,fg(2,nkldn4s+1,lat+1,ixt),
     |                fg(2,nkldn4s+1,lat,  ixt),
     |                fg(2,nkldn4s+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npdhn4s) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npdhn4s+k) = f(i+imax,npdhn4s+k)
          f(i+imaxp2,npdhn4s+k) = f(i+2,npdhn4s+k)
        enddo
      enddo
      f(:,npdhn4s+kmax) = 0.
!
! NOZ (use kldnoz):
      call lsqdsq(lat,fg(2,nkldnoz+1,lat+1,ixt),fg(2,nkldnoz+1,lat,ixt),
     |                fg(2,nkldnoz+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npdhnoz) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npdhnoz+k) = f(i+imax,npdhnoz+k)
          f(i+imaxp2,npdhnoz+k) = f(i+2,npdhnoz+k)
        enddo
      enddo
      f(:,npdhnoz+kmax) = 0.
!
! CO (use kldco):
      call lsqdsq(lat,fg(2,nkldco+1,lat+1,ixt),fg(2,nkldco+1,lat,ixt),
     |                fg(2,nkldco+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npdhco) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npdhco+k) = f(i+imax,npdhco+k)
          f(i+imaxp2,npdhco+k) = f(i+2,npdhco+k)
        enddo
      enddo
      f(:,npdhco+kmax) = 0.
!
! CO2 (use kldco2):
      call lsqdsq(lat,fg(2,nkldco2+1,lat+1,ixt),fg(2,nkldco2+1,lat,ixt),
     |                fg(2,nkldco2+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npdhco2) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npdhco2+k) = f(i+imax,npdhco2+k)
          f(i+imaxp2,npdhco2+k) = f(i+2,npdhco2+k)
        enddo
      enddo
      f(:,npdhco2+kmax) = 0.
!
! H2O (use kldh2o):
      call lsqdsq(lat,fg(2,nkldh2o+1,lat+1,ixt),fg(2,nkldh2o+1,lat,ixt),
     |                fg(2,nkldh2o+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npdhh2o) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npdhh2o+k) = f(i+imax,npdhh2o+k)
          f(i+imaxp2,npdhh2o+k) = f(i+2,npdhh2o+k)
        enddo
      enddo
      f(:,npdhh2o+kmax) = 0.
!
! H2 (use kldh2):
      call lsqdsq(lat,fg(2,nkldh2+1,lat+1,ixt),fg(2,nkldh2+1,lat,ixt),
     |                fg(2,nkldh2+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npdhh2) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npdhh2+k) = f(i+imax,npdhh2+k)
          f(i+imaxp2,npdhh2+k) = f(i+2,npdhh2+k)
        enddo
      enddo
      f(:,npdhh2+kmax) = 0.
!
! HOX (use kldhox):
      call lsqdsq(lat,fg(2,nkldhox+1,lat+1,ixt),fg(2,nkldhox+1,lat,ixt),
     |                fg(2,nkldhox+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npdhhox) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npdhhox+k) = f(i+imax,npdhhox+k)
          f(i+imaxp2,npdhhox+k) = f(i+2,npdhhox+k)
        enddo
      enddo
      f(:,npdhhox+kmax) = 0.
!
! CH4 (use kldch4):
      call lsqdsq(lat,fg(2,nkldch4+1,lat+1,ixt),fg(2,nkldch4+1,lat,ixt),
     |                fg(2,nkldch4+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npdhch4) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npdhch4+k) = f(i+imax,npdhch4+k)
          f(i+imaxp2,npdhch4+k) = f(i+2,npdhch4+k)
        enddo
      enddo
      f(:,npdhch4+kmax) = 0.
!
! AR (use kldar):
      call lsqdsq(lat,fg(2,nkldar+1,lat+1,ixt),fg(2,nkldar+1,lat,ixt),
     |                fg(2,nkldar+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npdhar) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npdhar+k) = f(i+imax,npdhar+k)
          f(i+imaxp2,npdhar+k) = f(i+2,npdhar+k)
        enddo
      enddo
      f(:,npdhar+kmax) = 0.
!
! HE (use kldhe):
      call lsqdsq(lat,fg(2,nkldhe+1,lat+1,ixt),fg(2,nkldhe+1,lat,ixt),
     |                fg(2,nkldhe+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npdhhe) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npdhhe+k) = f(i+imax,npdhhe+k)
          f(i+imaxp2,npdhhe+k) = f(i+2,npdhhe+k)
        enddo
      enddo
      f(:,npdhhe+kmax) = 0.
!
! NAT (use kldnat):
      call lsqdsq(lat,fg(2,nkldnat+1,lat+1,ixt),fg(2,nkldnat+1,lat,ixt),
     |                fg(2,nkldnat+1,lat-1,ixt),s1(3,1),long)
      do i=1,long
        f(i+2,npdhnat) = s1(i+2,1)*s2(i+2,1)
      enddo
      do i=1,2 ! periodic points
        do k=0,kmax
          f(i,npdhnat+k) = f(i+imax,npdhnat+k)
          f(i+imaxp2,npdhnat+k) = f(i+2,npdhnat+k)
        enddo
      enddo
      f(:,npdhnat+kmax) = 0.
! 
      return
      end