subroutine viper30(
     &     y0, z0, b0, R0, Gamma_oo, ac_speed, i_hw,
     &     nt, tempz, theta, rho, N_sq,
     &     ncw, cwz, cw, cw_spln, nhw, hwz, hw, nq, qz, q,
     &     t_initial, t_final, n_steps,
     &     nt_model, time_model,
     &     y_model_p, z_model_p, cir_model_p,
     &     y_model_s, z_model_s, cir_model_s )

      implicit none

      real y0, z0, b0, R0, Gamma_oo, ac_speed, t_initial, t_final,
     &     tempz(nt), theta(nt), rho(nt), N_sq(nt), cwz(ncw), cw(ncw),
     &     cw_spln(ncw), hwz(nhw), hw(nhw), qz(nq), q(nq),
     &     y_model_p(*), z_model_p(*), cir_model_p(*), time_model(*),
     &     y_model_s(*), z_model_s(*), cir_model_s(*)
      integer i_hw, nt, ncw, nhw, nq, n_steps, nt_model

      real y_p, z_p, circ_p, y_s, z_s, circ_s, baro_circ, v0,
     &     t_scale, t_scale_inv, time, dt, gamma_frac,Gamma_0, z_ige,
     &     z_avg, hw_avg, q_avg, theta_avg, N_sq_avg, cw_p, cw_s,
     &     cw_fit, dcw, ddcw, edr_star, N_sq_star, sg_star,
     &     theta_0_inv, t_run_ige, sol(7), r(7), r_old(7)
      integer i, n, nrk, k_hw, k_cwp, k_cws, k_q, k_t, k_n, ierr

      include 'vpr30.h'

      integer, parameter :: nptmax=160
      real dgam, delgam
      common/dgamb/ dgam(nptmax), delgam

      equivalence(       y_p, sol(1) )
      equivalence(       z_p, sol(2) )
      equivalence(    circ_p, sol(3) )
      equivalence(       y_s, sol(4) )
      equivalence(       z_s, sol(5) )
      equivalence(    circ_s, sol(6) )
      equivalence( baro_circ, sol(7) )

c   *** Define characteristic parameters and time step

      v0 = Gamma_oo/(2.0*pi*b0)
      t_scale = b0/v0
      t_scale_inv = 1.0/t_scale

      gamma_frac = R_c**2/(R0**2+R_c**2)
      Gamma_0 = Gamma_oo*gamma_frac

      dt = (t_final-t_initial)/float(n_steps)

      z_ige = z_ige_fac*b0

c   *** Initialize the solution variables

      time = t_initial
      y_p = y0 - 0.5*b0
      y_s = y0 + 0.5*b0
      z_p = z0
      z_s = z0
      circ_p = -Gamma_0
      circ_s =  Gamma_0
      baro_circ = 0.0
      sg_star = 0.0
      r_old = 0.0
      k_hw=0; k_cwp=0; k_cws=0; k_q=0; k_t=0; k_n=0

      nt_model = 1
      time_model( 1) = time
        y_model_p(1) = y_p
        z_model_p(1) = z_p
      cir_model_p(1) = circ_p
        y_model_s(1) = y_s
        z_model_s(1) = z_s
      cir_model_s(1) = circ_s

c   *** Advance the solution in time using the Runge-Kutta method.

      do i=2, n_steps+1

         do nrk=1, nrk_max

            z_avg = 0.5*( z_p + z_s )

            ierr = 0
            call interp(   hw,   hwz, nhw, k_hw, z_avg,    hw_avg, ierr)
            call interp(   cw,   cwz, ncw, k_cwp,z_p,      cw_p,   ierr)
            call interp(   cw,   cwz, ncw, k_cws,z_s,      cw_s,   ierr)
            call interp(    q,    qz, nq , k_q , z_avg,     q_avg, ierr)
            call interp(theta, tempz, nt , k_t , z_avg, theta_avg, ierr)
            call interp( N_sq, tempz, nt , k_n , z_avg,  N_sq_avg, ierr)
            if( ierr .ne. 0 ) then
               print *, ' '
               print *, 'vortex outside range of environmental data'
               print *, 'z_p, z_s = ', z_p, z_s
               return
            end if

            if( ncw .ge. 3 ) then
               call splint( cwz, cw, cw_spln, ncw, z_avg,
     &                      cw_fit, dcw, ddcw )
               sg_star = ddcw*b0**2/v0
            end if

            edr_star = ((q_avg*b0)**one_third)/v0
            N_sq_star = N_sq_avg*t_scale**2
            if( i .eq. 2 .and. nrk .eq. 1 ) theta_0_inv = 1.0/theta_avg

            call rhs( time, y_p, z_p, circ_p, y_s, z_s, circ_s,
     &                baro_circ, ac_speed, i_hw, hw_avg, cw_p, cw_s,
     &                Gamma_0, Gamma_oo, theta_0_inv, t_scale,
     &                t_scale_inv, t_initial,
     &                theta_avg, edr_star, N_sq_star, sg_star, r )

            do n=1, 7
               sol(n) = sol(n) + 
     &                  ( wt_rk1(nrk)*r(n) + wt_rk2(nrk)*r_old(n) )*dt
            end do

            time = time + (wt_rk1(nrk)+wt_rk2(nrk))*dt

            r_old = r

         end do

         nt_model = nt_model + 1
         time_model( i) = time
           y_model_p(i) = y_p
           z_model_p(i) = z_p
         cir_model_p(i) = circ_p
           y_model_s(i) = y_s
           z_model_s(i) = z_s
         cir_model_s(i) = circ_s

         z_avg = 0.5*( z_p + z_s )
         if( z_avg .lt. z_ige ) goto 90

      end do
      return

c   *** Compute the IGE portion via the viper_ige routine.

90    continue
      dgam(1) = r(3)
      dgam(2) = r(6)
      t_run_ige = t_final - time
      call viper_ige(t_run_ige,dt,b0,v0,ac_speed,i_hw,
     &               nt_model,time_model,y_model_p,z_model_p,
     &               cir_model_p,y_model_s,z_model_s,cir_model_s)

      if( abs(time_model(nt_model)-t_final) .gt. 1.0e-4 ) then
         print *, ' '
         print *, 'viper_ige halted'
      end if

      return
      end

      subroutine rhs( time, y_p, z_p, circ_p, y_s, z_s, circ_s, 
     &                baro_circ, ac_speed, i_hw, hw_avg, cw_p, cw_s,
     &                Gamma_0, Gamma_oo, theta_0_inv, t_scale,
     &                t_scale_inv, t_initial,
     &                theta_avg, edr_star, N_sq_star, sg_star, r )

      implicit none

      include 'vpr30.h'

      real time, y_p, z_p, circ_p, y_s, z_s, circ_s, baro_circ,
     &     ac_speed, hw_avg, cw_p, cw_s, Gamma_0, Gamma_oo, theta_0_inv,
     &     t_scale, t_scale_inv, t_initial, theta_avg, edr_star,
     &     N_sq_star, sg_star, r(7)
      integer i_hw

      real t_star, hw_fac, S_p, S_s, gam_shear_p, gam_shear_s, sum1,
     &     tau_dif, tau_1, tau_local, decay_rate, t_s_local_p,
     &     t_s_local_s, circ_avg, d, ann_fac, gam_halo_p, gam_halo_s,
     &     gam_baro_p, gam_baro_s, delt_y, delt_z, b, cos_fac, sin_fac,
     &     y_cent, z_cent, vv, wv, gv(8), yv(8), zv(8), gam_halo(6),
     &     gam_baro(6), gam_shear(6), halo_part(6), shear_part(6)
      integer i, j, m, n

      save decay_rate, t_s_local_p, t_s_local_s

      t_star = (time-t_initial)*t_scale_inv

      hw_fac = 1.0 - hw_avg/ac_speed

      if( sg_star .gt. 0.0 ) then
         S_p = c_sg1
         S_s = c_sg2
      else
         S_p = c_sg2
         S_s = c_sg1
      end if
      
      gam_shear_p = -0.5*(circ_p + circ_s)
      shear_part = 0.0
      do n=1, 6
         m = n
         if( gam_shear_p .lt. 0.0 ) m = iflip(n)
         sum1 = 0.0
         do j=1, n_fit
            sum1 = sum1 + sg_fac(m,j)*t_star**float(j-1)
         end do
         shear_part(n) = sum1
      end do

      do n=1, 2
         halo_part(n) = halo_fac(n,1) + 
     &                  0.5*(halo_fac(n,2)-halo_fac(n,1))*
     &    ( 1.0 + tanh(2.5*(t_star-halo_fac(n,3))/halo_fac(n,4)) )
      end do
      halo_part(3) = 1.0 - halo_part(1) - halo_part(2)

      halo_part(4) = halo_part(3)
      halo_part(5) = halo_part(2)
      halo_part(6) = halo_part(1)
      
      tau_dif = tau - tau_min
      tau_1 = tau_dif - tau_rate*edr_star
      tau_local = tau_min + max(tau_1,0.0) + 2.0*N_sq_star

      if( t_star .le. t_switch ) then
         decay_rate = Gamma_0/(tau_local*t_scale)
         t_s_local_p = tau_local - abs(circ_p)/decay_rate*t_scale_inv
         t_s_local_s = tau_local - abs(circ_s)/decay_rate*t_scale_inv
         r(3) =  decay_rate
         r(6) = -decay_rate
      else
         r(3) = -circ_p/( (tau_local-t_s_local_p)*t_scale )
         r(6) = -circ_s/( (tau_local-t_s_local_s)*t_scale )
      end if

      d = 1.0
      if( t_star .lt. t_roll ) d = 0.0
      r(3) = ( 1.0 - d*S_p*sg_star )*r(3)
      r(6) = ( 1.0 + d*S_s*sg_star )*r(6)

      circ_avg = 0.5*( circ_s - circ_p )
      ann_fac = A_ann + B_ann*t_star + C_ann*t_star**2
      gam_halo_p = -(Gamma_oo-circ_avg)*ann_fac
      gam_halo_s = -gam_halo_p
      gam_baro_p =  baro_circ*ann_fac
      gam_baro_s = -gam_baro_p
      gam_shear_s = gam_shear_p
      do i=1, 3
         gam_halo( i  ) = gam_halo_p
         gam_halo( i+3) = gam_halo_s
         gam_baro( i  ) = gam_baro_p
         gam_baro( i+3) = gam_baro_s
         gam_shear(i  ) = gam_shear_p
         gam_shear(i+3) = gam_shear_s
      end do

c   *** The sine and cosine factors are multiplied by the vortex separation b.
c   *** It is useful to do this since halo vortex positions (y1f, etc)
c   *** are normalized by b.  Thus we effectively undo the normalization
c   *** and multiply by the sine or cosine at the same time.

      delt_y = y_s - y_p
      delt_z = z_s - z_p
      b = sqrt( delt_y**2 + delt_z**2 )
      cos_fac = delt_y
      sin_fac = delt_z

      y_cent = 0.5*( y_s + y_p )
      z_cent = 0.5*( z_s + z_p )
      
      do i=1, 6
         gv(i) =  halo_part(i)*gam_halo( i) + baro_part(i)*gam_baro(i) +
     &           shear_part(i)*gam_shear(i)
         yv(i) = yhv(i)*cos_fac - zhv(i)*sin_fac + y_cent
         zv(i) = zhv(i)*cos_fac + yhv(i)*sin_fac + z_cent
      end do

      gv(7) = circ_p
      yv(7) = y_p
      zv(7) = z_p
      gv(8) = circ_s
      yv(8) = y_s
      zv(8) = z_s

      r(1) = cw_p
      r(2) = 0.0
      do n=1, 8
         call induction( gv(n), yv(n), zv(n), y_p, z_p, vv, wv )
         r(1) = r(1) + vv
         r(2) = r(2) + wv
      end do

      r(4) = cw_s
      r(5) = 0.0
      do n=1, 8
         call induction( gv(n), yv(n), zv(n), y_s, z_s, vv, wv )
         r(4) = r(4) + vv
         r(5) = r(5) + wv
      end do

      r(7) = grav*b*( 1.0 - theta_avg*theta_0_inv )

c   *** Account for headwind effects (only used when headwinds are light)

      if( i_hw .eq. 1 ) then
         r(2) = r(2) - hw_avg*tan3
         r(5) = r(5) - hw_avg*tan3
         do n=1, 7
            r(n) = r(n)*hw_fac
         end do
      end if

      return
      end

      subroutine induction( gv, yv, zv, yp, zp, v_ind, w_ind )

      implicit none

      real gv, yv, zv, yp, zp, v_ind, w_ind

      real gf, z_lower_boundary, zb, delta_z_pri, delta_z_img, delta_y,
     &     del_y(2), del_z(2), gam_fac(2), r_sq_inv(2)
      real, parameter :: two_pi_inv=0.15915494309189533577
      integer n

      gf = gv*two_pi_inv

      z_lower_boundary = 0.0

c   *** zb is the location of the mirror image below lower boundary

      zb = 2.0*z_lower_boundary - zv

      delta_z_pri = zp - zv
      delta_z_img = zp - zb

c   *** Primary vortex and its image below the lower boundary

      delta_y  = yp - yv
      del_y(1) = delta_y
      del_z(1) = delta_z_pri
      del_y(2) = delta_y
      del_z(2) = delta_z_img
      gam_fac(1) =  gf
      gam_fac(2) = -gf

      v_ind = 0.0
      w_ind = 0.0
      do n=1, 2
         if( abs(del_y(n))+abs(del_z(n)) .lt. 1.0e-12 ) then
            r_sq_inv(n) = 0.0
         else
            r_sq_inv(n) = 1.0/( del_y(n)**2 + del_z(n)**2 )
         end if
         v_ind = v_ind - gam_fac(n)*del_z(n)*r_sq_inv(n)
         w_ind = w_ind + gam_fac(n)*del_y(n)*r_sq_inv(n)
      end do

      return
      end

      subroutine interp( data, z, nz, k_interp, z_interp, dat_interp,
     &                   ierr )

      implicit none

      real data(nz), z(nz), z_interp, dat_interp
      integer nz, k_interp, ierr

      real w0, w1
      integer k

      if( k_interp .le. 0 .or. k_interp .gt. nz ) k_interp = nz

      if( z(k_interp) .ge. z_interp ) then
         do k=k_interp-1, 1, -1
            if( z(k) .lt. z_interp ) goto 10
         end do
         ierr = ierr + 1
         return
10       k_interp = k
      end if

      if( z(k_interp+1) .lt. z_interp ) then
         do k=k_interp+2, nz
            if( z(k) .gt. z_interp ) goto 20
         end do
         ierr = ierr + 1
         return
20       k_interp = k - 1
      end if

      w0 = ( z(k_interp+1) - z_interp )/( z(k_interp+1) - z(k_interp) )
      w1 = 1.0 - w0

      dat_interp = w0*data(k_interp) + w1*data(k_interp+1)

      return
      end