#include "dims.h"
!
      subroutine hdif1(lat,ixt)
      use cons_module,only: len1,kmax,imaxp4,re_inv,cs,t0,dlamda,dphi
      implicit none
!
! Save global kmh (eddy viscosity) and nrh (M/T) for use in hdif2
!   and hdif3. Global shared arrays fkmh and fnrh are in fgcom.h. 
!   These arrays include the latitude dimension to allow resolution
!   of dependencies when multi-tasking.
! This routine is called for lat indices = -2->jmax (kmh and nrh are
!   defined at those lat indices). On entry, fg contains fields for 
!   lat+1 and lat+2 corresponding to the old f-array indices njp1,njp2).
!
#include "params.h"
#include "fgcom.h"
#include "index.h"
!
! Args:
      integer,intent(in) :: lat,ixt
!
! Locals:
      real :: abcsj,abcsjp,con1,con2,con3,dt,ds
      integer :: i,ii,k,nunmk,nvnmk,ip1,ip1ku,ip1kv
      real :: cp2 = 0.2
!
! KMH = eddy viscosity = 2*K0*K0*SQRT(DS*DS+DT*DT)
! cs(lat) = cos(lat)
! c(1) = dlamda, c(2) = dphi
! c(52) = 1/r (r = earth radius)
! c(109) = 0.2
!
      abcsj = abs(cs(lat+1))
      abcsjp = abs(cs(lat+2))
      con1 = re_inv*.5/dlamda
      con2 = re_inv/(dphi*(abcsj+abcsjp))
      con3 = 2.*cp2*cp2
      ii = 0
      do k=1,kmax
        nunmk = nunm+k
        nvnmk = nvnm+k
        do i=1,imaxp4
          ii = ii+1
!
! ip1 = first dim i-index to fg corresponding to f(i+1..) in the original
!       code, which looped from 1,longp3
! ip1ku,v = 2nd dim k-indices to fg corresponding to f(i+1,unm and vnm) 
!       in the original code.
!
          if (mod(ii+1,imaxp4)==0) then
            ip1 = imaxp4
          else
            ip1 = mod(ii+1,imaxp4)
          endif 
          if (ip1==1) then
            ip1ku = nunmk+1
            ip1kv = nvnmk+1
          else
            ip1ku = nunmk
            ip1kv = nvnmk 
          endif
          dt = 
     | con1*((fg(ip1,ip1ku,lat+2,ixt)-fg(i,nunmk,lat+2,ixt))/cs(lat+2)+
     |       (fg(ip1,ip1ku,lat+1,ixt)-fg(i,nunmk,lat+1,ixt))/cs(lat+1))-
     | con2*((fg(ip1,ip1kv,lat+2,ixt)+fg(i,nvnmk,lat+2,ixt))*abcsjp -
     |       (fg(ip1,ip1kv,lat+1,ixt)+fg(i,nvnmk,lat+1,ixt))*abcsj)
          ds =
     | con1*((fg(ip1,ip1kv,lat+2,ixt)-fg(i,nvnmk,lat+2,ixt))/cs(lat+2)+
     |       (fg(ip1,ip1kv,lat+1,ixt)-fg(i,nvnmk,lat+1,ixt))/cs(lat+1))+
     | con2*((fg(ip1,ip1ku,lat+2,ixt)+fg(i,nunmk,lat+2,ixt))*abcsjp -
     |       (fg(ip1,ip1ku,lat+1,ixt)+fg(i,nunmk,lat+1,ixt))*abcsj)
          fkmh(i,k,lat) = con3*sqrt(ds*ds+dt*dt)
        enddo
      enddo
      fkmh(imaxp4,kmax,lat) = 0.
!
! NRH = ((nms(k)+nms(k+1))*0.5) / (tnm(k)+(t0(k)+t0(k+1))*0.5)
      do k=1,kmax
        do i=1,len1
          fnrh(i,k,lat) = 
     |      ((fg(i,nms+k,lat+1,ixt)+fg(i,nms+k+1,lat+1,ixt))*0.5) /
     |      (fg(i,ntnm+k,lat+1,ixt)+((t0(k)+t0(k+1))*.5))
        enddo
      enddo
      end subroutine hdif1