#include "dims.h"
      SUBROUTINE BDZDVB(PHI,DVB,H,J,ANS)
      use cons_module,only: len1,kmax,kmaxm1
      implicit none
C     ****
C     ****      EVALUATES  ANS = (BZ*D/(H*DZ)+DIVB)*PHI
C     ****
#include "params.h"
      real :: bxm,bx,bxp,by,byp,bz,bzp,bmod,bmodp
      COMMON/MAGFLD/BXM(76,2),BX(76,36),BXP(76,4),
     1  BY(76,36),BYP(76,4),BZ(76,36),BZP(76,4),
     2  BMOD(76,36),BMODP(76,2)
!
! Args:
      real,intent(in) :: PHI(ZIMXP,ZKMXP),DVB(ZIMXP),H(ZIMXP,ZKMXP)
      real,intent(out) :: ANS(ZIMXP,ZKMXP)
      integer,intent(in) :: j
!
! Local:
      integer :: k,i
!
      DO 1 K=2,KMAXM1
        DO 2 I=1,LEN1
          ANS(I,K) = BZ(I,J)/(2.*H(I,K)*dz)*(PHI(I,K+1)-PHI(I,K-1))
     A     +DVB(I)*PHI(I,K)
    2   CONTINUE
    1 CONTINUE

!     do i=1,len1
!       write(6,"(/,'bdzdvb: lat=',i3,' i=',i3,' dz=',f7.3,' dvb(i)=',
!    |    e12.4,' bz(i,lat)=', e12.4,' h(i,:)=',/,(6e12.4))") 
!    |    j,i,dz,dvb(i),bz(i,j),h(i,:)
!       write(6,"('phi(i,:)=',/,(6e12.4))") phi(i,:) 
!       write(6,"('ans(i,:)=',/,(6e12.4))") ans(i,:)
!     enddo
C     ****
C     ****      UPPER AND LOWER BOUNDARIES
C     ****
      DO 3 I=1,LEN1
        ANS(I,KMAX) = BZ(I,J)/(H(I,KMAX)*dz)*(PHI(I,KMAX)-
     A   PHI(I,KMAXM1))+DVB(I)*PHI(I,KMAX)
        ANS(I,1) = BZ(I,J)/(H(I,1)*dz)*(PHI(I,2)-
     A   PHI(I,1))+DVB(I)*PHI(I,1)
    3 CONTINUE
      RETURN
      END
C