C------------------------------------------------------------------
      SUBROUTINE RCM_FLTR(Mx,My,Mz,Gx,Gy,Gz,Gx_bx,Gy_bx,Gz_bx,Gx_by,
     &                    Gy_by,Gz_by,Gx_bz,Gy_bz,Gz_bz,Bx,By,Bz,Rr,Pp,
     &                    Np,Nt,Xlat,Xlon,Ftv,Brad,Bphi,Beq,Rrio,Ttio,
     &                    Ppio,Tmp_io1,Rcm_diptime,Tim1)
C------------------------------------------------------------------
      REAL Gx(*)
      REAL Gy(*)
      REAL Gz(*)
      REAL Gx_bx(*)
      REAL Gy_bx(*)
      REAL Gz_bx(*)
      REAL Gx_by(*)
      REAL Gy_by(*)
      REAL Gz_by(*)
      REAL Gx_bz(*)
      REAL Gy_bz(*)
      REAL Gz_bz(*)
      REAL Bx(Mx,My,Mz)
      REAL By(Mx,My,Mz)
      REAL Bz(Mx,My,Mz)
      REAL Rr(Mx,My,Mz)
      REAL Pp(Mx,My,Mz)
      REAL Xlat(Np,Nt)
      REAL Xlon(Np,Nt)
      REAL Ftv(Np,Nt)
      REAL Brad(Np,Nt)
      REAL Bphi(Np,Nt)
      REAL Beq(Np,Nt)
      REAL Rrio(Np,Nt)
      REAL Ttio(Np,Nt)
      REAL Ppio(Np,Nt)
      REAL Tmp_io1(Np,Nt)
      REAL eqx1(Np,Nt)
      REAL eqy1(Np,Nt)
      REAL eqz1(Np,Nt)
      REAL*8 Rcm_diptime
      CHARACTER*80 rec
      SAVE 
C
C..... the field line tracer, integrates flux tube volumes and other stuff
C      as needed by the RCM
C
C..... some constants
      pi = 4.0*ATAN(1.0)
      rad = pi/180.0
      deg = 1.0/rad
      riono = 1.02
      riono2 = (riono-0.005)**2
      dir = -0.01
      rflmax = 100.0
      max_nint = INT(rflmax/ABS(dir))
      rrmax2 = (20.0)**2
      dftv = 6372.0E3*ABS(dir)*1.0E-9
C
C..... convert real*8 time into year, month, ...
      CALL EPOCH1966(Rcm_diptime,iy,mo,id,ih,mi,se,1)
      WRITE(0,'(a,f15.4,5i5,f10.3)')'rcm_fltr:start:',Rcm_diptime,iy,mo,
     &                              id,ih,mi,se
C..... set up coordinate transforms
      CALL COTR_SET(iy,mo,id,ih,mi,se)
C..... loop over the ionosphere grid
C
C..... optimization opportunity:  by memorizing points with open fl we could
C      avoid tracing lines that are clearly in the polar cap.
C      (do later)
C
      DO ip = 1,Np
         lastcl = 0
         DO it = 1,Nt
C..... for open field lines these will stay negative
            Ftv(ip,it) = -1.001
            Brad(ip,it) = -1.001
            Bphi(ip,it) = -1.001
            Beq(ip,it) = -1.001
            Rrio(ip,it) = -1.001
            Ttio(ip,it) = -1.001
            Ppio(ip,it) = -1.001
            eqx1(ip,it) = -999.99
            eqy1(ip,it) = -999.99
            eqz1(ip,it) = -999.99
            Tmp_io1(ip,it) = 1.0
C
C.....  start point in SM coords
C
            xsm = riono*SIN(rad*Xlat(ip,it))*COS(rad*Xlon(ip,it))
            ysm = riono*SIN(rad*Xlat(ip,it))*SIN(rad*Xlon(ip,it))
            zsm = riono*COS(rad*Xlat(ip,it))
C.....  transform to MHD coords
            CALL COTR('sm ','mhd',xsm,ysm,zsm,xmhd,ymhd,zmhd)
            x1 = xmhd
            y1 = ymhd
            z1 = zmhd
C
C...... field line trace loop
            xlong = 0.0
            bb = 0.0
            brmax = 0.0
            xftv = 0.0
            axpp = 0.0
            axrr = 0.0
            axtt = 0.0
            nint = 0
 20         r2 = x1*x1+y1*y1+z1*z1
            IF(r2.LE.rrmax2)THEN
C
C..... in this case the field line is closed, done
               IF(r2.LT.riono2)THEN
                  Ftv(ip,it) = 1.0
                  phism = deg*ATAN2(y1s,x1s)
                  CALL DEGXYZ(1.0,phism,90.0,xmhd,ymhd,zmhd)
                  CALL COTR('mhd','sm ',xmhd,ymhd,zmhd,xsm,ysm,zsm)
                  CALL XYZDEG(xsm,ysm,zsm,rdum,pdum,tdum)
                  lastcl = 1
                  GOTO 40
               ENDIF
C
C...... get magnetic field values
               rbl1 = 4.5
               rbl2 = 6.5
               r2 = x1*x1+y1*y1+z1*z1
               r = SQRT(r2)
C
C...... dipole only
               IF(r.LE.rbl1)THEN
                  CALL COTR('mhd','sm ',x1,y1,z1,xsm,ysm,zsm)
                  xbbr = (30574.000)/(r2*r2*r)
                  xbbx = -3.0*xsm*zsm*xbbr
                  xbby = -3.0*ysm*zsm*xbbr
                  xbbz = (r2-3.0*zsm*zsm)*xbbr
                  CALL COTR('sm ','mhd',xbbx,xbby,xbbz,bbx,bby,bbz)
                  xrr = 0.0
                  xpp = 0.0
                  xtt = 0.0
C
C...... MHD only
               ELSEIF(r.GE.rbl2)THEN
                  iout = 0
                  CALL IPOL3ARP(Rr,Mx,My,Mz,xrr,Gx,Gy,Gz,x1,y1,z1,iout)
                  IF(iout.EQ.1)GOTO 40
                  iout = 0
                  CALL IPOL3ARP(Pp,Mx,My,Mz,xpp,Gx,Gy,Gz,x1,y1,z1,iout)
                  IF(iout.EQ.1)GOTO 40
                  IF(iout.EQ.1)GOTO 40
                  xtt = 72429.0*xpp/xrr/11600.0E3
                  iout = 0
                  CALL IPOL3AX(Bx,Mx,My,Mz,bbx,Gx_bx,Gy_bx,Gz_bx,x1,y1,
     &                         z1,iout)
                  IF(iout.EQ.1)GOTO 40
                  iout = 0
                  CALL IPOL3AY(By,Mx,My,Mz,bby,Gx_by,Gy_by,Gz_by,x1,y1,
     &                         z1,iout)
                  IF(iout.EQ.1)GOTO 40
                  iout = 0
                  CALL IPOL3AZ(Bz,Mx,My,Mz,bbz,Gx_bz,Gy_bz,Gz_bz,x1,y1,
     &                         z1,iout)
                  IF(iout.EQ.1)GOTO 40
C
C...... dipole-MHD blend
               ELSE
                  CALL COTR('mhd','sm ',x1,y1,z1,xsm,ysm,zsm)
                  xbbr = (30574.000)/(r2*r2*r)
                  xbbx = -3.0*xsm*zsm*xbbr
                  xbby = -3.0*ysm*zsm*xbbr
                  xbbz = (r2-3.0*zsm*zsm)*xbbr
                  CALL COTR('sm ','mhd',xbbx,xbby,xbbz,bbx2,bby2,bbz2)
                  s1 = MAX(0.0,MIN(1.0,(r-rbl1)/(rbl2-rbl1)))
                  s2 = 1.0-s1
                  iout = 0
                  CALL IPOL3AX(Bx,Mx,My,Mz,bbx1,Gx_bx,Gy_bx,Gz_bx,x1,y1,
     &                         z1,iout)
                  IF(iout.EQ.1)GOTO 40
                  iout = 0
                  CALL IPOL3AY(By,Mx,My,Mz,bby1,Gx_by,Gy_by,Gz_by,x1,y1,
     &                         z1,iout)
                  IF(iout.EQ.1)GOTO 40
                  iout = 0
                  CALL IPOL3AZ(Bz,Mx,My,Mz,bbz1,Gx_bz,Gy_bz,Gz_bz,x1,y1,
     &                         z1,iout)
                  IF(iout.EQ.1)GOTO 40
                  bbx = s1*bbx1+s2*bbx2
                  bby = s1*bby1+s2*bby2
                  bbz = s1*bbz1+s2*bbz2
                  iout = 0
                  CALL IPOL3ARP(Rr,Mx,My,Mz,xrr,Gx,Gy,Gz,x1,y1,z1,iout)
                  iout = 0
                  CALL IPOL3ARP(Pp,Mx,My,Mz,xpp,Gx,Gy,Gz,x1,y1,z1,iout)
                  IF(iout.EQ.1)GOTO 40
                  xtt = s1*72429.0*xpp/xrr/11600.0E3
                  xrr = s1*xrr
                  xpp = s1*xpp
               ENDIF
C
               bb = SQRT(bbx*bbx+bby*bby+bbz*bbz)
C
C..... bail out for zero field.  This should only happen close
C      to the boundary and should cause no harm
               IF(bb.GT.0.01)THEN
C
C.....  integrate 1/B
                  bbi = 1.0/bb
                  s = dir*bbi
                  xftv = xftv+bbi
                  nint = nint+1
C
C..... Field line is getting too long.  Either goes down the tail lobes
C      or is curled up.  In these cases the RCM would not care.
C      During nortward IMF this could also be a tail flank (LLBL) field
C      line.  Since these are convecting tailward there should be no
C      harm to the RCM.
C
                  IF(nint.LE.max_nint)THEN
C
                     axrr = axrr+xrr*bbi
                     axpp = axpp+xpp*bbi
                     axtt = axtt+xtt*bbi
                     x2 = x1+bbx*s
                     y2 = y1+bby*s
                     z2 = z1+bbz*s
                     xlong = xlong+SQRT((x2-x1)**2+(y2-y1)**2+(z2-z1)
     &                       **2)
C
C..... look for tip of the field line
                     IF(r.GT.brmax)THEN
                        x1s = x1
                        y1s = y1
                        z1s = z1
                        brmax = r
                        bbeq = bb
                        eqx1(ip,it) = x1s
                        eqy1(ip,it) = y1s
                        eqz1(ip,it) = z1s
                     ENDIF
C
                     x1 = x2
                     y1 = y2
                     z1 = z2
                     GOTO 20
                  ENDIF
               ENDIF
            ENDIF
C
C..... have a closed one
C
 40         IF(Ftv(ip,it).GT.0.0)THEN
               Ftv(ip,it) = ABS(xftv*dftv)
               axrr = axrr/xftv
               Rrio(ip,it) = axrr
               axpp = axpp/xftv
               Ppio(ip,it) = axpp
               axtt = axtt/xftv
               Ttio(ip,it) = axtt
               Brad(ip,it) = brmax
               Beq(ip,it) = bbeq
C
C..... need to convert the MLT into proper SM coords.
C
               phism = deg*ATAN2(y1s,x1s)
               CALL DEGXYZ(1.0,phism,90.0,xmhd,ymhd,zmhd)
               CALL COTR('mhd','sm ',xmhd,ymhd,zmhd,xsm,ysm,zsm)
               CALL XYZDEG(xsm,ysm,zsm,rdum,pdum,tdum)
               Bphi(ip,it) = pdum
            ENDIF
C
         ENDDO
      ENDDO
C
      END
C
C.... We keep 4 versions of ipol3a here so that each can
C     memorize it's own last indices for efficiency.
C
C--------------------------------------------------------------
      SUBROUTINE IPOL3AX(A,Nx,Ny,Nz,B,Gx,Gy,Gz,X,Y,Z,Iout)
C--------------------------------------------------------------
      REAL A(Nx,Ny,Nz),Gx(Nx),Gy(Ny),Gz(Nz)
      DATA ixl,iyl,izl/1,1,1/
      SAVE 
      ix = ixl
      IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
         ix = MIN0(ixl+1,Nx-1)
         IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
            ix = MAX0(ixl-1,1)
            IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
               DO i = 1,Nx-1
                  IF(X.GE.Gx(i).AND.X.LE.Gx(i+1))THEN
                     ix = i
                     GOTO 100
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 100  ixl = ix
      iy = iyl
      IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
         iy = MIN0(iyl+1,Ny-1)
         IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
            iy = MAX0(iyl-1,1)
            IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
               DO i = 1,Ny-1
                  IF(Y.GE.Gy(i).AND.Y.LE.Gy(i+1))THEN
                     iy = i
                     GOTO 200
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 200  iyl = iy
      iz = izl
      IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
         iz = MIN0(izl+1,Nz-1)
         IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
            iz = MAX0(izl-1,1)
            IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
               DO i = 1,Nz-1
                  IF(Z.GE.Gz(i).AND.Z.LE.Gz(i+1))THEN
                     iz = i
                     GOTO 300
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 300  izl = iz
      dxi = (X-Gx(ix))/(Gx(ix+1)-Gx(ix))
      dyi = (Y-Gy(iy))/(Gy(iy+1)-Gy(iy))
      dzi = (Z-Gz(iz))/(Gz(iz+1)-Gz(iz))
      x00 = A(ix,iy,iz)+dxi*(A(ix+1,iy,iz)-A(ix,iy,iz))
      x10 = A(ix,iy+1,iz)+dxi*(A(ix+1,iy+1,iz)-A(ix,iy+1,iz))
      x01 = A(ix,iy,iz+1)+dxi*(A(ix+1,iy,iz+1)-A(ix,iy,iz+1))
      x11 = A(ix,iy+1,iz+1)+dxi*(A(ix+1,iy+1,iz+1)-A(ix,iy+1,iz+1))
      y0 = x00+dyi*(x10-x00)
      y1 = x01+dyi*(x11-x01)
      B = y0+dzi*(y1-y0)
      END
C--------------------------------------------------------------
      SUBROUTINE IPOL3AY(A,Nx,Ny,Nz,B,Gx,Gy,Gz,X,Y,Z,Iout)
C--------------------------------------------------------------
      REAL A(Nx,Ny,Nz),Gx(Nx),Gy(Ny),Gz(Nz)
      DATA ixl,iyl,izl/1,1,1/
      SAVE 
      ix = ixl
      IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
         ix = MIN0(ixl+1,Nx-1)
         IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
            ix = MAX0(ixl-1,1)
            IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
               DO i = 1,Nx-1
                  IF(X.GE.Gx(i).AND.X.LE.Gx(i+1))THEN
                     ix = i
                     GOTO 100
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 100  ixl = ix
      iy = iyl
      IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
         iy = MIN0(iyl+1,Ny-1)
         IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
            iy = MAX0(iyl-1,1)
            IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
               DO i = 1,Ny-1
                  IF(Y.GE.Gy(i).AND.Y.LE.Gy(i+1))THEN
                     iy = i
                     GOTO 200
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 200  iyl = iy
      iz = izl
      IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
         iz = MIN0(izl+1,Nz-1)
         IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
            iz = MAX0(izl-1,1)
            IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
               DO i = 1,Nz-1
                  IF(Z.GE.Gz(i).AND.Z.LE.Gz(i+1))THEN
                     iz = i
                     GOTO 300
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 300  izl = iz
      dxi = (X-Gx(ix))/(Gx(ix+1)-Gx(ix))
      dyi = (Y-Gy(iy))/(Gy(iy+1)-Gy(iy))
      dzi = (Z-Gz(iz))/(Gz(iz+1)-Gz(iz))
      x00 = A(ix,iy,iz)+dxi*(A(ix+1,iy,iz)-A(ix,iy,iz))
      x10 = A(ix,iy+1,iz)+dxi*(A(ix+1,iy+1,iz)-A(ix,iy+1,iz))
      x01 = A(ix,iy,iz+1)+dxi*(A(ix+1,iy,iz+1)-A(ix,iy,iz+1))
      x11 = A(ix,iy+1,iz+1)+dxi*(A(ix+1,iy+1,iz+1)-A(ix,iy+1,iz+1))
      y0 = x00+dyi*(x10-x00)
      y1 = x01+dyi*(x11-x01)
      B = y0+dzi*(y1-y0)
      END
C--------------------------------------------------------------
      SUBROUTINE IPOL3AZ(A,Nx,Ny,Nz,B,Gx,Gy,Gz,X,Y,Z,Iout)
C--------------------------------------------------------------
      REAL A(Nx,Ny,Nz),Gx(Nx),Gy(Ny),Gz(Nz)
      DATA ixl,iyl,izl/1,1,1/
      SAVE 
      ix = ixl
      IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
         ix = MIN0(ixl+1,Nx-1)
         IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
            ix = MAX0(ixl-1,1)
            IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
               DO i = 1,Nx-1
                  IF(X.GE.Gx(i).AND.X.LE.Gx(i+1))THEN
                     ix = i
                     GOTO 100
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 100  ixl = ix
      iy = iyl
      IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
         iy = MIN0(iyl+1,Ny-1)
         IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
            iy = MAX0(iyl-1,1)
            IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
               DO i = 1,Ny-1
                  IF(Y.GE.Gy(i).AND.Y.LE.Gy(i+1))THEN
                     iy = i
                     GOTO 200
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 200  iyl = iy
      iz = izl
      IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
         iz = MIN0(izl+1,Nz-1)
         IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
            iz = MAX0(izl-1,1)
            IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
               DO i = 1,Nz-1
                  IF(Z.GE.Gz(i).AND.Z.LE.Gz(i+1))THEN
                     iz = i
                     GOTO 300
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 300  izl = iz
      dxi = (X-Gx(ix))/(Gx(ix+1)-Gx(ix))
      dyi = (Y-Gy(iy))/(Gy(iy+1)-Gy(iy))
      dzi = (Z-Gz(iz))/(Gz(iz+1)-Gz(iz))
      x00 = A(ix,iy,iz)+dxi*(A(ix+1,iy,iz)-A(ix,iy,iz))
      x10 = A(ix,iy+1,iz)+dxi*(A(ix+1,iy+1,iz)-A(ix,iy+1,iz))
      x01 = A(ix,iy,iz+1)+dxi*(A(ix+1,iy,iz+1)-A(ix,iy,iz+1))
      x11 = A(ix,iy+1,iz+1)+dxi*(A(ix+1,iy+1,iz+1)-A(ix,iy+1,iz+1))
      y0 = x00+dyi*(x10-x00)
      y1 = x01+dyi*(x11-x01)
      B = y0+dzi*(y1-y0)
      END
C--------------------------------------------------------------
      SUBROUTINE IPOL3ARP(A,Nx,Ny,Nz,B,Gx,Gy,Gz,X,Y,Z,Iout)
C--------------------------------------------------------------
      REAL A(Nx,Ny,Nz),Gx(Nx),Gy(Ny),Gz(Nz)
      DATA ixl,iyl,izl/1,1,1/
      SAVE 
      ix = ixl
      IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
         ix = MIN0(ixl+1,Nx-1)
         IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
            ix = MAX0(ixl-1,1)
            IF(X.LT.Gx(ix).OR.X.GT.Gx(ix+1))THEN
               DO i = 1,Nx-1
                  IF(X.GE.Gx(i).AND.X.LE.Gx(i+1))THEN
                     ix = i
                     GOTO 100
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 100  ixl = ix
      iy = iyl
      IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
         iy = MIN0(iyl+1,Ny-1)
         IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
            iy = MAX0(iyl-1,1)
            IF(Y.LT.Gy(iy).OR.Y.GT.Gy(iy+1))THEN
               DO i = 1,Ny-1
                  IF(Y.GE.Gy(i).AND.Y.LE.Gy(i+1))THEN
                     iy = i
                     GOTO 200
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 200  iyl = iy
      iz = izl
      IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
         iz = MIN0(izl+1,Nz-1)
         IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
            iz = MAX0(izl-1,1)
            IF(Z.LT.Gz(iz).OR.Z.GT.Gz(iz+1))THEN
               DO i = 1,Nz-1
                  IF(Z.GE.Gz(i).AND.Z.LE.Gz(i+1))THEN
                     iz = i
                     GOTO 300
                  ENDIF
               ENDDO
               Iout = 1
               RETURN
            ENDIF
         ENDIF
      ENDIF
 300  izl = iz
      dxi = (X-Gx(ix))/(Gx(ix+1)-Gx(ix))
      dyi = (Y-Gy(iy))/(Gy(iy+1)-Gy(iy))
      dzi = (Z-Gz(iz))/(Gz(iz+1)-Gz(iz))
      x00 = A(ix,iy,iz)+dxi*(A(ix+1,iy,iz)-A(ix,iy,iz))
      x10 = A(ix,iy+1,iz)+dxi*(A(ix+1,iy+1,iz)-A(ix,iy+1,iz))
      x01 = A(ix,iy,iz+1)+dxi*(A(ix+1,iy,iz+1)-A(ix,iy,iz+1))
      x11 = A(ix,iy+1,iz+1)+dxi*(A(ix+1,iy+1,iz+1)-A(ix,iy+1,iz+1))
      y0 = x00+dyi*(x10-x00)
      y1 = x01+dyi*(x11-x01)
      B = y0+dzi*(y1-y0)
      END
C------------------------------------------------------------------
      SUBROUTINE RCM_GETGRID(Mx,My,Mz,Gx,Gy,Gz,Gx_bx,Gy_bx,Gz_bx,Gx_by,
     &                       Gy_by,Gz_by,Gx_bz,Gy_bz,Gz_bz)
C------------------------------------------------------------------
      REAL Gx(*)
      REAL Gy(*)
      REAL Gz(*)
      REAL Gx_bx(*)
      REAL Gy_bx(*)
      REAL Gz_bx(*)
      REAL Gx_by(*)
      REAL Gy_by(*)
      REAL Gz_by(*)
      REAL Gx_bz(*)
      REAL Gy_bz(*)
      REAL Gz_bz(*)
      SAVE 
C...... centered grid
      DO i = 1,Mx+1
         CALL GRIDX(i,fff,fdd)
         Gx(i) = fff
         WRITE(0,*)'grid x ',i,Gx(i)
      ENDDO
      DO i = 1,My+1
         CALL GRIDY(i,fff,fdd)
         Gy(i) = fff
         WRITE(0,*)'grid y ',i,Gy(i)
      ENDDO
      DO i = 1,Mz+1
         CALL GRIDZ(i,fff,fdd)
         Gz(i) = fff
         WRITE(0,*)'grid z ',i,Gz(i)
      ENDDO
C...... B grid is staggered to the centers of cell faces
      DO i = 1,Mx
         CALL GRIDX(i+1,ff1,fdd)
         CALL GRIDX(i,ff2,fdd)
         Gx_bx(i) = 0.5*(ff1+ff2)
C      write(0,*)'grid gB_bA ',i,gB_bA(i)
      ENDDO
      DO i = 1,My
         Gy_bx(i) = Gy(i)
      ENDDO
      DO i = 1,Mz
         Gz_bx(i) = Gz(i)
      ENDDO
      DO i = 1,Mx
         Gx_by(i) = Gx(i)
      ENDDO
      DO i = 1,My
         CALL GRIDY(i+1,ff1,fdd)
         CALL GRIDY(i,ff2,fdd)
         Gy_by(i) = 0.5*(ff1+ff2)
C      write(0,*)'grid gB_bA ',i,gB_bA(i)
      ENDDO
      DO i = 1,Mz
         Gz_by(i) = Gz(i)
      ENDDO
      DO i = 1,Mx
         Gx_bz(i) = Gx(i)
      ENDDO
      DO i = 1,My
         Gy_bz(i) = Gy(i)
      ENDDO
      DO i = 1,Mz
         CALL GRIDZ(i+1,ff1,fdd)
         CALL GRIDZ(i,ff2,fdd)
         Gz_bz(i) = 0.5*(ff1+ff2)
C      write(0,*)'grid gB_bA ',i,gB_bA(i)
      ENDDO
      END