aboutsummaryrefslogtreecommitdiff
path: root/orrs/src/orrs.f
blob: 45722240c859988a9681d6a7a7a00e4f5e0cf25e (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
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677

      SUBROUTINE ORRS(LSTI,LREI,NI,YI,UI,UDI, REI, ITMAXI,
     &                ALPHAR,ALPHAI , OMEGAR,OMEGAI, 
     &                UTR,UTI, VTR,VTI, WTR, WTI, CTR, CTI, DELMAX)
      DIMENSION YI(NI), UI(NI), UDI(NI)
      DIMENSION UTR(NI), UTI(NI), VTR(NI), VTI(NI),
     &          WTR(NI), WTI(NI), CTR(NI), CTI(NI)
C---------------------------------------------------------------------
C     Routine for solving the Orr-Sommerfeld equation
C     in the spatial or temporal stability problems.
C
C     Input:
C     ------
C      LSTI    1: spatial  amplification problem
C              2: temporal amplification problem
C      LREI    1: Reynolds number fixed
C              2: Reynolds number variable 
C                 [ to obtain specified ai (LSTI=1), or wi (LSTI=2) ]
C      NI      total number of points in profiles
C      YI      normal BL coordinate array
C      UI      mean flow u(y) profile
C      UDI     mean flow du/dy profile
C      REI     Reynolds number
C      ITMAXI  max number of Newton iterations to seek eigenvalue
C      OMEGAR  real  part of temporal frequency (for initial guess)
C      OMEGAI  imag. part of temporal frequency (for initial guess)
C      ALPHAR  real  part of complex wavenumber (for initial guess)
C      ALPHAI  imag. part of complex wavenumber (for initial guess)
C
C     Output:
C     -------
C      OMEGAR  real  part of temporal frequency (if LSTI = 2)
C      OMEGAI  imag. part of temporal frequency (if LSTI = 2)
C      ALPHAR  real  part of complex wavenumber (if LSTI = 1)
C      ALPHAI  imag. part of complex wavenumber (if LSTI = 1)
C      UTR     real  part of perturbation x-velocity profile
C      UTI     imag. part of perturbation x-velocity profile
C      VTR     real  part of perturbation y-velocity profile
C      VTI     imag. part of perturbation y-velocity profile
C      WTR     real  part of perturbation vorticity profile
C      WTI     imag. part of perturbation vorticity profile
C      CTR     real  part of perturbation d(vorticity)/dy profile
C      CTI     imag. part of perturbation d(vorticity)/dy profile
C      DELMAX  max change in (UTR,UTI) in last iteration (~= 0 if converged)
C---------------------------------------------------------------------
C
      INCLUDE 'ORRS.INC'
C
C---- convergence tolerance
      DATA EPS / 1.0E-4 /
C
      IF(NI.GT.NMAX) STOP 'ORRS: Array overflow.'
C
C---- set initial BC flag
      IBC  = 1
C
C---- set du'/dy normalization constant (normally imposed at wall)
      FNORM = (1.0,-1.0)
ccc      FNORM = (1.0, 0.0)
C
C---- set input variables from parameter list
      LST = LSTI
      LRE = LREI
C
      ITMAX = ITMAXI
C
      N = NI
      DO I=1, N
        Y(I) = YI(I)
        U(I) = UI(I)
        UD(I) = UDI(I)
      ENDDO
C
      IRE   = CMPLX( 0.0  , REI  )
      ALPHA = CMPLX(ALPHAR,ALPHAI)
      OMEGA = CMPLX(OMEGAR,OMEGAI)
C
C---- save initial guess for restoration if normalization condition is relocated
      OMINIT = OMEGA
      ALINIT = ALPHA
C
C
C---- set number of righthand sides
      NRHS = 2
      IF(LRE .EQ. 2) NRHS = 3
C
      CALL OS_INIT
C
C---- Newton iteration loop
      DO 100 ITER=1, ITMAX
C
        CALL OS_SETUP
        CALL OS_SOLVE
        CALL OS_UPDATE
C
        CALL OS_BCCHEK
C
CCC        call newpen(1)
CCC        do 66 i=1, n
CCC          utr(i) = real(f1(i))
CCC   66   continue
CCC        call urplot(n,y,utr)
C
        DELMAX = DFMAX
C
        IF(ITMAX.EQ.1) GO TO 101
C
c        IF(LRE.EQ.1) THEN
c         IF(LST.EQ.1) 
c     &    WRITE(*,7011) ITER,DFMAX,REAL(ALPHA),IMAG(ALPHA)
c 7011     FORMAT(1X,I2,'  max =', E11.4,'  a =', 2F10.6)
c         IF(LST.EQ.2)
c     &    WRITE(*,7012) ITER,DFMAX,REAL(OMEGA),IMAG(OMEGA)
c 7012     FORMAT(1X,I2,'  max =', E11.4,'  w =', 2F10.6)
c        ELSE
c         IF(LST.EQ.1)
c     &    WRITE(*,7021) ITER,DFMAX,REAL(ALPHA),IMAG(ALPHA),IMAG(IRE)
c 7021    FORMAT(1X,I2,'  max =', E11.4,'  a =', 2F10.6,'  Re =',E11.4)
c         IF(LST.EQ.2)
c     &    WRITE(*,7022) ITER,DFMAX,REAL(OMEGA),IMAG(OMEGA),IMAG(IRE)
c 7022    FORMAT(1X,I2,'  max =', E11.4,'  w =', 2F10.6,'  Re =',E11.4)
c        ENDIF
C
        IF(ISOL.NE.0 .AND. DFMAX .LT. EPS) GO TO 101
  100 CONTINUE
      WRITE(*,*) 'ORRS: Convergence failed.  Continuing ...'
C
  101 CONTINUE
C
C---- save variables for passing back to calling routine
      ALPHAR = REAL(ALPHA)
      ALPHAI = IMAG(ALPHA)
      OMEGAR = REAL(OMEGA)
      OMEGAI = IMAG(OMEGA)
      REI = IMAG(IRE)
C
      DO 200 I=1, N
        UTR(I) =  REAL(F1(I))
        UTI(I) =  IMAG(F1(I))
        VTR(I) =  REAL((0.0,-1.0)*ALPHA*F0(I))
        VTI(I) =  IMAG((0.0,-1.0)*ALPHA*F0(I))
        WTR(I) = -REAL(F2(I))
        WTI(I) = -IMAG(F2(I))
        CTR(I) = -REAL(F3(I))
        CTI(I) = -IMAG(F3(I))
  200 CONTINUE
C
      RETURN
      END ! ORRS


      SUBROUTINE OS_INIT
      INCLUDE 'ORRS.INC'
C
      DO I=1, N
        F0(I) = 0.
        F1(I) = 0.
        F2(I) = 0.
        F3(I) = 0.
      ENDDO
C
      ISOL = 0
C
      RETURN
      END


      SUBROUTINE OS_BCCHEK
      INCLUDE 'ORRS.INC'
      COMPLEX FFACT
C
      FWALL = CABS(F2(1))
      FEDGE = CABS(F2(N))
C
      IF(IBC .EQ. 1  .AND.  FEDGE .GT. 2.0*FWALL) THEN
       WRITE(*,*) 'Switching normalizing condition to edge'
       IBC = 2
       FFACT = FNORM/F2(N)
      ELSE IF(IBC .EQ. 2  .AND.  FWALL .GT. 2.0*FEDGE) THEN
       WRITE(*,*) 'Switching normalizing condition to wall'
       IBC = 1
       FFACT = FNORM/F2(1)
      ELSE
       RETURN
      ENDIF
C
      DO I=1, N
        F0(I) = F0(I)*FFACT
        F1(I) = F1(I)*FFACT
        F2(I) = F2(I)*FFACT
        F3(I) = F3(I)*FFACT
      ENDDO
C
      ISOL = 0
      ITMAX = MIN0( ITMAX + 1 , 20 )
      IF(LST.EQ.1) ALPHA = ALINIT
      IF(LST.EQ.2) OMEGA = OMINIT
C
      RETURN
      END ! OS_BCCHEK


      SUBROUTINE OS_SETUP
      INCLUDE 'ORRS.INC'
C------------------------------------------------------
C     Sets up 4x4 block-tridiagnonal system 
C     for Orr-Sommerfeld equation solution.
C
C     The perturbation stream function has the form:
C
C         P(x,y,t;a,w,R,U)  =  p(y) exp[i(ax - wt)]
C
C     The four equations set up are:
C
C         p' = q    2
C         q' = r + a p
C         r' = s                     2
C         s' = iR[(aU-w)r - aU"p] + a p
C
C     p = streamfunction  = F0
C     q = velocity        = F1
C     r = vorticity       = F2
C     s = dr/dy           = F3
C
C------------------------------------------------------
C
C---- zero out A,B,C blocks and righthand sides R
      DO I=1, N
        DO J=1, 4
          DO K=1, 4
            A(J,K,I) = (0.0,0.0)
            B(J,K,I) = (0.0,0.0)
            C(J,K,I) = (0.0,0.0)
          ENDDO
          DO K=1, NRMAX
            R(J,K,I) = (0.0,0.0)
          ENDDO
        ENDDO
      ENDDO
C
      I = 1
C
C---- set 1st wall BC
      R(1,1,I) = F0(I)
      A(1,1,I) = 1.0
C
      IF(IBC.EQ.1) THEN
C
C----- set normalizing condition in lieu of 2nd wall BC (enforced in OS_UPDATE)
       R(2,1,I) = F2(I) - FNORM
       A(2,3,I) = 1.0
C
      ELSE
C
C----- set 2nd wall BC
       R(2,1,I) = F1(I)
       A(2,2,I) = 1.0
C
      ENDIF
C
C---- set interior equations
      DO 50 I=1,N-1
C
        DY = Y(I+1) - Y(I)
        UAV = 0.5*(U(I+1) + U(I))
        UDD = UD(I+1) - UD(I)
C---------------------------------------------------------------
C
        R(1,1,I+1) = F0(I+1) - F0(I) - 0.5*DY*(F1(I+1)+F1(I))
        B(1,1,I+1) = -1.0
        A(1,1,I+1) =  1.0
        B(1,2,I+1) = -0.5*DY
        A(1,2,I+1) = -0.5*DY
C---------------------------------------------------------------
C
        R(2,1,I+1) = F1(I+1) - F1(I)
     &             - 0.5*DY*(   F2(I+1)+F2(I)
     &                         + (F0(I+1)+F0(I))*ALPHA**2 )
        IF(LST.EQ.1)
     &  R(2,2,I+1) = -0.5*DY * (F0(I+1)+F0(I)) * 2.0*ALPHA
        B(2,1,I+1) = -0.5*DY*ALPHA**2
        A(2,1,I+1) = -0.5*DY*ALPHA**2
        B(2,2,I+1) = -1.0
        A(2,2,I+1) =  1.0
        B(2,3,I+1) = -0.5*DY
        A(2,3,I+1) = -0.5*DY
C---------------------------------------------------------------
C
        R(3,1,I) = F2(I+1) - F2(I) - 0.5*DY*(F3(I+1)+F3(I))
        A(3,3,I) = -1.0
        C(3,3,I) =  1.0
        A(3,4,I) = -0.5*DY
        C(3,4,I) = -0.5*DY
C---------------------------------------------------------------
C
        R(4,1,I) = F3(I+1) - F3(I)
     &             - 0.5*DY* (F2(I+1)+F2(I)) * ALPHA**2
     &             - IRE*( (ALPHA*UAV-OMEGA)*0.5*DY*(F2(I+1)+F2(I))
     &                    - ALPHA*UDD*0.5*(F0(I+1)+F0(I)) )
        IF(LST.EQ.1)
     &  R(4,2,I) = - 0.5*DY* (F2(I+1)+F2(I)) * 2.0*ALPHA
     &             - IRE*(        UAV       *0.5*DY*(F2(I+1)+F2(I))
     &                    -       UDD*0.5*(F0(I+1)+F0(I)) )
        IF(LST.EQ.2)
     &  R(4,2,I) = - IRE*( (         -1.0  )*0.5*DY*(F2(I+1)+F2(I)) )
        R(4,3,I) = 
     &        -(0.0,1.0)*( (ALPHA*UAV-OMEGA)*0.5*DY*(F2(I+1)+F2(I))
     &                    - ALPHA*UDD*0.5*(F0(I+1)+F0(I)) )
        A(4,1,I) =                      IRE* ALPHA*UDD*0.5
        C(4,1,I) =                      IRE* ALPHA*UDD*0.5
        A(4,3,I) = -0.5*DY*ALPHA**2 - IRE*(ALPHA*UAV-OMEGA)*0.5*DY
        C(4,3,I) = -0.5*DY*ALPHA**2 - IRE*(ALPHA*UAV-OMEGA)*0.5*DY
        A(4,4,I) = -1.0
        C(4,4,I) =  1.0
C---------------------------------------------------------------
C
   50 CONTINUE
C
C---- set asymptotic regularity conditions at outer edge
C
      FACSQ  = ALPHA**2 + IRE*(ALPHA*U(N)-OMEGA)
      FAC    = CSQRT(FACSQ)
      IF(REAL(FACSQ) .LT. 0.0  .AND.  IMAG(FACSQ) .LT. 0.0) THEN
CCC       WRITE(*,*) 'ORRS: Overdamped mode.'
       FAC = -FAC
      ENDIF
      FAC_AL = (2.0*ALPHA + IRE*U(N)) * 0.5/FAC
      FAC_OM = (          - IRE     ) * 0.5/FAC
      FAC_RE = (0.0,1.0)*(ALPHA*U(N)-OMEGA) * 0.5/FAC
C
      IF(IBC.EQ.2) THEN
C
C----- set normalization condition in lieu of asymptotic regularity condition
       R(3,1,N) = F2(N) - FNORM
       A(3,3,N) = 1.0
C
      ELSE
C
       R(3,1,N) = (ALPHA + FAC   )*(F1(N) + F0(N)*ALPHA) + F2(N)
       IF(LST.EQ.1)
     & R(3,2,N) = (1.0   + FAC_AL)*(F1(N) + F0(N)*ALPHA)
     &          + (ALPHA + FAC   )*(        F0(N)      )
       IF(LST.EQ.2)
     & R(3,2,N) = (        FAC_OM)*(F1(N) + F0(N)*ALPHA)
       R(3,3,N) = (        FAC_RE)*(F1(N) + F0(N)*ALPHA)
       A(3,1,N) = (ALPHA + FAC   )*(              ALPHA)
       A(3,2,N) = (ALPHA + FAC   )
       A(3,3,N) = 1.0
C
      ENDIF
C
      R(4,1,N) = F3(N) + F2(N)*FAC
      IF(LST.EQ.1)
     &R(4,2,N) =         F2(N)*FAC_AL
      IF(LST.EQ.2)
     &R(4,2,N) =         F2(N)*FAC_OM
      R(4,3,N) =         F2(N)*FAC_RE
      A(4,3,N) =               FAC
      A(4,4,N) = 1.0
C
      RETURN
      END ! OS_SETUP


      SUBROUTINE OS_SOLVE
      INCLUDE 'ORRS.INC'
      COMPLEX PIVOT, TEMP
C---------------------------------------------------
C     4x4 complex tridiagonal block solver.
C     Customized for Orr-Sommerfeld equation system,
C     with certain entries assumed to be zero.
C     (Gives large CPU speedup).
C
C      Assumed initial structure for a block row:
C
C         p q r s   p q r s   p q r s
C        |* * # 0| |* * 0 0| |0 0 0 0|  <--  p' = q    2
C        |* * * 0| |* * * 0| |0 0 0 0|  <--  q' = r + a p
C        |# # # 0| |* * * *| |0 0 * *|  <--  r' = s                     2
C        |# # # 0| |* * * *| |* 0 * *|  <--  s' = iR[(au-w)r - au"p] + a p
C
C         B block   A block   C block
C
C        *  assumed nonzero in initial system
C        #  assumed zero in initial system, becoming nonzero due to fill-in
C        0  assumed zero always
C---------------------------------------------------
C
CCC** Backward sweep: Elimination of upper block diagonal (C's).
      DO 1 I=N, 1, -1
C
        IP = I+1
C
C------ don't eliminate Cn block because it doesn't exist
        IF(I.EQ.N) GO TO 12
C
C------ eliminate Ci block, thus modifying Ai and Ri blocks
        DO 111 L=1, 3
          K = 3
          A(K,L,I) = A(K,L,I)
     &             - C(K,3,I)*B(3,L,IP)
     &             - C(K,4,I)*B(4,L,IP)
          K = 4
          A(K,L,I) = A(K,L,I)
     &             - C(K,1,I)*B(1,L,IP)
     &             - C(K,3,I)*B(3,L,IP)
     &             - C(K,4,I)*B(4,L,IP)
  111   CONTINUE
        DO 112 L=1, NRHS
          K = 3
          R(K,L,I) = R(K,L,I)
     &             - C(K,3,I)*R(3,L,IP)
     &             - C(K,4,I)*R(4,L,IP)
          K = 4
          R(K,L,I) = R(K,L,I)
     &             - C(K,1,I)*R(1,L,IP)
     &             - C(K,3,I)*R(3,L,IP)
     &             - C(K,4,I)*R(4,L,IP)
  112   CONTINUE
C
C                                                              -1
CCC---- multiply Bi block and righthand side Ri vectors by (Ai)
C       using Gaussian elimination.
C
   12   CONTINUE
C        
        DO 13 KPIV=4, 2, -1
C
          KM1 = KPIV-1
C
          PIVOT = 1.0/A(KPIV,KPIV,I)
C
C-------- normalize pivot row
          DO 132 L=1, KM1
            A(KPIV,L,I) = A(KPIV,L,I)*PIVOT
  132     CONTINUE
C
          B(KPIV,1,I) = B(KPIV,1,I)*PIVOT
          B(KPIV,2,I) = B(KPIV,2,I)*PIVOT
          B(KPIV,3,I) = B(KPIV,3,I)*PIVOT
C
          DO 134 L=1, NRHS
            R(KPIV,L,I) = R(KPIV,L,I)*PIVOT
  134     CONTINUE
C
C-------- eliminate upper off-diagonal element in Ai block
          K = KM1
            TEMP = A(K,KPIV,I)
            DO 1351 L=KM1, 1, -1
              A(K,L,I) = A(K,L,I) - TEMP*A(KPIV,L,I)
 1351       CONTINUE
            B(K,1,I) = B(K,1,I) - TEMP*B(KPIV,1,I)
            B(K,2,I) = B(K,2,I) - TEMP*B(KPIV,2,I)
            B(K,3,I) = B(K,3,I) - TEMP*B(KPIV,3,I)
            DO 1352 L=1, NRHS
              R(K,L,I) = R(K,L,I) - TEMP*R(KPIV,L,I)
 1352       CONTINUE
C
   13   CONTINUE
C
C
C------ solve for first row
        PIVOT = 1.0/A(1,1,I)
        B(1,1,I) = B(1,1,I)*PIVOT
        B(1,2,I) = B(1,2,I)*PIVOT
        B(1,3,I) = B(1,3,I)*PIVOT
        DO 14 L=1, NRHS
          R(1,L,I) = R(1,L,I)*PIVOT
   14   CONTINUE
C
C------ back substitute (eliminate everything below diagonal in Ai block)
        DO 15 L=1, 3
          B(2,L,I) = B(2,L,I) - A(2,1,I)*B(1,L,I)
          B(3,L,I) = B(3,L,I) - A(3,1,I)*B(1,L,I)
     &                        - A(3,2,I)*B(2,L,I)
          B(4,L,I) = B(4,L,I) - A(4,1,I)*B(1,L,I)
     &                        - A(4,2,I)*B(2,L,I)
     &                        - A(4,3,I)*B(3,L,I)
   15   CONTINUE
C
        DO 16 L=1, NRHS
          R(2,L,I) = R(2,L,I) - A(2,1,I)*R(1,L,I)
          R(3,L,I) = R(3,L,I) - A(3,1,I)*R(1,L,I)
     &                        - A(3,2,I)*R(2,L,I)
          R(4,L,I) = R(4,L,I) - A(4,1,I)*R(1,L,I)
     &                        - A(4,2,I)*R(2,L,I)
     &                        - A(4,3,I)*R(3,L,I)
  16    CONTINUE
C
    1 CONTINUE
C
CCC** Forward sweep: Back substitution using lower block diagonal (Bi's).
      DO 2 I=2, N
        IM = I-1
        DO 21 L=1, NRHS
          DO 211 K=1, 4
            R(K,L,I) = R(K,L,I)
     &        - (  R(1,L,IM)*B(K,1,I)
     &           + R(2,L,IM)*B(K,2,I)
     &           + R(3,L,IM)*B(K,3,I) )
  211     CONTINUE
   21   CONTINUE
    2 CONTINUE
C
      RETURN
      END ! OS_SOLVE


      SUBROUTINE OS_UPDATE
      INCLUDE 'ORRS.INC'
      COMPLEX DF0,DF1,DF2,DF3
      COMPLEX DAW
      COMPLEX RES, RES_AL, RES_OM, RES_RE, RES_F0, RES_F1, RES_F2,
     &             RES_AW
C
      IF(ISOL.EQ.0) THEN
C
C----- no mode solution yet -- don't try to converge on eigenvalue
       DAW = (0.0,0.0)
       DRE = 0.0
C
      ELSE
C
C----- drive eigenvalue (alpha or omega) to satisfy dropped BC at wall or edge
       IF(IBC.EQ.1) THEN
C
C------ wall BC was dropped -- enforce it here
        I = 1
        DAW = (F1(I) - R(2,1,I)) / R(2,2,I)
        DRE = 0.0
C
       ELSE
C
C------ edge BC was dropped -- enforce it here
        RES    = (ALPHA + FAC   )*(F1(N) + F0(N)*ALPHA) + F2(N)
        RES_AL = (1.0   + FAC_AL)*(F1(N) + F0(N)*ALPHA)
     &         + (ALPHA + FAC   )*(        F0(N)      )
        RES_OM = (        FAC_OM)*(F1(N) + F0(N)*ALPHA)
        RES_RE = (        FAC_RE)*(F1(N) + F0(N)*ALPHA)
        RES_F0 = (ALPHA + FAC   )*(              ALPHA)
        RES_F1 = (ALPHA + FAC   )
        RES_F2 = 1.0
C
        IF(LST.EQ.1) RES_AW = RES_AL
        IF(LST.EQ.2) RES_AW = RES_OM
C
        DAW =-(RES   -RES_F0*R(1,1,N)-RES_F1*R(2,1,N)-RES_F2*R(3,1,N))
     &      / (RES_AW-RES_F0*R(1,2,N)-RES_F1*R(2,2,N)-RES_F2*R(3,2,N))
C
       ENDIF
C
      ENDIF
C
C---- set either alpha or omega change (spatial or temporal problem)
      IF(LST.EQ.1) THEN
       DALPHA = DAW
       DOMEGA = (0.0,0.0)
      ELSE
       DALPHA = (0.0,0.0)
       DOMEGA = DAW
      ENDIF
C
C
      RLX = 1.0
C
      DALF = REAL(DALPHA)/ABS(ALPHA)
      IF(RLX*DALF .LT. -.1) RLX = -.1/DALF
      IF(RLX*DALF .GT. 0.1) RLX = 0.1/DALF
C
      DALF = IMAG(DALPHA)/ABS(ALPHA)
      IF(RLX*DALF .LT. -.1) RLX = -.1/DALF
      IF(RLX*DALF .GT. 0.1) RLX = 0.1/DALF
C
      DOMF = REAL(DOMEGA)/ABS(OMEGA)
      IF(RLX*DOMF .LT. -.1) RLX = -.1/DOMF
      IF(RLX*DOMF .GT. 0.1) RLX = 0.1/DOMF
C
      DOMF = IMAG(DOMEGA)/ABS(OMEGA)
      IF(RLX*DOMF .LT. -.1) RLX = -.1/DOMF
      IF(RLX*DOMF .GT. 0.1) RLX = 0.1/DOMF
C
C      DREF = DRE / IMAG(IRE)
C      IF(RLX*DREF .LT. -.2) RLX = -.2/DREF
C      IF(RLX*DREF .GT. 0.3) RLX = 0.3/DREF
CC
C
C==== see if normalizing condition position needs to be changed
C
cC---- predicted wall and edge f" values at next iteration level
c      FWALL = CABS(F2(1) - R(3,1,1) - DAW*R(3,2,1) - DRE*R(3,3,1))
c      FEDGE = CABS(F2(N) - R(3,1,N) - DAW*R(3,2,N) - DRE*R(3,3,N))
cC
cC---- set flag to normalize whatever is bigger by factor of 2
cC
c      IF(IBC .EQ. 1  .AND.  FEDGE .GT. 2.0*FWALL) THEN
cC
c       WRITE(*,*) 'Switching normalizing condition to edge'
c       IBC = 2
c       ITMAX = MIN0( ITMAX+1 , 20 )
c       IF(LST.EQ.1) ALPHA = ALINIT
c       IF(LST.EQ.2) OMEGA = OMINIT
c       RETURN
cC
c      ELSE IF(IBC .EQ. 2  .AND.  FWALL .GT. 2.0*FEDGE) THEN
cC
c       WRITE(*,*) 'Switching normalizing condition to wall'
c       IBC = 1
c       ITMAX = MIN0( ITMAX+1 , 20 )
c       IF(LST.EQ.1) ALPHA = ALINIT
c       IF(LST.EQ.2) OMEGA = OMINIT
c       RETURN
cC
c      ENDIF
C
C
      DFMAX = 0.0
      DFRMS = 0.0
C
C---- perform Newton update on modes
      DO 50 I=1, N
        DF0 = -R(1,1,I) - DAW*R(1,2,I) - DRE*R(1,3,I)
        DF1 = -R(2,1,I) - DAW*R(2,2,I) - DRE*R(2,3,I)
        DF2 = -R(3,1,I) - DAW*R(3,2,I) - DRE*R(3,3,I)
        DF3 = -R(4,1,I) - DAW*R(4,2,I) - DRE*R(4,3,I)
C
        F0(I) = F0(I) + RLX*DF0
        F1(I) = F1(I) + RLX*DF1
        F2(I) = F2(I) + RLX*DF2
        F3(I) = F3(I) + RLX*DF3
C
        D0SQ = (REAL(DF0)**2 + IMAG(DF0)**2)
        D1SQ = (REAL(DF1)**2 + IMAG(DF1)**2)
        D2SQ = (REAL(DF2)**2 + IMAG(DF2)**2)
        D3SQ = (REAL(DF3)**2 + IMAG(DF3)**2)
C
C       IF(D0SQ .GT. DFMAX) THEN
C        KVMAX = 0
C        IVMAX = I
C        DFMAX = D0SQ
C       ENDIF
C
C       IF(D1SQ .GT. DFMAX) THEN
C        KVMAX = 1
C        IVMAX = I
C        DFMAX = D1SQ
C       ENDIF
C
C       IF(D2SQ .GT. DFMAX) THEN
C        KVMAX = 2
C        IVMAX = I
C        DFMAX = D2SQ
C       ENDIF
C
C       IF(D3SQ .GT. DFMAX) THEN
C        KVMAX = 3
C        IVMAX = I
C        DFMAX = D3SQ
C       ENDIF
C
        DFMAX = MAX( DFMAX , D0SQ , D1SQ , D2SQ , D3SQ )
        DFRMS =      DFRMS + D0SQ + D1SQ + D2SQ + D3SQ
   50 CONTINUE
C
      DFMAX = SQRT( DFMAX )
      DFRMS = SQRT( DFRMS / (4.0*FLOAT(N)) )
C
C---- perform Newton update on eigenvalues
      ALPHA = ALPHA + RLX*DALPHA
      OMEGA = OMEGA + RLX*DOMEGA
      IRE   = IRE   + RLX*CMPLX(0.0,DRE)
C
C---- modes are now available 
      ISOL = 1
C
      RETURN
      END ! OS_UPDATE