-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathDivstab.for
105 lines (103 loc) · 4.25 KB
/
Divstab.for
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
SUBROUTINE DIVSTAB
INCLUDE 'Soldiv.fi'
C DIVERTOR THERMAL STABILITY ANALYSIS
C SETS UP DISPERSION RELATION FOR N-V-T LINEARIZED EQUATIONS
C FINDS COMPLEX EIGENVALUES BY CALLING GVCCG
INTEGER LDA,LDB,LDEVEC,NN
PARAMETER (NN=3, LDA=NN, LDB=NN, LDEVEC=NN)
INTEGER II,NOUT
REAL GPICG,PI,ZTT,ZTN,ZTV,ZVV,ZVN,ZVT,ZNN,ZNV,ZNT,ZMACH2,
2 RNN,RNV,RNT,RVN,RVV,RVT,RTN,RTV,RTT
COMPLEX ADIV(LDA,NN),OMEGA(NN),BDIV(LDB,NN),BETA(NN),
2 EVAL(NN), EVEC(LDEVEC,NN)
EXTERNAL GPICG, GVCCG, UMACH
C EVALUATE TERMS FOR STABILITY ANALYSIS
c KAPPA = KAPPA*REDKAPPA2/REDKAPPA
DO 500 M = 1,1
DXLDIV = XLPAR-XLPERP-DELLT
XKLL = M/(XLPAR)
xkll = 1./dxldiv
XMDIV = 0.5*(1.+XMSOL)
CONDUCT = XKAPPA*((TD**2.5)*DELLT/XND + (TDIV**2.5)*DXLDIV/XNDIV
2 + (TSOL**2.5)*XLPERP/XNSOL)*(XKLL**2)/XK
VM2AV = CSD*(DELLT + (XMDIV**3)*DXLDIV + (XMSOL**3)*XLPERP)
VAV = CSD*(DELLT + (XMDIV)*DXLDIV + (XMSOL)*XLPERP)
ZMACH2 = (DELLT + (XMDIV**3)*DXLDIV + (XMSOL**3)*XLPERP)
CONDUCT1 = XKLL*(3.*VAV + 2.*XKAPPA*((TSEP**3.5)-(TD**3.5))/
2 (XK*XNSEP*TSEP))
DNDPSID = 0.25*CSD*LOG(XND/XNSEP)
HEATIN = FLUXHEAT*XLPERP/(XNSEP*XK*TSEP*DELEA)
FREQIONRECD = FREQIOND - FREQRECD
FREQIONRECM = FREQIONM - FREQRECD
C EVALUATE REAL & IMAGINARY MATRIX ELEMENTS
RNN = FLUXPART*XLPERP/(XNSEP*DELNV)
RNV = RNN + FREQIONREC
RNT = -1.*(DFREQIONDT - DFREQRECDT)
RVN = DFREQIONDN - (XMSOL**2)*RNN
RVV = -2.*DNDPSID + 3.*FREQIONM -2.*FREQRECD + FREQATM
RVT = 0.25*CSD*LOG(XND*TD/(XNSEP*TSEP)) + DFREQIONDTM + DFREQATDTM
RTN = -2.*DNDPSID + 2.*(1.-3.*(XMSOL**2))*RNN + 2.*RAD1 +
2 (1.5*FREQAT - 5.*FREQATM) +
3 (5.*FREQIONREC + FREQIONE - 7.*FREQIONM + 2.*FREQRECD) +
4 ((2.+EIOND/TD)*DFREQIONDN + 4.*DFREQRECDN)
RTV = -3.*DNDPSID +(2.-8.*(XMSOL**2))*RNN -7.*(FREQATM + FREQRECD)
2 + (2.*FREQIONREC - 8.*FREQIONRECM)
RTT = CONDUCT + 3.5*HEATIN - (16.5-10.5*(XMSOL**2))*RNN +
1 5.*DNDPSID - RADDIV -
2 (3.75*FREQAT - 14.5*FREQATM) +(1.5*DFREQATDT-2.*DFREQATDTM)-
3 (16.5*FREQIONREC + 3.5*FREQIONE - 10.5*FREQIONM-4.*FREQRECD)
4 +(3.*DFREQIONDT-DFREQIONDTM-3.5*DFREQIONDTE) - 4.*DFREQRECDT
ZTT = XKLL*(CONDUCT1 + 3.*VAV)
ZTN = 0.
ZTV = 2.*XKLL*VAV
ZVT = XKLL*VAV
ZVV = XKLL*VM2AV
ZVN = XKLL*VAV
ZNT = 0.
ZNV = XKLL*VAV
ZNN = XKLL*VAV
C SET UP MATRICES
ADIV(1,1) = CMPLX(RNN,ZNN)
ADIV(1,2) = CMPLX(RNV,ZNV)
ADIV(1,3) = CMPLX(RNT,ZNT)
ADIV(2,1) = CMPLX(RVN,ZVN)
ADIV(2,2) = CMPLX(RVV,ZVV)
ADIV(2,3) = CMPLX(RVT,ZVT)
ADIV(3,1) = CMPLX(RTN,ZTN)
ADIV(3,2) = CMPLX(RTV,ZTV)
ADIV(3,3) = CMPLX(RTT,ZTT)
BDIV(1,1) = -1.0*CMPLX(1.0,0.0)
BDIV(1,2) = -1.0*CMPLX(0.0,0.0)
BDIV(1,3) = -1.0*CMPLX(0.0,0.0)
BDIV(2,1) = -1.0*CMPLX(0.0,0.0)
BDIV(2,2) = -1.0*CMPLX(ZMACH2,0.0)
BDIV(2,3) = -1.0*CMPLX(0.0,0.0)
BDIV(3,1) = -1.0*CMPLX(0.0,0.0)
BDIV(3,2) = -1.0*CMPLX(0.0,0.0)
BDIV(3,3) = -1.0*CMPLX(3.0,0.0)
C CALCULATE EIGENVALUES AND EIGENVECTORS
CALL GVCCG(NN,ADIV,LDA,BDIV,LDB,OMEGA,BETA,EVEC,LDEVEC)
DO 400 II=1,NN
EVAL(II) = OMEGA(II)/BETA(II)
400 CONTINUE
C PRINT RESULTS
c WRITE(6,999) -0.3333*RTT,-1.0*RNN,EVAL
WRITE(6,997) EVAL
OPEN(111,FILE='SOLDIV.TXT',STATUS='UNKNOWN')
WRITE(111,'(1X,29A)') 'DIVERTOR INSTABILITY RESULTS'
WRITE(111,998) -0.3333*RTT,-1.0*RNN,-1.0*RVV/ZMACH2
WRITE(111,997) EVAL
500 CONTINUE
WRITE(111,'(1X,A)')'DENSITY-FLOW-TEMPERATURE PERTS EIGENFUNCTIONS'
WRITE(111,1000) EVEC
C CALL UMACH (2,NOUT)
C CALL WRCRN ('EVAL',1,NN,EVAL,1,0)
C CALL WRCRN ('EVEC',NN,NN,EVEC,LDEVEC,0)
997 FORMAT (1X,'REAL/IMAG OMEG:#1=',2E9.3,1X,'#2=',2E9.3,1X,
2'#3=',2E9.3)
998 FORMAT (1X,'TEMP GR=',E9.3,1X,'DENSITY GR=',E9.3,1X,
2'FLOW GR=',E9.3)
999 FORMAT (1X,'DIV OM=',8E9.3)
1000 FORMAT (6E9.3)
RETURN
END