!
!  IGDRVS - modified to reorder YZ data. 
!  1st dimension of YZ() et al. set to 2
!
      SUBROUTINE IGDRVS(YZ,DYZDT)
      PARAMETER (NPTMAX=160,NZMAX=401)
      DIMENSION YZ(2,NPTMAX),DYZDT(2,NPTMAX)
      COMMON /NGAM/GAM(NPTMAX),CGAM(NPTMAX),NPTS
      COMMON /KFLAGB/KUCR,KFLAG
      COMMON /RHOB/BVFSQ(NZMAX),RHO(NZMAX)
      COMMON /ATMB/TDEGZ(NZMAX),TDEG(NZMAX)
     X            ,QZ(NZMAX),Q(NZMAX)
     X            ,UCRZ(NZMAX),UCR(NZMAX)
     X            ,UHWZ(NZMAX),UHW(NZMAX)
     &            ,NZTDEG,NZQ,NZUCR,NZUHW
      COMMON /IZTOPB/ZLAST,QB,QB1,QB4
     &              ,IZTOPT,IZTOPQ,IZTOPU,IZTOPHW,NUTTOP
      COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2,
     X            THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2,
     X            DELTH1,DELTH2,FACMIN,FACMAX,DELFAC
      COMMON /AIRPLANE/AIRPORTALT,RHOAIR,ACSPEEDI,TAN3
      DIMENSION Y(NPTMAX),Z(NPTMAX),DYDT(NPTMAX),DZDT(NPTMAX)
      SMALL = 1.0e-30
C
CC      PRINT 995, T,(YZ(1,I),I=1,6),(YZ(2,I),I=1,6)
CC   &                      ,(GAM(I),I=1,6),R2PI,NPTS
CC  995 FORMAT('T,YZ = ',25E13.4/'GAM,R2PI,NPTS = ',14E13.4,I5)

!  set up crosswind conditional so it's not in the loop below
      IF(KUCR.EQ.1 .AND. KFLAG.EQ.1) THEN
        FUCR=0.0
      ELSE
        FUCR=1.0
      ENDIF

      DO 10 I=1,NPTS
      Y(I)=YZ(1,I)
      Z(I)=YZ(2,I)
   10 CONTINUE
!  for the Jth vortex ...
      DO 30 J=1,NPTS
      SUMDY =0.
      SUMDZ =0.
!  compute crosswind and headwind at vortex altitude
      UCRK = FK(Z(J),UCRZ,UCR,NZUCR,IZTOPU)
      UHWK = FK(Z(J),UHWZ,UHW,NZUHW,IZTOPU)

!  compute velocity contribution from all vortices and their images
      DO 20 I=1,NPTS
!!!   IF(I.EQ.J) GO TO 20
      DELY=Y(J)-Y(I)
      DELZ=Z(J)-Z(I)
      SUMZ=Z(J)+Z(I)
      DSQ=DELY*DELY+DELZ*DELZ+SMALL
      DSI=DELY*DELY+SUMZ*SUMZ
CC      PRINT 993, J,I,DELY,DELZ,DSQ
CC  993 FORMAT(/' J,I,DELY,DELZ,DSQ = ',2I10,3E13.4)
      SUMDY=SUMDY+CGAM(I)*(SUMZ/DSI-DELZ/DSQ)
      SUMDZ=SUMDZ+CGAM(I)*(DELY/DSQ-DELY/DSI)
CC      PRINT 994, SUMDY,SUMDZ
CC  994 FORMAT(/' SUMDY,SUMDZ = ',2E13.4)
   20 CONTINUE
!  total advected velocity is given by ...
      DYDT(J) = ( R2PI*SUMDY + UCRK*FUCR ) * ( 1 - UHWK*ACSPEEDI )
      DZDT(J) = ( R2PI*SUMDZ - UHWK*TAN3 ) * ( 1 - UHWK*ACSPEEDI )
!     DZDT(J) = ( R2PI*SUMDZ - UHWK* 0   ) * ( 1 - UHWK*ACSPEEDI )
!!!
!!!      IF(J.EQ.1 .OR. J.EQ.2 .OR. J.EQ.5 .OR. J.EQ.6 .OR.
!!!     X                           J.EQ.9 .OR. J.EQ.10) THEN
!!!         UCRK = FK(Z(J),UCRZ,UCR,NZUCR,IZTOPU)
!!!         UHWK = FK(Z(J),UHWZ,UHW,NZUHW,IZTOPU)
!!!         IF(KUCR.EQ.1 .AND. KFLAG.EQ.1) UCRK = 0.
!!!CC         PRINT *, "J,Z(J),UCRK = ",J,Z(J),UCRK
!!!         DYDT(J)=-R2PI*SUMDY + UCRK
!!!      ELSEIF(J.EQ.3) THEN
!!!         DYDT(J) = DYDT(1)
!!!      ELSEIF(J.EQ.4) THEN
!!!         DYDT(J) = DYDT(2)
!!!      ELSEIF(J.EQ.7) THEN
!!!         DYDT(J) = DYDT(5)
!!!      ELSEIF(J.EQ.8) THEN
!!!         DYDT(J) = DYDT(6)
!!!      ELSEIF(J.EQ.11) THEN
!!!         DYDT(J) = DYDT(9)
!!!      ELSEIF(J.EQ.12) THEN
!!!         DYDT(J) = DYDT(10)
!!!      ENDIF
!!!      DZDT(J)=R2PI*SUMDZ
!!!

!wj   write(2,*)"hi joe: J,DYDT,DZDT =",J,DYDT(J),DZDT(J)

   30 CONTINUE
C
      DO 90 J=1,NPTS
      DYZDT(1,J)=DYDT(J)
      DYZDT(2,J)=DZDT(J)
!     write(6,*)"IGDRVS - J,Y,Z,CGAM =",J,Y(J),Z(J),CGAM(J)
   90 CONTINUE
      RETURN
      END
!
!  RKQC - this is a general-use routine with no common blocks; there
!  is no need to reveal details of the variable dimension structure.
!
!  added variable GAM to disregard zero circulation vortices from 
!  error calculation
!
      SUBROUTINE RKQC(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT
     &                                  ,GAM,GAMSMALL,ZMIN)
      PARAMETER (NPTMAX=160,NMAX=2*NPTMAX,FCOR=.0666666667,
     &    ONE=1.,SAFETY=0.9)
      DIMENSION Y(N),DYDX(N),YSCAL(N),YTEMP(NMAX),YSAV(NMAX),DYSAV(NMAX)
      DIMENSION GAM(N/2),ZMIN(N/2)
      PGROW=-0.10
      PSHRNK=-0.125
      ERRCON = (4./SAFETY)**(1./PGROW)
      HMIN = 0.01
      HMAX = 0.4
      HTRY = AMAX1(HMIN,HTRY)
      HTRY = AMIN1(HMAX,HTRY)
      XSAV=X
      do i=1,n/2
        y(i*2) = amax1( y(i*2) , zmin(i) )
        yscal(i*2) = amax1( yscal(i*2) , zmin(i) )
      enddo
      DO 11 I=1,N
        YSAV( I)=Y(   I)
        DYSAV(I)=DYDX(I)
11    CONTINUE
      H=HTRY
1     HH=0.5*H
CC      print 96
CC   96 format(' rk4old try1')
      CALL RK4OLD(YSAV,DYSAV,N,XSAV,HH,YTEMP,ZMIN)
      do i=1,n/2
        ytemp(i*2) = amax1( ytemp(i*2) , zmin(i) )
      enddo
CC      print 99
CC   99 format(' rk4old ok1')
      X=XSAV+HH
      CALL IGDRVS(YTEMP,DYDX)
      CALL RK4OLD(YTEMP,DYDX,N,X,HH,Y,ZMIN)
      do i=1,n/2
        y(i*2) = amax1( y(i*2) , zmin(i) )
      enddo
CC      print 98
CC   98 format(' rk4old ok2')
      X=XSAV+H
C
CCC      IF(X.EQ.XSAV)PAUSE 'Stepsize not significant in RKQC.'
      IF(X.EQ.XSAV)THEN
         WRITE(2,991) 
 991     FORMAT(/'Stepsize not significant in RKQC.')
         CLOSE(2)
CCC         CLOSE(3)
         STOP
      ENDIF
C
      CALL RK4OLD(YSAV,DYSAV,N,XSAV,H,YTEMP,ZMIN)
      do i=1,n/2
        ytemp(i*2) = amax1( ytemp(i*2) , zmin(i) )
      enddo
CC      print 97
CC   97 format(' rk4old ok3')
      ERRMAX=0.
      DO 12 I=1,N
        YTEMP(I)=Y(I)-YTEMP(I)
        I2 = (I+1)/2
        if(abs(GAM(I2)).ge.GAMSMALL) then
          ERRMAX=MAX(ERRMAX,ABS(YTEMP(I)/YSCAL(I)))
        endif
12    CONTINUE
      ERRMAX=ERRMAX/EPS
CC      PRINT 997, SAFETY,H,ERRMAX,PSHRNK
CC  997 FORMAT('SAFETY,H,ERRMAX,PSHRNK = ',4E15.5)
CC      PRINT *, Y
!JOE  IF(ERRMAX.GT.ONE) THEN
      IF(ERRMAX.GT.ONE .AND. H.GT.HMIN) THEN
        H=SAFETY*H*(ERRMAX**PSHRNK)
        GOTO 1
      ELSE
        HDID=H
        IF(ERRMAX.GT.ERRCON)THEN
          HNEXT=SAFETY*H*(ERRMAX**PGROW)
        ELSE
          HNEXT=amax1(4.*H,HMAX)
        ENDIF
      ENDIF
      DO 13 I=1,N
        Y(I)=Y(I)+YTEMP(I)*FCOR
13    CONTINUE
      do i=1,n/2
        y(i*2) = amax1( y(i*2) , zmin(i) )
      enddo
      RETURN
      END
C
!  RK4OLD - this is a general-use routine with no common blocks; there
!  is no need to reveal details of the variable dimension structure.
!
      SUBROUTINE RK4OLD(Y,DYDX,N,X,H,YOUT,ZMIN)
      PARAMETER (NPTMAX=160,NMAX=2*NPTMAX)
      DIMENSION Y(N),DYDX(N),YOUT(N),YT(NMAX),DYT(NMAX),DYM(NMAX)
     &         ,ZMIN(N/2)
      HH=H*0.5
      H6=H/6.
      XH=X+HH
      DO 11 I=1,N
        YT(I)=Y(I)+HH*DYDX(I)
11    CONTINUE
      do i=1,n/2
        yt(2*i) = amax1(yt(2*i),zmin(i))
      enddo
CC      print 91
CC   91 format(' igdrvs try1')
      CALL IGDRVS(YT,DYT)
CC      print 92
CC   92 format(' igdrvs ok1')
      DO 12 I=1,N
        YT(I)=Y(I)+HH*DYT(I)
12    CONTINUE
      do i=1,n/2
        yt(2*i) = amax1(yt(2*i),zmin(i))
      enddo
CC      print 93
CC   93 format(' igdrvs try2')
      CALL IGDRVS(YT,DYM)
CC      print 94
CC   94 format(' igdrvs ok2')
      DO 13 I=1,N
        YT(I)=Y(I)+H*DYM(I)
13    CONTINUE
      do i=1,n/2
        yt(2*i) = amax1(yt(2*i),zmin(i))
      enddo
      do I=1,N
        DYM(I)=DYT(I)+DYM(I)
      enddo
CC      print 95
CC   95 format(' igdrvs try3')
      CALL IGDRVS(YT,DYT)
CC      print 96
CC   96 format(' igdrvs ok3')
      DO 14 I=1,N
        YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I))
14    CONTINUE
      do i=1,n/2
        yout(2*i) = amax1(yout(2*i),zmin(i))
      enddo
      RETURN
      END
C
!  GTOGL - this is a general-use routine with no common blocks; there
!  is no need to reveal details of the variable dimension structure.
!  1D interpolation, that's all that's going on here..
!
      SUBROUTINE GTOGL(NX,X,Y,NXI,XI,YI)
      DIMENSION X(*),Y(*),XI(*),YI(*)
      INTEGER IPT,IXI,IX1,IX,IXP1
      REAL FPT,PT
C
      IPT=1
      DO 20 IXI=1,NXI
      FPT=0.
      PT=XI(IXI)
      IF(PT.LT.X(1) .OR. PT.GT.X(NX)) GO TO 19
C
      IX1=IPT
      DO 10 IX=IX1,NX
      IXP1=IX+1
      IF(PT.LT.X(IX) .OR. PT.GT.X(IXP1)) GO TO 10
      FPT=Y(IX) + (Y(IXP1)-Y(IX))*(PT-X(IX))/(X(IXP1)-X(IX))
      IPT=IX
      GO TO 19
   10 CONTINUE
C
   19 CONTINUE
      YI(IXI)=FPT
   20 CONTINUE
      RETURN
      END
C
      REAL FUNCTION ANGLE(DY,DZ)
C
      COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2,
     X            THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2,
     X            DELTH1,DELTH2,FACMIN,FACMAX,DELFAC
C
      HYP = SQRT(DY*DY+DZ*DZ)
C
      IF (DY .GE. 0. .AND. DZ .GE. 0.) THEN
         ANGLE = ASIN(DZ/HYP)
      ELSEIF (DY .LT. 0. .AND. DZ .GE. 0.) THEN
         ANGLE = PI - ASIN(DZ/HYP)
      ELSEIF (DY .LE. 0. .AND. DZ .LT. 0.) THEN
         ANGLE = PI + ASIN(-DZ/HYP)
      ELSEIF (DY .GT. 0. .AND. DZ .LT. 0.) THEN
         ANGLE = TWOPI - ASIN(-DZ/HYP)
      ENDIF
C
      RETURN
      END
C
      REAL FUNCTION THFACF(THIN)
C
      COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2,
     X            THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2,
     X            DELTH1,DELTH2,FACMIN,FACMAX,DELFAC
C
      TH = THIN
      IF (TH.LT.THBOT) TH = TH + TWOPI
      IF (TH.GT.THTOP) TH = TH - TWOPI
C
      IF(THMAX1.LE.TH .AND. TH.LE.THMAX2) THEN
         THFACF = FACMAX
      ELSEIF(THMAX2.LT.TH .AND. TH.LT.THMIN1) THEN
         THFACF = FACMIN + DELFAC*(THMIN1-TH)/DELTH1
      ELSEIF(THMIN1.LE.TH .AND. TH.LE.THMIN2) THEN
         THFACF = FACMIN
      ELSE
         THFACF = FACMIN + DELFAC*(TH-THMIN2)/DELTH2
      ENDIF
C
      RETURN
      END
C
      REAL FUNCTION FREDDG(TARG)
C
      COMMON /SARPB/C0,BZ,VZIN,R3,F0,VSB,VVSB,RVSB2,AFAC,DBFAC,
     X              DIMT,BRATIO,GAMAVB,ARGTB,RHOINIT,TDEGINIT
      COMMON /FREDB/C1S3,C2S3,C4S15,GAM0,DGAM0,TNDSAV,
     X              DGAM25,TG,TGG,TNDB,Z1B,Z2B
C
      EARG = -0.4 * TARG**C2S3
      EFAC = C4S15 / TARG**C1S3
      DGAMND = EFAC * EXP(EARG)
      FREDDG = GAM0 * VSB * DGAMND

      RETURN
      END
