#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
!     if (lat > 0)
!    |  call addfsech('AVKMH',' ',' ',avkmh,zimxp,zkmxp,zkmx,lat)
!
! S2=RHO*AVG(KMH)
      do k=1,kmax
        do i=1,imaxp4
          rhoavkmh(i,k) = avkmh(i,k)*fnrh(i,k,lat)
        enddo
      enddo
!     if (lat > 0)
!    |  call addfsech('RHOKMH',' ',' ',rhoavkmh,zimxp,zkmxp,zkmx,lat)
!
! CALC RHO*KMH*(L*L(D*D)(PSI)) AT J+1 AND N-1
!
! FOR PSI = U
      longp2 = long+2
      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
!     if (lat > 0)
!    |  call addfsech('FKLDU',' ',' ',fg(1,nkldu+1,lat+1,ixt),
!    |    zimxp,zkmxp,zkmx,lat)
!
! 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
!     if (lat > 0)
!    |  call addfsech('FKLDV',' ',' ',fg(1,nkldv+1,lat+1,ixt),
!    |    zimxp,zkmxp,zkmx,lat)
!
! 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
!     if (lat > 0)
!    |  call addfsech('FKLDT',' ',' ',fg(1,nkldt+1,lat+1,ixt),
!    |    zimxp,zkmxp,zkmx,lat)
!
! 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
!     if (lat > 0)
!    |  call addfsech('FKLDO2',' ',' ',fg(1,nkldps+1,lat+1,ixt),
!    |    zimxp,zkmxp,zkmx,lat)
!
! 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
!     if (lat > 0)
!    |  call addfsech('FKLDO1',' ',' ',fg(1,nkldp2+1,lat+1,ixt),
!    |    zimxp,zkmxp,zkmx,lat)
      return
      end