! 
      module hdif_module
      use params_module,only: nlevp1,nlonp4,nlat
      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 btf: changed lon dim from nlonp4+1 to 0:nlonp4+1 
      real,dimension(nlevp1,0:nlonp4+1,-2:nlat) ::
     |  fnrh,  ! eddy viscosity
     |  fkmh   ! M/T
!
! 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,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
!
! 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 btf: 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
!
! 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
      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 ,o2_nm  ,ox_nm  ,n4s_nm ,noz_nm ,co_nm  ,
     |  co2_nm ,h2o_nm ,h2_nm ,hox_nm ,ch4_nm ,ar_nm  ,he_nm  ,nat_nm ,
     |  fkldt  ,fkldu  ,fkldv ,fkldo2 ,fkldox ,fkldn4s,fkldnoz,fkldco ,
     |  fkldco2,fkldh2o,fkldh2,fkldhox,fkldch4,fkldar ,fkldhe ,fkldnat, 
     |  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)
     |  o2_nm , ! molecular oxygen at time n-1 (mmr)
     |  ox_nm , ! atomic oxygen at time n-1 (mmr)
     |  n4s_nm, ! minor species follow (time-gcm only)...
     |  noz_nm, ! 
     |  co_nm , ! 
     |  co2_nm, ! 
     |  h2o_nm, ! 
     |  h2_nm , ! 
     |  hox_nm, ! 
     |  ch4_nm, ! 
     |  ar_nm , ! 
     |  he_nm , ! 
     |  nat_nm  ! 

!
! 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
     |  fkldox , ! o1 horizontal diffusion term
     |  fkldo2 , ! o2 horizontal diffusion term
     |  fkldn4s, ! minor species follow (time-gcm only)...
     |  fkldnoz, !
     |  fkldco , !
     |  fkldco2, !
     |  fkldh2o, !
     |  fkldh2 , !
     |  fkldhox, !
     |  fkldch4, !
     |  fkldar , !
     |  fkldhe , !
     |  fkldnat  !
!
! Local:
      integer :: k,i,ip1,ip1k,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
!
! avkmh = average fkmh:
! Allocation of fkmh in sub init_kmh_nrh (fields_module.F):
!     allocate(fkmh%data3d(nlevp1,nlonp4+1,-2:nlat),stat=ier)
!
      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
      enddo
!     if (lat >= lat0)
!    |  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 >= lat0)
!    |  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
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDU',' ',' ',fkldu(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! 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
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDV',' ',' ',fkldv(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! 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
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDT',' ',' ',fkldt(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! FOR PSI = O2
      call lsqdsq(o2_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            o2_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            o2_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
          fkldo2(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDO2',' ',' ',fkldo2(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! FOR PSI = OX
      call lsqdsq(ox_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            ox_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            ox_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
          fkldox(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDO1',' ',' ',fkldox(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! Minor species (time-gcm only):
!
! N4S
      call lsqdsq(n4s_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            n4s_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            n4s_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
          fkldn4s(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDN4S',' ',' ',fkldn4s(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! NOZ
      call lsqdsq(noz_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            noz_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            noz_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
          fkldnoz(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDNOZ',' ',' ',fkldnoz(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! CO
      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
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDCO',' ',' ',fkldco(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! CO2
      call lsqdsq(co2_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            co2_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            co2_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
          fkldco2(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDCO2',' ',' ',fkldco2(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! H2O:
      call lsqdsq(h2o_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            h2o_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            h2o_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
          fkldh2o(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDH2O',' ',' ',fkldh2o(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! H2
      call lsqdsq(h2_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            h2_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            h2_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
          fkldh2(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDH2',' ',' ',fkldh2(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! HOX:
      call lsqdsq(hox_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            hox_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            hox_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
          fkldhox(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDHOX',' ',' ',fkldhox(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! CH4:
      call lsqdsq(ch4_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            ch4_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            ch4_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
          fkldch4(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDCH4',' ',' ',fkldch4(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! AR:
      call lsqdsq(ar_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            ar_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            ar_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
          fkldar(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDAR',' ',' ',fkldar(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! HE:
      call lsqdsq(he_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            he_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            he_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
          fkldhe(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDHE',' ',' ',fkldhe(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
! NAT:
      call lsqdsq(nat_nm(lev0:lev1,lonbeg-1:lonend+1,lat+2),
     |            nat_nm(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            nat_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
          fkldnat(k,i,lat+1) = avkmh(k,i)*rhokmh(k,i)
        enddo
      enddo
!     if (lat >= lat0-1 .and. lat <= lat1-1)
!    |  call addfld('FKLDNAT',' ',' ',fkldnat(:,lon0:lon1,lat+1),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,lat+1)
!
      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,kldo2,kldox,kldn4s,kldnoz,
     |  kldco,kldco2,kldh2o,kldh2,kldhox,kldch4,kldar,kldhe,kldnat,
     |  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.
!
      use fields_module,only: f3d
#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,kldo2,kldox,kldn4s,kldnoz,
     |  kldco,kldco2,kldh2o,kldh2,kldhox,kldch4,kldar,kldhe,kldnat
!
! Local:
      integer :: i,nlevs
      real :: f(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2,16)
!
#ifdef MPI
!
      nlevs = lev1-lev0+1
      do i=1,16
        f(:,:,:,i) = f3d(i)%data(:,:,:)
      enddo

      call mp_bndlons_f3d(f,nlevs,lon0,lon1,lat0,lat1,16,0)

      do i=1,16
        f3d(i)%data(:,:,:) = f(:,:,:,i)
      enddo
#endif
      end subroutine hdif_bndlons
!-----------------------------------------------------------------------
      subroutine hdif3(cp,
     |  kldt   ,kldu   ,kldv  ,kldo2  ,kldox  ,kldn4s ,kldnoz ,kldco  ,
     |  kldco2 ,kldh2o ,kldh2 ,kldhox ,kldch4 ,kldar  ,kldhe  ,kldnat ,
     |  hdt    ,hdu    ,hdv   ,hdo2   ,hdox   ,hdn4s  ,hdnoz  ,hdco   ,
     |  hdco2  ,hdh2o  ,hdh2  ,hdhox  ,hdch4  ,hdar   ,hdhe   ,hdnat  ,
     |  lev0,lev1,lon0,lon1,lat)
!
! Calculate horizontal diffusion terms for t,u,v,o2,o1 and time-gcm minor 
! species 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,o2,o1 and minor species.
!
      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 (input from hdif2):
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in):: 
     |  kldt   ,kldu   ,kldv  ,kldo2  ,kldox  ,kldn4s ,kldnoz ,kldco ,
     |  kldco2 ,kldh2o ,kldh2 ,kldhox ,kldch4 ,kldar  ,kldhe  ,kldnat 
!
! 2d output for current latitude (output):
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) :: 
     |  hdt    ,hdu    ,hdv   ,hdo2   ,hdox   ,hdn4s  ,hdnoz  ,hdco   ,
     |  hdco2  ,hdh2o  ,hdh2  ,hdhox  ,hdch4  ,hdar   ,hdhe   ,hdnat
!
! Local:
      integer :: k,i,lonbeg,lonend
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  fnrh_inv, hdout, cpi 
      integer :: i0,i1
!
      i0=lon0 ; i1=lon1
!
! 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

!
! 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
!
! 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
!
! Make hdo2 from kldo2 (was npsdh):
      call lsqdsq(kldo2(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldo2(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldo2(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
          hdo2(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdo2(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!
! Make hdox from kldox (was npsdh2):
      call lsqdsq(kldox(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldox(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldox(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
          hdox(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdox(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!
! Make hdn4s from kldn4s (was npdhn4s):
      call lsqdsq(kldn4s(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldn4s(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldn4s(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
          hdn4s(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdn4s(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!
! Make hdnoz from kldnoz (was npdhnoz):
      call lsqdsq(kldnoz(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldnoz(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldnoz(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
          hdnoz(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdnoz(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!
! Make hdco from kldco (was npdhco):
      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
!
! Make hdco2 from kldco2 (was npdhco2):
      call lsqdsq(kldco2(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldco2(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldco2(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
          hdco2(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdco2(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!
! Make hdh2o from kldh2o (was npdhh2o):
      call lsqdsq(kldh2o(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldh2o(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldh2o(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
          hdh2o(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdh2o(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!
! Make hdh2 from kldh2 (was npdhh2):
      call lsqdsq(kldh2(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldh2(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldh2(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
          hdh2(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdh2(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!
! Make hdhox from kldhox (was npdhhox):
      call lsqdsq(kldhox(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldhox(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldhox(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
          hdhox(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdhox(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!
! Make hdch4 from kldch4 (was npdhch4):
      call lsqdsq(kldch4(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldch4(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldch4(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
          hdch4(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdch4(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!
! Make hdar from kldar (was npdhar):
      call lsqdsq(kldar(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldar(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldar(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
          hdar(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdar(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!
! Make hdhe from kldhe (was npdhhe):
      call lsqdsq(kldhe(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldhe(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldhe(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
          hdhe(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdhe(lev1,i) = 0.
      enddo ! i=lonbeg,lonend
!
! Make hdnat from kldnat (was npdhnat):
      call lsqdsq(kldnat(lev0:lev1,lonbeg-1:lonend+1,lat+1),
     |            kldnat(lev0:lev1,lonbeg-1:lonend+1,lat  ),
     |            kldnat(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
          hdnat(k,i) = hdout(k,i)*fnrh_inv(k,i)
        enddo ! k=lev0,lev1-1
        hdnat(lev1,i) = 0.
      enddo ! i=lonbeg,lonend

!     call addfld('HDT'  ,' ',' ',hdt  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDU'  ,' ',' ',hdu  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDV'  ,' ',' ',hdv  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDO2' ,' ',' ',hdo2 (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDOX' ,' ',' ',hdox (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDN4S',' ',' ',hdn4s(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDNOZ',' ',' ',hdnoz(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDCO' ,' ',' ',hdco (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDCO2',' ',' ',hdco2(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDH2O',' ',' ',hdh2o(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDH2' ,' ',' ',hdh2 (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDHOX',' ',' ',hdhox(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDCH4',' ',' ',hdch4(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDAR' ,' ',' ',hdar (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDHE' ,' ',' ',hdhe (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HDNAT',' ',' ',hdnat(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      end subroutine hdif3
!-----------------------------------------------------------------------
      subroutine hdif_periodic(hdt,hdu,hdv,hdo2,hdox,
     |  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,hdo2,hdox
!
! Local:
      real :: ftmp(lev0:lev1,lon0:lon1,lat0:lat1,5)
#ifdef MPI
      ftmp(:,:,:,1) = hdt (:,lon0:lon1,lat0:lat1)
      ftmp(:,:,:,2) = hdu (:,lon0:lon1,lat0:lat1)
      ftmp(:,:,:,3) = hdv (:,lon0:lon1,lat0:lat1)
      ftmp(:,:,:,4) = hdo2(:,lon0:lon1,lat0:lat1)
      ftmp(:,:,:,5) = hdox(:,lon0:lon1,lat0:lat1)

      call mp_periodic_f3d(ftmp,lev0,lev1,lon0,lon1,lat0,lat1,5)

      hdt (:,lon0:lon1,lat0:lat1) = ftmp(:,:,:,1)
      hdu (:,lon0:lon1,lat0:lat1) = ftmp(:,:,:,2)
      hdv (:,lon0:lon1,lat0:lat1) = ftmp(:,:,:,3)
      hdo2(:,lon0:lon1,lat0:lat1) = ftmp(:,:,:,4)
      hdox(:,lon0:lon1,lat0:lat1) = ftmp(:,:,:,5)
#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,
     |  mxlon,mpi_timing,time_bndlats_kmh,total_bndlats_kmh
      use timing_module,only: timer
      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,save :: sndbuf(:,:,:),rcvbuf(:,:,:)
      integer :: ier,nlons,len,jprev,jnext,ireqrecv,ireqsend
      real :: time0,time1
#ifdef VT
!     call vtsymdef(111, 'mp_bndlats_kmh','Communication',ier)
      call vtbegin(111,ier)
#endif
!
      if (mpi_timing) call timer(time0,time1,"mp_bndlats_kmh",0,0)
      len = mxlon*nlevp1*3
      if (.not.allocated(sndbuf)) then
        allocate(sndbuf(nlevp1,mxlon,3),stat=ier)
        if (ier /= 0) write(6,"('>>> mp_bndlats_kmh: error allocating',
     |    ' sndbuf: mxlon=',i3,' nlevp1=',i3)") mxlon,nlevp1
        sndbuf = 0.
!       write(6,"('hdif mp_bndlats_kmh: allocated sndbuf: mxlon=',
!    |    i3,' len=',i5)") mxlon,len
      endif
      if (.not.allocated(rcvbuf)) then
        allocate(rcvbuf(nlevp1,mxlon,3),stat=ier)
        if (ier /= 0) write(6,"('>>> mp_bndlats_kmh: error allocating',
     |    ' rcvbuf: mxlon=',i3,' nlevp1=',i3)") mxlon,nlevp1
        rcvbuf = 0.
!       write(6,"('hdif mp_bndlats_kmh: allocated rcvbuf: mxlon=',
!    |    i3,' len=',i5)") mxlon,len
      endif
      nlons = lon1-lon0+1
!
      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(:,1:nlons,:) = 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(:,1:nlons,:)
!
      if (mpi_timing) then
        call timer(time0,time1,"mp_bndlats_kmh",1,0)
        time_bndlats_kmh = time_bndlats_kmh+time1
      endif
#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,mxlat
      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,save :: 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:
      if (.not.allocated(rcvbuf)) then
        allocate(rcvbuf(nlevp1,mxlat+3),stat=ier)
        if (ier /= 0) write(6,"('>>> mp_bndlons_kmh: error allocating',
     |    ' rcvbuf: nlevp1=',i4,' mxlat+3=',i3)") nlevp1,mxlat+3
        rcvbuf = 0.
!       write(6,"('hdif mp_bndlons_kmh: allocated rcvbuf: nlevp1=',i4,
!    |    ' mxlat+3=',i4)") nlevp1,mxlat+3
      endif
      if (.not.allocated(sndbuf)) then
        allocate(sndbuf(nlevp1,mxlat+3),stat=ier)
        if (ier /= 0) write(6,"('>>> mp_bndlons_kmh: error allocating',
     |    ' sndbuf: nlevp1=',i4,' mxlat+3=',i4)") nlevp1,mxlat+3
        sndbuf = 0.
!       write(6,"('hdif mp_bndlons_kmh: allocated sndbuf: nlevp1=',i4,
!    |    ' mxlat+3=',i4)") nlevp1,mxlat+3
      endif
      len = nlevp1*(mxlat+3)
!
! 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
!
#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