GTEDGE/Distr.for at master · gt-frc/GTEDGE · GitHub

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
      SUBROUTINE DISTR
      INCLUDE 'Soldiv.fi'
C           DISTRIBUTES N & T

      AP = 39.44*RMAJOR*AMINOR*SQRT((1.+ELONG**2)/2.)
c     ap = 39.44*rmajor*aminor*sqrt(elong)
      VP = 19.72*RMAJOR*(AMINOR**2)*0.5*(1+ELONG**2)
c     VP = 19.72*RMAJOR*(AMINOR**2)*ELONG
      PSEP = XK*XNSEP*TSEP
      CBAL = CBALLOON
      FAC = CBAL*((BFIELD**2)/(2.52E-6))*DELTB/((Q95**2)*RMAJOR*PSEP)
c     XNAVPED = (XNAV/XNEL)/APED
c      XNCTRPED = (XNAVPED*(1.+ALPHAN)-ALPHAN)

C     XN0 = (1.+ALPHAN)*XNAV/(1.0 + ALPHAN/XNCTRPED)
C     ZT = 1.0/TCTRPED
C     ZN = 1.0/XNCTRPED
C     A1 = ((1.0+ALPHAN*ZN)/(1.0+ALPHAN))
C     A2 = ((1.0-ZT)*(1.0-ZN)/(ALPHAT2+ALPHAN+1.0))
C     A3 = ZT*(1.0-ZN)/(ALPHAN+1.0)
C     A4 = ZN*(1.0-ZT)/(ALPHAT2+1.0)
C     B1 = A1 + A2 + A3 + ZN*ZT
C     T0 = A1*TAV/B1
C     X = (1.+ALPHAN+ALPHAT2)*XNAV*TAV/(XNSEP*TSEP*(1.+FAC))
C     Y = (ALPHAT2/(1.+ALPHAN))*(XNCTRPED+ALPHAN*(2.+ALPHAN+ALPHAT2)/
C     2    (1.+ALPHAT2))
C     Z = XNCTRPED + ALPHAN/(1.+ALPHAT2)

c     RATNLINEAV = XNEL/XNAV
c     this expression valid for -1<alphn<1, including limits 1
      if(jjoptped.eq.10) xnctrped = xn0/xnped
      if(jjoptped.eq.10) tctrped = t0/tped
      y = xnctrped - 1.

      ratnlineav = (1.+y*(1.-alphan/3. + alphan*(alphan-1.)/10. -
     2            alphan*(alphan-1.)*(alphan-2.)/42.))/
     3           (1. + y/(1.+alphan))
      IF(JJOPTPED.NE.9) GOTO 10
      XNPED = APED*XNAV*RATNLINEAV
      write(*,*) 'ratnlineav =', ratnlineav
c     ALPHAN = (XNCTRPED*(APED*RATNLINEAV)-1.)/(1.-APED*RATNLINEAV)
      XN0 = (1.+ALPHAN)*XNAV - ALPHAN*XNPED
      TPED = XNSOL*TSOL*(1.+FAC)/XNPED
      alphat = alphat2
      X1 = (TAV/(APED*RATNLINEAV)) - ALPHAT*TPED*(1./(1.+ALPHAT) +
     2      (XNCTRPED-1.)/((1.+ALPHAN)*(1.+ALPHAN+ALPHAT)))
      X2 = 1./(1.+ALPHAT) + (XNCTRPED-1.)/(1.+ALPHAN+ALPHAT)
      T0 = X1/X2
      GOTO 15
 10   IF(JJOPTPED.NE.10) GOTO 15
C     XNPED = sqrt(1.+fac)*XNSOL
c     ALPHAN = (XNCTRPED*(APED*RATNLINEAV)-1.)/(1.-APED*RATNLINEAV)
C     IF(ITERN.GT.15)
C     2           ALPHAN = (XNCTRPED*(XNPED/XNAV)-1.)/(1.-XNPED/XNAV)
      XN0 = (1.+ALPHAN)*XNAV - ALPHAN*XNPED
C     TPED = sqrt(1.+fac)*TSOL
      alphat = alphat2
c     X1 = (TAV/(APED*RATNLINEAV)) - ALPHAT*TPED*(1./(1.+ALPHAT) +
c     2      (XNCTRPED-1.)/((1.+ALPHAN)*(1.+ALPHAN+ALPHAT)))
c           X2 = 1./(1.+ALPHAT) + (XNCTRPED-1.)/(1.+ALPHAN+ALPHAT)
      x1 = xnav*tav - (xn0-xnped)*tped*((1./(1.+alphan))-
     2     (1./(1.+alphan+alphat))) - xnped*tped*(1.- (1./(1.+alphat)))
      x2 = ((xn0-xnped)/(1.+alphan+alphat)) + xnped/(1.+alphat)
      T0 = X1/X2

 15   CONTINUE
C     T0 = TPED*(X - Y)/Z
      DELC = 1./56.

C     SE = SQRT(0.5*(1.+ELONG**2))
C     DELC = (1.0-DELTB/(AMINOR*SE))/56.0
      RHO(1) = 0.5*DELC
      TE(1) = TPED + A2*(T0-TPED)*((1.0 - (RHO(1)**2))**ALPHAT2) +
     2                    A1*(T0-TPED)*((1.0 - (RHO(1)))**ALPHAT1)
      XNC(1) = XNPED + (XN0-XNPED)*((1.0 - (RHO(1)**2))**ALPHAN)
      XNEL = XNC(1)*DELC
      if (ITERN.eq.100) then
        write(9016,'(A1,A10,A7,A10)') 'n','Rho','T','ne'
        write(9016,'(I2,3E10.3)') 1,RHO(1),TE(1),XNC(1)
      endif
c         write(*,*) 'TE(1) =', TE(1)
      DO 35 I = 2,55
        RHO(I) = RHO(I-1) + DELC
        TE(I) = TPED + A2*(T0-TPED)*((1.0 - (RHO(I)**2))**ALPHAT2) +
     2                      A1*(T0-TPED)*((1.0 - (RHO(I)))**ALPHAT1)
        XNC(I) = XNPED + (XN0-XNPED)*((1.0 - (RHO(I)**2))**ALPHAN)
        XNEL = XNEL + XNC(I)*DELC
        IF (ITERN.eq.100) THEN
          write(9016,'(I2,3E10.3)') I,RHO(I),TE(I),XNC(I)
        endif
c        write(*,*) 'I =', I
c        write(*,*) 'TE(I) =', TE(I)
 35   CONTINUE
      TE(56) = TPED
      XNC(56) = XNPED
      TE(57) = TBAR
      XNC(57) = XNBAR
C     AVG TEMP IN SOL = TSEP(1. - EXP(-DELN/DELT))
      TE(58) = TSOL*(1. - EXP(-1.*DELRATNT))
C     AVG DENSITY IN SOL =NSEP( 1. - EXP(-1.))= 0.632
      XNC(58) = 0.632*XNSOL
      XNEL = XNEL + XNC(56)*DELC
c     RELATING XNEL TO XNAV
c     XNEL = (1.333/(1.0+1.0/XNCTRPED))*XNAV
c     XNEL = XNAV

      IF (ITERN.eq.100) THEN
        write(9016,'(I3,3E10.3)') 56,RHO(55)+DELC,TE(56),XNC(56)
        write(9016,'(I3,3E10.3)') 57,RHO(55)+2*DELC,TE(57),XNC(57)
        write(9016,'(I3,3E10.3)') 58,RHO(55)+3*DELC,TE(58),XNC(58)
      ENDIF


      xnel = ratNlineav*xnav

      JIT = 0
      DO 50 I = 1,55
        IF(TE(I).LE.0.) TE(I) = TPED
        IF(XNC(I).LE.0.) JIT = 1
        IF(XNC(I).LE.0.) XNC(I) = 1.E17
 50   CONTINUE

      RETURN
      END