      PROGRAM INTAI
      PARAMETER (NH=14, NF=20, NR=100)
      REAL H(NH), A(NR,NF,NH), R(NR,NH)
      REAL M(NH), L(NH)
      REAL F(NF,NH)
      REAL FEN(NR)
      REAL RTN(7)
      INTEGER IH1(7),IH2(7)
C
      CHARACTER*1 ANS
C
      DATA H / 2.2,    2.3,    2.4,    2.5,    2.6,    2.8,    3.0, 
     &         3.4,    4.0,    5.0,    7.0,   10.0,   15.0,   20.0 /
C
      DATA M / 1.4240, 0.3225, 0.1231, 0.0408, -.0031, -.0475, -.0683, 
     &         -.0852, -.0904, -.0868, -.0724, -.0568, -.0415, -.0331 /
C
      DATA L / 0.0637, 0.1876, 0.2902, 0.3756, 0.4471, 0.5569, 0.6330, 
     &         0.7168, 0.7540, 0.7153, 0.5731, 0.4083, 0.2582, 0.1786 /
C
      RTN(1) = 10000.0
      RTN(2) = 5000.0
      RTN(3) = 2000.0
      RTN(4) = 500.0
      RTN(5) = 200.0
      RTN(6) = 100.0
      RTN(7) = 50.0
C
      IH1(1) = 1
      IH2(1) = 3
C
      IH1(2) = 3
      IH2(2) = 5
C
      IH1(3) = 5
      IH2(3) = 7
C
      IH1(4) = 7
      IH2(4) = 9
C
      IH1(5) = 9
      IH2(5) = 11
C
      IH1(6) = 11
      IH2(6) = 13
C
      IH1(7) = 13
      IH2(7) = 14
C
      NRANN = 10
C
      ANN = 20.0
      NANN = 10
C
ccc   LMASK = -32640
      LMASK = -30584
ccc   LMASK = -21846
C
      IDEV = 1
      IDEVRP = 2
C
      SIZE = 8.0
      IPSLU = 0
      SCRNFR = 0.85
C
      CALL PLINITIALIZE
C
      PAR = 0.8
      CH = 0.02
C
      DO 5 IH=1, NH
        write(*,*) ih
        CALL NCALC(H(IH),M(IH),L(IH), NR,R(1,IH), NF,F(1,IH),A(1,1,IH))
 5    CONTINUE
C
C
      DO 100 IPLOT=1, 7
        IF(IPLOT.GT.1) CALL PLOT(0.0,0.0,-999)
C
        CALL PLOPEN(SCRNFR,IPSLU,IDEV)
        CALL NEWFACTOR(SIZE)
        CALL PLOT(5.0*CH,5.0*CH,-3)
C
        DELR = RTN(IPLOT)/FLOAT(NRANN)
        RWT = 1.0/RTN(IPLOT)
        CALL XAXIS(0.0,0.0,1.0,RWT*DELR,0.0,DELR,CH,-1)
C
        DA = ANN/FLOAT(NANN)
        AWT = PAR/ANN
        CALL YAXIS(0.0,0.0,PAR,AWT*DA,0.0,DA,CH,-1)
C
        CALL PLGRID(0.0,0.0,NRANN,RWT*DELR,NANN,AWT*DA,LMASK)
C
C
        DO 10 IH=IH1(IPLOT), IH2(IPLOT)
          DO 102 IR=1, NR
            CALL DAMPL(H(IH),R(IR,IH),FEN(IR))
 102      CONTINUE
C
          CALL XYPLOT(NR,R(1,IH),FEN,0.0,RWT,0.0,AWT,2,0.0,0)
C
          DO 105 IR=2, NR
            IFMAX = 1
            AFMAX = 0.0
            DO 1052 IF=1, NF
              IF(A(IR,IF,IH) .GT. AFMAX) THEN
               AFMAX = A(IR,IF,IH)
               IFMAX = IF
              ENDIF
 1052       CONTINUE
            IF(AFMAX.EQ.0.0) GO TO 105
C
ccc            DO 1055 IF=IFMAX, IFMAX
            DO 1055 IF=1, NF
              XPLT1 = RWT*R(IR-1,IH)
              YPLT1 = AWT*A(IR-1,IF,IH)
              XPLT2 = RWT*R(IR,IH)
              YPLT2 = AWT*A(IR,IF,IH)
              IF(YPLT2 .LE. YPLT1) GO TO 1055
              IF(YPLT2 .GT.   PAR) GO TO 1055
              IF(XPLT2 .GT.   1.0) GO TO 10
C
              CALL PLOT(XPLT1,YPLT1,3)
              CALL PLOT(XPLT2,YPLT2,2)
 1055       CONTINUE
 105      CONTINUE
 10     CONTINUE
C
        CALL PLFLUSH
C
        WRITE(*,*) 'Hardcopy ?   N'
        READ(*,8000) ANS
 8000   FORMAT(A)
        IF(INDEX('Yy',ANS).NE.0) CALL REPLOT(IDEVRP)
C
        CALL PLEND
C
 100  CONTINUE
C
      CALL PLCLOSE
      STOP
      END



      SUBROUTINE NCALC(HK,AM,AL, NR,RT, NF,F, A)
C---------------------------------------------------------------
C     Computes N factor for a range of frequencies
C     and Reynolds numbers by integrating growth rates.
C
C     Input:  HK     shape parameter
C             AM     x/Ue dUe/dx
C             AL     theta^2 / nu  dUe/dx
C             NR     number of Rthetas
C             NF     number of frequencies
C
C     Output: RT(.)  Rtheta values
C             F(.)   frequency values
C             A(..)  TS wave amplitudes
C---------------------------------------------------------------
      REAL RT(NR), F(NF), A(NR,NF)
      LOGICAL OK
C
      DW = -2.00/FLOAT(NF-1)
C
      DO 10 IR=1, NR
        DO 105 IF=1, NF
          A(IR,IF) = 0.0
  105   CONTINUE
   10 CONTINUE
C
      HKB = 1.0 / (HK - 1.0)
      RDLC = 2.23 + 1.35*HKB + 0.85*TANH(10.4*HKB - 7.07) - 0.1
      RDC = 10.0**RDLC
      RTC = RDC/HK
C
      WRITE(*,*) 'H Rcr =', HK, RTC
C
      IF(HK.LE.2.21) THEN
       RTN = 3.0*RTC
       DW = -0.20/FLOAT(NF-1)
       WL1 = -1.7
      ELSE IF(HK.LE.2.31) THEN
       RTN = 4.0*RTC
       DW = -0.30/FLOAT(NF-1)
       WL1 = -1.6
      ELSE IF(HK.LE.2.41) THEN
       RTN = 8.0*RTC
       DW = -0.7/FLOAT(NF-1)
       WL1 = -1.5
      ELSE IF(HK.LE.2.51) THEN
       RTN = 12.0*RTC
       DW = -1.20/FLOAT(NF-1)
       WL1 = -1.4
      ELSE IF(HK.LE.2.61) THEN
       RTN = 20.0*RTC
       DW = -1.75/FLOAT(NF-1)
       WL1 = -1.2
      ELSE IF(HK.LE.2.81) THEN
       RTN = 30.0*RTC
       DW = -2.00/FLOAT(NF-1)
       WL1 = -1.0
      ELSE
       RTN = 50.0*RTC
       DW = -2.25/FLOAT(NF-1)
       WL1 = -0.7
      ENDIF
C
ccc      DW = -2.00/FLOAT(NF-1)
C
C
      GEO = (RTN/RTC)**(1.0/FLOAT(NR-1))
      RT(1) = RTC
      DO 20 IR=2, NR
        RT(IR) = RT(IR-1)*GEO
   20 CONTINUE
C
   21 ISTART = 1
C
      REXP = (1.0 - 3.0*AM)/(1.0 + AM)
      AFAC = 0.5*(1.0 + AM) * AL
c
ccc      write(*,*) rexp, afac
C
      IR = ISTART
      UOT1 = RT(IR)**REXP
C
      DO 30 IF=1, NF
        WLOG = WL1 + DW*FLOAT(IF-1)
        F(IF) = 10.0 ** WLOG
   30 CONTINUE
C
      DO 40 IR=ISTART+1, NR
        IRM = IR-1
C
        DRT = RT(IR) - RT(IRM)
        RSP = 0.5*(RT(IR) + RT(IRM))
        HSP = HK
        HSP = AMIN1( HSP , 19.999 )
C
        DO 405 IF=1, NF
          UOT = RSP**REXP
          FSP = F(IF) * (UOT/UOT1)
          CALL OSMAP(RSP,FSP,HSP,
     &                AR,
     &                AR_R, AR_F, AR_H,
     &                ARF_R,ARF_F,ARF_H ,
     &                AI,
     &                AI_R, AI_F, AI_H,
     &                AIF_R,AIF_F,AIF_H , OK)
C
          IF(IR .EQ. ISTART+1) THEN
           IF(AI.LT.0.0) WRITE(*,*) 'Rcrit too high.  H =', HSP
          ENDIF
C
          IF(OK) THEN
           DNDRT = -AI/AFAC
          ELSE
           DNDRT = 0.
          ENDIF
C
          A(IR,IF) = A(IRM,IF) + DNDRT * DRT
          A(IR,IF) = MAX( A(IR,IF) , 0.0 )
  405   CONTINUE
   40 CONTINUE
C
      RETURN
      END


      SUBROUTINE DAMPL(H,RT,AN)
C------------------------------------------------------
C     Returns envelope amplitude for a similar flow.
C
C     Input:  H   shape parameter
C             RT  Rtheta
C
C     Output: AN  n-factor (envelope amplitude)
C------------------------------------------------------
      HMI = H - 1.0
C
      RLCRIT = 2.492/HMI**0.43 + 0.7*(1.0 + TANH(14.0/HMI - 9.24))
      RCRIT = 10.0**RLCRIT
C
      AN = 0.0
      IF(RT .LE. RCRIT) RETURN
C
      DNDR = 0.028*HMI - 0.0345*EXP(-(3.87/HMI - 2.52)**2)
C
      AN = DNDR*(RT - RCRIT)
      RETURN
      END
