#include "dims.h"
!
      subroutine hdif1(lat,ixt)
      use cons_module,only: len1,imaxp4,kmax,t0,cs,dlamda,dphi,re_inv,
     |  kmaxp1
      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,fmin,fmax
      integer :: i,ii,k,nunmk,nvnmk,ip1,ip1ku,ip1kv
      real :: cp2 = 0.2
!
!     call fminmax(fg(1,ntnm+1,lat+1,ixt),imaxp4*kmaxp1,
!    |  fmin,fmax)
!     write(6,"('hdif1: tn_nm at lat+1=',i2,' min,max=',2e12.4)")
!    |  lat+1,fmin,fmax
!     call fminmax(fg(1,ntnm+1,lat+2,ixt),imaxp4*kmaxp1,
!    |  fmin,fmax)
!     write(6,"('hdif1: tn_nm at lat+2=',i2,' min,max=',2e12.4)")
!    |  lat+2,fmin,fmax
!
!     call fminmax(fg(1,nunm+1,lat+1,ixt),imaxp4*kmaxp1,
!    |  fmin,fmax)
!     write(6,"('hdif1: un_nm at lat+1=',i2,' min,max=',2e12.4)")
!    |  lat+1,fmin,fmax
!     call fminmax(fg(1,nunm+1,lat+2,ixt),imaxp4*kmaxp1,
!    |  fmin,fmax)
!     write(6,"('hdif1: un_nm at lat+2=',i2,' min,max=',2e12.4)")
!    |  lat+2,fmin,fmax
!
!     call fminmax(fg(1,nvnm+1,lat+1,ixt),imaxp4*kmaxp1,
!    |  fmin,fmax)
!     write(6,"('hdif1: vn_nm at lat+1=',i2,' min,max=',2e12.4)")
!    |  lat+1,fmin,fmax
!     call fminmax(fg(1,nvnm+1,lat+2,ixt),imaxp4*kmaxp1,
!    |  fmin,fmax)
!     write(6,"('hdif1: vn_nm at lat+2=',i2,' min,max=',2e12.4)")
!    |  lat+2,fmin,fmax
!
!     call fminmax(fg(1,nms+1,lat+1,ixt),imaxp4*kmaxp1,
!    |  fmin,fmax)
!     write(6,"('hdif1: fnms  at lat+1=',i2,' min,max=',2e12.4)")
!    |  lat+1,fmin,fmax
!     call fminmax(fg(1,nms+1,lat+2,ixt),imaxp4*kmaxp1,
!    |  fmin,fmax)
!     write(6,"('hdif1: fnms  at lat+2=',i2,' min,max=',2e12.4)")
!    |  lat+2,fmin,fmax
!
! KMH = eddy viscosity = 2*K0*K0*SQRT(DS*DS+DT*DT)
! cs(lat) = cos(lat)
!
      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

!         write(6,"('lat=',i2,' k=',i2,' i=',i2,' ip1=',i2,
!    |      /,' un_nm=',4e12.4,/,' vn_nm=',4e12.4)") lat,k,i,ip1,
!    |      fg(ip1,ip1ku,lat+2,ixt),fg(i,nunmk,lat+2,ixt),
!    |      fg(ip1,ip1ku,lat+1,ixt),fg(i,nunmk,lat+1,ixt),
!    |      fg(ip1,ip1kv,lat+2,ixt),fg(i,nvnmk,lat+2,ixt),
!    |      fg(ip1,ip1kv,lat+1,ixt),fg(i,nvnmk,lat+1,ixt)

          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 ! i=1,imaxp4

!         if (lat >= zjmx-1) then ! look at top 2 lats
!           write(6,"(/,'hdif1: lat=',i2,' k=',i2,
!    |        ' fg(1:imaxp4,ip1ku,lat+2,ixt)=',/,(6e12.4))")
!    |        lat,k,fg(1:imaxp4,ip1ku,lat+2,ixt)
!
!           write(6,"(/,'hdif1: lat=',i2,' k=',i2)") lat,k
!           write(6,"('  fkmh(1:imaxp4,k,lat)=',/,(6e12.4))") 
!    |        fkmh(1:38,k,lat)
!           write(6,"(6e12.4)") fkmh(39:imaxp4,k,lat)
!         endif

      enddo ! k=1,kmax

!     call fminmax(fkmh(:,1:kmax,lat),imaxp4*kmax,fmin,fmax)
!     write(6,"('hdif1: fkmh  at lat=',i2,' min,max=',2e12.4)")
!    |  lat,fmin,fmax
!     if (lat > 0)
!    |  call addfsech('FKMH',' ',' ',fkmh(1,1,lat),
!    |    zimxp,zkmxp,zkmx,lat)

      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

!       if (lat >= zjmx-1) then ! look at top 2 lats
!         write(6,"(/,'hdif1: lat=',i2,' k=',i2,
!    |      ' fg(1:imaxp4,k,lat+1,ixt)=',/,(6e12.4))")
!    |      lat,k,fg(1:38,k,lat+1,ixt)
!         write(6,"(6e12.4)") fg(39:imaxp4,k,lat+1,ixt)

!         write(6,"(/,'hdif1: lat=',i2,' k=',i2)") lat,k
!         write(6,"('  fnrh(1:imaxp4,k,lat)=',/,(6e12.4))")
!    |      fnrh(1:38,k,lat)
!         write(6,"(6e12.4)") fnrh(39:imaxp4,k,lat)
!       endif

      enddo

!     call fminmax(fnrh(:,1:kmax,lat),imaxp4*kmax,fmin,fmax)
!     write(6,"('hdif1: fnrh  at lat=',i2,' min,max=',2e12.4)")
!    |  lat,fmin,fmax
!     if (lat > 0)
!    |  call addfsech('FNRH',' ',' ',fnrh(1,1,lat),
!    |    zimxp,zkmxp,zkmx,lat)

      end subroutine hdif1