-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathNeutral3.for
817 lines (767 loc) · 30.3 KB
/
Neutral3.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
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
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
SUBROUTINE NEUTRAL3
INCLUDE 'Soldiv.fi'
DIMENSION AMAT(12,12),SOURCE (12), IPIV(12)
real svelsol,svcxsol,alphar,vodiv,ascatd,yionrec1,
1 factorrec,factorion,zxpt
C TEP NEUTRAL MODEL 7/97
C SET DELN = DELNV
DELNREAL = DELN
DELN = DELNV
RATNENI2D=(1.+ IZINJECT*FZINJECT+4.*FBE+2.*FHE+
2 IZINTRIN*FZINTRIN)**2
CNEUT = 1.
IF(IOPTELN.EQ.0) CNEUT = 0.
C HIGH RECYCLING REGION IN FRONT OF DIVERTOR PLATE
TDD = TD
IF(TDD.LT.1.E-1) TDD = 1.05E-1
IF(TDD.GT.1E3) TDD = .95E3
XNDD = XND
IF(XND.GT.1E22) XNDD = 0.95E22
IF(XND.LT.1E16) XNDD = 1.1E16
TND = TDD
IF(TND.GE.1000) TND = 995.
CALL INTERP(TDD,TDD,TND,XNDD)
SVELD = SEL(1)
SVELDN= SELN(1)
SVCXD = SCX(1)
SVATD = SEL(1) + CNEUT*SELN(1)*XNOD/XND + SCX(1)
C RADIAL INTEGRATION OF IONIZATION & RECOMBINATION INTEGRALS,
C LEADING TO EFFECTIVE CROSS SECTIONS FOR RECYCLING REGION 4/22/98
5 SVIOND = SION(1)
SVRECD = RECOM(1)
IF(TD.LT.5.0) GOTO 20
FACTORREC = 0.25
FACTORION = 0.25
TDD = TD*EXP(-1.*DELRATNT/2.)
IF(TDD.LT.1.E-1) TDD = 1.05E-1
IF(TDD.GT.1E3) TDD = .95E3
XNDD = XND*EXP(-1./2.)
IF(XND.GT.1E22) XNDD = 0.95E22
IF(XND.LT.1E16) XNDD = 1.1E16
TND = TDD
IF(TND.GE.100) TND = 95.
CALL INTERP(TDD,TDD,TND,XNDD)
FACTORION = FACTORION + 0.5*SION(1)/SVIOND
FACTORREC = FACTORREC + 0.5*RECOM(1)/SVRECD
TDD = TD*EXP(-1.*DELRATNT)
IF(TDD.LT.1.E-1) TDD = 1.05E-1
IF(TDD.GT.1E3) TDD = .95E3
XNDD = XND*EXP(-1.)
IF(XND.GT.1E22) XNDD = 0.95E22
IF(XND.LT.1E16) XNDD = 1.1E16
TND = TDD
IF(TND.GE.100) TND = 95.
CALL INTERP(TDD,TDD,TND,XNDD)
FACTORION = FACTORION + 0.25*SION(1)/SVIOND
FACTORREC = FACTORREC + 0.25*RECOM(1)/SVRECD
GOTO 30
C MORE DETAILED RADIAL INTEGRATION FOR TD < 5 eV
20 FACTORREC = 0.0
FACTORION = 0.0
DL = -0.5*DELN/10.
DO 25 J=1,10
DL = DL + DELN/10.
TDD = TD*EXP(-1.*DELRATNT*DL/DELN)
IF(TDD.LT.1.E-1) TDD = 1.05E-1
IF(TDD.GT.1E3) TDD = .95E3
XNDD = XND*EXP(-1.*DL/DELN)
IF(XND.GT.1E22) XNDD = 0.95E22
IF(XND.LT.1E16) XNDD = 1.1E16
TND = TDD
IF(TND.GE.1000) TND = 995.
CALL INTERP(TDD,TDD,TND,XNDD)
FACTORION = FACTORION + 0.1*SION(1)/SVIOND
FACTORREC = FACTORREC + 0.1*RECOM(1)/SVRECD
25 CONTINUE
30 CONTINUE
C FACTORION = 1.0
C FACTORREC = 1.0
SVIOND = SVIOND*FACTORION
SVRECD = SVRECD*FACTORREC
YATD = SVATD
SVTOTD = SVIOND + YATD
RSCATD = YATD/(YATD + SVIOND)
ASCATD = YATD/(YATD + SVIOND)
C ASSUME REFLECTED ATOMS HAVE 1/2 INCIDENT ION ENERGY, NEUTRALS
C FROM REEMITTED MOLECULES HAVE 2 eV, IONS RE-EMITTED AS MOLECULES
CALL MOLECULE
EMOL = 2.
SFRAC = GAMZEROM/(GAMZERO/RWDP)
C ASSUMPTION--INCIDENT IONS RECYCLE AS MOLECULES--DEEP PENETRATION
TODD = (SFRAC*(0.5*(1.-FMOL)*TD + FMOL*EMOL) + (1.-SFRAC)*EMOL)
IF(SFRAC.LE.0.0) TODD = 2.0
VOD = SQRT(XK*TODD/XMASS)
IF(XND.LE.1.E21)
2 EIOND = 17.5 + (5.+37.5/TD)*LOG10(1.E21/(XND))
IF(XND.GT.1.E21)
2 EIOND = (30.6 - 16.4*EXP(-5.E19/XND))*
3 EXP(5.45/(TD*EXP((XND/1.37E20)**0.26)))
C TEP TRANSMISSION & ESCAPE PARAMETERS
C FIRST-FLIGHT ESCAPE PROBABILITY
C ******OUTBOARD*****************
C NEUTRALS FROM PLATE OR PL OR PF
YMFP = 1./(XND*SVTOTD/VOD)
C NEUTRALS FROM DIVERTOR PLASMA
VODIV = SQRT(XK*TDIV/XMASS)
YMFP1 = 1./(XND*SVTOTD/VODIV)
DELNS1 = EPDIV1*DELN1/SIN(THETA1)
X = 2.*(DELLT1*BETAG + DELNS1)
RATIOVOLSURF = EPDIV1*DELN1*DELLT1*BETAG/X
XX = 4.*RATIOVOLSURF/YMFP
POD1 = (1. - (1.+XX/2.09)**(-2.09))/XX
C TOTAL ESCAPE PROBABILITY
PESCAPER1 = POD1/(1. - ASCATD*(1. - POD1))
C FIRST-FLIGHT TRANSMISSION PROBABILITIES
Z = BETAG*DELLT1
YPERP = BETAG*DELLT1*SIN(THETA1)
YPAR = EPDIV1*DELN1
PHI = 3.1416 - THETA1
CALL TRANSM(1,0,YMFP,DELNS1,Z,THETA1,0.0,0.0,DELNS1,TOD2PF1)
CALL TRANSM(2,0,YMFP,DELNS1,DELNS1,THETA1,0.0,YPERP,DELNS1,TOD2P1)
CALL TRANSM(1,0,YMFP,DELNS1,Z,PHI,0.0,0.0,DELNS1,TOD2PL1)
CALL TRANSM(1,0,YMFP,Z,DELNS1,THETA1,0.0,0.0,0.0,TOPF2D1)
CALL TRANSM(2,0,YMFP,Z,Z,PHI,0.0,YPAR,0.0,TOPF2PL1)
TOPL2PF1 = TOPF2PL1
TOPL2P1 = TOPF2D1
CALL TRANSM(1,0,YMFP,Z,DELNS1,PHI,0.0,0.0,0.0,TOPF2P1)
TOPL2D1 = TOPF2P1
CALL TRANSM(1,0,YMFP1,DELNS1,Z,THETA1,0.0,0.0,0.0,TOP2PL1)
CALL TRANSM(1,0,YMFP1,DELNS1,Z,PHI,0.0,0.0,0.0,TOP2PF1)
CALL TRANSM(2,0,YMFP1,DELNS1,DELNS1,THETA1,0.0,YPERP,0.0,TOP2D1)
C TOTAL FIRST-FLIGHT TRANSMISSION PROBABILITIES
TOD1 = TOD2PL1 + TOD2PF1 + TOD2P1
TOP1 = TOP2PL1 + TOP2PF1 + TOP2D1
TOPF1 = TOPF2P1 + TOPF2D1 + TOPF2PL1
TOPL1 = TOPL2P1 + TOPL2D1 + TOPL2PF1
C DIRECTIONAL ESCAPE FACTORS
X = (DELLT1*BETAG + DELNS1)
FPF = 0.5
FPD = 0.5
XLAMPF1 = DELLT1*BETAG*FPF/X
XLAMPL1 = DELLT1*BETAG*(1.-FPF)/X
XLAMD1 = DELNS1*FPD/X
XLAMP1 = DELNS1*(1.-FPD)/X
C DIRECTIONAL ESCAPE PROBABILITIES
XX = ASCATD*PESCAPER1
GPF1 = XX*XLAMPF1
GPL1 = XX*XLAMPL1
GD1 = XX*XLAMD1
GP1 = XX*XLAMP1
C TRANSMISSION PROBABILITIES INCLUDING REFLECTION
TPF2P1 = TOPF2P1 + (1.-TOPF2P1-TOPF2PL1-TOPF2D1)*GP1
TPF2D1 = TOPF2D1 + (1.-TOPF2P1-TOPF2PL1-TOPF2D1)*GD1
TPF2PL1 = TOPF2PL1 + (1.-TOPF2P1-TOPF2PL1-TOPF2D1)*GPL1
TPF2PF1 = (1.-TOPF2P1-TOPF2PL1-TOPF2D1)*GPF1
TPL2P1 = TOPL2P1 + (1.-TOPL2P1-TOPL2PF1-TOPL2D1)*GP1
TPL2D1 = TOPL2D1 + (1.-TOPL2P1-TOPL2D1-TOPL2PF1)*GD1
TPL2PF1 = TOPL2PF1 + (1.-TOPL2P1-TOPL2D1-TOPL2PF1)*GPF1
TPL2PL1 = (1.-TOPL2P1-TOPL2D1-TOPL2PF1)*GPL1
TP2PF1 = TOP2PF1 + (1.-TOP2PF1-TOP2PL1-TOP2D1)*GPF1
TP2PL1 = TOP2PL1 + (1.-TOP2PF1-TOP2PL1-TOP2D1)*GPL1
TP2D1 = TOP2D1 + (1.-TOP2PF1-TOP2PL1-TOP2D1)*GD1
TP2P1 = (1.-TOP2PF1-TOP2PL1-TOP2D1)*GP1
TD2PF1 = TOD2PF1 + (1.-TOD2PF1-TOD2PL1-TOD2P1)*GPF1
TD2PL1 = TOD2PL1 + (1.-TOD2PF1-TOD2PL1-TOD2P1)*GPL1
TD2P1 = TOD2P1 + (1.-TOD2PF1-TOD2PL1-TOD2P1)*GP1
TD2D1 = (1.-TOD2PF1-TOD2PL1-TOD2P1)*GD1
C ************INBOARD****************************
C NEUTRALS FROM PLATE OR PL OR PF
YMFP = 1./(XND*SVTOTD/VOD)
C NEUTRALS FROM DIVERTOR PLASMA
VODIV = SQRT(XK*TDIV/XMASS)
YMFP1 = 1./(XND*SVTOTD/VODIV)
DELNS2= EPDIV2*DELN2/SIN(THETA2)
X = 2.*(DELLT2*BETAG + DELNS2)
RATIOVOLSURF = EPDIV2*DELN2*DELLT2*BETAG/X
XX = 4.*RATIOVOLSURF/YMFP
POD2 = (1. - (1.+XX/2.09)**(-2.09))/XX
C TOTAL ESCAPE PROBABILITY
PESCAPER2 = POD2/(1. - ASCATD*(1. - POD2))
C FIRST-FLIGHT TRANSMISSION PROBABILITIES
Z = BETAG*DELLT1
YPERP = BETAG*DELLT2*SIN(THETA2)
YPAR = EPDIV2*DELN2
PHI = 3.1416 - THETA2
CALL TRANSM(1,0,YMFP,DELNS2,Z,THETA2,0.0,0.0,DELNS2,TOD2PF2)
CALL TRANSM(2,0,YMFP,DELNS2,DELNS2,THETA2,0.0,YPERP,DELNS2,TOD2P2)
CALL TRANSM(1,0,YMFP,DELNS2,Z,PHI,0.0,0.0,DELNS2,TOD2PL2)
CALL TRANSM(1,0,YMFP,Z,DELNS2,THETA2,0.0,0.0,0.0,TOPF2D2)
CALL TRANSM(2,0,YMFP,Z,Z,PHI,0.0,YPAR,0.0,TOPF2PL2)
TOPL2PF2 = TOPF2PL2
TOPL2P2 = TOPF2D2
CALL TRANSM(1,0,YMFP,Z,DELNS2,PHI,0.0,0.0,0.0,TOPF2P2)
TOPL2D2 = TOPF2P2
CALL TRANSM(1,0,YMFP1,DELNS2,Z,THETA2,0.0,0.0,0.0,TOP2PL2)
CALL TRANSM(1,0,YMFP1,DELNS2,Z,PHI,0.0,0.0,0.0,TOP2PF2)
CALL TRANSM(2,0,YMFP1,DELNS2,DELNS2,THETA2,0.0,YPERP,0.0,TOP2D2)
C TOTAL FIRST-FLIGHT TRANSMISSION PROBABILITIES
TOD2 = TOD2PL2 + TOD2PF2 + TOD2P2
TOP2 = TOP2PL2 + TOP2PF2 + TOP2D2
TOPF2 = TOPF2P2 + TOPF2D2 + TOPF2PL2
TOPL2 = TOPL2P2 + TOPL2D2 + TOPL2PF2
C DIRECTIONAL ESCAPE FACTORS
X = (DELLT2*BETAG + DELNS2)
FPF = 0.5
FPD = 0.5
XLAMPF2 = DELLT2*BETAG*FPF/X
XLAMPL2 = DELLT2*BETAG*(1.-FPF)/X
XLAMD2 = DELNS2*FPD/X
XLAMP2 = DELNS2*(1.-FPD)/X
C DIRECTIONAL ESCAPE PROBABILITIES
XX = ASCATD*PESCAPER2
GPF2 = XX*XLAMPF2
GPL2 = XX*XLAMPL2
GD2 = XX*XLAMD2
GP2 = XX*XLAMP2
C TRANSMISSION PROBABILITIES INCLUDING REFLECTION
TPF2P2 = TOPF2P2 + (1.-TOPF2P2-TOPF2PL2-TOPF2D2)*GP2
TPF2D2 = TOPF2D2 + (1.-TOPF2P2-TOPF2PL2-TOPF2D2)*GD2
TPF2PL2 = TOPF2PL2 + (1.-TOPF2P2-TOPF2PL2-TOPF2D2)*GPL2
TPF2PF2 = (1.-TOPF2P2-TOPF2PL2-TOPF2D2)*GPF2
TPL2P2 = TOPL2P2 + (1.-TOPL2P2-TOPL2PF2-TOPL2D2)*GP2
TPL2D2 = TOPL2D2 + (1.-TOPL2P2-TOPL2D2-TOPL2PF2)*GD2
TPL2PF2 = TOPL2PF2 + (1.-TOPL2P2-TOPL2D2-TOPL2PF2)*GPF2
TPL2PL2 = (1.-TOPL2P2-TOPL2D2-TOPL2PF2)*GPL2
TP2PF2 = TOP2PF2 + (1.-TOP2PF2-TOP2PL2-TOP2D2)*GPF2
TP2PL2 = TOP2PL2 + (1.-TOP2PF2-TOP2PL2-TOP2D2)*GPL2
TP2D2 = TOP2D2 + (1.-TOP2PF2-TOP2PL2-TOP2D2)*GD2
TP2P2 = (1.-TOP2PF2-TOP2PL2-TOP2D2)*GP2
TD2PF2 = TOD2PF2 + (1.-TOD2PF2-TOD2PL2-TOD2P2)*GPF2
TD2PL2 = TOD2PL2 + (1.-TOD2PF2-TOD2PL2-TOD2P2)*GPL2
TD2P2 = TOD2P2 + (1.-TOD2PF2-TOD2PL2-TOD2P2)*GP2
TD2D2 = (1.-TOD2PF2-TOD2PL2-TOD2P2)*GD2
C SCRAPE-OFF LAYER, TRANSPORT BARRIER & PEDESTAL
C SCRAPE-OFF LAYER
C FIRST-FLIGHT ESCAPE PROBABILITY
TSOL = TE(58)
VSOL = HVSOL*0.0
XNSOL = XNC(58)
TSOLL = TSOL
IF(TSOL.LT.1.E-1) TSOLL = 1.05E-1
IF(TSOL.GT.1E3) TSOLL = .95E3
XNSOLL = XNSOL
IF(XNSOL.GT.1E22) XNSOLL = 0.95E22
IF(XNSOL.LT.1E16) XNSOLL = 1.1E16
TNSOL = TSOLL
IF(TNSOL.GE.1000) TNSOL = 995.
CALL INTERP(TSOLL,TSOLL,TNSOL,XNSOLL)
SVELSOL = SEL(1)
SVELSOLN= SELN(1)
SVCXSOL = SCX(1)
SVATSOL = SEL(1) + SCX(1)
SVIONSOL = SION(1)
SVTOTSOL = SION(1) + SEL(1) + SCX(1) + CNEUT*SELN(1)*XNOSOL/XNSOL
Y = SVCXSOL + SVELSOL + CNEUT*SELN(1)*XNOSOL/XNSOL
SVRECSOL = RECOM(1)
RSCATSOL = Y/(Y + SVIONSOL)
ASCATSOL = Y/(Y + SVIONSOL)
C AVG TSOL = TSOL*DELT*(1-e(-DELN/DELT))/(1-e(-1))
AVG = DELT*(1.-EXP(-1.*DELN/DELT))/((DELN+DELT)*0.632)
VOSOL = SQRT(XK*TSOL*AVG/XMASS)
EIONSOL = 17.5
IF(XNSOL.LE.1.E21)
2 EIONSOL = 17.5 + (5.+37.5/TSOL)*LOG10(1.E21/XNSOL)
IF(XNSOL.GT.1.E21)
2 EIONSOL = (30.6 - 16.4*EXP(-5.E19/XNSOL))*
3 EXP(5.45/(TSOL*EXP((XNSOL/1.37E20)**0.26)))
C TRANSPORT BARRIER
TBARR = TBAR
IF(TBAR.LT.1.E-1) TBARR = 1.05E-1
IF(TBAR.GT.1E3) TBARR = .95E3
XNBARR = XNBAR
IF(XNBAR.GT.1E22) XNBARR = 0.95E22
IF(XNBAR.LT.1E16) XNBARR = 1.1E16
TNBAR = TBAR
IF(TNBAR.GE.1000) TNBAR = 995.
CALL INTERP(TBARR,TBARR,TNBAR,XNBARR)
SVELBAR = SEL(1)
SVELBARN= SELN(1)
SVCXBAR = SCX(1)
SVATBAR = SEL(1) + SCX(1)
SVIONBAR = SION(1)
SVTOTBAR = SION(1) + SEL(1) + SCX(1) + CNEUT*SELN(1)*XNOBAR/XNBAR
Y = SVCXBAR + SVELBAR + CNEUT*SVELBARN*XNOBAR/XNBAR
ASCAT = Y/(Y + SVIONBAR)
RSCATTB = ASCAT
VOBAR = SQRT(XK*TBAR/XMASS)
EIONBAR = 17.5
IF(XNBAR.LE.1.E21)
2 EIONBAR = 17.5 + (5.+37.5/TBAR)*LOG10(1.E21/XNBAR)
IF(XNBAR.GT.1.E21)
2 EIONBAR = (30.6 - 16.4*EXP(-5.E19/XNBAR))*
3 EXP(5.45/(TBAR*EXP((XNBAR/1.37E20)**0.26)))
C PEDESTAL
Tepedd = 70.
TPEDD = TPED
IF(TPED.LT.1.E-1) TPEDD = 1.05E-1
IF(TPED.GT.1E3) TPEDD = .95E3
XNPEDD = XNPED
IF(XNPED.GT.1E22) XNPEDD = 0.95E22
IF(XNPED.LT.1E16) XNPEDD = 1.1E16
TNPED = TPED
IF(TNPED.GE.1000) TNPED = 995.
CALL INTERP(TePEDD,TPEDD,TNPED,XNPEDD)
SVELPED = SEL(1)
SVCXPED = SCX(1)
SVATPED = SEL(1) + SCX(1)
SVIONPED = SION(1)
SVTOTPED = SION(1) + SEL(1) + SCX(1)
Y = SVCXPED + SVELPED
ASCAT = Y/(Y + SVIONPED)
RSCATPED = ASCAT
VOPED = SQRT(XK*TPED/XMASS)
EIONPED = 17.5
IF(XNPED.LE.1.E21)
2 EIONPED = 17.5 + (5.+37.5/TPED)*LOG10(1.E21/XNPED)
IF(XNPED.GT.1.E21)
2 EIONPED = (30.6 - 16.4*EXP(-5.E19/XNPED))*
3 EXP(5.45/(TPED*EXP((XNPED/1.37E20)**0.26)))
C CORE ALBEDO
IF(ITERN.GT.1) TPLAS = TPED
IF(ITERN.GT.1) XNPLAS = XNPED
c IF(ITERN.GT.1) TPLAS = TPEDEXI
c IF(ITERN.GT.1) XNPLAS = XNPEDEX
TPLASD = TPLAS
IF(TPLASD.LT.1.E-1) TPLASD = 1.05E-1
IF(TPLASD.GT.1E3) TPLASD = .95E3
XNPLASD = XNPLAS
IF(XNPLASD.GT.1E22) XNPLASD = 0.95E22
IF(XNPLASD.LT.1E16) XNPLASD = 1.1E16
TNPLAS = TPLAS
IF(TNPLAS.GE.1000) TNPLAS = 995.
CALL INTERP(TPLAS,TPLASD,TNPLAS,XNPLASD)
SVELPLAS = SEL(1)
SVCXPLAS = SCX(1)
SVATPLAS = SEL(1) + SCX(1)
SVIONPLAS = SION(1)
SVTOTPLAS = SION(1) + SEL(1) + SCX(1)
Y = SVCXPLAS + SVELPLAS
ASCAT = Y/(Y + SVIONPLAS)
VOPLAS = SQRT(XK*TPLAS/XMASS)
EIONPLAS = 13.5
IF(XNPLAS.LE.1.E21)
2 EIONPLAS = 17.5 + (5.+37.5/TPLASD)*LOG10(1.E21/XNPLAS)
IF(XNPLAS.GT.1.E21)
2 EIONPLAS = (30.6 - 16.4*EXP(-5.E19/XNPLAS))*
3 EXP(5.45/(TPLAS*EXP((XNPLAS/1.37E20)**0.26)))
X = SQRT((1./ASCAT)-1.)
ALPHACORE = (1.- 1.155*X)/(1.+1.155*X)
C SOLUTIONS FOR THE SOL-TRANS BARRIER--PEDESTAL
CALL NEUTFRAC
CALL NEUTXPT
C DIVERTOR PLASMA
C FIRST-FLIGHT ESCAPE PROBABILITY
TDIV = HTDIV*TSOL + (1.-HTDIV)*TD
VDIV = HVDIV*0.0 + (1.-HVDIV)*CSD
XNDIV = HNDIV*XNSOL + (1.-HNDIV)*XND
TDIVV = TDIV
IF(TDIV.LT.1.E-1) TDIVV = 1.05E-1
IF(TDIV.GT.1E3) TDIVV = .95E3
XNDIVV = XNDIV
IF(XNDIV.GT.1E22) XNDIVV = 0.95E22
IF(XNDIV.LT.1E16) XNDIVV = 1.1E16
TNDIV = TDIVV
IF(TNDIV.GE.1000) TNDIV = 995.
CALL INTERP(TDIVV,TDIVV,TNDIV,XNDIVV)
SVELDIV = SEL(1)
SVELDIVN= SELN(1)
SVCXDIV = SCX(1)
SVATDIV = SEL(1) + SCX(1)
SVIONDIV = SION(1)
SVTOTDIV = SION(1) + SEL(1) + SCX(1) + CNEUT*SELN(1)*XNODIV/
2 (XNDIV)
Y = SVCXDIV + SVELDIV + CNEUT*SELN(1)*XNODIV/(XNDIV)
RSCATDIV = Y/(Y + SVIONDIV)
ASCATDIV = Y/(Y + SVIONDIV)
SVRECDIV = RECOM(1)
VODIV = SQRT(XK*TDIV/XMASS)
EIONDIV = 17.5
IF(XNDIV.LE.1.E21)
2 EIONDIV = 17.5 + (5.+37.5/TDIV)*LOG10(1.E21/(XNDIV))
IF(XNDIV.GT.1.E21)
2 EIONDIV = (30.6 - 16.4*EXP(-5.E19/XNDIV))*
3 EXP(5.45/(TDIV*EXP((XNDIV/1.37E20)**0.26)))
ZDIV = EPDIV*DELN*XNDIV*SVTOTDIV/VODIV
CALL EXPINT(ZDIV,E1,E2,E3,E4)
PO = (0.5 - E3)/ZDIV
C TOTAL ESCAPE PROBABILITY
PESCAPE = PO/(1. - ASCATDIV*(1. - PO))
C FIRST-FLIGHT TRANSMISSION PROBABILITY
TODIV = E2
C REFLECTION PROBABILITY
RDIV = 0.5*ASCATDIV*PESCAPE*(1. - TODIV)
C TOTAL TRANSMISSION PROBABILITY
TRDIV = TODIV + RDIV
C X-POINT PLASMA
C FIRST-FLIGHT ESCAPE PROBABILITY
TXPT = TSOL
XNXPT = XNSOL
TDXPT = TXPT
IF(TXPT.LT.1.E-1) TDXPT = 1.05E-1
IF(TXPT.GT.1E3) TDXPT = .95E3
XNDXPT = XNXPT
IF(XNDXPT.GT.1E22) XNDXPT = 0.95E22
IF(XNDXPT.LT.1E16) XNDXPT = 1.1E16
TNDXPT = TDXPT
IF(TNDXPT.GE.1000) TNDXPT = 995.
CALL INTERP(TDXPT,TDXPT,TNDXPT,XNDXPT)
SVELXPT = SEL(1)
SVELXPTN= SELN(1)
SVCXXPT = SCX(1)
SVATXPT = SEL(1) + SCX(1)
SVIONXPT = SION(1)
SVTOTXPT = SION(1) + SEL(1) + SCX(1) + CNEUT*SELN(1)*XNOXPT(58)/
2 (XNXPT)
Y = SVCXXPT + SVELXPT + CNEUT*SELN(1)*XNOXPT(58)/(XNXPT)
RSCATXPT = Y/(Y + SVIONXPT)
ASCATXPT = Y/(Y + SVIONXPT)
SVRECXPT = RECOM(1)
VOXPT = VOSOL
EIONXPT = 17.5
IF(XNXPT.LE.1.E21)
2 EIONXPT = 17.5 + (5.+37.5/TXPT)*LOG10(1.E21/(XNXPT))
IF(XNXPT.GT.1.E21)
2 EIONXPT = (30.6 - 16.4*EXP(-5.E19/XNXPT))*
3 EXP(5.45/(TXPT*EXP((XNXPT/1.37E20)**0.26)))
YXPT = DELN
ZXPT = YXPT*XNXPT*SVTOTXPT/VOXPT
CALL EXPINT(ZXPT,E1,E2,E3,E4)
PO = (0.5 - E3)/ZXPT
C TOTAL ESCAPE PROBABILITY
PESCAPE = PO/(1. - ASCATXPT*(1. - PO))
C FIRST-FLIGHT TRANSMISSION PROBABILITY
TOXPT = E2
C REFLECTION PROBABILITY
RXPT = 0.5*ASCATXPT*PESCAPE*(1. - TOXPT)
C TOTAL TRANSMISSION PROBABILITY
TRXPT = TOXPT + RXPT
C ALBEDO
ALPHAXPT = RXPT + ALPHAX*TRXPT/(1.-ALPHAX*RXPT)
C DIFFUSION THEORY ALBEDO FOR NEUTRALS ENTERING FROM
C RECYCLING REGION
C ***********OUTBOARD*******************
X = SQRT((1./ASCATDIV)-1.)
ALPHAR = (1.- 1.155*X)/(1.+1.155*X)
PHI = 3.1416 - THETA1
X = (XLPAR1 - XLPERP1 - DELLT1)*BETAG
ADIV = 0.5*ASCATDIV*PESCAPE*(1.- ALPHAR)*EPDIV1*DELN1/X
PDIV = PESCAPE
ZDIV = (1.-ALPHAR)*EPDIV1*DELN1*(1.-PDIV)
C FIRST-FLIGHT TRANSMISSION OF NEUTRALS ENTERING FROM RECYCLING REGIONS
VODIV2 = SQRT(XK*TD/XMASS)
YMFPDIV = 1./(XND*SVTOTDIV/VODIV2)
Z = BETAG*(XLPAR1-XLPERP1-DELLT1)
CALL TRANSM(1,0,YMFPDIV,DELNS1,Z,THETA1,0.0,0.0,0.0,TOR2PF1)
CALL TRANSM(1,0,YMFPDIV,DELNS1,Z,PHI,0.0,0.0,0.0,TOR2PL1)
c REDUCE ALBEDO TO ACCOUNT FOR NON-INFINITE MEDIUM
TOR = TOR2PF1 + TOR2PL1
TRANS = TOR + (1.-TOR)*ASCATDIV*PDIV
ALPHAREFF = ALPHAR*(1. - TRANS)/(1. - ALPHAR*TRANS)
C ESCAPE QUANTITITES
TDIVPF = (TOR2PF1 + 0.5*(1.-TOR)*ASCATDIV*PDIV)
TDIVPL = (TOR2PL1 + 0.5*(1.-TOR)*ASCATDIV*PDIV)
ESCPF = (1.-ALPHAREFF)*(TOR2PF1 + 0.5*(1.-TOR)*ASCATDIV*PDIV)
ESCPL = (1.-ALPHAREFF)*(TOR2PL1 + 0.5*(1.-TOR)*ASCATDIV*PDIV)
ESCPL1 = ESCPL
ESCPL2 = ESCPL
ESCPF1 = ESCPF
ESCPF2 = ESCPF
C PUMP GEOMETRY
N1 = 1
CALL GEOMETRY(N1)
C SET UP SOURCE VECTOR & TEP MATRIX
C SOURCE
C OUTBOARD
100 AD = BETAG*DELLT1*DELN1*EPDIV1
SAT = (XND)*((XND)*SVRECD)/SQRT(RATNENI2D)
CSD = SQRT(2.*TD*XK/XMASS)
SOURCE(1) = SAT*PESCAPER1*XLAMPF1*AD
SOURCE(2) = SAT*PESCAPER1*XLAMPL1*AD
SOURCE(3) = SAT*PESCAPER1*XLAMP1*AD
SOURCE(4) = SAT*PESCAPER1*XLAMD1*AD +
2 XND*CSD*DELN1*EPDIV1*FAC6*BETAG*FSHEATH
Y =(BETAG*XLPERP1-0.5*DELXPT)*(1.-ALPHASOL) +
2 XLWSPL*(1.-RWSOL)+DELX
SOURCE(5) = 0.5*FUELPF/(6.2832*RMAJOR)
SOURCE(6)=(FUELPLOUT+DELX*FUELMP/Y)/(6.2832*RMAJOR) +
2 betag*DNRAD*RWPL
C INBOARD
AD = BETAG*DELLT2*DELN2*EPDIV2
SAT = (XND)*((XND)*SVRECD)/SQRT(RATNENI2D)
CSD = SQRT(2.*TD*XK/XMASS)
SOURCE(7) = SAT*PESCAPER2*XLAMPF2*AD
SOURCE(8) = SAT*PESCAPER2*XLAMPL2*AD
SOURCE(9) = SAT*PESCAPER2*XLAMP2*AD
SOURCE(10) = SAT*PESCAPER2*XLAMD2*AD +
2 XND*CSD*DELN2*EPDIV2*FAC6*BETAG*FSHEATH
SOURCE(11) = 0.5*FUELPF/(6.2832*RMAJOR)
SOURCE(12)=FUELPLIN/(6.2832*RMAJOR) +
2 betag*DNRAD*RWPL
C THE ADD-ON TREATS LEAKAGE INTO PLASMA AT XPT
C TRANSPORT MATRIX ELEMENTS
C *******OUTBOARD*********
Z = BETAG*DELLT1
ZQ = BETAG*(XLPAR1 - XLPERP1 - DELLT1)
ZW = BETAG*(XLPAR1-XLPERP1)
ZZ = BETAG*XLUP1 - 0.5*DELXPT
DPL1 = ZQ + XLWPL1*(1.-RWPL) + DELX*(1.-DELX/Y) + DELPUMP1 +
2 ZZ*(1. - ALPHASOL) - (1.-GEOM21)*ZQ*RDIV*EPOUT +
3 YPUMP1*(1.-RPUMP1)*APL1
DPF1 = ZW - (RDIV+DOUT22*TRDIV*TRDIV)*ZQ + XLWPF1*(1.-RWPF) +
2 0.5*DELXPT*(1. - ALPHASOLXPT) + DELPF +
3 YPUMP1*(1.-RPUMP1)*APF1
C
AMAT(1,1) = Z
AMAT(1,2) = -1.*Z*DOUT11*TPL2PF1
AMAT(1,3) = -1.*DELNS1*(ALPHAREFF*TP2PF1 + DOUT21*ESCPL1*TPL2PF1)
AMAT(1,4) = -1.*(TD2PF1 + HDIS0*GPF1)*DELNS1
AMAT(1,5) = -1.*(TPF2PF1*Z + ZQ*TRDIV*TPL2PF1*DOUT21)
AMAT(1,6) = -1.*(DELPUMP1 + ZQ*DOUT21*RDIV)*TPL2PF1
AMAT(2,1) = 0.
AMAT(2,2) = Z*(1.-DOUT11*TPL2PL1)
AMAT(2,3) = -1.*DELNS*(ALPHAREFF*TP2PL1 + DOUT21*TPL2PL1*ESCPL1)
AMAT(2,4) = -1.*(TD2PL1+HDIS0*GPL1)*DELNS1
AMAT(2,5) = -1.*(Z*TPF2PL1 + DOUT21*ZQ*TRDIV*TPL2PL1)
AMAT(2,6) = -1.*(DELPUMP1 + DOUT21*ZQ*RDIV)*TPL2PL1
AMAT(3,1) = 0.
AMAT(3,2) = -1.*Z*DOUT11*TPL2P1
AMAT(3,3) = DELNS1*(1.-ALPHAREFF*TP2P1 - DOUT21*TPL2P1*ESCPL1)
AMAT(3,4) = -1.*(TD2P1+HDIS0*GP1)*DELNS1
AMAT(3,5) = -1.*(Z*TPF2P1 + DOUT21*ZQ*TRDIV*TPL2P1)
AMAT(3,6) = -1.*(DELPUMP1 + DOUT21*ZQ*RDIV*TPL2P1)
AMAT(4,1) = 0.
AMAT(4,2) = -1.*Z*DOUT11*TPL2D1
AMAT(4,3) = -1.*DELNS1*(ALPHAREFF*TP2D1 + DOUT21*TPL2D1*ESCPL1)
AMAT(4,4) = ((1./(RNN*RWDP)) - TD2D1 - HDIS0*GD1)*DELNS1
AMAT(4,5) = -1.*(Z*TPF2D1 + DOUT21*ZQ*TRDIV*TPL2D1)
AMAT(4,6) = -1.*(DELPUMP1 + DOUT21*ZQ*RDIV)*TPL2D1
AMAT(5,1) = -1.*Z
AMAT(5,2) = -1.*Z*CPUMP21*GEOM11*(TRDIV+DOUT22*RDIV*TRDIV)
AMAT(5,3) = -1.*DELNS1*(ESCPF1 + DOUT22*TRDIV*ESCPL1)
AMAT(5,4) = 0.
AMAT(5,5) = DPF1
AMAT(5,6) = -1.*ZQ*(TRDIV + DOUT22*TRDIV*RDIV)
AMAT(5,11) = -1.*DELPF
AMAT(6,1) = 0.
AMAT(6,2) = -1.*Z*((1.-GEOM11) + (1.-GEOM21)*DOUT12*RDIV)
AMAT(6,3) = -1.*(1.-GEOM21)*ESCPL1*DELNS1*EPOUT
AMAT(6,4) = 0.
AMAT(6,5) = -1.*(1.-GEOM21)*ZQ*TRDIV*EPOUT
AMAT(6,6) = DPL1
C ******INBOARD*******
Z = BETAG*DELLT2
ZQ = BETAG*(XLPAR2 - XLPERP2 - DELLT2)
ZW = BETAG*(XLPAR2-XLPERP2)
ZZ = BETAG*XLPERP2
DPL2 = ZW + XLWPL2*(1.-RWPL) + ZZ*(1. - ALPHASOL)
2 - ZQ*(RDIV + DIN22*TRDIV*TRDIV) +
3 YPUMP2*(1.-RPUMP2)*APL2
DPF2 = ZW - (1. - GEOM22)*RDIV*ZQ*EPIN + DELPUMP2 +
2 0.5*DELXPT*(1. - ALPHASOLXPT) + DELPF + XLWPF2*(1.-RWPF)
AMAT(7,7) = Z*(1.- DIN11*TPF2PF2)
AMAT(7,8) = 0.
AMAT(7,9) = -1.*DELNS2*(ALPHAREFF*TP2PF2 + DIN21*ESCPF2*TPF2PF2)
AMAT(7,10) = -1.*(TD2PF2 + HDIS0*GPF2)*DELNS2
AMAT(7,11) = -1.*(DELPUMP2 + ZQ*RDIV*DIN21)*TPF2PF2
AMAT(7,12) = -1.*(Z*TPL2PF2 + ZQ*DIN21*TRDIV*TPL2PF2)
AMAT(8,7) = -1.*Z*DIN11*TPF2PL2
AMAT(8,8) = Z
AMAT(8,9) = -1.*DELNS2*(ALPHAREFF*TP2PL2 + DIN21*TPF2PL2*ESCPF2)
AMAT(8,10) = -1.*(TD2PL2+HDIS0*GPL2)*DELNS2
AMAT(8,11) = -1.*(DELPUMP2 + DIN21*ZQ*RDIV)*TPF2PL2
AMAT(8,12) = -1.*(Z*TPL2PL2 + DIN21*ZQ*TRDIV*TPF2PL2)
AMAT(9,7) = -1.*DIN11*Z*TPF2P2
AMAT(9,8) = 0.
AMAT(9,9) = DELNS2*((1.-ALPHAREFF*TP2P2) - DIN21*TPF2P2*ESCPF2)
AMAT(9,10) = -1.*(TD2P2+HDIS0*GP2)*DELNS2
AMAT(9,11) = -1.*(DELPUMP2 + DIN21*ZQ*RDIV)*TPF2P2
AMAT(9,12) = -1.*(Z*TPF2P1 + DIN21*ZQ*TRDIV*TPF2P2)
AMAT(10,7) = -1.*DIN11*Z*TPF2D2
AMAT(10,8) = 0.
AMAT(10,9) = -1.*DELNS2*(ALPHAREFF*TP2D2 + DIN21*TPF2D2*ESCPF2)
AMAT(10,10) = ((1./(RNN*RWDP)) - TD2D2 - HDIS0*GD2)*DELNS2
AMAT(10,11) = -1.*(DELPUMP2 + DIN21*ZQ*RDIV)*TPF2D2
AMAT(10,12) = -1.*(Z*TPL2D2 + DIN21*ZQ*TRDIV*TPF2D2)
AMAT(11,5) = -1.*DELPF
AMAT(11,7) = -1.*Z*((1.-GEOM12)+(1.-GEOM22)*DIN12*RDIV)
AMAT(11,8) = 0.
AMAT(11,9) = -1.*(1.-GEOM22)*DELNS2*ESCPF2*EPOUT
AMAT(11,10) = 0.
AMAT(11,11) = DPF2
AMAT(11,12) = -1.*ZQ*TRDIV*EPOUT
AMAT(12,7) = -1.*Z*DIN22*TRDIV
AMAT(12,8) = -1.*Z
AMAT(12,9) = -1.*DELNS2*(ESCPL2 + DIN22*TRDIV*ESCPF2)
AMAT(12,10) = 0.
AMAT(12,11) = -1.*ZQ*EPIN*TRDIV
AMAT(12,12) = DPL2
C INVERT TRANSPORT MATRIX & CONSTRUCT FLUX SOLUTION
CALL sgesv_NR (12,1,AMAT,12,IPIV,SOURCE,12,INFO)
C THIS RETURNS THE FLUXES IN THE SOURCE VECTOR
GAMR2PF1 = SOURCE(1)
GAMR2PL1 = SOURCE(2)
GAMR2P1 = SOURCE(3)
GAMZERO1 = SOURCE(4)
GAMOUTPF1= SOURCE(5)
GAMOUTPL1= SOURCE(6)
GAMR2PF2 = SOURCE(7)
GAMR2PL2 = SOURCE(8)
GAMR2P2 = SOURCE(9)
GAMZERO2 = SOURCE(10)
GAMOUTPF2= SOURCE(11)
GAMOUTPL2= SOURCE(12)
C CONSTRUCT IONIZATION TERMS & TOTAL NEUTRAL DENSITIES
C THE 'IONIZATION RATES' ARE THE IONIZATION RATES/2*PI*R IN THE POLOIDAL
C PLANE. TRUE IONIZATION RATES IN INCLINED PLASMA STRIP OBTAIN BY
C DIVIDING BY BETAG--FOR USE IN PLASMA BALANCE.
c *************OUTSIDE******************
C GEOMETRIC
GAMP2R1 = ALPHAREFF*GAMR2P1
DXLDIV = XLPAR1-XLPERP1-DELLT1
DLDIV = XLPAR1-XLPERP1
Z = BETAG*DELLT1
ZQ = BETAG*(XLPAR1 - XLPERP1 - DELLT1)
ZZ = BETAG*XLUP1 - 0.5*DELXPT
C2P = DOUT22*TDIV
CPUMP1 = CPUMP11*YPUMP1
CPUMP2 = CPUMP21*YPUMP1
C NEUTRAL FLUX TO SOL IN MIDPLANE
GAMOUTSPL =(FUELMP+DELX*GAMOUTPL1)/
1 (BETAG*XLPERP1*(1.-ALPHASOL)+XLWSPL*(1.-RWSOL)+DELX)
C NEUTRAL DENSITY IN PRIVATE FLUX REGIONS
XNOPF1 = GAMOUTPF1/VODIV
XNOPF2 = GAMOUTPF2/VODIV
C NEUTRAL DENSITY IN PLENUM REGIONS
XNOPL1 = GAMOUTPL1/VODIV
XNOPL2 = GAMOUTPL2/VODIV
XNOSPL = GAMOUTSPL/VOSOL
C DIVERTOR PLASMA CHANNELS
XNODIV1 = (GAMDIV2PL1+GAMDIV2PF1)/VODIV
XNODIV2 = (GAMDIV2PL2+GAMDIV2PF2)/VODIV
C PUMP FLUXES
Z = BETAG*DELLT1
ZQ = BETAG*(XLPAR1 - XLPERP1 - DELLT1)
GAMPUMP1 = (GEOM11*Z*GAMR2PL1 + GEOM21*ZQ*GAMDIV2PL1 +
2 (APL1*GAMOUTPL1 + APF1*GAMOUTPF1)*YPUMP1)/YPUMP1
Z = BETAG*DELLT2
ZQ = BETAG*(XLPAR2 - XLPERP2 - DELLT2)
GAMPUMP2 = (GEOM12*Z*GAMR2PF2 + GEOM22*ZQ*GAMDIV2PF2 +
2 (APL2*GAMOUTPL2)*YPUMP2)/YPUMP2
C FLUXES OUT OF DIVERTOR PLASMA CHANNELS
GAMDIV2PL1 = (BETAG*DELX1*(GAMOUTPL1*RDIV + GAMOUTPF1*TRDIV) +
1 CPUMP21*GAMPUMP1*YPUMP1*RDIV + DELNS1*ESCPL*GAMR2P1)/
2 (BETAG*DELX1)
GAMDIV2PF1 = (BETAG*DELX1*(GAMOUTPL1*TRDIV + GAMOUTPF1*RDIV) +
1 CPUMP21*GAMPUMP1*YPUMP1*TRDIV +DELNS1*ESCPF*GAMR2P1)/
2 (BETAG*DELX1)
GAMDIV2PL2 = (BETAG*DELX2*(GAMOUTPL2*RDIV + GAMOUTPF2*TRDIV) +
1 CPUMP22*GAMPUMP2*YPUMP2*TRDIV + DELNS2*ESCPL*GAMR2P2)/
2 (BETAG*DELX2)
GAMDIV2PF2 = (BETAG*DELX2*(GAMOUTPL2*TRDIV + GAMOUTPF2*RDIV) +
1 CPUMP22*GAMPUMP2*YPUMP2*RDIV +DELNS2*ESCPF*GAMR2P2)/
2 (BETAG*DELX2)
C DIVERTOR PLASMA CHANNELS BY IONIZATION FLUX BALANCE
YIONDIV1 = ((GAMOUTPL1+GAMOUTPF1)*BETAG*DELX1 + DELNS1*GAMR2P1 +
2 CPUMP21*YPUMP1*GAMPUMP1) -
3 ((GAMDIV2PL1+GAMDIV2PF1)*BETAG*DELX1 +
4 ALPHAREFF*DELNS1*GAMR2P1)
xiondiv = (betag*delx1*(gamoutpl1+gamoutpf1)+
2 cpump21*ypump1*gampump1)*(1.-rdiv-trdiv)+
3 delns1*gamr2p1*(1.-alphareff-escpf-escpl)
XNODIV1 = SQRT(RATNENI2D)*(YIONDIV1)/
2 (XNDIV*SVIONDIV*EPDIV1*DELN1*BETAG*DELX1)
YIONDIV2 = ((GAMOUTPL2+GAMOUTPF2)*BETAG*DELX2 + DELNS2*GAMR2P2 +
2 CPUMP22*YPUMP2*GAMPUMP2) -
3 ((GAMDIV2PL2+GAMDIV2PF2)*BETAG*DELX2 +
4 ALPHAREFF*DELNS2*GAMR2P2)
XNODIV2 = SQRT(RATNENI2D)*(YIONDIV2)/
2 (XNDIV*SVIONDIV*EPDIV2*DELN2*BETAG*DELX2)
C RECYCLING REGIONS BY IONIZATION FLUX BALANCES
Z = BETAG*DELLT1
GAMINDD1 = GAMOUTPF1*Z + CPUMP11*GAMPUMP1*YPUMP1 +
2 DELPUMP1*GAMOUTPL1 +
2 (GAMZERO1/(1.-FMOL)+ALPHAREFF*GAMR2P1)*DELNS1
GAMZEROM1 = ((GAMZERO1/((1.-FMOL)*RWDP)) -
2 XND*CSD*SIN(THETA1)*betag*FAC6)
GAMOUTD1 = (GAMR2PF1+GAMR2PL1)*Z + (GAMZEROM1+GAMR2P1)*DELNS1
YIONREC1 = (GAMINDD1 - GAMOUTD1) + XND*SVRECD*XND*AD -
2 GAMZERO1*DELNS1*HDISP
XNOD1 = SQRT(RATNENI2D)*(YIONREC1)*FACTORION/(XND*SVIOND*AD)
YIONREC1 = (GAMINDD1 - GAMOUTD1)
YALPHD1 = YIONREC1/(XNOD1*SVTOTD*XND*AD)
C (IONIZATION - RECOMBINATION - DISSOCIATION) NEUTRAL RATE
Z = BETAG*DELLT2
GAMINDD2 = GAMOUTPF2*Z + CPUMP12*GAMPUMP2*YPUMP2 +
2 DELPUMP2*GAMOUTPF2 +
2 (GAMZERO2/(1.-FMOL)+ALPHAREFF*GAMR2P2)*DELNS2
GAMZEROM2 = ((GAMZERO2/((1.-FMOL)*RWDP)) -
2 XND*CSD*SIN(THETA2)*betag*FAC6)
GAMOUTD2 = (GAMR2PF2+GAMR2PL2)*Z + (GAMZEROM2+GAMR2P2)*DELNS2
YIONREC2 = (GAMINDD2 - GAMOUTD2) + XND*SVRECD*XND*AD -
2 GAMZERO2*DELNS2*HDISP
XNOD2 = SQRT(RATNENI2D)*(YIONREC2)*FACTORION/(XND*SVIOND*AD)
YIONREC2 = (GAMINDD2 - GAMOUTD2)
YALPHD2 = YIONREC2/(XNOD2*SVTOTD*XND*AD)
C IONIZATION & NEUTRAL DENSITIES IN SOL, TB & PEDESTAL
CALL NEUTFRAC
CALL NEUTXPT
C CONSTRUCT COLD NEUTRAL DENSITIES
XNCOLDD =((1.-TOD1)*GAMZERO1*DELNS1 +
2 ((1.-TOPF1)*GAMOUTPF+(1.-TOPL1)*GAMOUTPL1)*BETAG*DELLT1)/
3 (XND*SVTOTD*EPDIV1*DELN1*BETAG*DELLT1)
TODIV = RDIV + TRDIV
XNCOLDDIV = (GAMOUTPF1 + GAMOUTPL1)*(1.-TODIV)/
2 (XNDIV*SVTOTDIV*EPDIV1*DELN1)
C FIRST-COLLISION SCATTERING RATES
FCRATDOLD = FCRATD
FCRATDIVOLD = FCRATDIV
FCRATSOLOLD = FCRATSOL
FCRATD = RSCATD*((1.-TOD1)*GAMZERO1*EPDIV1*DELN1)
C 2 + (1.-TOPF)*GAMOUTPF*BETAG*DELLT +
C 3 (1.-TOPL)*GAMOUTPL*BETAG*DELLT)
IF(ITERN.GT.2) FCRATD = FCRATDOLD + CHANGE*(FCRATD-FCRATDOLD)
X = (XLPAR1 - XLPERP1 - DELLT1)*BETAG
FCRATDIV = RSCATDIV*(GAMOUTPF1 + GAMOUTPL1)*(1.-TODIV)*X
IF(ITERN.GT.2) FCRATDIV=FCRATDIVOLD+CHANGE*(FCRATDIV-FCRATDIVOLD)
C TRUE FIRST COLLISION RATES FOR USE IN PLASMA BALANCE OBTAINED BY
C DIVIDING BY BETAG
C FRACTION SOURCE PARTICLES PUMPED--SHOULD BE UNITY
PUMPRATE=(GAMPUMP1*YPUMP1*(1.-RPUMP1) +
2 GAMPUMP2*YPUMP2*(1.-RPUMP2))*6.2832*RMAJOR
FUELPL = FUELPLIN + FUELPLOUT
SOURCEPART = (FUELMP+FUELPF+FUELPL+SPELLET)
PUMPFRAC = PUMPRATE/SOURCEPART
C RESET DELN = DELN
DELN = DELNREAL
IF(ITERN.LT.66) GOTO 500
C DO SOL/DIV PLASMA CALCULATION FOR OUTBOARD LEG.
C SET NEUTRAL PARAMETERS TO OUTBOARD PARAMETERS
500 XLPAR = XLPAR1
XLPERP = XLPERP1
DELLT = DELLT1
EPDIV = EPDIV1
YIONREC = YIONREC1
YIONDIV = YIONDIV1
GEOM1 = GEOM11
GEOM2 = GEOM21
GAMR2PF = GAMR2PF1
GAMR2PL = GAMR2PL1
GAMR2P = GAMR2P1
GAMZERO = GAMZERO1
GAMZEROM = GAMZEROM1
GAMOUTPL = GAMOUTPL1
GAMOUTPF = GAMOUTPF1
GAMDIV2PL = GAMDIV2PL1
GAMDIV2PF = GAMDIV2PF1
XNOD = XNOD1
XNODIV = XNODIV1
YALPHD = YALPHD1
THETA = THETA1
gpl = gpl1
tod2pf = tod2pf1
TOD = TOD1
TOD2PL = TOD2PL1
GPF = GPF1
XLPF = XLWPF1
RETURN
END