aboutsummaryrefslogtreecommitdiff
path: root/orrs/src/mapgen.f
blob: 2af02e078fe50ba737cf3112364331f57ea32d29 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
      PROGRAM MAPGEN
      PARAMETER (NMAX=257,NRX=101,NWX=101)
      REAL ETA(NMAX), F(NMAX), U(NMAX), S(NMAX)
      REAL UTR(NMAX), UTI(NMAX), VTR(NMAX), VTI(NMAX)
      CHARACTER*48 FNAME
      REAL AR(NRX,NWX), AI(NRX,NWX)
      REAL RT(NRX),RTL(NRX), WS(NWX),WSL(NWX)
      LOGICAL CONV(NRX,NWX)
C
      LST = 1
      LRE = 1
C
      WRMAX = 0.25
C
      RESMAX = 0.1
C
C---- default profile parameters
      N = 256
      GEO = 1.02
      ETAE = 14.0
C
      DO 5 IR=1, NRX
        DO 4 IW=1, NWX
          CONV(IR,IW) = .FALSE.
    4   CONTINUE
    5 CONTINUE
C
C---- generate or read in profile
      CALL PFLGET(N,GEO,ETAE,ETA,F,U,S,H)
C
C
      CALL ASKR('Enter lower log10(Rtheta)^',RT1L)
      CALL ASKR('Enter upper log10(Rtheta)^',RT2L)
      CALL ASKI('Enter number of log10(Rtheta )intervals^',NR)
C
      CALL ASKR('Enter lower log10(Wr*sqrt(Rtheta))^',WS1L)
      CALL ASKR('Enter upper log10(Wr*sqrt(Rtheta))^',WS2L)
      CALL ASKI('Enter number of log10(Wr) intervals^',NW)
C
      NRP = NR + 1
      NWP = NW + 1
C
      IF(NRP.GT.NRX) STOP 'Array overflow'
      IF(NWP.GT.NWX) STOP 'Array overflow'
C
      RT1 = 10.0 ** RT1L
      RT2 = 10.0 ** RT2L
      DO 10 IR=1, NRP
        RTL(IR) = RT1L + (RT2L-RT1L)*FLOAT(IR-1)/FLOAT(NR)
        RT(IR) = 10.0 ** RTL(IR)
   10 CONTINUE
C
      WS1 = 10.0 ** WS1L
      WS2 = 10.0 ** WS2L
      DO 15 IW=1, NWP
        WSL(IW) = WS1L + (WS2L-WS1L)*FLOAT(IW-1)/FLOAT(NW)
        WS(IW) = 10.0 ** WSL(IW)
   15 CONTINUE
C
C
      CALL ASKR('Enter initial  ar  for lower Rtheta, upper Wr^',AR0)
      CALL ASKR('Enter initial  ai  for lower Rtheta, upper Wr^',AI0)
C
C
      CALL ASKS('Enter map output filename^',FNAME)
      OPEN(19,FILE=FNAME,STATUS='NEW',FORM='UNFORMATTED')
      WRITE(19) N, H
      WRITE(19) (ETA(I),I=1, N)
      WRITE(19) (U(I)  ,I=1, N)
      WRITE(19) (S(I)  ,I=1, N)
      WRITE(19) NRP, NWP
      WRITE(19) (RTL(IR),IR=1,NRP)
      WRITE(19) (WSL(IW),IW=1,NWP)
C
      IR1 = NRP
      IR2 = 1
      IRD = -1
C
      DO 100 IW=1, NWP
        WRITE(6,2010)
 2010   FORMAT(/1X,'--------------------')
        DO 90 IR=IR1, IR2, IRD
C
          WR = WS(IW)/SQRT(RT(IR))
C
          WRITE(6,2020) IR,IW, RT(IR), WR
 2020     FORMAT(/1X,2I4,'    Rth =', E12.4, '    Wr =', E12.4)
C
          WR0 = WR
          WI0 = 0.0
C
C-------- set initial wavenumber guess
          IRM1 = IR -   IRINCR
          IRM2 = IR - 2*IRINCR
          IRM3 = IR - 3*IRINCR
C
          IWM1 = IW -   IWINCR
          IWM2 = IW - 2*IWINCR
          IWM3 = IW - 3*IWINCR
C
          IF(IRM2.GE.1 .AND. IRM2.LE.NRP  .AND.
     &       IWM1.GE.1 .AND. IWM1.LE.NWP        ) THEN
            AR0 =                 2.0*AR(IRM1,IW  ) - AR(IRM2,IW  )
     &          + AR(IR  ,IWM1) - 2.0*AR(IRM1,IWM1) + AR(IRM2,IWM1)
            AI0 =                 2.0*AI(IRM1,IW  ) - AI(IRM2,IW  )
     &          + AI(IR  ,IWM1) - 2.0*AI(IRM1,IWM1) + AI(IRM2,IWM1)
          ELSE IF(IRM1.GE.1 .AND. IRM1.LE.NRP  .AND.
     &            IWM2.GE.1 .AND. IWM2.LE.NWP        ) THEN
            AR0 =                         AR(IRM1,IW  )
     &          + 2.0*AR(IR  ,IWM1) - 2.0*AR(IRM1,IWM1)
     &          -     AR(IR  ,IWM2) +     AR(IRM1,IWM2)
            AI0 =                         AI(IRM1,IW  )
     &          + 2.0*AI(IR  ,IWM1) - 2.0*AI(IRM1,IWM1)
     &          -     AI(IR  ,IWM2) +     AI(IRM1,IWM2)
          ELSE IF(IRM1.GE.1 .AND. IRM1.LE.NRP  .AND.
     &            IWM1.GE.1 .AND. IWM1.LE.NWP        ) THEN
            AR0 =                 AR(IRM1,IW  )
     &          + AR(IR  ,IWM1) - AR(IRM1,IWM1)
            AI0 =                 AI(IRM1,IW  )
     &          + AI(IR  ,IWM1) - AI(IRM1,IWM1)
          ELSE IF(IRM2.GE.1 .AND. IRM2.LE.NRP) THEN
            AR0 = 2.0*AR(IRM1,IW) - AR(IRM2,IW)
            AI0 = 2.0*AI(IRM1,IW) - AI(IRM2,IW)
          ELSE IF(IWM2.GE.1 .AND. IWM2.LE.NWP) THEN
            AR0 = 2.0*AR(IR,IWM1) - AR(IR,IWM2)
            AI0 = 2.0*AI(IR,IWM1) - AI(IR,IWM2)
          ELSE IF(IRM1.GE.1 .AND. IRM1.LE.NRP) THEN
            AR0 = AR(IRM1,IW)
            AI0 = AI(IRM1,IW)
          ELSE IF(IWM1.GE.1 .AND. IWM1.LE.NWP) THEN
            AR0 = AR(IR,IWM1)
            AI0 = AI(IR,IWM1)
CCC          ELSE
CCC            STOP 'Cannot start in corner and go in'
          ENDIF
c
          AR(IR,IW) = AR0
          AI(IR,IW) = AI0
C
C-------- don't bother with absurdly high frequency
          IF(WR .GE. WRMAX) THEN
           DELMAX = 0.0
           GO TO 89
          ENDIF
C
          ITMAX = 10
          CALL ORRS(LST,LRE,N,ETA,U,S, RT(IR), ITMAX,
     &              AR0,AI0, WR0,WI0, UTR,UTI,VTR,VTI,DELMAX)
C
   89     IF(DELMAX.LT.RESMAX) CONV(IR,IW) = .TRUE.
C
          AR(IR,IW) = AR0
          AI(IR,IW) = AI0
C
   90   CONTINUE
C
        WRITE(19) (AR(IR,IW),IR=1,NRP)
        WRITE(19) (AI(IR,IW),IR=1,NRP)
C
  100 CONTINUE
C
      CLOSE(19)
C
      STOP
      END


      SUBROUTINE PFLGET(N,GEO,ETAE,ETA,F,U,S,H)
      DIMENSION ETA(N),F(N),U(N),S(N)
      CHARACTER*48 FNAME
C
C---- eta coordinate normalized with momentum thickness
      INORM = 3
C
      WRITE(6,*) ' '
      WRITE(6,*) '  1   Falkner-Skan parameter m = x/U dU/dx'
      WRITE(6,*) '  2   Falkner-Skan parameter beta = 2m/(m+1)'
      WRITE(6,*) '  3   Falkner-Skan shape parameter H'
      WRITE(6,*) '  4   General profile input file'
      WRITE(6,*) ' '
      CALL ASKI('Select profile option^',IOPT)
C
      IF(IOPT.NE.4) THEN
       CALL ASKI('Enter number of BL points^',N)
       CALL ASKR('Enter geometric stretching factor^',GEO)
       CALL ASKR('Enter edge eta value^',ETAE)
      ENDIF
C
C
      IF(IOPT.EQ.1) THEN
C
       CALL ASKR('Enter m^',BU)
       CALL FS(INORM,1,BU,H,N,ETAE,GEO,ETA,F,U,S)
C
      ELSE IF(IOPT.EQ.2) THEN
C
       CALL ASKR('Enter beta^',BETA)
       BU = BETA/(2.0-BETA)
       CALL FS(INORM,1,BU,H,N,ETAE,GEO,ETA,F,U,S)
C
      ELSE IF(IOPT.EQ.3) THEN
C
       CALL ASKR('Enter H^',H)
       CALL FS(INORM,2,BU,H,N,ETAE,GEO,ETA,F,U,S)
C
      ELSE
C
       CALL ASKS('Enter profile filename^',FNAME)
       OPEN(1,FILE=FNAME,STATUS='OLD')
       READ(1,*) N, H
       DO 5 I=1, N
         READ(1,*) ETA(I), U(I), S(I)
    5  CONTINUE
       CLOSE(1)
C
       GEO = (ETA(3)-ETA(2)) / (ETA(2)-ETA(1))
      ENDIF
C
      WRITE(6,1050) N, H, ETA(N), GEO
 1050 FORMAT(/' n =', I4,'   H =', F7.3,
     &                   '   Ye =', F7.3,
     &                   '   dYi+1/dYi =',F6.3 /)
C
      RETURN
      END