#include "dims.h"
!
      subroutine hdif2(lat,ixt)
      use cons_module,only: long,kmax,imaxp4
      implicit none
!
! This routine defines KLD terms in fg, for use by hdif3.
! These are 16 (nfkld) fields NKLDU->NKLDNAT, intermediate
!   to the final horizontal diffusion fields.
! This uses fkmh (lat-1) and fnrh(lat), which were defined 
!   at all lats by hdif1. 
! There are no references to f-array here, only fg, fkmh, and fnrh.
!
#include "params.h"
#include "fgcom.h"
#include "index.h"
!
! Args:
      integer,intent(in) :: 
     |  lat, ! current latitude index
     |  ixt  ! index to 4th dimension of fg
!
! Local:
      integer :: ii,i,k,longp2,ip1,ip1k,n,nn
      real :: 
     |  avkmh(zimxp,zkmxp),    ! average kmh (s1 in old code)
     |  rhoavkmh(zimxp,zkmxp)  ! nrh*avkmh (s2 in old code)
!
! CALC AVG KMH VALUE AT EACH GRID PT ON J+1
      ii = 0
      do k=1,kmax
        do i=1,imaxp4
          ii = ii+1
          if (mod(ii+1,imaxp4)==0) then
            ip1 = imaxp4
          else
            ip1 = mod(ii+1,imaxp4)
          endif
          if (ip1==1) then
            ip1k = k+1
          else
            ip1k = k
          endif
          avkmh(ip1,ip1k) = ((fkmh(i,k,lat-1)+fkmh(ip1,ip1k,lat-1)+
     |                   fkmh(i,k,lat))+fkmh(ip1,ip1k,lat))*0.25
        enddo
      enddo
!
! S2=RHO*AVG(KMH)
      do k=1,kmax
        do i=1,imaxp4
          rhoavkmh(i,k) = avkmh(i,k)*fnrh(i,k,lat)
        enddo
      enddo
!
! From con.f: LONG=LEN2-4, LEN2=IMAXP4*KMAX
      longp2 = long+2
!
! CALC RHO*KMH*(L*L(D*D)(PSI)) AT J+1 AND N-1
!
! FOR PSI = U
      call lsqdsq(lat+1,fg(1,nunm+1,lat+2,ixt),fg(1,nunm+1,lat+1,ixt),
     |            fg(1,nunm+1,lat,ixt),avkmh(2,1),longp2)
      do i=1,longp2
        fg(i+1,nkldu+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! FOR PSI = V
      call lsqdsq(lat+1,fg(1,nvnm+1,lat+2,ixt),fg(1,nvnm+1,lat+1,ixt),
     |            fg(1,nvnm+1,lat,ixt),avkmh(2,1),longp2)
      do i=1,longp2
        fg(i+1,nkldv+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! FOR PSI=T
      call lsqdsq(lat+1,fg(1,ntnm+1,lat+2,ixt),fg(1,ntnm+1,lat+1,ixt),
     |            fg(1,ntnm+1,lat,ixt),avkmh(2,1),longp2)
      do i=1,longp2
        fg(i+1,nkldt+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! FOR PSI=PSI (O2)
      call lsqdsq(lat+1,fg(1,npsnm+1,lat+2,ixt),
     |  fg(1,npsnm+1,lat+1,ixt),fg(1,npsnm+1,lat,ixt),avkmh(2,1),longp2)
      do i=1,longp2
        fg(i+1,nkldps+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! FOR PSI=PSI (O)
      call lsqdsq(lat+1,fg(1,nps2nm+1,lat+2,ixt),
     |  fg(1,nps2nm+1,lat+1,ixt),fg(1,nps2nm+1,lat,ixt),avkmh(2,1),
     |  longp2)
      do i=1,longp2
        fg(i+1,nkldp2+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! Following 11 minor species are for time-gcm only (not tiegcm):
!
! N4S:
      call lsqdsq(lat+1,fg(1,nn4snm+1,lat+2,ixt),
     |  fg(1,nn4snm+1,lat+1,ixt),fg(1,nn4snm+1,lat,ixt),
     |    avkmh(2,1),longp2)
      do i=1,longp2
        fg(i+1,nkldn4s+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! NOZ:
      call lsqdsq(lat+1,fg(1,npnoznm+1,lat+2,ixt),
     |  fg(1,npnoznm+1,lat+1,ixt),fg(1,npnoznm+1,lat,ixt),avkmh(2,1),
     |  longp2)
      do i=1,longp2
        fg(i+1,nkldnoz+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! CO:
      call lsqdsq(lat+1,fg(1,npconm+1,lat+2,ixt),
     |  fg(1,npconm+1,lat+1,ixt),fg(1,npconm+1,lat,ixt),avkmh(2,1),
     |  longp2)
      do i=1,longp2
        fg(i+1,nkldco+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! CO2:
      call lsqdsq(lat+1,fg(1,npco2nm+1,lat+2,ixt),
     |  fg(1,npco2nm+1,lat+1,ixt),fg(1,npco2nm+1,lat,ixt),avkmh(2,1),
     |  longp2)
      do i=1,longp2
        fg(i+1,nkldco2+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! H2O:
      call lsqdsq(lat+1,fg(1,nph2onm+1,lat+2,ixt),
     |  fg(1,nph2onm+1,lat+1,ixt),fg(1,nph2onm+1,lat,ixt),avkmh(2,1),
     |  longp2)
      do i=1,longp2
        fg(i+1,nkldh2o+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! H2:
      call lsqdsq(lat+1,fg(1,nph2nm+1,lat+2,ixt),
     |  fg(1,nph2nm+1,lat+1,ixt),fg(1,nph2nm+1,lat,ixt),avkmh(2,1),
     |  longp2)
      do i=1,longp2
        fg(i+1,nkldh2+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! HOX:
      call lsqdsq(lat+1,fg(1,nphoxnm+1,lat+2,ixt),
     |  fg(1,nphoxnm+1,lat+1,ixt),fg(1,nphoxnm+1,lat,ixt),avkmh(2,1),
     |  longp2)
      do i=1,longp2
        fg(i+1,nkldhox+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! CH4:
      call lsqdsq(lat+1,fg(1,npch4nm+1,lat+2,ixt),
     |  fg(1,npch4nm+1,lat+1,ixt),fg(1,npch4nm+1,lat,ixt),avkmh(2,1),
     |  longp2)
      do i=1,longp2
        fg(i+1,nkldch4+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! AR:
      call lsqdsq(lat+1,fg(1,nparnm+1,lat+2,ixt),
     |  fg(1,nparnm+1,lat+1,ixt),fg(1,nparnm+1,lat,ixt),avkmh(2,1),
     |  longp2)
      do i=1,longp2
        fg(i+1,nkldar+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! HE:
      call lsqdsq(lat+1,fg(1,nphenm+1,lat+2,ixt),
     |  fg(1,nphenm+1,lat+1,ixt),fg(1,nphenm+1,lat,ixt),avkmh(2,1),
     |  longp2)
      do i=1,longp2
        fg(i+1,nkldhe+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
! NAT:
      call lsqdsq(lat+1,fg(1,npnatnm+1,lat+2,ixt),
     |  fg(1,npnatnm+1,lat+1,ixt),fg(1,npnatnm+1,lat,ixt),avkmh(2,1),
     |  longp2)
      do i=1,longp2
        fg(i+1,nkldnat+1,lat+1,ixt) = avkmh(i+1,1)*rhoavkmh(i+1,1) 
      enddo
!
      return
      end