V3FIT
add_fluxes.f90
1  SUBROUTINE add_fluxes_par(overg, bsupu, bsupv, lcurrent)
2  USE vmec_main
3  USE realspace, ONLY: pwint, pguu, pguv, pchip, pphip
4  USE vmec_input, ONLY: nzeta
5  USE vmec_dim, ONLY: ntheta3
6  USE parallel_include_module
7  IMPLICIT NONE
8 !-----------------------------------------------
9 ! D u m m y A r g u m e n t s
10 !-----------------------------------------------
11  REAL(rprec), DIMENSION(nznt,ns), INTENT(in) :: overg
12  REAL(rprec), DIMENSION(nznt,ns), INTENT(inout) :: bsupu, bsupv
13  LOGICAL, INTENT(in) :: lcurrent
14 !-----------------------------------------------
15 ! L o c a l P a r a m e t e r s
16 !-----------------------------------------------
17  REAL(rprec), PARAMETER :: p5=0.5_dp, c1p5=1.5_dp
18  REAL(rprec), PARAMETER :: iotaped = 0.10
19 !-----------------------------------------------
20 ! L o c a l V a r i a b l e s
21 !-----------------------------------------------
22  INTEGER :: js, l
23  REAL(rprec) :: top, bot
24 
25  INTEGER :: i, j, k, nsmin, nsmax, lnsnum, istat
26 !-----------------------------------------------
27 !
28 ! ADD MAGNETIC FLUX (CHIP, PHIP) TERMS TO BSUPU=-OVERG*LAM_V, BSUPV=OVERG*LAM_U
29 ! COMPUTE FLUX FROM ITOR = <B_u>, ITOR(s) = integrated toroidal current (icurv)
30 ! IF ncurr == 1, AND ictrl_prec2d != 0, COMPUTE FORCE IN TOMNSP TO UPDATE chips
31 !
32 
33  IF (.NOT.lcurrent .OR. ncurr.EQ.0) GOTO 100
34 
35  nsmin=max(2,t1lglob); nsmax=t1rglob
36  DO js = nsmin, nsmax
37  top = icurv(js)
38  bot = 0
39  DO j=1,nznt
40  top = top - pwint(j,js)*(pguu(j,js)*bsupu(j,js) &
41  + pguv(j,js)*bsupv(j,js))
42  bot = bot + pwint(j,js)*overg(j,js)*pguu(j,js)
43  END DO
44  IF (bot.ne.zero) chips(js) = top/bot
45  IF (phips(js).ne.zero) iotas(js) = chips(js)/phips(js)
46  END DO
47 
48  100 CONTINUE
49 
50  nsmin=max(2,t1lglob); nsmax=t1rglob
51 ! CHANGE THIS FOR lRFP = T (solve for phips?)
52  IF (ncurr .EQ. 0) THEN
53  chips(nsmin:nsmax) = iotas(nsmin:nsmax)*phips(nsmin:nsmax)
54  ELSE IF (.NOT.lcurrent) THEN
55  WHERE (phips(nsmin:nsmax) .NE. zero) &
56  iotas(nsmin:nsmax) = chips(nsmin:nsmax)/phips(nsmin:nsmax)
57  END IF
58 
59  DO js = nsmin, nsmax
60  pchip(:,js) = chips(js)
61  END DO
62 
63  nsmin=max(2,t1lglob); nsmax=min(ns-1,trglob)
64  IF (t1lglob .eq. 1 .and. trglob .gt. 2) THEN
65  chipf(1) = c1p5*chips(2) - p5*chips(3)
66  ELSE IF (t1lglob .eq. 1) THEN
67  chipf(1) = chips(2)
68  END IF
69  chipf(nsmin:nsmax) = (chips(nsmin:nsmax) + chips(nsmin+1:nsmax+1))/2
70  IF (nsmax.EQ.ns) chipf(ns) = c1p5*chips(ns)- p5*chips(ns-1)
71 
72 ! Do not compute iota too near origin
73  IF(trglob_arr(1).LE.2) THEN
74 #if defined(MPI_OPT)
75  CALL mpi_bcast(iotas(3),1,mpi_real8,1,ns_comm,mpi_err)
76 #endif
77  END IF
78  IF (lrfp) THEN
79  IF (nsmin.EQ.1) iotaf(1) = one/(c1p5/iotas(2) - p5/iotas(3))
80  IF (nsmax.EQ.ns) iotaf(ns)=one/(c1p5/iotas(ns)-p5/iotas(ns-1))
81  DO js = max(2,t1lglob), min(ns-1,t1rglob)
82  iotaf(js) = 2.0_dp/(one/iotas(js) + one/iotas(js+1))
83  END DO
84  ELSE
85  IF (nsmin.EQ.1) iotaf(1) = c1p5*iotas(2) - p5*iotas(3)
86  IF (nsmax.EQ.ns) iotaf(ns)=c1p5*iotas(ns) - p5*iotas(ns-1)
87  DO js = max(2,t1lglob), min(ns-1,trglob)
88  iotaf(js) = p5*(iotas(js) + iotas(js+1))
89  END DO
90  END IF
91 
92  nsmin=max(1,t1lglob); nsmax=min(ns,t1rglob)
93  bsupu(:,nsmin:nsmax) = bsupu(:,nsmin:nsmax)+pchip(:,nsmin:nsmax)*overg(:,nsmin:nsmax)
94 
95  END SUBROUTINE add_fluxes_par
96 
97  SUBROUTINE add_fluxes(overg, bsupu, bsupv, lcurrent)
98  USE vmec_main
99  USE realspace, ONLY: wint, guu, guv, chip, phip
100 
101  IMPLICIT NONE
102 !-----------------------------------------------
103 ! D u m m y A r g u m e n t s
104 !-----------------------------------------------
105  REAL(rprec), DIMENSION(nrzt), INTENT(in) :: overg
106  REAL(rprec), DIMENSION(nrzt), INTENT(inout) :: bsupu, bsupv
107  LOGICAL, INTENT(in) :: lcurrent
108 !-----------------------------------------------
109 ! L o c a l P a r a m e t e r s
110 !-----------------------------------------------
111  REAL(rprec), PARAMETER :: p5=0.5_dp, c1p5=1.5_dp
112  REAL(rprec), PARAMETER :: iotaped = 0.10_dp
113 !-----------------------------------------------
114 ! L o c a l V a r i a b l e s
115 !-----------------------------------------------
116  INTEGER :: js, l
117  REAL(rprec) :: top, bot
118 
119 !-----------------------------------------------
120 !
121 ! ADD MAGNETIC FLUX (CHIP, PHIP) TERMS TO BSUPU=-OVERG*LAM_V, BSUPV=OVERG*LAM_U
122 ! COMPUTE FLUX FROM ITOR = <B_u>, ITOR(s) = integrated toroidal current (icurv)
123 ! IF ncurr == 1
124 !
125  IF (.not.lcurrent .or. ncurr.eq.0) GOTO 100
126 
127 ! nsmin=MAX(2,tlglob); nsmax=trglob
128 
129  DO js = 2, ns
130  top = icurv(js)
131  bot = 0
132  DO l = js, nrzt, ns
133  top = top - wint(l)*(guu(l)*bsupu(l) + guv(l)*bsupv(l))
134  bot = bot + wint(l)*overg(l)*guu(l)
135  END DO
136  IF (bot .ne. zero) chips(js) = top/bot
137  IF (phips(js) .ne. zero) iotas(js) = chips(js)/phips(js)
138  END DO
139 
140  100 CONTINUE
141 
142 ! CHANGE THIS FOR lRFP = T (solve for phips?)
143  IF (ncurr .eq. 0) THEN
144  chips = iotas*phips
145  ELSE IF (.not.lcurrent) THEN
146  WHERE (phips .ne. zero) iotas = chips/phips
147  END IF
148 
149  DO js = 2, ns
150  chip(js:nrzt:ns) = chips(js)
151  END DO
152 
153  chipf(1) = c1p5*chips(2) - p5*chips(3) !SPH ADDED THIS 4-8-16
154  chipf(2:ns1) = (chips(2:ns1) + chips(3:ns1+1))/2
155  chipf(ns) = c1p5*chips(ns)- p5*chips(ns1) !SPH FIXED THIS 4-8-16
156 
157 ! Do not compute iota too near origin
158  IF (lrfp) THEN
159  iotaf(1) = one/(c1p5/iotas(2) - p5/iotas(3))
160  iotaf(ns) = one/(c1p5/iotas(ns) - p5/iotas(ns1))
161  DO js = 2, ns-1
162  iotaf(js) = 2.0_dp/(one/iotas(js) + one/iotas(js+1))
163  END DO
164 
165  ELSE
166  iotaf(1) = c1p5*iotas(2) - p5*iotas(3) !zero gradient near axis
167  iotaf(ns) = c1p5*iotas(ns) - p5*iotas(ns-1)
168  DO js = 2, ns-1
169  iotaf(js) = p5*(iotas(js) + iotas(js+1))
170  END DO
171  END IF
172 
173  bsupu(:nrzt) = bsupu(:nrzt)+chip(:nrzt)*overg(:nrzt)
174 
175  END SUBROUTINE add_fluxes
Generated on Thu Mar 5 2020 15:49:26 for V3FIT by   doxygen 1.8.17