c********************************************************************
c
c                     CINTI.FOR
c
c   ROUTINE TO INTERPOLATE THE MEASURED VELOCITY DISTRIBUTION
c 
c********************************************************************
      SUBROUTINE CINTI (F,XZ,YZ,ZZ,MVZ,MAZZ,MEZ)

      IMPLICIT REAL*8 (A-H,O-Z)
      REAL PI,PIH,UMRF
      DIMENSION F(33,9,9),S(4),A(3,3,3),X(33),Y(9),Z(9),VDX(3,3,3),
     1IK(3),JK(3),KK(3),AC2(2,2),BC2(2),AC4V(2,2,2,2),AC4(4,4),
     2BC4V(2,2),BC4(4),AC8V(2,2,2,2,2,2),AC8(8,8),BC8(8),FV2(3),FV3(3),
     3FV4(2,2,2),AQA(3,3,3,16,4,4),AQM(3,3,3),XZ(33),YZ(9),ZZ(9)
      COMMON /KONST/ PI,UMRF,LH1,LH2,LH3,PIH
      COMMON/A7/MEN,MAZ,MEL,LC1,LC2,LC3
      LOGICAL IPL
      
	write(6,*) '### entry of cinti'
	 write(6,*)'f',( f(k,5,3),k=17,27)

	read(5,*)

	it=lc1
	jt=lc2
	kt=lc3

      DO 101 I=1,it
  101 X(I) = XZ(I)
      DO 102 J=1,jt
  102 Y(J) = YZ(J)
      DO 103 K=1,kt
  103 Z(K) = ZZ(K)
      MV = MVZ
      MA = MAZZ
      ME = MEZ
      IK(1)=1
      IK(2)=0
      IK(3)=0
      JK(1)=0
      JK(2)=1
      JK(3)=0
      KK(1)=0
      KK(2)=0
      KK(3)=1
	 
c      WRITE(6,1113) (((F(I,J,K),I,J,K,I=1,LC1),J=1,LC2),K=1,LC3)
c 1113 FORMAT(4H  F=D20.8,4H  I=I5,4H  J=I5,4H  K=I5)
c      WRITE(6,1114) LC1,LC2,LC3
c 1114 FORMAT(6H  LC1=I5,6H  LC2=I5,6H  LC3=I5)
c      WRITE(6,1115)   (   X(I),I,I=1,LC1)
c 1115 FORMAT(7H  DGES=D20.8,4H  I=I5)
c      WRITE(6,1116) (  Y(J),J,J=1,LC2)
c 1116 FORMAT(6H  DAZ=D20.8,4H  J=I5)
c      WRITE(6,1117) (  Z(K),K,K=1,LC3)
c 1117 FORMAT(6H  DEL=D20.8,4H  K=I5)
C 
      DO 601 IQ=1,LH1
      I=2*IQ+1
      DO 601 JQ=1,LH2
      J=2*JQ+1
      DO 601 KQ=1,LH3
      K=2*KQ+1
      A(1,1,1) = F(I-2,J-2,K-2)
      VDX(1,1,1) =0.
      VDX(2,1,1)= X(I-1)-X(I-2)
      VDX(3,1,1) = X(I)-X(I-2)
      VDX(1,2,1) = Y(J-1)-Y(J-2)
      VDX(1,3,1) = Y(J)-Y(J-2)
      VDX(1,1,2) = Z(K-1) - Z(K-2)
      VDX(1,1,3) = Z(K) - Z(K-2)
C 
      DO 1 I1=1,3
      DO 2 L1=1,2
      IA=L1*IK(I1)
      JA=L1*JK(I1)
      KA=L1*KK(I1)
      DO 3 L2=1,2
    3 AC2(L1,L2)=(VDX(1+IA,1+JA,1+KA))**L2
    2 BC2(L1)=(F(I-2+IA,J-2+JA,K-2+KA)-F(I-2,J-2,K-2))
C
C     Write(6,*) 'call DGELG 1'
      CALL DGELG(BC2,AC2,2,1,1.0D-15,IER)
C     Write(6,*) 'ier 1', IER
C
      DO 4 L1=1,2
      IB=L1*IK(I1)+1
      JB=L1*JK(I1)+1
      KB=L1*KK(I1)+1
      IG = IB+JB+KB-3
    4 A(IB,JB,KB)=BC2(IG)
    1 CONTINUE
C
      DO 5 I1=1,3
      L1=1
      DO6 L12=1,2
      DO6 L11=1,2
      JO=L11*IK(I1)
      KO=L11*JK(I1)
      IO=L11*KK(I1)
      IT=L12*IK(I1)
      JT=L12*JK(I1)
      KT=L12*KK(I1)
      L2=1
      DO 7 L22=1,2
      DO7 L21=1,2
      AC4V(L11,L12,L21,L22)=(VDX(1+IO,1+JO,1+KO))**L21*(VDX(1+IT,1+JT,
     11+KT))**L22
      AC4(L1,L2)=AC4V(L11,L12,L21,L22)
    7 L2=L2+1
      FV1=F(I-2,J-2,K-2)
      FV2(1)=F(I-2+IO,J-2+JO,K-2+KO)-F(I-2,J-2,K-2)
      FV2(2)=F(I-2+IT,J-2+JT,K-2+KT)-F(I-2,J-2,K-2)
      BC4V(L11,L12)=F(I-2+IO+IT,J-2+JO+JT,K-2+KO+KT)-FV2(1)-FV2(2)-FV1
      BC4(L1)=BC4V(L11,L12)
    6 L1=L1+1
C 
C     Write(6,*) 'call DGELG 2'
      CALL DGELG(BC4,AC4,4,1,1.0D-15,IER)   
C     Write(6,*) 'ier 2', IER
C

      LZ=1
      DO 8 L12=1,2
      DO 8 L11=1,2
      ID= L11*KK(I1)+L12*IK(I1)+1
      JD= L11*IK(I1)+L12*JK(I1)+1
      KD= L11*JK(I1)+L12*KK(I1)+1
      A(ID,JD,KD)=BC4(LZ)
    8 LZ = LZ+1
    5 CONTINUE
      L = 1
      DO 9 L1=1,2
      DO 9 L2=1,2
      DO 9 L3=1,2
      FV2(1)=F(I-2+L1,J-2,K-2)-F(I-2,J-2,K-2)
      FV2(2)=F(I-2,J-2+L2,K-2)-F(I-2,J-2,K-2)
      FV2(3)=F(I-2,J-2,K-2+L3)-F(I-2,J-2,K-2)
      FV3(1)=F(I-2+L1,J-2+L2,K-2)-FV2(1)-FV2(2)-FV1
      FV3(2)=F(I-2+L1,J-2,K-2+L3)-FV2(1)-FV2(3)-FV1
      FV3(3)=F(I-2,J-2+L2,K-2+L3)-FV2(2)-FV2(3)-FV1
      FV4(L1,L2,L3)=F(I-2+L1,J-2+L2,K-2+L3)-FV3(1)-FV3(2)-FV3(3)-FV2(1)-
     1FV2(2)-FV2(3)-FV1
      BC8(L)=FV4(L1,L2,L3)
    9 L=L+1
      L1=1
      DO 10 L11=1,2
      DO 10 L12=1,2
      DO 10 L13=1,2
      L2=1
      DO 11 L21=1,2
      DO 11 L22=1,2
      DO 11 L23=1,2
      AC8V(L11,L12,L13,L21,L22,L23)=(VDX(1+L11,1,1))**L21*(VDX(1,1+L12,1
     1))**L22*(VDX(1,1,1+L13))**L23
      AC8(L1,L2) = AC8V(L11,L12,L13,L21,L22,L23)
   11 L2=L2+1
   10 L1=L1+1
C 
C     Write(6,*) 'call DGELG 3' 
      CALL DGELG(BC8,AC8,8,1,1.0D-15,IER)
C     Write(6,*) 'ier 3', IER
C
      LZ=1
      DO 13 I13=2,3
      DO 13 J13=2,3
      DO 13 K13=2,3
      A(I13,J13,K13)=BC8(LZ)
   13 LZ=LZ+1
      DO 501 IQ1=1,3
      DO 501 JQ1=1,3
      DO 501 KQ1=1,3
      AQA(IQ1,JQ1,KQ1,IQ,JQ,KQ)=A(IQ1,JQ1,KQ1)
C      Write(6,*) 'cinti 501', AQA(IQ1,JQ1,KQ1,IQ,JQ,KQ)
  501 CONTINUE
C
  601 CONTINUE
      IF(MV.EQ.1.OR.MV.EQ.LC1) GO TO 270
      IF(MA.EQ.1.OR.MA.EQ.LC2) GO TO 270
      IF(ME.EQ.1.OR.ME.EQ.LC3) GO TO 270
      I = MV+1
      J = MA+1
      K = ME+1
      A(1,1,1) = F(I-2,J-2,K-2)
      VDX(1,1,1) =0.
      VDX(2,1,1)= X(I-1)-X(I-2)
      VDX(3,1,1) = X(I)-X(I-2)
      VDX(1,2,1) = Y(J-1)-Y(J-2)
      VDX(1,3,1) = Y(J)-Y(J-2)
      VDX(1,1,2) = Z(K-1) - Z(K-2)
      VDX(1,1,3) = Z(K) - Z(K-2)
      DO 701 I1=1,3
      DO 702 L1=1,2
      IA=L1*IK(I1)
      JA=L1*JK(I1)
      KA=L1*KK(I1)
      DO703 L2=1,2
  703 AC2(L1,L2)=(VDX(1+IA,1+JA,1+KA))**L2
  702 BC2(L1)=(F(I-2+IA,J-2+JA,K-2+KA)-F(I-2,J-2,K-2))
C 
C     Write(6,*) 'call DGELG 4' 
      CALL DGELG(BC2,AC2,2,1,1.0D-15,IER)
C     Write(6,*) 'ier 4', IER
C
      DO 704 L1=1,2
      IB=L1*IK(I1)+1
      JB=L1*JK(I1)+1
      KB=L1*KK(I1)+1
      IG = IB+JB+KB-3
  704 A(IB,JB,KB)=BC2(IG)
  701 CONTINUE
      DO 705 I1=1,3
      L1=1
      DO706 L12=1,2
      DO706 L11=1,2
      JO=L11*IK(I1)
      KO=L11*JK(I1)
      IO=L11*KK(I1)
      IT=L12*IK(I1)
      JT=L12*JK(I1)
      KT=L12*KK(I1)
      L2=1
      DO 707 L22=1,2
      DO707 L21=1,2
      AC4V(L11,L12,L21,L22)=(VDX(1+IO,1+JO,1+KO))**L21*(VDX(1+IT,1+JT,
     11+KT))**L22
      AC4(L1,L2)=AC4V(L11,L12,L21,L22)
  707 L2=L2+1
      FV1=F(I-2,J-2,K-2)
      FV2(1)=F(I-2+IO,J-2+JO,K-2+KO)-F(I-2,J-2,K-2)
      FV2(2)=F(I-2+IT,J-2+JT,K-2+KT)-F(I-2,J-2,K-2)
      BC4V(L11,L12)=F(I-2+IO+IT,J-2+JO+JT,K-2+KO+KT)-FV2(1)-FV2(2)-FV1
      BC4(L1)=BC4V(L11,L12)
  706 L1=L1+1
C
C     Write(6,*) 'call DGELG 5'
      CALL DGELG(BC4,AC4,4,1,1.0D-15,IER)
C     Write(6,*) 'ier 5', IER
C
      LZ=1
      DO 708 L12=1,2
      DO 708 L11=1,2
      ID= L11*KK(I1)+L12*IK(I1)+1
      JD= L11*IK(I1)+L12*JK(I1)+1
      KD= L11*JK(I1)+L12*KK(I1)+1
      A(ID,JD,KD)=BC4(LZ)
  708 LZ = LZ+1
  705 CONTINUE
      L = 1
      DO 709 L1=1,2
      DO 709 L2=1,2
      DO 709 L3=1,2
      FV2(1)=F(I-2+L1,J-2,K-2)-F(I-2,J-2,K-2)
      FV2(2)=F(I-2,J-2+L2,K-2)-F(I-2,J-2,K-2)
      FV2(3)=F(I-2,J-2,K-2+L3)-F(I-2,J-2,K-2)
      FV3(1)=F(I-2+L1,J-2+L2,K-2)-FV2(1)-FV2(2)-FV1
      FV3(2)=F(I-2+L1,J-2,K-2+L3)-FV2(1)-FV2(3)-FV1
      FV3(3)=F(I-2,J-2+L2,K-2+L3)-FV2(2)-FV2(3)-FV1
      FV4(L1,L2,L3)=F(I-2+L1,J-2+L2,K-2+L3)-FV3(1)-FV3(2)-FV3(3)-FV2(1)-
     1FV2(2)-FV2(3)-FV1
      BC8(L)=FV4(L1,L2,L3)
  709 L=L+1
      L1=1
      DO710 L11=1,2
      DO710 L12=1,2
      DO710 L13=1,2
      L2=1
      DO711 L21=1,2
      DO711 L22=1,2
      DO711 L23=1,2
      AC8V(L11,L12,L13,L21,L22,L23)=(VDX(1+L11,1,1))**L21*(VDX(1,1+L12,1
     1))**L22*(VDX(1,1,1+L13))**L23
      AC8(L1,L2) = AC8V(L11,L12,L13,L21,L22,L23)
  711 L2=L2+1
  710 L1=L1+1
C
C     Write(6,*) 'call DGELG 6'
      CALL DGELG(BC8,AC8,8,1,1.0D-15,IER)
C     Write(6,*) 'ier 6', IER
C
      LZ=1
      DO713 I=2,3
      DO713 J=2,3
      DO713 K=2,3
      A(I,J,K)=BC8(LZ)
  713 LZ=LZ+1
      DO 505 I=1,3
      DO 505 J=1,3
      DO 505 K=1,3
      AQM(I,J,K)=A(I,J,K)
  505 CONTINUE
C
  270 RETURN

C***************************************************************
C        ENTRY IN SUBROUTINE CINTI 
C***************************************************************
	ENTRY CINT(S,IPL,IERR)
  130 XS=S(1)
      YS=S(2)
      ZS=S(3)
      IF( IPL ) GO TO 1270
      IF(MV.EQ.1.OR.MV.EQ.LC1) GO TO 1270
      IF(MA.EQ.1.OR.MA.EQ.LC2) GO TO 1270
      IF(ME.EQ.1.OR.ME.EQ.LC3) GO TO 1270
      IF(XS.LT.X(MV-1).OR.YS.LT.Y(MA-1).OR.ZS.LT.Z(ME-1)) GO TO 1270
      IF(XS.GT.X(MV+1).OR.YS.GT.Y(MA+1).OR.ZS.GT.Z(ME+1)) GO TO 1270
      I = MV+1
      J = MA+1
      K = ME+1
      XO=X(I-2)
      YO=Y(J-2)
      ZO=Z(K-2)
      DO 1501 I=1,3
      DO 1501 J=1,3
      DO 1501 K=1,3
 1501 A(I,J,K)=AQM(I,J,K)
      GO TO 193
 1270 IF(XS.LT.X(1).OR.YS.LT.Y(1).OR.ZS.LT.Z(1)) GO TO 170
      DO 140 I=3,LC1,2
      IF(XS.LE.X(I)) GO TO 150
  140 CONTINUE
      GO TO 170
  150 XO=X(I-2)
      DO 160 J=3,LC2,2
      IF(YS.LE.Y(J)) GO TO 151
  160 CONTINUE
      GO TO 170
  151 YO=Y(J-2)
      DO 180 K=3,LC3,2
      IF(ZS.LE.Z(K)) GO TO 200
  180 CONTINUE
  170 IERR = 1
      RETURN
C 
  200 ZO = Z(K-2)
      IQ=(I-1)/2
      JQ=(J-1)/2
      KQ=(K-1)/2
      DO 502 I=1,3
      DO 502 J=1,3
      DO 502 K=1,3
  502 A(I,J,K)=AQA(I,J,K,IQ,JQ,KQ)
  193 SU = 0.
      DO 14 I=2,3
      DO 14 J=2,3
      DO 14 K=2,3
   14 SU = SU + A(I,J,K)*(S(1)-XO)**(I-1)*(S(2)-YO)**(J-1)*(S(3)-ZO)**
     1(K-1)
      DO15 J=2,3
      DO15 K=2,3
   15 SU=SU+A(1,J,K)*(S(2)-YO)**(J-1)*(S(3)-ZO)**(K-1)
      DO16 I=2,3
      DO 16 K=2,3
   16 SU=SU+A(I,1,K)*(S(1)-XO)**(I-1)*(S(3)-ZO)**(K-1)
      DO17 I=2,3
      DO17 J=2,3
   17 SU=SU+A(I,J,1)*(S(1)-XO)**(I-1)*(S(2)-YO)**(J-1)
      DO 18 I=2,3
   18 SU=SU+A(I,1,1)*(S(1)-XO)**(I-1)
      DO 19 J=2,3
   19 SU=SU+A(1,J,1)*(S(2)-YO)**(J-1)
      DO 20 K=2,3
   20 SU=SU+A(1,1,K)*(S(3)-ZO)**(K-1)
      SU = SU+A(1,1,1)
      S(4)=SU
      IERR=0
      RETURN
      END
C
C
C**********************************************************************
C        SUBROUTINE DGELG
C
C        PURPOSE
C           TO SOLVE A GENERAL SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS.
C
C        USAGE
C           CALL DGELG(R,A,M,N,EPS,IER)
C
C        DESCRIPTION OF PARAMETERS
C           R      - DOUBLE PRECISION M BY N RIGHT HAND SIDE MATRIX
C                    (DESTROYED). ON RETURN R CONTAINS THE SOLUTIONS
C                    OF THE EQUATIONS.
C           A      - DOUBLE PRECISION M BY M COEFFICIENT MATRIX
C                    (DESTROYED).
C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.
C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.
C           EPS    - SINGLE PRECISION INPUT CONSTANT WHICH IS USED AS
C                    RELATIVE TOLERANCE FOR TEST ON LOSS OF
C                    SIGNIFICANCE.
C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
C                    IER=0  - NO ERROR,
C                    IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR
C                             PIVOT ELEMENT AT ANY ELIMINATION STEP
C                             EQUAL TO 0,
C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
C                             CANCE INDICATED AT ELIMINATION STEP K+1,
C                             WHERE PIVOT ELEMENT WAS LESS THAN OR
C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
C                             ABSOLUTELY GREATEST ELEMENT OF MATRIX A.
C
C        REMARKS
C           INPUT MATRICES R AND A ARE ASSUMED TO BE STORED COLUMNWISE
C           IN M*N RESP. M*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN
C           SOLUTION MATRIX R IS STORED COLUMNWISE TOO.
C           THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS
C           GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS
C           ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN -
C           INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL
C           SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE
C           INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS
C           GIVEN IN CASE M=1.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH
C           COMPLETE PIVOTING.
C
C     ..................................................................
C
      SUBROUTINE DGELG(R,A,M,N,EPS,IER)

      DIMENSION A(1),R(1)
      DOUBLE PRECISION R,A,PIV,TB,TOL,PIVI
      IF(M)23,23,1

C     SEARCH FOR GREATEST ELEMENT IN MATRIX A
    1 IER=0
      PIV=0.D0
      MM=M*M
      NM=N*M
      DO 3 L=1,MM
      TB=DABS(A(L))
      IF(TB-PIV)3,3,2
    2 PIV=TB
      I=L
    3 CONTINUE
      TOL=EPS*PIV

C     A(I) IS PIVOT ELEMENT. PIV CONTAINS THE ABSOLUTE VALUE OF A(I).
C
C     START ELIMINATION LOOP
      LST=1
      DO 17 K=1,M
C	
C     TEST ON SINGULARITY
      IF(PIV)23,23,4
    4 IF(IER)7,5,7
    5 IF(PIV-TOL)6,6,7
    6 IER=K-1
    7 if(abs(A(I)).lt.1.0e-30)then 
       A(I)=sign(A(I),A(I))
	endif
      PIVI=1.D0/A(I)
      J=(I-1)/M
      I=I-J*M-K
      J=J+1-K

C     I+K IS ROW-INDEX, J+K COLUMN-INDEX OF PIVOT ELEMENT
C
C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R
      DO 8 L=K,NM,M
      LL=L+I
      TB=PIVI*R(LL)
      R(LL)=R(L)
    8 R(L)=TB
C
C     IS ELIMINATION TERMINATED
      IF(K-M)9,18,18
C
C     COLUMN INTERCHANGE IN MATRIX A
    9 LEND=LST+M-K
      IF(J)12,12,10
   10 II=J*M
      DO 11 L=LST,LEND
      TB=A(L)
      LL=L+II
      A(L)=A(LL)
   11 A(LL)=TB
C
C     ROW INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A
   12 DO 13 L=LST,MM,M
      LL=L+I
      TB=PIVI*A(LL)
      A(LL)=A(L)
   13 A(L)=TB
C
C     SAVE COLUMN INTERCHANGE INFORMATION
      A(LST)=J
C
C     ELEMENT REDUCTION AND NEXT PIVOT SEARCH
      PIV=0.D0
      LST=LST+1
      J=0
      DO 16 II=LST,LEND
      PIVI=-A(II)
      IST=II+M
      J=J+1
      DO 15 L=IST,MM,M
      LL=L-J
      A(L)=A(L)+PIVI*A(LL)
      TB=DABS(A(L))
      IF(TB-PIV)15,15,14
   14 PIV=TB
      I=L
   15 CONTINUE
      DO 16 L=K,NM,M
      LL=L+J
   16 R(LL)=R(LL)+PIVI*R(L)
   17 LST=LST+M

C     END OF ELIMINATION LOOP
C
C     BACK SUBSTITUTION AND BACK INTERCHANGE
   18 IF(M-1)23,22,19
   19 IST=MM+M
      LST=M+1
      DO 21 I=2,M
      II=LST-I
      IST=IST-LST
      L=IST-M
      L=A(L)+.5D0
      DO 21 J=II,NM,M
      TB=R(J)
      LL=J
      DO 20 K=IST,MM,M
      LL=LL+1
   20 TB=TB-A(K)*R(LL)
      K=J+L
      R(J)=R(K)
   21 R(K)=TB
   22 print*,'return'
      RETURN
C
C     ERROR RETURN
   23 IER=-1
      print*,'error return'
      RETURN
      END

C***************************************************************
C     MAXIMUMSSUCHROUTINE
      SUBROUTINE MAXSU(VXS,VYS,VZS,FMAX,VL,IPL)
      LOGICAL IPL
      COMMON/KONST/PI,UMRF,LH1,LH2,LH3,PIH
      DIMENSION Z(3)
      DOUBLE PRECISION S(4)
      NZAHL = 0
      FMAX = EXP(FMAX)
      VL1 = VL
      VXS1 = VXS
      VYS1 = VYS
      VZS1 = VZS
      VXS2 = VXS1
      VYS2 = VYS1
      VZS2 = VZS1
      FMAX1 = FMAX
      FMAX2 = FMAX1
  100 CONTINUE
      NZAHL = NZAHL + 1
C     WRITE(6,101) NZAHL
  101 FORMAT(7H NZAHL=I1)
      Z(1) = -VL1
      Z(2) = 0.
      Z(3) = +VL1
      DO 1  I=1,3
      DO 1  J=1,3
      DO 1  K=1,3
      VX = VXS1 + Z(I)
      VY = VYS1 + Z(J)
      VZ = VZS1 + Z(K)
      R = VX**2 + VY**2 + VZ**2
      R1 = SQRT(VX**2 + VY**2)
      S(1) = 1.D0*SQRT(R)
      S(2) = 1.D0*ATAN2(VY,VX)/UMRF
      S(3) = 1.D0*ATAN2(VZ,R1)/UMRF
      CALL CINT(S,IPL,IERR)
      FFMAX = 0.
      IF ( IERR .EQ. 0. ) FFMAX = SNGL(DEXP(S(4)))
C     IF ( IERR .EQ. 1 ) FFMAX= 0.
C     WRITE(6,102) VX,VY,VZ,FFMAX
  102 FORMAT(5X,3E10.3,6HFFMAX=E10.3)
      IF ( FFMAX .LE.  FMAX2 ) GO TO 1
      FMAX2 = FFMAX
      VXS2 = VX
      VYS2 = VY
      VZS2 = VZ
    1 CONTINUE
      IF ( FMAX2 .EQ. FMAX1 ) GO TO 120
      EPS = ABS( (FMAX2-FMAX1)/FMAX1 )
C     WENN DIE ABWEICHUNG DES NEUEN MAXIMUMS VOM ALTEN <0.2% SPRING ANS
C     ENDE
      IF ( EPS .LT. 0.02 ) GO TO 200
      GO TO 121
  120 IF ( FMAX2 .NE. FMAX ) GO TO 200
  121 VL1 = VL1/2.
      FMAX1 = FMAX2
      VXS1 = VXS2
      VYS1 = VYS2
      VZS1 = VZS2
      IF ( NZAHL .GT. 3 ) GO TO 200
      GO TO 100
  200 VXS = VXS2
      VYS = VYS2
      VZS = VZS2
      FMAX = ALOG(FMAX2)
  300 RETURN
      END