! 
      module hdif_module
      use params_module,only: nlevp1,nlonp4,nlat,spval
      use mpi_module,only: lon0,lon1,lat0,lat1
      use addfld_module,only: addfld
      implicit none
!
! An extra latitude is needed in fnrh and fkmh: 
! 10/17/02 bf: changed lon dim from nlonp4+1 to 0:nlonp4+1 
! 9/13/06 btf: Changed ch4 to o1 and h2 to co for venus vtgcm.
!
      real,dimension(nlevp1,0:nlonp4+1,-2:nlat) ::
     |  fnrh,  ! eddy viscosity
     |  fkmh   ! M/T
! -----------------------------------------------------------------
! Titan TGCM modifications. Latest version: 9/15/04 S. W. BOUGHER
! -----------------------------------------------------------------
!
! VT vampir tracing:
!
#ifdef VT
#include "VT.inc"
#endif
!
      contains
!-----------------------------------------------------------------------
!
      subroutine hdif1(tn_nm,un_nm,vn_nm,barm,
     |  lev0,lev1,lon0,lon1,lat0,lat1)
!
! Save global kmh (eddy viscosity) and nrh (M/T) for use in hdif2
!   and hdif3. This routine is called from advance.
!
      use cons_module,only: t0,cs,dlamda,dphi,re_inv
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in)::
     |  tn_nm, ! neutral temperature at time n-1
     |  un_nm, ! zonal velocity at time n-1
     |  vn_nm, ! meridional velocity at time n-1
     |  barm   ! mean molecular weight
!
! Local:
      integer :: k,ii,i,j,ip1,ip1k,nlons,nlevs,lat,ier
      real :: fmin,fmax
      real :: abcsj,abcsjp,con1,con2,con3,delt,dels
      real :: cp2 = 0.2
!
#ifdef VT
!     code = 119 ; state = 'hdif1' ; activity='ModelCode'
      call vtbegin(119,ier)
#endif
!
      nlons = lon1-lon0+1
      nlevs = lev1-lev0+1
!
! Latitude loop starts at lat0-3. Inputs are referenced only through
! lat+2, but kmh and nrh are defined from lat0-3 to lat1.
!
      do lat=lat0-3,lat1
!
! 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 = 1 -> nlev*nlonp4
!   i = 1 -> 75,76
! ip1 = 2 -> 76, 1
!
        do i=lon0,lon1
          do k=lev0,lev1-1
            ip1 = i+1

!           if (i==nlonp4) ip1 = 1
!           if (ip1==1) then

! 10/17/02 bf: if i==nlonp4, then current task does not have i==1, so try
!   using ip1==nlonp4-3 instead (it should be same as i==1 because of
!   periodic points). 
!
            if (i==nlonp4) ip1 = 9999
            if (ip1==9999) then
              ip1k = k+1
              ip1 = nlonp4-3
            else
              ip1k = k
            endif

!           write(6,"('lat=',i2,' i=',i2,' k=',i2,' ip1k=',i2,' ip1=',
!    |        i2,' lat0,1=',2i3,' lon0,1=',2i3,' lev0,1=',2i3)")
!    |        lat,i,k,ip1k,ip1,lat0,lat1,lon0,lon1,lev0,lev1

            delt = 
     |        con1*((un_nm(ip1k,ip1,lat+2)-un_nm(k,i,lat+2))/cs(lat+2)+
     |              (un_nm(ip1k,ip1,lat+1)-un_nm(k,i,lat+1))/cs(lat+1))-
     |        con2*((vn_nm(ip1k,ip1,lat+2)+vn_nm(k,i,lat+2))*abcsjp -
     |              (vn_nm(ip1k,ip1,lat+1)+vn_nm(k,i,lat+1))*abcsj)
            dels = 
     |        con1*((vn_nm(ip1k,ip1,lat+2)-vn_nm(k,i,lat+2))/cs(lat+2)+
     |              (vn_nm(ip1k,ip1,lat+1)-vn_nm(k,i,lat+1))/cs(lat+1))+
     |        con2*((un_nm(ip1k,ip1,lat+2)+un_nm(k,i,lat+2))*abcsjp -
     |              (un_nm(ip1k,ip1,lat+1)+un_nm(k,i,lat+1))*abcsj)
            fkmh(k,i,lat) = con3*sqrt(dels*dels+delt*delt)
          enddo ! i=lon0,lon1
        enddo ! k=lev0,lev1-1
        fkmh(lev1,:,lat) = spval

!       if (lat >= lat0)
!    |    call addfld('FKMH',' ',' ',fkmh(:,lon0:lon1,lat),
!    |      'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! NRH = ((barm(k)+barm(k+1))*0.5) / (tnm(k)+(t0(k)+t0(k+1))*0.5)
        do i=lon0,lon1
          do k=lev0,lev1-1
            fnrh(k,i,lat) = ((barm(k,i,lat+1)+barm(k+1,i,lat+1))*0.5) /
     |                      (tn_nm(k,i,lat+1)+((t0(k)+t0(k+1))*.5))
          enddo
        enddo
        fnrh(lev1,:,lat) = spval
 
!       if (lat >= lat0)
!    |    call addfld('FNRH',' ',' ',fnrh(:,lon0:lon1,lat),
!    |      'lev',lev0,lev1,'lon',lon0,lon1,lat)

      enddo ! lat=lat0-3,lat1
#ifdef MPI
!
! Exchange boundary lats and lons in fkmh:
      call mp_bndlats_kmh(lev0,lev1,lon0,lon1,lat0,lat1)
      call mp_bndlons_kmh(lev0,lev1,lon0,lon1,lat0,lat1)
#endif
!     do j=lat0,lat1
!         call addfld('FNRH',' ',' ',fnrh(:,lon0:lon1,j),
!    |      'lev',lev0,lev1,'lon',lon0,lon1,j)
!         call addfld('FKMH',' ',' ',fkmh(:,lon0:lon1,j),
!    |      'lev',lev0,lev1,'lon',lon0,lon1,j)
!     enddo ! j=lat0,lat1
!
#ifdef VT
!     code = 119 ; state = 'hdif1' ; activity='ModelCode'
      call vtend(119,ier)
#endif
      end subroutine hdif1
!-----------------------------------------------------------------------
!
      subroutine hdif2(tn_nm,un_nm,vn_nm,o1_nm,co_nm,
     |  fkldt,fkldu,fkldv,fkldo1,fkldco,lev0,lev1,lon0,lon1,lat0,lat1)
!
! Save fkldx horizontal diffusion terms for hdif3. fkmh and fnrh are 
! module data, and were calculated by hdif1.
!
      implicit none
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in):: 
     |  tn_nm,  ! neutral temperature at time n-1 (deg K)
     |  un_nm,  ! zonal velocity at time n-1 (cm/s)
     |  vn_nm,  ! meridional velocity at time n-1 (cm/s)
     |  o1_nm,  ! atomic oxygen at time n-1 (mmr)
     |  co_nm   ! carbon monoxide at time n-1 (mmr)
!
! Ouput args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  fkldt,  ! tn horizontal diffusion term
     |  fkldu,  ! un horizontal diffusion term
     |  fkldv,  ! vn horizontal diffusion term
     |  fkldo1, ! o horizontal diffusion term
     |  fkldco  ! co horizontal diffusion term
!
! Local:
      integer :: k,i,ip1,ip1k,nlevs,lonbeg,lonend,lat,ier
      real :: 
     |  avkmh (lev0:lev1,lon0:lon1), ! average kmh
     |  rhokmh(lev0:lev1,lon0:lon1)  ! nrh*avkmh
      real :: fmin,fmax
!
#ifdef VT
!     code = 120 ; state = 'hdif2' ; activity='ModelCode'
      call vtbegin(120,ier)
#endif
!
! Latitude scan:
      do lat=lat0-2,lat1
!       write(6,"('hdif2: lat=',i3)") lat
!
! avkmh = average fkmh:
! Allocation of fkmh in sub init_kmh_nrh (fields_module.F):
!     allocate(fkmh%data3d(nlevp1,nlonp4+1,-2:nlat),stat=ier)
!
      nlevs = lev1-lev0+1 ! for addfsech calls
      do i=lon0,lon1
        do k=lev0,lev1-1
          avkmh(k,i) = ((fkmh(k,i-1,lat-1)+fkmh(k,i,lat-1)+
     |                   fkmh(k,i-1,lat)) +fkmh(k,i,lat))*0.25
        enddo
        avkmh(lev1,i) = spval
      enddo
!     if (lat > 0)
!    |  call addfld('AVKMH',' ',' ',avkmh(lon0:lon1,:),
!    |      'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! rhokmh = rho*avg(kmh)
      do i=lon0,lon1
        do k=lev0,lev1-1
          rhokmh(k,i) = avkmh(k,i)*fnrh(k,i,lat)
        enddo
        rhokmh(lev1,i) = 0. ! added to prevent NaNS init fpe
      enddo
!
!     if (lat > 0)
!    |  call addfld('RHOKMH',' ',' ',rhokmh(lon0:lon1,:),
!    |      'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! CALC RHO*KMH*(L*L(D*D)(PSI)) AT J+1 AND N-1
!
! Define fkldx at i=2,nlonp4-1 and lat+1.
! (note pronostic inputs at time n-1 have lon0-2->lon1+2 from bndlons 
!  call in advance)
!
      lonbeg = lon0
      if (lon0==1) lonbeg = 2
      lonend = lon1
      if (lon1==nlonp4) lonend = nlonp4-1
!
! FOR PSI = U
      call lsqdsq(un_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            un_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            un_nm(lev0:lev1,lonbeg-1:lonend+1,lat),
     |            avkmh(lev0:lev1,lonbeg:lonend),lonbeg,lonend,
     |            lev0,lev1,lat+1)
      do i=lonbeg,lonend
        do k=lev0,lev1
          fkldu(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
        fkldu(lev1,i,lat+1) = spval
      enddo
!     if (lat > 0)
!    |  call addfld('FKLDU',' ',' ',fkldu(:,lon0:lon1,lat+1),
!    |      'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! FOR PSI = V
      call lsqdsq(vn_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            vn_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            vn_nm(lev0:lev1,lonbeg-1:lonend+1,lat),
     |            avkmh(lev0:lev1,lonbeg:lonend),lonbeg,lonend,
     |            lev0,lev1,lat+1)
      do i=lonbeg,lonend
        do k=lev0,lev1
          fkldv(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
        fkldv(lev1,i,lat+1) = spval
      enddo
!     if (lat > 0)
!    |  call addfld('FKLDV',' ',' ',fkldv(:,lon0:lon1,lat+1),
!    |      'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! FOR PSI = T
      call lsqdsq(tn_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            tn_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            tn_nm(lev0:lev1,lonbeg-1:lonend+1,lat),
     |            avkmh(lev0:lev1,lonbeg:lonend),lonbeg,lonend,
     |            lev0,lev1,lat+1)
      do i=lonbeg,lonend
        do k=lev0,lev1
          fkldt(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
        fkldt(lev1,i,lat+1) = spval
      enddo
!     if (lat > 0)
!    |  call addfld('FKLDT',' ',' ',fkldt(:,lon0:lon1,lat+1),
!    |      'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! FOR PSI = CH4
      call lsqdsq(o1_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            o1_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            o1_nm(lev0:lev1,lonbeg-1:lonend+1,lat),
     |            avkmh(lev0:lev1,lonbeg:lonend),lonbeg,lonend,
     |            lev0,lev1,lat+1)
      do i=lonbeg,lonend
        do k=lev0,lev1
          fkldo1(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
        fkldo1(lev1,i,lat+1) = spval
      enddo
!     if (lat > 0)
!    |  call addfld('FKLDO1',' ',' ',fkldo1(:,lon0:lon1,lat+1),
!    |      'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! FOR PSI = H2
      call lsqdsq(co_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            co_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            co_nm(lev0:lev1,lonbeg-1:lonend+1,lat),
     |            avkmh(lev0:lev1,lonbeg:lonend),lonbeg,lonend,
     |            lev0,lev1,lat+1)
      do i=lonbeg,lonend
        do k=lev0,lev1
          fkldco(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
        fkldco(lev1,i,lat+1) = spval
      enddo
!     if (lat > 0)
!    |  call addfld('FKLDH2',' ',' ',fkldco(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
      enddo ! lat=lat0-2,lat1
!
#ifdef VT
!     code = 120 ; state = 'hdif2' ; activity='ModelCode'
      call vtend(120,ier)
#endif
      end subroutine hdif2
!-----------------------------------------------------------------------
      subroutine hdif_bndlons(kldt,kldu,kldv,kldo1,kldco,
     |  lev0,lev1,lon0,lon1,lat0,lat1)
!
! Exchange boundary longitudes of horizontal diffusion coeffs
! calculated by hdif2. This is called from dynamics, in preparation 
! for hdif3.
!
#ifdef MPI
      use mpi_module,only: mp_bndlons_f3d
#endif
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(inout) :: kldt,kldu,kldv,kldo1,kldco
!
! Local:
      integer :: nlevs
!
#ifdef MPI
!
! subroutine mp_bndlons_f3d(f,id1,i0,i1,j0,j1,nf)
! From oplus: call mp_bndlons_f3d(dj,nlevs,lon0,lon1,lat0-2,lat1+2,1)
! Could also try transferring to a local 5-field array, and make only
!   one call to mp_bndlons_f3d.
!
      nlevs = lev1-lev0+1
      call mp_bndlons_f3d(kldt ,nlevs,lon0,lon1,lat0-2,lat1+2,1)
      call mp_bndlons_f3d(kldu ,nlevs,lon0,lon1,lat0-2,lat1+2,1)
      call mp_bndlons_f3d(kldv ,nlevs,lon0,lon1,lat0-2,lat1+2,1)
      call mp_bndlons_f3d(kldo1, nlevs,lon0,lon1,lat0-2,lat1+2,1)
      call mp_bndlons_f3d(kldco,nlevs,lon0,lon1,lat0-2,lat1+2,1)
#endif
      end subroutine hdif_bndlons
!-----------------------------------------------------------------------
      subroutine hdif3(cp,
     |  kldt,kldu,kldv,kldo1,kldco,     ! input
     |  hdt ,hdu ,hdv ,hdo1 ,hdco,      ! output
     |  lev0,lev1,lon0,lon1,lat)
!
! Calculate horizontal diffusion terms for t,u,v,o1,co at current latitude, 
! using coefficients that were output by hdif2. (hdif2 is called from advance, 
! hdif3 is called from dynamics). Sub kld_bndlons has been called prior to
! this routine so tasks have boundary longitudes for kldt,u,v,o1,co.
!
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
!
! 2d input:
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: cp
!
! 3d input at full task subdomain:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in):: 
     |  kldt,kldu,kldv,kldo1,kldco ! input from hdif2
!
! 2d output for current latitude:
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out):: 
     |  hdt ,hdu ,hdv ,hdo1 ,hdco  ! output
!
! Local:
      integer :: k,i,lonbeg,lonend,nlevs
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  fnrh_inv, hdout, cpi 
!
      nlevs = lev1-lev0+1
!
! fnrh is eddy viscosity (hdif module data), output by hdif1:
      do i=lon0,lon1
        do k=lev0,lev1-1
          fnrh_inv(k,i) = -1. / fnrh(k,i,lat-1)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
! At 5 deg res (nlonp4=76), lonbeg,end will be 3,74
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = nlonp4-2
!
! Make hdu from kldu (hdu was nflh):
      call lsqdsq(kldu(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldu(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldu(lev0:lev1,lonbeg-1:lonend+1,lat-1),
     |            hdout(lev0:lev1,lonbeg:lonend),lonbeg,lonend,
     |            lev0,lev1,lat)
      do i=lonbeg,lonend
        do k=lev0,lev1-1
          hdu(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdu(lev1,i) = 0.
      enddo ! i=lonbeg,lonend

!     call addfld('HDU',' ',' ',hdu(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Make hdv from kldv (hdv was nfph):
      call lsqdsq(kldv(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldv(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldv(lev0:lev1,lonbeg-1:lonend+1,lat-1),
     |            hdout(lev0:lev1,lonbeg:lonend),lonbeg,lonend,
     |            lev0,lev1,lat)
      do i=lonbeg,lonend
        do k=lev0,lev1-1
          hdv(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdv(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!     call addfld('HDV',' ',' ',hdv(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Make hdt from kldt (was nqdh):
      call lsqdsq(kldt(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldt(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldt(lev0:lev1,lonbeg-1:lonend+1,lat-1),
     |            hdout(lev0:lev1,lonbeg:lonend),lonbeg,lonend,
     |            lev0,lev1,lat)
      do i=lon0,lon1
        do k=lev0,lev1-1
          cpi(k,i) = .5*(cp(k,i)+cp(k+1,i))
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
      do i=lonbeg,lonend
        do k=lev0,lev1-1
          hdt(k,i) = hdout(k,i)*fnrh_inv(k,i)*cpi(k,i)
        enddo ! k=lev0,lev1-1
        hdt(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!     call addfld('HDT',' ',' ',hdt(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Make hdo1 from kldco (was npsdh):
      call lsqdsq(kldo1(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldo1(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldo1(lev0:lev1,lonbeg-1:lonend+1,lat-1),
     |            hdout(lev0:lev1,lonbeg:lonend),lonbeg,lonend,
     |            lev0,lev1,lat)
      do i=lonbeg,lonend
        do k=lev0,lev1-1
          hdo1(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdo1(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!     call addfld('HDCH4',' ',' ',hdo1(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Make hdco from kldco
      call lsqdsq(kldco(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldco(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldco(lev0:lev1,lonbeg-1:lonend+1,lat-1),
     |            hdout(lev0:lev1,lonbeg:lonend),lonbeg,lonend,
     |            lev0,lev1,lat)
      do i=lonbeg,lonend
        do k=lev0,lev1-1
          hdco(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdco(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!     call addfld('HDH2',' ',' ',hdco(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
      end subroutine hdif3
!-----------------------------------------------------------------------
      subroutine hdif_periodic(hdt,hdu,hdv,hdo1,hdco,
     |  lev0,lev1,lon0,lon1,lat0,lat1)
#ifdef MPI
      use mpi_module,only: mp_periodic_f3d
#endif
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(inout) :: hdt,hdu,hdv,hdo1,hdco
#ifdef MPI
      call mp_periodic_f3d(hdt(:,lon0:lon1,lat0-1:lat1+1),lev0,lev1,
     |  lon0,lon1,lat0,lat1)
      call mp_periodic_f3d(hdu(:,lon0:lon1,lat0-1:lat1+1),lev0,lev1,
     |  lon0,lon1,lat0,lat1)
      call mp_periodic_f3d(hdv(:,lon0:lon1,lat0-1:lat1+1),lev0,lev1,
     |  lon0,lon1,lat0,lat1)
      call mp_periodic_f3d(hdo1(:,lon0:lon1,lat0-1:lat1+1),lev0,lev1,
     |  lon0,lon1,lat0,lat1)
      call mp_periodic_f3d(hdco(:,lon0:lon1,lat0-1:lat1+1),lev0,lev1,
     |  lon0,lon1,lat0,lat1)
#endif
      end subroutine hdif_periodic
!-----------------------------------------------------------------------
#ifdef MPI
      subroutine mp_bndlats_kmh(lev0,lev1,lon0,lon1,lat0,lat1)
      use mpi_module,only: mytidi,mytidj,itask_table,handle_mpi_err
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
#include "mpif.h"
!
! Local:
      integer ::
     |  irstat(MPI_STATUS_SIZE)      ! mpi receive status
      real,allocatable :: sndbuf(:,:,:),rcvbuf(:,:,:)
      integer :: ier,nlons,len,jprev,jnext,ireqrecv,ireqsend
#ifdef VT
!     call vtsymdef(111, 'mp_bndlats_kmh','Communication',ier)
      call vtbegin(111,ier)
#endif
!
      nlons = lon1-lon0+1
      allocate(sndbuf(nlevp1,nlons,3),stat=ier)
      if (ier /= 0) write(6,"('>>> mp_bndlats_kmh: error allocating',
     |  ' sndbuf: nlons=',i3,' nlevp1=',i3)") nlons,nlevp1
      allocate(rcvbuf(nlevp1,nlons,3),stat=ier)
      if (ier /= 0) write(6,"('>>> mp_bndlats_kmh: error allocating',
     |  ' rcvbuf: nlons=',i3,' nlevp1=',i3)") nlons,nlevp1
      len = nlons*nlevp1*3
!
      jprev = itask_table(mytidi,mytidj-1) ! task to south
      jnext = itask_table(mytidi,mytidj+1) ! task to north
!
! Send lat1-1,lat1 to jnext:
      sndbuf(:,:,:) = fkmh(:,lon0:lon1,lat1-2:lat1)
      call mpi_isend(sndbuf,len,MPI_REAL8,jnext,1,MPI_COMM_WORLD,
     |  ireqsend,ier)      
      if (ier /= 0) 
     |  call handle_mpi_err(ier,'mp_bndlats_kmh send to jnext')
!
! Receive lat0-2,lat0-1 from jprev:
      call mpi_irecv(rcvbuf,len,MPI_REAL8,jprev,1,MPI_COMM_WORLD,
     |  ireqrecv,ier)
      if (ier /= 0) 
     |  call handle_mpi_err(ier,'mp_bndlats_kmh recv from jprev')
!
! Wait for completions:
      call mpi_wait(ireqsend,irstat,ier)
      if (ier /= 0)
     |  call handle_mpi_err(ier,'mp_bndlats_kmh wait for send')
      call mpi_wait(ireqrecv,irstat,ier)
      if (ier /= 0)
     |  call handle_mpi_err(ier,'mp_bndlats_kmh wait for recv')
!
! Define array from receive buffer:
      if (lat0 > 1) ! i.e., mytaskj > 0
     |  fkmh(:,lon0:lon1,lat0-3:lat0-1) = rcvbuf(:,:,:)
!
! Release local buffer space:
      deallocate(sndbuf)
      deallocate(rcvbuf)
#ifdef VT
!     call vtsymdef(111, 'mp_bndlats_kmh','Communication',ier)
      call vtend(111,ier)
#endif
      end subroutine mp_bndlats_kmh
!-----------------------------------------------------------------------
      subroutine mp_bndlons_kmh(lev0,lev1,lon0,lon1,lat0,lat1)
      use mpi_module,only: mytidi,mytidj,itask_table,ntaski,
     |  handle_mpi_err
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
#include "mpif.h"
!
! Each task receives fkmh(:,lon0-1,lat0-3:lat1) from task iprev, 
!   and sends fkmh%data(:,lon1,lat0-3:lat1) to task inext.
! lon1==nlonp4 is sent to task mytidi==0, and lon0-1==0 is received 
!   from task ntaski-1
!
! Local:
      integer ::
     |  irstat(MPI_STATUS_SIZE)      ! mpi receive status
      integer :: j,jj,ier,len,iprev,inext,isend,irecv,nlats
      real,allocatable :: rcvbuf(:,:),sndbuf(:,:)
      real :: fmin,fmax
#ifdef VT
!     call vtsymdef(112, 'mp_bndlons_kmh','Communication',ier)
      call vtbegin(112,ier)
#endif
!
! Allocate send and receive buffers:
      nlats = lat1-lat0+1+3
      allocate(rcvbuf(nlevp1,nlats),stat=ier)
      if (ier /= 0) write(6,"('>>> mp_bndlons_kmh: error allocating',
     |  ' rcvbuf: nlevp1=',i3,' nlats=',i3)") nlevp1,nlats
      allocate(sndbuf(nlevp1,nlats),stat=ier)
      if (ier /= 0) write(6,"('>>> mp_bndlons_kmh: error allocating',
     |  ' sndbuf: nlevp1=',i3,' nlats=',i3)") nlevp1,nlats
      len = nlevp1*nlats
!
! Receive lon0-1 from iprev:
! (if lon0==1, receive lon0-1 from task with mytidi==ntask-1)
      iprev = itask_table(mytidi-1,mytidj)
      if (lon0==1) iprev = itask_table(ntaski-1,mytidj) 
      call mpi_irecv(rcvbuf,len,MPI_REAL8,iprev,1,MPI_COMM_WORLD,
     |  irecv,ier)
      if (ier /= 0) 
     |  call handle_mpi_err(ier,'mp_bndlons_kmh recv from iprev')
!
! Send lon1 to inext:
! (if lon1==nlonp4, send lon1 to task with mytidi==0)
      jj = 0
      do j=lat0-3,lat1
        jj = jj+1
        sndbuf(:,jj) = fkmh(:,lon1,j)
      enddo
      inext = itask_table(mytidi+1,mytidj)
      if (lon1==nlonp4) inext = itask_table(0,mytidj)

!     call fminmax(sndbuf(:,:),nlevp1*nlats,fmin,fmax)
!     write(6,"('Send lon1=',i2,' to inext=',i2,' sndbuf min,max=',
!    |  2e12.4)") lon1,inext,fmin,fmax

      call mpi_isend(sndbuf,len,MPI_REAL8,inext,1,MPI_COMM_WORLD,
     |  isend,ier)
      if (ier /= 0) 
     |  call handle_mpi_err(ier,'mp_bndlons_kmh send to inext')
!
! Wait for completions:
      call mpi_wait(isend,irstat,ier)
      if (ier /= 0) 
     |  call handle_mpi_err(ier,'mp_bndlons_kmh wait for send')
      call mpi_wait(irecv,irstat,ier)
      if (ier /= 0) 
     |  call handle_mpi_err(ier,'mp_bndlons_kmh wait for recv')
!
! Copy lon0-1 from rcvbuf:
      call fminmax(rcvbuf(:,:),nlevp1*nlats,fmin,fmax)
!     write(6,"('Recv lon0-1=',i2,' from iprev=',i2,' rcvbuf min,max=',
!    |  2e12.4)") lon0-1,iprev,fmin,fmax

      jj = 0
      do j=lat0-3,lat1
        jj = jj+1
        fkmh(:,lon0-1,j) = rcvbuf(:,jj)
      enddo
!
! Release local buffer space:
      deallocate(sndbuf)
      deallocate(rcvbuf)
#ifdef VT
!     call vtsymdef(112, 'mp_bndlons_kmh','Communication',ier)
      call vtend(112,ier)
#endif
      end subroutine mp_bndlons_kmh
!
! #endif MPI:
#endif
      end module hdif_module