pro mtraj_runge,newlat,newlon,curlat,curlon,dayni,err,$
    forward=forward,timestep=dt,$
    uout=uout,vout=vout,hout=hout,tout=tout,$
    dayfracs1=dayfracs1,u1=u1,v1=v1,q1=q1,t1=t1x1,$
    dayfracs2=dayfracs2,u2=u2,v2=v2,q2=q2,t2=t2x,$
    dir=dir
;+
; NAME:
;    mraj_runge
; PURPOSE:
;    given positions curlat,curlon we want the positions
;    at dayni which is at time t+dt (curlat and curlon are at time t)
;    This code does that computation or new positions using 4th order runge kutta
; CATEGORY:
;    trajectory
; CALLING SEQUENCE:
;    pro ftraj_runge,newlat,newlon,curlat,curlon,dayni,err,$
;    newtheta=newtheta,curtheta=curtheta,$
;    diabatic=diabatic,forward=forward,dt=dt,$
;    uout=uout,vout=vout,hout=hout,tout=tout,data_available=data_available
; INPUTS:
;    curlat = current latitude array (degrees)
;    curlon = current longitude array
; OPTIONAL INPUT PARAMETERS:
; KEYWORD PARAMETERS:
;    newtheta = new value of theta
; OUTPUTS:
;    newlat,newlon = new postion arrays
; OPTIONAL OUTPUT PARAMETERS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; REQUIRED ROUTINES:
; MODIFICATION HISTORY:
;    mrs
;    $Header: /science/trajectory/programs/traj_runge.pro,v 1.1 1995/03/08 19:04:23 morris Exp $
;    GM 6/3/93 -- fixed pole-crossing parcel bugs.
;-
       dtd=dt*86400.
       c1=!pi/180.
       c2=1./(6.37e6*c1)


       dlat1 = forward*dtd*vout*c2

       t1=curlat+dlat1*0.5 ; polar fixup
       n1=where(t1 gt 90.)
       if n1(0) ge 0 then t1(n1) = 180.-t1(n1)
       n2=where(t1 lt -90.)
       if n2(0) ge 0 then t1(n2)=-t1(n2)-180.
       ctemp = cos(t1*c1)

       dlon1=dlat1*0.
       n = where(ctemp ne 0.0)
       dlon1(n) = forward*dtd*uout(n)*c2/(ctemp(n))

       t2=dlon1*0.5+curlon
       if n1(0) ge 0 then t2(n1)=t2(n1)+180.
       if n2(0) ge 0 then t2(n2)=t2(n2)+180.
       t2 = t2 mod 360.

       ; dt is used instead of dtd because heating is in deg/day

        mtraj_fixvals,t1,t2

      mtraj_getdata,dayni-forward*dt*.5,t1,t2,uout,vout,tout,qout,$
      err,forward=forward,$
      dayfracs1=dayfracs1,u1=u1,v1=v1,q1=q1,t1=t1x1,$
      dayfracs2=dayfracs2,u2=u2,v2=v2,q2=q2,t2=t2x,dir=dir

    if err ne 0 then begin
       print,' ftraj_getdata error in ftraj_runge -1'
       return
       endif

       dlat2 = forward*dtd*vout*c2

       t1=curlat+dlat2*0.5
       n1=where(t1 gt 90.)
       if n1(0) ge 0 then t1(n1) = 180.-t1(n1)
       n2=where(t1 lt -90.)
       if n2(0) ge 0 then t1(n2)=-t1(n2)-180.
       ctemp = cos(t1*c1)
       dlon2=dlat2*0.

       n = where(ctemp ne 0.0)

       dlon2(n) = forward*dtd*uout(n)*c2/(ctemp(n))

       t2=dlon2*0.5+curlon
       if n1(0) ge 0 then t2(n1)=t2(n1)+180.
       if n2(0) ge 0 then t2(n2)=t2(n2)+180.
       t2 = t2 mod 360.


       mtraj_fixvals,t1,t2

    mtraj_getdata,dayni-forward*dt*.5,t1,t2,uout,vout,tout,qout,err,$
      forward=forward,$
      dayfracs1=dayfracs1,u1=u1,v1=v1,q1=q1,t1=t1x1,$
      dayfracs2=dayfracs2,u2=u2,v2=v2,q2=q2,t2=t2x,dir=dir

       dlat3 = forward*dtd*vout*c2

       t1=curlat+dlat3*.5
       n1=where(t1 gt 90.)
       if n1(0) ge 0 then t1(n1) = 180.-t1(n1)
       n2=where(t1 lt -90.)
       if n2(0) ge 0 then t1(n2)=-t1(n2)-180.
       ctemp = cos(t1*c1)
       dlon3=dlat3*0.
       n = where(ctemp ne 0.0)


       dlon3(n) = forward*dtd*uout(n)*c2/(ctemp(n))


       t1=curlat+dlat3
       n1=where(t1 gt 90.)
       if n1(0) ge 0 then t1(n1) = 180.-t1(n1)
       n2=where(t1 lt -90.)
       if n2(0) ge 0 then t1(n2) = -t1(n2)-180.

       t2=dlon3+curlon
       if n1(0) ge 0 then t2(n1)=t2(n1)+180.
       if n2(0) ge 0 then t2(n2)=t2(n2)+180.
       t2 = t2 mod 360.
       mtraj_fixvals,t1,t2

    mtraj_getdata,dayni,t1,t2,uout,vout,tout,qout,err,$
      forward=forward,$
      dayfracs1=dayfracs1,u1=u1,v1=v1,q1=q1,t1=t1x1,$
      dayfracs2=dayfracs2,u2=u2,v2=v2,q2=q2,t2=t2x,dir=dir


dlat4 = forward*dtd*vout*c2

       t1=curlat+dlat4*.5
       n1=where(t1 gt 90.)
       if n1(0) ge 0 then t1(n1) = 180.-t1(n1)
       n2=where(t1 lt -90.)
       if n2(0) ge 0 then t1(n2)=-t1(n2)-180.
       ctemp = cos(t1*c1)
       dlon4=dlat4*0.
       n = where(ctemp ne 0.0)

dlon4(n) = forward*dtd*uout(n)*c2/(ctemp(n))


;  *** Finish step ***

       newlat = curlat + (dlat1 + 2.*dlat2 + 2.*dlat3 + dlat4)/6.
       newlon = curlon + (dlon1 + 2.*dlon2 + 2.*dlon3 + dlon4)/6.
       mtraj_fixvals,newlat,newlon


     return
     end