SUBROUTINE Btrace_geopack (iyear, iday, ihour, imin, isec, &
                                  R_i, R_e, colat, aloct, pss, V_flag, &
                                  parmod, imax, jmax, ipdim, &
                                  vm, bmin, xmin, ymin, &
                                  exname, inname, n_bf)
       IMPLICIT NONE
       INTEGER, INTENT (IN) :: iyear, iday, ihour, imin, isec, &
                               imax, jmax, ipdim, n_bf
       REAL, INTENT (IN) :: R_i, R_e, pss, parmod(ipdim), V_flag
       REAL, DIMENSION (imax,jmax), INTENT (IN) :: colat, aloct
       REAL, DIMENSION (Imax,jmax), INTENT (OUT):: vm, bmin, xmin, ymin
       EXTERNAL :: exname, inname
!
!
       LOGICAL :: L_mgnp
       COMMON /stanislav_1/ L_mgnp
!
       REAL :: aa(10), sps, cps, bb(3), ps, cc(19)
       COMMON /geopack1/ aa,sps,cps,bb,ps,cc
!
       INTEGER :: iopt, i, j, m, n_open, n
       REAL :: xx(1000), yy(1000), zz(1000), btot(1000), pi, &
               dir, r0, rlim, gphi, gthet, grr, glat, xf, yf, zf, &
               bx, by, bz, bbx, bby, bbz, sum, h, xgsm, ygsm, zgsm
!
!
!
! * * field line tracing from ionospheric grid specified through colat
!     and aloct to compute flux tube volumes.
!
!
! * * paramerers:  colat, aloct   =  ionospheric grip point coordinates
!
!                  parmod         = input vector
!                   parmod(1)     = solar wind pressure p_dyn (nPa)
!                   parmod(2)     = Dst (nT)
!                   parmod(3)     = imf_by (nT)
!                   parmod(4)     = imf_bz (nT)
!                   parmod(5)     = g_2 index
!                   parmod(6)     = g_3 index
! (see Tsyganenko, Singer, and Kasper [2003] for def. of g_2 and g_3)
!
!                  ps             = dipole tilt angle
!
!                  latdim = idim  =  dimension of latitudinal points
!                  ltdim  = jdim  =  dimension of longitudinal points
!                  ipdim          =  dimension of parmod
!
!
!                  vm             =  resulting fluxtube volumes
!                  bmin           =  minimum b value along a field line
!                  xmin,ymin,zmin =  gsm position of bmin value.
!
!
      pi = acos(-1.000000)
!
      CALL Recalc (iyear,iday,ihour,imin,isec)
!
      ps = pss; sps = SIN(ps); cps = COS(ps)
!
!
      vm   = -1.0
      bmin = v_flag
      xmin = v_flag
      ymin = v_flag
!
      dir  =  1.0  ! start trace from northern hemisphere
      rlim = 40.0  ! limit the tracing region within 40 Re
      r0   =  1.0  ! landing point on Earth's surface
      iopt =    0  ! dummy unused parameter
!
      do_j_loop: DO j = 1, jmax, 1
!
!        For each j-line, start at the equat. boundary, and count open field lines:
!
         n_open = 0
!
      do_i_loop: DO i = imax, 1, -1
!
         gphi = aloct(i,j) 
         gthet = colat(i,j) 
         glat  = pi/2. - gthet
         grr   = R_i / R_e
! 
         xgsm = grr*sin(gthet)*cos(gphi)
         ygsm = grr*sin(gthet)*sin(gphi)
         zgsm = grr*cos(gthet)
!
         L_mgnp = .TRUE.   ! Flag if a point is inside magnetosphere
!                            false means either outside magnetopause
!                            or in the "boundary layer"
!
!        write (*,'(a,3I5)') '  tracing:',n_bf,i,j
         CALL Trace (xgsm, ygsm,zgsm, dir, rlim, r0, iopt, parmod,  &
                     exname, inname, xf, yf, zf, xx, yy, zz, m)
!
         IF (m>998 .OR. .not.L_mgnp .OR. SQRT(xf**2+yf**2+zf**2)>1.6) THEN !OPEN FIELD 
            n_open = n_open + 1
            IF (n_open > 4)  EXIT do_i_loop
!
         ELSE  ! closed field line, proceed with integration
!
            DO n=1,m   !   Get B-field values along a field line:
               CALL Exname (iopt,parmod,ps,xx(n),yy(n),zz(n),bx,by,bz)
               CALL Inname (xx(n), yy(n), zz(n), bbx, bby, bbz)
               btot(n) = SQRT ((bx+bbx)**2 + (by+bby)**2 + (bz+bbz)**2)
            END DO
!
            sum = 0.0
            DO n = 1, m-1   !  Integrate flux-tube volume:
               h = SQRT((xx(n+1)-xx(n))**2 + &
                        (yy(n+1)-yy(n))**2 + &
                        (zz(n+1)-zz(n))**2)
               sum = sum + 2.*h/(btot(n+1)+btot(n))
            END DO
            vm(i,j) = sum**(-2.0/3.0)
!
!           Find equatorial crossings (note: cannot use min B, because
!           B minimizes off the z=0 plane):
!
            do_find_eq_crossing: DO n = 1, m-1
               IF (zz(n) > 0.0 .AND. zz(n+1) < 0.0) THEN
                   xmin (i,j) = (xx(n)+xx(n+1))/2.0
                   ymin (i,j) = (yy(n)+yy(n+1))/2.0
                   EXIT do_find_eq_crossing
               END IF
            END DO do_find_eq_crossing
!
         END IF
!
      END DO do_i_loop
      END DO do_j_loop
!
      RETURN
      END SUBROUTINE Btrace_geopack