#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(zimxp,2),BX(zimxp,zjmx),BXP(zimxp,4),
     1  BY(zimxp,zjmx),BYP(zimxp,4),BZ(zimxp,zjmx),BZP(zimxp,4),
     2  BMOD(zimxp,zjmx),BMODP(zimxp,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
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