aboutsummaryrefslogtreecommitdiff
path: root/orrs/old
diff options
context:
space:
mode:
authorDimitri Sokolyuk <demon@dim13.org>2009-05-11 00:27:49 +0000
committerDimitri Sokolyuk <demon@dim13.org>2009-05-11 00:27:49 +0000
commit0d4f43d355de79178b1142e9735902cf641670b6 (patch)
tree2ced2323f6351db2a51090b3fd13eb11f69ff53f /orrs/old
Xfoil 6.97
Diffstat (limited to 'orrs/old')
-rw-r--r--orrs/old/aigen.f421
-rwxr-xr-xorrs/old/oshai.f338
2 files changed, 759 insertions, 0 deletions
diff --git a/orrs/old/aigen.f b/orrs/old/aigen.f
new file mode 100644
index 0000000..b584cc2
--- /dev/null
+++ b/orrs/old/aigen.f
@@ -0,0 +1,421 @@
+
+ PROGRAM AIGEN
+C-----------------------------------------------------------------------
+C Reads OS amplification data ai(R,w) stored in separate files,
+C one file for each H value.
+C
+C Distills this data into arrays which define a tri-cubic spline
+C which can be efficiently interrogated to return the ai(R,W,H)
+C function and its derivatives.
+C
+C The tri-cubic spline data is written out as a binary file,
+C to be read and used in SUBROUTINE OSHAI.
+C-----------------------------------------------------------------------
+C
+ PARAMETER (NMAX=257,NRX=111,NWX=91,NHX=21)
+ REAL ETA(NMAX,NHX), U(NMAX,NHX), S(NMAX,NHX)
+C
+ REAL ATMP(NRX+NWX+NHX), ADTMP(NRX+NWX+NHX)
+ REAL AC(NRX,NWX,NHX,2),
+ & AC_R(NRX,NWX,NHX,2), AC_W(NRX,NWX,NHX,2), AC_H(NRX,NWX,NHX,2),
+ & AC_RW(NRX,NWX,NHX,2),AC_RH(NRX,NWX,NHX,2),AC_WH(NRX,NWX,NHX,2),
+ & AC_RFH(NRX,NWX,NHX,2)
+ REAL RTL(NRX,NHX), WSL(NWX,NHX), HHL(NHX)
+ INTEGER N(NHX), NRP(NHX), NWP(NHX)
+ INTEGER IRP1(NHX),IRP2(NHX),IWP1(NHX),IWP2(NHX)
+C
+ PARAMETER (NRZ=31, NWZ=41, NHZ=21)
+ INTEGER IW1(NHZ), IW2(NHZ), IR1(NHZ), IR2(NHZ)
+ REAL RL(NRZ), WL(NWZ), HL(NHZ),
+ & A(NRZ,NWZ,NHZ),
+ & AR(NRZ,NWZ,NHZ), AW(NRZ,NWZ,NHZ), AH(NRZ,NWZ,NHZ),
+ & ARW(NRZ,NWZ,NHZ),ARH(NRZ,NWZ,NHZ),AWH(NRZ,NWZ,NHZ),
+ & ARWH(NRZ,NWZ,NHZ)
+C
+ CHARACTER ARGP1
+ LOGICAL LSPLINE
+C
+C---- if T, use splines to compute derivatives, otherwise use finite-diff.
+ LSPLINE = .TRUE.
+C
+C---- strides in R and W file values selected for storage in binary table
+C- (i.e. binary table can be less dense than the source storage files)
+ IRINC = 4
+ IWINC = 2
+C
+C---- set expeced format of source files
+ IFORM = 0 ! binary
+ccc IFORM = 1 ! ascii
+C
+ CALL READOS(ARGP1,IFORM,
+ & N,NMAX,ETA,U,S,
+ & NRP,NWP,NHP,
+ & RTL,WSL,HHL, AC(1,1,1,1), AC(1,1,1,2),
+ & NRX,NWX,NHX)
+C
+C
+ RTLMIN = RTL(1,1)
+ WSLMIN = WSL(1,1)
+ RTLMAX = RTL(1,1)
+ WSLMAX = WSL(1,1)
+ DO 10 IHP=1, NHP
+ RTLMIN = MIN( RTLMIN , RTL(1,IHP) )
+ WSLMIN = MIN( WSLMIN , WSL(1,IHP) )
+ RTLMAX = MAX( RTLMAX , RTL(NRP(IHP),IHP) )
+ WSLMAX = MAX( WSLMAX , WSL(NWP(IHP),IHP) )
+ 10 CONTINUE
+C
+ DRTL = RTL(2,1) - RTL(1,1)
+ DWSL = WSL(2,1) - WSL(1,1)
+C
+ NRPTOT = INT( (RTLMAX - RTLMIN)/DRTL + 1.001 )
+ NWPTOT = INT( (WSLMAX - WSLMIN)/DWSL + 1.001 )
+C
+ IF(NRPTOT .GT. NRX) STOP 'AIGEN: R index overflow'
+ IF(NWPTOT .GT. NWX) STOP 'AIGEN: W index overflow'
+C
+C---- move ai array for each H to a common origin for splining
+ DO 20 IHP=1, NHP
+ IROFF = INT( (RTL(1,IHP) - RTLMIN)/DRTL + 0.001 )
+ IWOFF = INT( (WSL(1,IHP) - WSLMIN)/DWSL + 0.001 )
+ IF(IROFF.EQ.0 .AND. IWOFF.EQ.0) GO TO 19
+C
+ DO IC = 1, 2
+ DO IRP=NRP(IHP), 1, -1
+ DO IWP=NWP(IHP), 1, -1
+ AC(IRP+IROFF,IWP+IWOFF,IHP) = AC(IRP,IWP,IHP)
+ AC(IRP,IWP,IHP) = 0.0
+ ENDDO
+ ENDDO
+ ENDDO
+C
+ IF(IROFF.GT.0) THEN
+ DO IRP=NRP(IHP), 1, -1
+ RTL(IRP+IROFF,IHP) = RTL(IRP,IHP)
+ RTL(IRP,IHP) = 0.0
+ ENDDO
+ ENDIF
+C
+ IF(IWOFF.GT.0) THEN
+ DO IWP=NWP(IHP), 1, -1
+ WSL(IWP+IWOFF,IHP) = WSL(IWP,IHP)
+ WSL(IWP,IHP) = 0.0
+ ENDDO
+ ENDIF
+C
+ 19 IRP1(IHP) = IROFF + 1
+ IWP1(IHP) = IWOFF + 1
+ IRP2(IHP) = IROFF + NRP(IHP)
+ IWP2(IHP) = IWOFF + NWP(IHP)
+C
+C------ set newly-defined R and W coordinate values
+ DO IRP=1, IRP1(IHP)-1
+ RTL(IRP,IHP) = RTL(IRP1(IHP),IHP) + DRTL*FLOAT(IRP-IRP1(IHP))
+ ENDDO
+ DO IRP=IRP2(IHP)+1, NRPTOT
+ RTL(IRP,IHP) = RTL(IRP2(IHP),IHP) + DRTL*FLOAT(IRP-IRP2(IHP))
+ ENDDO
+C
+ DO IWP=1, IWP1(IHP)-1
+ WSL(IWP,IHP) = WSL(IWP1(IHP),IHP) + DWSL*FLOAT(IWP-IWP1(IHP))
+ ENDDO
+ DO IWP=IWP2(IHP)+1, NWPTOT
+ WSL(IWP,IHP) = WSL(IWP2(IHP),IHP) + DWSL*FLOAT(IWP-IWP2(IHP))
+ ENDDO
+C
+ 20 CONTINUE
+C
+C---- differentiate in H with spline routine to get AC_H
+ DO 40 IRP=1, NRPTOT
+ DO 401 IWP=1, NWPTOT
+C
+C-------- find first H index at this R,w
+ DO IHP=1, NHP
+ IF(IRP.GE.IRP1(IHP) .AND. IRP.LE.IRP2(IHP) .AND.
+ & IWP.GE.IWP1(IHP) .AND. IWP.LE.IWP2(IHP) ) GO TO 4012
+ ENDDO
+ GO TO 401
+ 4012 IHP1 = IHP
+C
+C-------- find last H index at this R,w
+ DO IHP=NHP, 1, -1
+ IF(IRP.GE.IRP1(IHP) .AND. IRP.LE.IRP2(IHP) .AND.
+ & IWP.GE.IWP1(IHP) .AND. IWP.LE.IWP2(IHP) ) GO TO 4022
+ ENDDO
+ GO TO 401
+ 4022 IHP2 = IHP
+C
+ DO IC = 1, 2
+ DO IHP=IHP1, IHP2
+ ATMP(IHP) = AC(IRP,IWP,IHP,IC)
+ ENDDO
+C
+ IHPNUM = IHP2 - IHP1 + 1
+ CALL SPLINE(ATMP(IHP1),ADTMP(IHP1),HHL(IHP1),IHPNUM)
+C
+ DO IHP=IHP1, IHP2
+ AC_H(IRP,IWP,IHP,IC) = ADTMP(IHP)
+ ENDDO
+ ENDDO
+C
+ 401 CONTINUE
+ 40 CONTINUE
+C
+C
+ DO IC = 1, 2
+C
+ IF(LSPLINE) THEN
+C----- calculate AC_R and AC_W arrays from spline coefficients
+ CALL RDIFFS(IRP1,IRP2,IWP1,IWP2,NHP,RTL,NRX,NWX,
+ & AC(1,1,1,1), AC_R(1,1,1,1) )
+ CALL WDIFFS(IRP1,IRP2,IWP1,IWP2,NHP,WSL,NRX,NWX,
+ & AC(1,1,1,1), AC_W(1,1,1,1) )
+C
+ CALL RDIFFS(IRP1,IRP2,IWP1,IWP2,NHP,RTL,NRX,NWX,
+ & AC_H(1,1,1,1), AC_RH(1,1,1,1) )
+ CALL WDIFFS(IRP1,IRP2,IWP1,IWP2,NHP,WSL,NRX,NWX,
+ & AC_H(1,1,1,1), AC_WH(1,1,1,1) )
+ CALL WDIFFS(IRP1,IRP2,IWP1,IWP2,NHP,WSL,NRX,NWX,
+ & AC_R(1,1,1,1), AC_RF(1,1,1,1) )
+C
+ CALL WDIFFS(IRP1,IRP2,IWP1,IWP2,NHP,WSL,NRX,NWX,
+ & AC_RH(1,1,1,1), AC_RFH(1,1,1,1) )
+C
+ ELSE
+C----- calculate AC_R and AC_W arrays by finite-differencing
+ CALL RDIFF(IRP1,IRP2,IWP1,IWP2,NHP,RTL,NRX,NWX,
+ & AC(1,1,1,1), AC_R(1,1,1,1) )
+ CALL WDIFF(IRP1,IRP2,IWP1,IWP2,NHP,WSL,NRX,NWX,
+ & AC(1,1,1,1), AC_W(1,1,1,1) )
+C
+ CALL RDIFF(IRP1,IRP2,IWP1,IWP2,NHP,RTL,NRX,NWX,
+ & AC_H(1,1,1,1), AC_RH(1,1,1,1) )
+ CALL WDIFF(IRP1,IRP2,IWP1,IWP2,NHP,WSL,NRX,NWX,
+ & AC_H(1,1,1,1), AC_WH(1,1,1,1) )
+ CALL WDIFF(IRP1,IRP2,IWP1,IWP2,NHP,WSL,NRX,NWX,
+ & AC_R(1,1,1,1), AC_RF(1,1,1,1) )
+C
+ CALL WDIFF(IRP1,IRP2,IWP1,IWP2,NHP,WSL,NRX,NWX,
+ & AC_RH(1,1,1,1), AC_RFH(1,1,1,1) )
+ ENDIF
+C
+ ENDDO
+C
+C
+ NR = (NRPTOT-1)/IRINC + 1
+ NW = (NWPTOT-1)/IWINC + 1
+ NH = NHP
+C
+ DO 60 IHP=1, NHP
+ IH = IHP
+ IR1(IH) = (IRP1(IHP)-1)/IRINC + 1
+ IR2(IH) = (IRP2(IHP)-1)/IRINC + 1
+ IW1(IH) = (IWP1(IHP)-1)/IWINC + 1
+ IW2(IH) = (IWP2(IHP)-1)/IWINC + 1
+C
+ IC = 2
+ DO 605 IR=1, NR
+ IRP = IRINC*(IR-1) + 1
+ DO 6055 IW=1, NW
+ IWP = IWINC*(IW-1) + 1
+ A (IR,IW,IH) = AC (IRP,IWP,IHP,IC)
+ AR (IR,IW,IH) = AC_R (IRP,IWP,IHP,IC)
+ AW (IR,IW,IH) = AC_W (IRP,IWP,IHP,IC)
+ AH (IR,IW,IH) = AC_H (IRP,IWP,IHP,IC)
+ ARW (IR,IW,IH) = AC_RW (IRP,IWP,IHP,IC)
+ ARH (IR,IW,IH) = AC_RH (IRP,IWP,IHP,IC)
+ AWH (IR,IW,IH) = AC_WH (IRP,IWP,IHP,IC)
+ ARWH(IR,IW,IH) = AC_RWH(IRP,IWP,IHP,IC)
+ 6055 CONTINUE
+ 605 CONTINUE
+C
+ 60 CONTINUE
+C
+C
+ IHP = 1
+C
+ DO IR=1, NR
+ IRP = IRINC*(IR-1) + 1
+ RL(IR) = RTL(IRP,IHP)
+ ENDDO
+C
+ DO IW=1, NW
+ IWP = IWINC*(IW-1) + 1
+ WL(IW) = WSL(IWP,IHP)
+ ENDDO
+C
+ DO IH=1, NH
+ IHP = IH
+ HL(IH) = HHL(IHP)
+ ENDDO
+C
+C
+ OPEN(UNIT=31,FILE='oshai.dat',STATUS='NEW',FORM='UNFORMATTED')
+ WRITE(31) NR, NW, NH
+ WRITE(31) (RL(IR), IR=1,NR)
+ WRITE(31) (WL(IW), IW=1,NW)
+ WRITE(31) (HL(IH), IH=1,NH)
+ WRITE(31) (IR1(IH),IR2(IH),IW1(IH),IW2(IH), IH=1,NH)
+ DO 3 IH=1, NH
+ DO 2 IW=IW1(IH), IW2(IH)
+ WRITE(31) ( A(IR,IW,IH), IR=IR1(IH),IR2(IH))
+ WRITE(31) ( AR(IR,IW,IH), IR=IR1(IH),IR2(IH))
+ WRITE(31) ( AW(IR,IW,IH), IR=IR1(IH),IR2(IH))
+ WRITE(31) ( AH(IR,IW,IH), IR=IR1(IH),IR2(IH))
+ WRITE(31) ( ARW(IR,IW,IH), IR=IR1(IH),IR2(IH))
+ WRITE(31) ( ARH(IR,IW,IH), IR=IR1(IH),IR2(IH))
+ WRITE(31) ( AWH(IR,IW,IH), IR=IR1(IH),IR2(IH))
+ WRITE(31) (ARWH(IR,IW,IH), IR=IR1(IH),IR2(IH))
+ 2 CONTINUE
+ 3 CONTINUE
+ CLOSE(31)
+C
+ STOP
+ END
+
+
+
+
+ SUBROUTINE RDIFF(IRP1,IRP2,IWP1,IWP2,NHP,RTL,NRX,NWX,
+ & AC, AC_R )
+ COMPLEX AC(NRX,NWX,*),AC_R(NRX,NWX,*)
+ REAL RTL(NRX,*)
+ INTEGER IRP1(*),IRP2(*),IWP1(*),IWP2(*)
+C
+ DO 1 IHP=1, NHP
+C
+C------ differentiate in R with finite differences
+ DO 10 IWP=IWP1(IHP), IWP2(IHP)
+ IRP = IRP1(IHP)
+ DELR = RTL(IRP+1,IHP) - RTL(IRP,IHP)
+ AC_R(IRP,IWP,IHP) = (-3.0*AC(IRP ,IWP,IHP)
+ & + 4.0*AC(IRP+1,IWP,IHP)
+ & - AC(IRP+2,IWP,IHP) )/DELR
+ IRP = IRP2(IHP)
+ DELR = RTL(IRP,IHP) - RTL(IRP-1,IHP)
+ AC_R(IRP,IWP,IHP) = ( 3.0*AC(IRP ,IWP,IHP)
+ & - 4.0*AC(IRP-1,IWP,IHP)
+ & + AC(IRP-2,IWP,IHP) )/DELR
+ DO 101 IRP=IRP1(IHP)+1, IRP2(IHP)-1
+ DELR = RTL(IRP+1,IHP) - RTL(IRP-1,IHP)
+ AC_R(IRP,IWP,IHP) = ( AC(IRP+1,IWP,IHP)
+ & - AC(IRP-1,IWP,IHP) )/DELR
+ 101 CONTINUE
+ 10 CONTINUE
+C
+ 1 CONTINUE
+C
+ RETURN
+ END
+
+
+
+ SUBROUTINE WDIFF(IRP1,IRP2,IWP1,IWP2,NHP,WSL,NRX,NWX,
+ & AC, AC_W)
+ COMPLEX AC(NRX,NWX,*),AC_W(NRX,NWX,*)
+ REAL WSL(NWX,*)
+ INTEGER IRP1(*),IRP2(*),IWP1(*),IWP2(*)
+C
+ DO 1 IHP=1, NHP
+C
+C------ differentiate in F with finite differences
+ DO 10 IRP=IRP1(IHP), IRP2(IHP)
+ IWP = IWP1(IHP)
+ DELF = WSL(IWP+1,IHP) - WSL(IWP,IHP)
+ AC_W(IRP,IWP,IHP) = (-3.0*AC(IRP,IWP ,IHP)
+ & + 4.0*AC(IRP,IWP+1,IHP)
+ & - AC(IRP,IWP+2,IHP) )/DELF
+ IWP = IWP2(IHP)
+ DELF = WSL(IWP,IHP) - WSL(IWP-1,IHP)
+ AC_W(IRP,IWP,IHP) = ( 3.0*AC(IRP,IWP ,IHP)
+ & - 4.0*AC(IRP,IWP-1,IHP)
+ & + AC(IRP,IWP-2,IHP) )/DELF
+ DO 101 IWP=IWP1(IHP)+1, IWP2(IHP)-1
+ DELF = WSL(IWP+1,IHP) - WSL(IWP-1,IHP)
+ AC_W(IRP,IWP,IHP) = ( AC(IRP,IWP+1,IHP)
+ & - AC(IRP,IWP-1,IHP) )/DELF
+ 101 CONTINUE
+ 10 CONTINUE
+C
+ 1 CONTINUE
+C
+ RETURN
+ END
+
+
+
+ SUBROUTINE RDIFFS(IRP1,IRP2,IWP1,IWP2,NHP,RTL,NRX,NWX,
+ & AC, AC_R )
+ COMPLEX AC(NRX,NWX,*),AC_R(NRX,NWX,*)
+ REAL RTL(NRX,*)
+ INTEGER IRP1(*),IRP2(*),IWP1(*),IWP2(*)
+C
+ PARAMETER (NDIM=500)
+ REAL ATMP(NDIM), ADTMP(NDIM)
+C
+ DO 1 IHP=1, NHP
+ IF(IRP2(IHP).GT.NDIM) THEN
+ WRITE(*,*) 'RDIFFS: Array overflow. Increase NDIM to',IRP2(IHP)
+ STOP
+ ENDIF
+C
+C------ differentiate in R with spline
+ DO 10 IWP=IWP1(IHP), IWP2(IHP)
+C
+ DO 101 IRP=IRP1(IHP), IRP2(IHP)
+ ATMP(IRP) = AC(IRP,IWP,IHP)
+ 101 CONTINUE
+C
+ IRP = IRP1(IHP)
+ NUM = IRP2(IHP) - IRP1(IHP) + 1
+ CALL SPLINE(ATMP(IRP),ADTMP(IRP),RTL(IRP,IHP),NUM)
+C
+ DO 102 IRP=IRP1(IHP), IRP2(IHP)
+ AC_R(IRP,IWP,IHP) = ADTMP(IRP)
+ 102 CONTINUE
+C
+ 10 CONTINUE
+C
+ 1 CONTINUE
+C
+ RETURN
+ END
+
+
+
+ SUBROUTINE WDIFFS(IRP1,IRP2,IWP1,IWP2,NHP,WSL,NRX,NWX,
+ & AC, AC_W)
+ COMPLEX AC(NRX,NWX,*),AC_W(NRX,NWX,*)
+ REAL WSL(NWX,*)
+ INTEGER IRP1(*),IRP2(*),IWP1(*),IWP2(*)
+
+ PARAMETER (NDIM=500)
+ REAL ATMP(NDIM), ADTMP(NDIM)
+C
+ DO 1 IHP=1, NHP
+ IF(IWP2(IHP).GT.NDIM) THEN
+ WRITE(*,*) 'WDIFFS: Array overflow. Increase NDIM to',IWP2(IHP)
+ STOP
+ ENDIF
+C
+C------ differentiate in F with spline
+ DO 10 IRP=IRP1(IHP), IRP2(IHP)
+C
+ DO 101 IWP=IWP1(IHP), IWP2(IHP)
+ ATMP(IWP) = AC(IRP,IWP,IHP)
+ 101 CONTINUE
+C
+ IWP = IWP1(IHP)
+ NUM = IWP2(IHP) - IWP1(IHP) + 1
+ CALL SPLINE(ATMP(IWP),ADTMP(IWP),WSL(IWP,IHP),NUM)
+C
+ DO 102 IWP=IWP1(IHP), IWP2(IHP)
+ AC_W(IRP,IWP,IHP) = ADTMP(IWP)
+ 102 CONTINUE
+C
+ 10 CONTINUE
+C
+ 1 CONTINUE
+C
+ RETURN
+ END
diff --git a/orrs/old/oshai.f b/orrs/old/oshai.f
new file mode 100755
index 0000000..2029c59
--- /dev/null
+++ b/orrs/old/oshai.f
@@ -0,0 +1,338 @@
+ SUBROUTINE OSHAI(RSP,FSP,HSP, AI,
+ & AI_R, AI_F, AI_H,
+ & AIF_R,AIF_F,AIF_H , OK)
+C---------------------------------------------------------------------
+C
+C Returns imaginary part of complex wavenumber (Ai) eigenvalue
+C from Orr-Sommerfeld solution with mean profiles characterized
+C by shape parameter H. Also returns the sensitivities of Ai
+C with respect to the input parameters.
+C
+C The eigenvalue Ai(Rtheta,f,H) is stored as a 3-D array at
+C discrete points, which is then interpolated to any (Rtheta,f,H)
+C via a tricubic spline. The spline coordinates actually used are:
+C
+C RL = log10(Rtheta)
+C FL = log10(f) + 0.5 log10(Rtheta)
+C HL = H
+C
+C
+C Input:
+C ------
+C RSP momentum thickness Reynolds number Rtheta = Theta Ue / v
+C FSP normalized disturbance frequency f = w Theta/Ue
+C HSP shape parameter of mean profile H = Dstar/Theta
+C
+C Output:
+C -------
+C AI imaginary part of complex wavenumber * Theta
+C AI_R d(AI)/dRtheta
+C AI_F d(AI)/df
+C AI_H d(AI)/dH
+C AIF_R d(dAI/df)/dRtheta
+C AIF_F d(dAI/df)/df
+C AIF_H d(dAI/df)/dH
+C OK T if look up was successful; all values returned are valid
+C F if point fell outside (RL,FL) spline domain limits;
+C all values (AI, AI_R, etc.) are returned as zero.
+C Exception: If points only falls outside HL spline limits,
+C then the HL limit is used and an AI value is calculated,
+C but OK is still returned as F.
+C
+C---------------------------------------------------------------------
+C
+ REAL B(2,2), BR(2,2), BF(2,2), BH(2,2),
+ & BRF(2,2),BRH(2,2),BFH(2,2),BRFH(2,2)
+ REAL C(2) , CR(2) , CF(2) , CH(2) ,
+ & CRF(2) ,CRH(2) ,CFH(2) ,CRFH(2)
+C
+ PARAMETER (NRX=31, NFX=41, NHX=21)
+ COMMON /AICOM/ NR, NF, NH,
+ & IF1(NHX), IF2(NHX), IR1(NHX),IR2(NHX),
+ & RINCR, FINCR, RL(NRX), FL(NFX), HL(NHX),
+ & A(NRX,NFX,NHX),
+ & AR(NRX,NFX,NHX),
+ & AF(NRX,NFX,NHX),
+ & AH(NRX,NFX,NHX),
+ & ARF(NRX,NFX,NHX),
+ & ARH(NRX,NFX,NHX),
+ & AFH(NRX,NFX,NHX),
+ & ARFH(NRX,NFX,NHX)
+ LOGICAL LOADED, OK
+ SAVE LOADED
+C
+ DATA LOADED /.FALSE./
+C
+C---- set ln(10) for derivatives of log10 function
+ DATA AL10 /2.302585093/
+C
+C
+ IF(LOADED) GO TO 9
+C
+C---- first time OSHAI is called ... load in 3-D spline data
+ OPEN(UNIT=31,FILE='~/codes/mses/orrs/oshai.dat',
+ & STATUS='OLD',FORM='UNFORMATTED')
+ WRITE(*,*) 'Loading Orr-Sommerfeld maps...'
+ READ(31) NR, NF, NH
+ IF(NR.GT.NRX) STOP 'OSHAI: R array limit overflow.'
+ IF(NF.GT.NFX) STOP 'OSHAI: F array limit overflow.'
+ IF(NH.GT.NHX) STOP 'OSHAI: H array limit overflow.'
+ READ(31) (RL(IR), IR=1,NR)
+ READ(31) (FL(IF), IF=1,NF)
+ READ(31) (HL(IH), IH=1,NH)
+ READ(31) (IR1(IH),IR2(IH),IF1(IH),IF2(IH), IH=1,NH)
+ DO 3 IH=1, NH
+ DO 2 IF=IF1(IH), IF2(IH)
+ READ(31) ( A(IR,IF,IH), IR=IR1(IH),IR2(IH))
+ READ(31) ( AR(IR,IF,IH), IR=IR1(IH),IR2(IH))
+ READ(31) ( AF(IR,IF,IH), IR=IR1(IH),IR2(IH))
+ READ(31) ( AH(IR,IF,IH), IR=IR1(IH),IR2(IH))
+ READ(31) ( ARF(IR,IF,IH), IR=IR1(IH),IR2(IH))
+ READ(31) ( ARH(IR,IF,IH), IR=IR1(IH),IR2(IH))
+ READ(31) ( AFH(IR,IF,IH), IR=IR1(IH),IR2(IH))
+ READ(31) (ARFH(IR,IF,IH), IR=IR1(IH),IR2(IH))
+ 2 CONTINUE
+ 3 CONTINUE
+ CLOSE(31)
+C
+ RINCR = (RL(NR) - RL(1))/FLOAT(NR-1)
+ FINCR = (FL(NF) - FL(1))/FLOAT(NF-1)
+ LOADED = .TRUE.
+C
+ 9 CONTINUE
+C
+C---- set returned variables in case of out-of-limits error
+ AI = 0.0
+ AI_R = 0.0
+ AI_F = 0.0
+ AI_H = 0.0
+ AIF_R = 0.0
+ AIF_F = 0.0
+ AIF_H = 0.0
+C
+C---- define specified spline coordinates
+ RLSP = ALOG10(RSP)
+ FLSP = ALOG10(FSP) + 0.5*RLSP
+ HLSP = HSP
+C
+ OK = .TRUE.
+C
+C---- find H interval
+ DO 10 IH=2, NH
+ IF(HL(IH) .GE. HLSP) GO TO 11
+ 10 CONTINUE
+ IH = NH
+ 11 CONTINUE
+C
+ IF(HSP.LT.HL(1) .OR. HSP.GT.HL(NH)) THEN
+ OK = .FALSE.
+CCC WRITE(6,*) 'Over H limits. R w H:', RSP,FSP,HSP
+CCC RETURN
+ ENDIF
+C
+C---- find R interval
+ IR = INT((RLSP-RL(1))/RINCR + 2.001)
+ IR1MAX = MAX0( IR1(IH) , IR1(IH-1) )
+ IR2MIN = MIN0( IR2(IH) , IR2(IH-1) )
+ IF(IR-1.LT.IR1MAX .OR. IR.GT.IR2MIN) THEN
+ OK = .FALSE.
+CCC WRITE(6,*) 'Over R limits. R w H:', RSP,FSP,HSP
+CCC RETURN
+ ENDIF
+C
+C---- find F interval
+ IF = INT((FLSP-FL(1))/FINCR + 2.001)
+ IF1MAX = MAX0( IF1(IH) , IF1(IH-1) )
+ IF2MIN = MIN0( IF2(IH) , IF2(IH-1) )
+ IF(IF-1.LT.IF1MAX .OR. IF.GT.IF2MIN) THEN
+ OK = .FALSE.
+CCC WRITE(6,*) 'Over w limits. R w H:', RSP,FSP,HSP
+CCC RETURN
+ ENDIF
+C
+C
+ DRL = RL(IR) - RL(IR-1)
+ DFL = FL(IF) - FL(IF-1)
+ DHL = HL(IH) - HL(IH-1)
+ TR = (RLSP - RL(IR-1)) / DRL
+ TF = (FLSP - FL(IF-1)) / DFL
+ TH = (HLSP - HL(IH-1)) / DHL
+C
+C---- evaluate spline in Rtheta at the corners of HL,FL cell
+ DO 20 KH=1, 2
+ JH = IH + KH-2
+ DO 205 KF=1, 2
+ JF = IF + KF-2
+ A1 = A (IR-1,JF,JH)
+ AR1 = AR (IR-1,JF,JH)
+ AF1 = AF (IR-1,JF,JH)
+ AH1 = AH (IR-1,JF,JH)
+ ARF1 = ARF (IR-1,JF,JH)
+ ARH1 = ARH (IR-1,JF,JH)
+ AFH1 = AFH (IR-1,JF,JH)
+ ARFH1 = ARFH(IR-1,JF,JH)
+C
+ A2 = A (IR ,JF,JH)
+ AR2 = AR (IR ,JF,JH)
+ AF2 = AF (IR ,JF,JH)
+ AH2 = AH (IR ,JF,JH)
+ ARF2 = ARF (IR ,JF,JH)
+ ARH2 = ARH (IR ,JF,JH)
+ AFH2 = AFH (IR ,JF,JH)
+ ARFH2 = ARFH(IR ,JF,JH)
+C
+ DA1 = DRL*AR1 - A2 + A1
+ DA2 = DRL*AR2 - A2 + A1
+ DAF1 = DRL*ARF1 - AF2 + AF1
+ DAF2 = DRL*ARF2 - AF2 + AF1
+ DAH1 = DRL*ARH1 - AH2 + AH1
+ DAH2 = DRL*ARH2 - AH2 + AH1
+ DAFH1 = DRL*ARFH1 - AFH2 + AFH1
+ DAFH2 = DRL*ARFH2 - AFH2 + AFH1
+C
+C-------- set AI, dAI/dFL, dAI/dHL, d2AI/dHLdFL
+ B(KF,KH) = (1.0-TR)* A1 + TR* A2
+ & + ((1.0-TR)*DA1 - TR*DA2 )*(TR-TR*TR)
+ BF(KF,KH) = (1.0-TR)* AF1 + TR* AF2
+ & + ((1.0-TR)*DAF1 - TR*DAF2 )*(TR-TR*TR)
+ BH(KF,KH) = (1.0-TR)* AH1 + TR* AH2
+ & + ((1.0-TR)*DAH1 - TR*DAH2 )*(TR-TR*TR)
+ BFH(KF,KH) = (1.0-TR)* AFH1 + TR* AFH2
+ & + ((1.0-TR)*DAFH1 - TR*DAFH2)*(TR-TR*TR)
+C
+C-------- also, the RL derivatives of the quantities above
+ BR(KF,KH) = (A2 - A1
+ & + (1.0-4.0*TR+3.0*TR*TR)*DA1 + (3.0*TR-2.0)*TR*DA2 )/DRL
+ BRF(KF,KH) = (AF2 - AF1
+ & + (1.0-4.0*TR+3.0*TR*TR)*DAF1 + (3.0*TR-2.0)*TR*DAF2 )/DRL
+ BRH(KF,KH) = (AH2 - AH1
+ & + (1.0-4.0*TR+3.0*TR*TR)*DAH1 + (3.0*TR-2.0)*TR*DAH2 )/DRL
+ BRFH(KF,KH) = (AFH2 - AFH1
+ & + (1.0-4.0*TR+3.0*TR*TR)*DAFH1 + (3.0*TR-2.0)*TR*DAFH2)/DRL
+C
+ 205 CONTINUE
+ 20 CONTINUE
+C
+C---- evaluate spline in HL at the two FL-interval endpoints
+ DO 30 KF=1, 2
+ B1 = B (KF,1)
+ BR1 = BR (KF,1)
+ BF1 = BF (KF,1)
+ BH1 = BH (KF,1)
+ BRF1 = BRF (KF,1)
+ BRH1 = BRH (KF,1)
+ BFH1 = BFH (KF,1)
+ BRFH1 = BRFH(KF,1)
+C
+ B2 = B (KF,2)
+ BR2 = BR (KF,2)
+ BF2 = BF (KF,2)
+ BH2 = BH (KF,2)
+ BRF2 = BRF (KF,2)
+ BRH2 = BRH (KF,2)
+ BFH2 = BFH (KF,2)
+ BRFH2 = BRFH(KF,2)
+C
+ DB1 = DHL*BH1 - B2 + B1
+ DB2 = DHL*BH2 - B2 + B1
+ DBR1 = DHL*BRH1 - BR2 + BR1
+ DBR2 = DHL*BRH2 - BR2 + BR1
+ DBF1 = DHL*BFH1 - BF2 + BF1
+ DBF2 = DHL*BFH2 - BF2 + BF1
+ DBRF1 = DHL*BRFH1 - BRF2 + BRF1
+ DBRF2 = DHL*BRFH2 - BRF2 + BRF1
+C
+C------ set AI, dAI/dRL, dAI/dFL
+ C(KF) = (1.0-TH)* B1 + TH* B2
+ & + ((1.0-TH)*DB1 - TH*DB2 )*(TH-TH*TH)
+ CR(KF) = (1.0-TH)* BR1 + TH* BR2
+ & + ((1.0-TH)*DBR1 - TH*DBR2 )*(TH-TH*TH)
+ CF(KF) = (1.0-TH)* BF1 + TH* BF2
+ & + ((1.0-TH)*DBF1 - TH*DBF2 )*(TH-TH*TH)
+ CRF(KF) = (1.0-TH)* BRF1 + TH* BRF2
+ & + ((1.0-TH)*DBRF1 - TH*DBRF2)*(TH-TH*TH)
+C
+C------ also, the HL derivatives of the quantities above
+ CH(KF) = (B2 - B1
+ & + (1.0-4.0*TH+3.0*TH*TH)*DB1 + (3.0*TH-2.0)*TH*DB2 )/DHL
+ CRH(KF) = (BR2 - BR1
+ & + (1.0-4.0*TH+3.0*TH*TH)*DBR1 + (3.0*TH-2.0)*TH*DBR2 )/DHL
+ CFH(KF) = (BF2 - BF1
+ & + (1.0-4.0*TH+3.0*TH*TH)*DBF1 + (3.0*TH-2.0)*TH*DBF2 )/DHL
+ CRFH(KF) = (BRF2 - BRF1
+ & + (1.0-4.0*TH+3.0*TH*TH)*DBRF1 + (3.0*TH-2.0)*TH*DBRF2)/DHL
+C
+ 30 CONTINUE
+C
+C---- evaluate cubic in FL
+ C1 = C (1)
+ CR1 = CR (1)
+ CF1 = CF (1)
+ CH1 = CH (1)
+ CRF1 = CRF (1)
+ CRH1 = CRH (1)
+ CFH1 = CFH (1)
+ CRFH1 = CRFH(1)
+C
+ C2 = C (2)
+ CR2 = CR (2)
+ CF2 = CF (2)
+ CH2 = CH (2)
+ CRF2 = CRF (2)
+ CRH2 = CRH (2)
+ CFH2 = CFH (2)
+ CRFH2 = CRFH(2)
+C
+ DC1 = DFL*CF1 - C2 + C1
+ DC2 = DFL*CF2 - C2 + C1
+ DCH1 = DFL*CFH1 - CH2 + CH1
+ DCH2 = DFL*CFH2 - CH2 + CH1
+ DCR1 = DFL*CRF1 - CR2 + CR1
+ DCR2 = DFL*CRF2 - CR2 + CR1
+ DCRH1 = DFL*CRFH1 - CRH2 + CRH1
+ DCRH2 = DFL*CRFH2 - CRH2 + CRH1
+C
+C---- set AI, dAI/dRL, dAI/dHL
+ AI = (1.0-TF)* C1 + TF* C2
+ & + ((1.0-TF)*DC1 - TF*DC2 )*(TF-TF*TF)
+ AI_RL = (1.0-TF)* CR1 + TF* CR2
+ & + ((1.0-TF)*DCR1 - TF*DCR2 )*(TF-TF*TF)
+ AI_HL = (1.0-TF)* CH1 + TF* CH2
+ & + ((1.0-TF)*DCH1 - TF*DCH2 )*(TF-TF*TF)
+C
+C---- also, the FL derivatives of the quantities above
+ AI_FL = (C2 - C1
+ & + (1.0-4.0*TF+3.0*TF*TF)*DC1 + (3.0*TF-2.0)*TF*DC2 )/DFL
+ AIF_RL = (CR2 - CR1
+ & + (1.0-4.0*TF+3.0*TF*TF)*DCR1 + (3.0*TF-2.0)*TF*DCR2 )/DFL
+ AIF_HL = (CH2 - CH1
+ & + (1.0-4.0*TF+3.0*TF*TF)*DCH1 + (3.0*TF-2.0)*TF*DCH2 )/DFL
+C
+ AIF_FL = ((6.0*TF-4.0)*DC1 + (6.0*TF-2.0)*DC2 )/DFL**2
+C
+C
+C---- convert derivatives wrt to spline coordinates (RL,FL,HL) into
+C- derivatives wrt input variables (Rtheta,f,H)
+ AI_R = (AI_RL + 0.5*AI_FL) / (AL10 * RSP)
+ AI_F = (AI_FL ) / (AL10 * FSP)
+ AI_H = AI_HL
+C
+ AIF_R = (AIF_RL + 0.5*AIF_FL) / (AL10**2 * FSP*RSP)
+ AIF_F = (AIF_FL - AL10*AI_FL) / (AL10**2 * FSP*FSP)
+ AIF_H = AIF_HL / (AL10 * FSP )
+C
+C---- if we're within the spline data space, the derivatives are valid
+ IF(OK) RETURN
+C
+C---- if not, the ai value is clamped, and its derivatives are zero
+ AI_R = 0.0
+ AI_F = 0.0
+ AI_H = 0.0
+C
+ AIF_R = 0.0
+ AIF_F = 0.0
+ AIF_H = 0.0
+C
+ RETURN
+ END
+