ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
•All Classes Files Functions Variables Pages
emeq_equil.f
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
7 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
8  module emeq_equil
9 
10  implicit none
11 
12  private
13  public :: alloc_emeq_equil, dealloc_emeq_equil, e3astr
14 
15  integer, save :: np, iter
16  double precision, save :: wsacc
17  double precision, save ::
18  1 wgbj, wgb, wgb0, wgbst, wgbd,
19  2 wgmj, wgmjex, wsqc, wbbs0, wbj0,
20  3 wbr00, wbr0, wsac1, wsac2, wsac3
21  double precision, allocatable, save ::
22  1 wsa(:), wsaa(:),
23  4 wsd1(:), wgl(:), wdgl(:),
24  5 wdsd1(:), wd2sd1(:),wbs1(:),
25  6 wd2gl(:), wbs2(:),
26  7 wsd3(:), wdsd3(:), wd2sd3(:),wbs3(:),
27  8 wbk02(:), wbk0(:), wdbk00(:),wbd02(:),
28  9 wbd0(:), wbg02(:), wbg0(:), wbd12(:),
29  & wbg332(:),wbg33(:), wbg222(:),wbg22(:)
30  double precision, allocatable, save ::
31  1 wsl0(:), wsv0(:), wsu0(:), wsw0(:),
32  2 wgmc(:), wdgmc(:), wba(:), wdba(:),
33  3 wbb(:), wdbb(:),
34  4 wscj1(:), wsci1(:), wscj3(:), wsci3(:),
35  5 wscj5(:), wsci5(:), wscj7(:), wsci7(:),
36  6 wbk10(:), wbk11(:), wbk20(:), wbk22(:),
37  7 wbk13(:), wbk30(:), wbk31(:), wbk33(:),
38  8 wsl1(:), wsv1(:), wsu1(:),
39  9 wsl22(:), wsl2(:), wsv2(:), wsu2(:),
40  & wsl3(:), wsv3(:), wsu3(:)
41  double precision, allocatable, save ::
42  1 wsp(:), wdsp(:), wsjp(:), wdsjp(:),
43  2 wsj(:), wdsj(:), wsq(:), wdsq(:),
44  3 wdsqrq(:),wbf(:), wbff(:), wsjsl(:),
45  4 wsjsr(:), wgp(:)
46  double precision, allocatable, save ::
47  1 skdr(:),skga(:),sqg22r(:)
48 
49 C Total length 76*NP (+21*NP)
50 
51  contains
52 
53  subroutine alloc_emeq_equil(NP_in)
54  implicit none
55  integer np_in
56  np=np_in
57  allocate(
58  1 wsa(np), wsaa(np),
59  4 wsd1(np), wgl(np), wdgl(np),
60  5 wdsd1(np), wd2sd1(np),wbs1(np),
61  6 wd2gl(np), wbs2(np),
62  7 wsd3(np), wdsd3(np), wd2sd3(np),wbs3(np),
63  8 wbk02(np), wbk0(np), wdbk00(np),wbd02(np),
64  9 wbd0(np), wbg02(np), wbg0(np), wbd12(np),
65  & wbg332(np),wbg33(np), wbg222(np),wbg22(np))
66  allocate(
67  1 wsl0(np), wsv0(np), wsu0(np), wsw0(np),
68  2 wgmc(np), wdgmc(np), wba(np), wdba(np),
69  3 wbb(np), wdbb(np),
70  4 wscj1(np), wsci1(np), wscj3(np), wsci3(np),
71  5 wscj5(np), wsci5(np), wscj7(np), wsci7(np),
72  6 wbk10(np), wbk11(np), wbk20(np), wbk22(np),
73  7 wbk13(np), wbk30(np), wbk31(np), wbk33(np),
74  8 wsl1(np), wsv1(np), wsu1(np),
75  9 wsl22(np), wsl2(np), wsv2(np), wsu2(np),
76  & wsl3(np), wsv3(np), wsu3(np))
77  allocate(
78  1 wsp(np), wdsp(np), wsjp(np), wdsjp(np),
79  2 wsj(np), wdsj(np), wsq(np), wdsq(np),
80  3 wdsqrq(np),wbf(np), wbff(np), wsjsl(np),
81  4 wsjsr(np), wgp(np))
82  allocate(
83  1 skdr(np),skga(np),sqg22r(np))
84  wsd1=0.d0
85  wdsd1=0.d0
86  wdgl=0.d0
87  wbs1=0.d0
88  wbs2=0.d0
89  wbs3=0.d0
90 
91  end subroutine alloc_emeq_equil
92 
93  subroutine dealloc_emeq_equil
94 
95  deallocate(
96  1 wsa, wsaa,
97  4 wsd1, wgl, wdgl,
98  5 wdsd1, wd2sd1,wbs1,
99  6 wd2gl, wbs2,
100  7 wsd3, wdsd3, wd2sd3,wbs3,
101  8 wbk02, wbk0, wdbk00,wbd02,
102  9 wbd0, wbg02, wbg0, wbd12,
103  & wbg332,wbg33, wbg222,wbg22)
104  deallocate(
105  1 wsl0, wsv0, wsu0, wsw0,
106  2 wgmc, wdgmc, wba, wdba,
107  3 wbb, wdbb,
108  4 wscj1, wsci1, wscj3, wsci3,
109  5 wscj5, wsci5, wscj7, wsci7,
110  6 wbk10, wbk11, wbk20, wbk22,
111  7 wbk13, wbk30, wbk31, wbk33,
112  8 wsl1, wsv1, wsu1,
113  9 wsl22, wsl2, wsv2, wsu2,
114  & wsl3, wsv3, wsu3)
115  deallocate(
116  1 wsp, wdsp, wsjp, wdsjp,
117  2 wsj, wdsj, wsq, wdsq,
118  3 wdsqrq,wbf, wbff, wsjsl,
119  4 wsjsr, wgp)
120  deallocate(
121  1 skdr,skga,sqg22r)
122  end subroutine dealloc_emeq_equil
123 
124 C======================================================================|
125 C Input:
126 C BA(1:NEQL)[MA/m^2] - flux function in front of R/r on the rhs of G-Sh. eq.
127 C BB(1:NEQL)[MA/m^2] - flux function in front of (r/R-R/r) - " - " - " - " -
128 C BR00 [m] - R the geom. center of the edge magnetic surface
129 C SA0 [m] - a_{edge} - "a" for the edge magnetic surface
130 C GL0 [d/l] - elongation of the edge magnetic surface
131 C GD30 [m] - triangularity of the edge magnetic surface (\delta(a_{edge}))
132 C NEQL [d/l]- # of grid points
133 C B0T [T] - vacuum toroidal magnetic field at the point r=R
134 C PLCUR[MA] - plasma current
135 C
136 C Output:
137 C GR(1:NEQL) [m] - \rho(a)
138 C GBD(1:NEQL) [m] - \Delta(a)
139 C GL(1:NEQL) [d/l] - \lambda(a)
140 C GSD(1:NEQL) [d/l] - \delta(a)/[a_{edge}(a/a_{edge})^2]
141 C GRA(1:NEQL) [d/l] - <(\nabla a)^2>
142 C GRAR(1:NEQL) [1/m^2] - <[(\nabla a)/r]^2>
143 C AVR2(1:NEQL) [1/m^2] - <[1/r]^2>
144 C AI0(1:NEQL) [mT] - I = r B_{tor}
145 C Everywhere above
146 C (i) the same equidistant grid with respect to the "radial"
147 C variable "a" is assumed to be used both for input and output.
148 C (ii) average <f> is an [0.5/\pi] integral over the poloidal angle \tau
149 C of the quantity [f\sqrt{g}], where g_{ij} is the metric tensor
150 C of {a,\tau,\zeta}.
151 C (iii) r is the current polar radius (the distance to the major torus axis)
152 C
153 
154 C - BR00,SA0 - MAJOR & MINOR RADII /METER/
155 C - GLO,GD3O - ELONGATION AND TRIANGULARITY /BOUNDARY VALUE/
156 C - GSD = TRIANGULARITY
157 
158 C=======================================================================
159 C Input: BA [MA/m^2] Zakharov's input function on the rhs current
160 C BB [MA/m^2] --"--"--"--"--"--"--"--"--"--- presure
161 C BR00 [m] major radius of the edge magnetic surface
162 C SA0 [m] minor radius in the equatorial plane
163 C GL0 [d/l] edge elongation
164 C GD30 [m] edge triangularity
165 C NA1 number of radial points
166 C B0T [T] vacuum magnetic field at the position BR00
167 C PLCUR [MA] total plasma current
168 C Output: GR - rho
169 C GBD - shift
170 C GL - elongation
171 C GSD - triangularity
172 C GRA - <G11>
173 C SQGRA - <sqrt(G11)>
174 C GRAR - <G11/R#2>
175 C AVR2 - <1/R#2>
176 C AI0 - IPOL*RTOR*BTOR
177 C======================================================================|
178  subroutine e3astr
179 C Input:
180  & (ba,bb ! j_zeta = BA*(R_0/r) + BB*(r/R_0-R_0/r)
181  & ,br00 ! R_0+\Delta_edge ! RTOR+SHIFT
182  & ,sa0 ! a_edge ! ABC
183  & ,gl0 ! \lambda_edge ! ELONG
184  & ,gd30 ! \delta_edge ! TRIAN*ABC
185  & ,na1 ! radial grid point No. ! NEQL
186  & ,acc ! accuracy required ! ACEQLB
187  & ,itermx ! Maximum number of iterations
188  & ,b0t ! B_0 at the magnetic axis ! BTOR*RTOR/(RTOR+SHIFT)
189  & ,plcur ! Total plasma current ! IPL
190 C Output:
191  & ,gr ! \Rho(a)
192  & ,gbd ! \Delta(a)
193  & ,gl ! \lambda(a)
194  & ,gsd ! \delta(a)
195  & ,gra ! <g^{11}> = <[nabla(a)]^2>
196  & ,sqgra ! <\sqrt{|g^{11}|}> = <|nabla(a)|>
197  & ,grar ! <g^{11}g^{33}> = <[nabla(a)/r]^2>
198  & ,avr2 ! <g^{33}> = <1/r^2>=G33/RTOR**2
199  & ,ai0 ! I -> IPOL*RTOR*BTOR
200  & ,dgrda ! \prti{\rho}{a}
201  & ,avsqg ! \prti{V}{a}/(4\pi^2)
202  & ,vol ! V(a)
203  & ,b2b0eq ! <B**2/B0**2>
204  & ,b0b2eq ! <B0**2/B**2>
205  & ,bmaxeq ! BMAXT
206  & ,bmineq ! BMINT
207  & ,bmodeq ! <B/BTOR>
208  & ,fofbeq ! <(BTOR/B)**2*(1.-SQRT(1-B/Bmax)*(1+.5B/Bmax))>
209  & ,grdaeq ! <grad a>
210  & ,time)
211 
212 
213 
214 C - BR00,SA0 - MAJOR & MINOR RADII /METER/
215 C - GLO,GD3O - ELONGATION AND TRIANGULARITY /BOUNDARY VALUE/
216 C - GSD = TRIANGULARITY
217 C \Delta' =WDSD1(I)*WSA(I)
218 C \lambda'=WDGL(I)*WGL(I)*WSAA(I)
219 C \delta' =WDSD3(I)*WSA(I)
220 
221  implicit none
222  integer na1,naold,na,nt,nt1,i,i1,j,k,niter,itermx
223  double precision br00,sa0,gl0,gd30,acc,b0t,plcur,time
224  double precision ba(na1),bb(na1),gr(na1),gbd(na1),gl(na1),
225  & gsd(na1),gra(na1),sqgra(na1),grar(na1),avr2(na1),ai0(na1),
226  & dgrda(na1),avsqg(na1),vol(na1),gr2aux(na1)
227  double precision b2b0eq(na1),b0b2eq(na1),bmaxeq(na1),bmineq(na1),
228  1 bmodeq(na1),fofbeq(na1),grdaeq(na1),
229  2 a,aa,c,cc,s,ss,sr,sx,sx1,t,y,y1,sdt,sdt0,
230  3 drda,dzda,drdt,dzdt,dmetr,da2,dgr2,fi,fj,d0,cgp,
231  4 gp,gp2,gr2,glold,g3dold,aold,g22a2,skggg,
232  5 sqg,ylin,yvol,ymin,ymax
233  save aold,glold,g3dold,naold,cgp,niter
234  data aold/0.d0/glold/0.d0/g3dold/0.d0/naold/1/
235  data cgp/3.14159265359d0/ niter/60/
236 
237 *************************************************
238 C if (TIME .gt. .1) then
239 C write(*,*)BR00 ! R_0+\Delta_edge ! RTOR+SHIFT
240 C write(*,*)SA0 ! a_edge ! ABC
241 C write(*,*)GL0 ! \lambda_edge ! ELONG
242 C write(*,*)GD30 ! \delta_edge ! TRIAN*ABC0
243 C write(*,*)NA1
244 C write(*,*)B0T ! B_0 at the magnetic axis ! BTOR*RTOR/(RTOR+SHIFT)
245 C write(*,*)PLCUR ! Total plasma current ! IPL
246 C write(*,'(2(E14.6))')(BA(j),BB(j),j=1,NA1)
247 C stop
248 C endif
249 **************************************************
250 
251  wbbs0=b0t
252  wbj0=0.2*plcur
253  na=na1-1
254  nt=12*max(1,na1/8)
255 C Set initial conditions / zero iteration
256 ! if ( NAOLD .eq. NA
257 ! & .and. abs(AOLD-SA0) .lt. 1.E-6
258 ! & .and. abs(GLOLD-GL0) .lt. 1.E-6
259 ! & .and. abs(G3DOLD-GD30) .lt. 1.E-6
260 ! * ) goto 4
261  call eqgb3(br00,sa0,gl0,gd30,na)
262  aold=sa0
263  glold=gl0
264  g3dold=gd30
265  naold=na
266  wbr00=br00
267  4 continue
268  do i=1,na1
269  wsjp(i)=ba(i)
270  wsp(i)=bb(i)
271  enddo
272  niter=itermx ! Use 60 for the 1st entry only
273  call eqab3(na,nt,niter,acc) ! Call MEM equil solver
274  if (na .le. 1) then
275  na1 = na
276  return
277  endif
278  call eqppab(na)
279  d0=wsd1(na1)*wsaa(na1)
280  gr(1)=0.
281  vol(1)=0.
282  gr2=0.
283  fi=0.
284  s=4.*cgp*cgp
285  do i=1,na1
286  vol(i)=s*wsaa(i)*(wbr0*wsl0(i)+wsaa(i)*wsl1(i))
287  j=i-1
288  fj=fi
289  ai0(i)=wbf(i)*wbr0
290  sqg=wbg0(i)*wgl(i)*wbr0
291  avsqg(i)=wsa(i)*sqg
292  gra(i) =skga(i)/sqg
293  sqgra(i)=sqg22r(i)/sqg
294  grar(i)=wbg22(i)/(wgl(i)*wbr0*sqg)
295  avr2(i)=wbg33(i)*wgl(i)/(wbr0*sqg)
296  gl(i)=wgl(i)
297  gbd(i)=d0-wsd1(i)*wsaa(i)
298  gsd(i)=wsd3(i)*wsaa(i)
299  fi=wbg33(i)*wgl(i)*ai0(i)/(wbr0*wbbs0)
300  if (i.gt.1) then
301  dgr2=fj*wscj1(j)+fi*wsci1(j)
302  gr2=gr2+dgr2
303  if (dgr2.le.0.) then
304  write(*,*)"NA1 = 0",gr2,dgr2
305  na1 = 0
306  return
307  endif
308  gr(i)=sqrt(gr2*2.)
309  dgrda(i)=fi*wsa(i)/gr(i)
310  else
311  dgrda(i)=sqrt(fi)
312  endif
313  enddo
314 CMR extra quantities
315  gr2aux(1)=0.
316  gp=3.1415926
317  gp2=2.*gp
318  nt1=nt+1
319  sdt0=1./nt
320  ylin=0.
321  yvol=0.
322  do 10 i=1,na1
323  a=wsa(i)
324  aa=a*a
325  b2b0eq(i) =0.
326  b0b2eq(i) =0.
327  bmodeq(i) =0.
328  fofbeq(i) =0.
329  grdaeq(i) =0.
330  skggg =0.
331  ymin =99999.
332  ymax =0.
333  ai0(i)=wbf(i)*wbr0
334  IF(i.EQ.1) goto 3
335  k=i-1
336  gr2aux(i)=gr2aux(k)+(skdr(k)*ai0(k)*wscj1(k)+
337  + skdr(i)*ai0(i)*wsci1(k))*2./wbbs0
338  3 CONTINUE
339 C Need to calculate Ymax and Ymin, before doing flux surface averages
340 C ( Y=B^2 )
341  DO 50 k=1,nt1
342  t=sdt0*gp*(k-1)
343  sdt=sdt0
344  IF(k.EQ.1.OR.k.EQ.nt1) sdt=sdt0*0.5
345  c=cos(t)
346  s=sin(t)
347  ss=s*s
348  cc=1.-ss
349  sx1=-wsd1(i)-wsd3(i)*ss
350  sx=c+a*sx1
351  sr=wbr0+a*sx
352 C R,Z derivatives
353  drda=-wdsd1(i)*a+c-wdsd3(i)*a*ss
354  dzda=s*wgl(i)*(aa*wdgl(i)+1)
355  drdt=-a*s-2.*aa*wsd3(i)*c*s
356  dzdt=wgl(i)*a*c
357 C metric tensor components
358  dmetr=drda*dzdt-drdt*dzda
359  skggg=skggg+dmetr*sr*sdt
360  if(i.NE.1)then
361  i1=i-1
362  y=(ai0(i)/sr/wbbs0)**2 + (drdt**2+dzdt**2)*
363  * ((gr2aux(i)-gr2aux(i1))/(wsa(i)-wsa(i1))/
364  / (wsq(i)+wsq(i1))/sr/dmetr)**2
365  else
366  y=(ai0(i)/sr/wbbs0)**2
367  endif
368  ymax=max(y,ymax)
369  ymin=min(y,ymin)
370  50 continue
371  DO 40 k=1,nt1
372  t=sdt0*gp*(k-1)
373  sdt=sdt0
374  IF(k.EQ.1.OR.k.EQ.nt1) sdt=sdt0*0.5
375  c=cos(t)
376  s=sin(t)
377  ss=s*s
378  cc=1.-ss
379  sx1=-wsd1(i)-wsd3(i)*ss
380  sx=c+a*sx1
381  sr=wbr0+a*sx
382 C R,Z derivatives
383  drda=-wdsd1(i)*a+c-wdsd3(i)*a*ss
384  dzda=s*wgl(i)*(aa*wdgl(i)+1)
385  drdt=-a*s-2.*aa*wsd3(i)*c*s
386  dzdt=wgl(i)*a*c
387 c G22/A**2
388  g22a2=ss+4.*a*wsd3(i)*ss*c+(2.*a*wsd3(i)*s*c)**2+
389  + (wgl(i)*c)**2
390 C D**2/A**2
391  da2=(c-a*(wdsd1(i)+wdsd3(i)*ss))*wgl(i)*c+
392  + wgl(i)*ss*(wdgl(i)*aa+1.)*(1.+2.*a*c*wsd3(i))
393  da2=da2*da2
394 C metric tensor components
395  dmetr =drda*dzdt-drdt*dzda
396  if(i.NE.1)then
397  i1=i-1
398  y=(ai0(i)/sr/wbbs0)**2 + (drdt**2+dzdt**2)*
399  * ((gr2aux(i)-gr2aux(i1))/(wsa(i)-wsa(i1))/
400  / (wsq(i)+wsq(i1))/sr/dmetr)**2
401  b2b0eq(i)=b2b0eq(i)+ y*dmetr*sr*sdt
402  b0b2eq(i)=b0b2eq(i)+ dmetr*sr*sdt/y
403  bmodeq(i)=bmodeq(i)+ sqrt(y)*dmetr*sr*sdt
404  IF (y.GT.ymax) ymax=y
405  y1 = sqrt(y/ymax)
406  fofbeq(i)=fofbeq(i)+1./y*(1.-sqrt(abs(1.-y1))
407  + *(1.+.5*y1))*dmetr*sr*sdt
408  grdaeq(i)=grdaeq(i)+ sqrt(g22a2/da2)*dmetr*sr*sdt
409  ylin =ylin+(drdt**2+dzdt**2)*
410  * ((gr2aux(i)-gr2aux(i1))/(wsa(i)-wsa(i1))/
411  / (wsq(i)+wsq(i1))/sr/dmetr)**2*dmetr*sr*sdt
412  yvol = yvol + dmetr*sr*sdt
413  else
414  y=(ai0(i)/sr/wbbs0)**2
415  b2b0eq(i)=b2b0eq(i)+y*sdt
416  b0b2eq(i)=b0b2eq(i)+ sdt/y
417  bmodeq(i)=bmodeq(i)+ sqrt(y)*sdt
418  IF (y.GT.ymax) ymax=y
419  y1 = sqrt(y/ymax)
420  y1 = 1.d0/y*(1.d0-sqrt(abs(1.d0-y1))*(1.d0+.5d0*y1))*sdt
421  fofbeq(i)=fofbeq(i)+y1
422  grdaeq(i)=grdaeq(i)+sqrt(g22a2/da2)*sdt
423  endif
424  40 continue
425  bmaxeq(i)=sqrt(ymax)*wbbs0
426  bmineq(i)=sqrt(ymin)*wbbs0
427  if(i.NE.1)then
428  b0b2eq(i)=b0b2eq(i)/skggg
429  b2b0eq(i)=b2b0eq(i)/skggg
430  bmodeq(i)=bmodeq(i)/skggg
431  fofbeq(i)=fofbeq(i)/skggg
432  grdaeq(i)=grdaeq(i)/skggg
433  endif
434  10 continue
435 CMR
436  ylin=ylin*wbbs0*wbbs0*a*2./(.4*cgp*plcur)**2/br00
437 C open(33,file='dat/lin3')
438 C write(33,*) YLIN
439 !,yvol*A
440 C close(33)
441  30 format(5(f10.4))
442  END subroutine e3astr
443 
444 C - 3 MOMENT EQUILIBRIUM SOLVER AUGUST 17,1988
445 C NITER - max number of iterations
446 C ACC - relative tolerance parameter
447  subroutine eqab3(NA,NT,NITER,ACC)
448  implicit none
449  integer na,nt,niter,na1,i,j,i1,jj
450  double precision
451  1 acc,aa,aai,a4ai,a6ai,s,sc2,w0,w1,w2,w3,bi0,bp1,bp3,
452  2 fu0,fu0i,fu0j,fu1i,fu1j,fu2i,fu2j,fu3i,fu3j,
453  3 fv0i,fv0j,fv1i,fv1j,fv2i,fv2j,fv3i,fv3j
454  na1= na+1
455  iter= 0
456  200 iter= iter+1
457  wsac1= wdsd1(na1)
458  wsac2= wdgl(na1)
459  wsac3= wdsd3(na1)
460  wbr0= wbr00+wsd1(na1)*wsaa(na1)
461  CALL eqk3(na,nt)
462  CALL eqlvu3(na)
463  CALL eqc1(na)
464 C GMC,SW0,SDD1,GDL
465  wgmc(1)=wba(1)*wsl0(1)/wbk0(1)
466  bi0= 0.
467  w0= 0.
468  fu0i= wdba(1)*wsl0(1)
469  fv0i= wdbb(1)*wsv0(1)
470  fu0= wgmc(1)*wsu0(1)
471  wsw0(1)=0.25*(fu0-fu0i)
472  w2= 0.
473  fu2i= wgmc(1)*wsu2(1)-wdba(1)*wsl2(1)
474  fv2i= wdbb(1)*wsv2(1)
475  sc2= 0.25*(wgl(1)**2-1.)
476  wbs2(1)=((wba(1)*wsl22(1)+wbb(1)*(wsv2(1)+sc2*wsv0(1))+
477  , fu2i/6.+sc2*wsw0(1))/wgmc(1)-wbk20(1))/wbk22(1)
478  wdgl(1)=0.
479  w1= 0.
480  fu1i= wgmc(1)*wsu1(1)
481  fv1i= wdba(1)*wsl1(1)+wdbb(1)*wsv1(1)
482  bp1= (wba(1)*wsl1(1)+wbb(1)*wsv1(1)+
483  * 0.25*fu1i)/wgmc(1)-wbk10(1)
484  w3= 0.
485  fu3i= wgmc(1)*wsu3(1)
486  fv3i= wdba(1)*wsl3(1)+wdbb(1)*wsv3(1)
487  bp3= (wba(1)*wsl3(1)+wbb(1)*wsv3(1)+
488  * fu3i/6.)/wgmc(1)-wbk30(1)
489  s= 1./(wbk11(1)*wbk33(1)-wbk13(1)*wbk31(1))
490  wbs1(1)=(bp1*wbk33(1)-bp3*wbk13(1))*s
491  wbs3(1)=(bp3*wbk11(1)-bp1*wbk31(1))*s
492  wdsd1(1)=wbs1(1)
493  DO 1 i=2,na1
494  j= i-1
495  aa= wsaa(i)
496  aai= 1./aa
497  a4ai= aai*aai
498  a6ai= a4ai*aai
499  fu0j= fu0i
500  fu0i= wdba(i)*wsl0(i)
501  fv0j= fv0i
502  fv0i= wdbb(i)*wsv0(i)
503  bi0= bi0+wscj3(j)*fu0j+wsci3(j)*fu0i+
504  * wscj5(j)*fv0j+wsci5(j)*fv0i
505  w0= w0+wscj3(j)*fu0
506  wgmc(i)=(wba(i)*wsl0(i)+wbb(i)*aa*wsv0(i)+(w0-bi0)*aai)
507  , /(wbk0(i)-wsci3(j)*wsu0(i)*aai)
508  fu0= wgmc(i)*wsu0(i)
509  w0= w0+wsci3(j)*fu0
510  wsw0(i)=(w0-bi0)*a4ai
511  fu2j= fu2i
512  fu2i= wgmc(i)*wsu2(i)-wdba(i)*wsl2(i)
513  fv2j= fv2i
514  fv2i= wdbb(i)*wsv2(i)
515  w2= w2+fu2j*wscj5(j)+fu2i*wsci5(j)-
516  * fv2j*wscj7(j)-fv2i*wsci7(j)
517  sc2= 0.25*(wgl(i)**2-1.)
518  wbs2(i)=((wba(i)*wsl22(i)+wbb(i)*(wsv2(i)+sc2*wsv0(i))+
519  , w2*a6ai+sc2*wsw0(i))/wgmc(i)-wbk20(i))/wbk22(i)
520  wdgl(i)=wdgl(j)+(wbs2(j)+wbs2(i))*wscj1(j)
521  fv1j= fv1i
522  fv1i= wdba(i)*wsl1(i)+wdbb(i)*wsv1(i)
523  fu1j= fu1i
524  fu1i= wgmc(i)*wsu1(i)
525  w1= w1+fu1j*wscj3(j)+fu1i*wsci3(j)-
526  * fv1j*wscj5(j)-fv1i*wsci5(j)
527  bp1= (wba(i)*wsl1(i)+wbb(i)*wsv1(i)+
528  * w1*a4ai)/wgmc(i)-wbk10(i)
529  fu3j= fu3i
530  fu3i= wgmc(i)*wsu3(i)
531  fv3j= fv3i
532  fv3i= wdba(i)*wsl3(i)+wdbb(i)*wsv3(i)
533  w3= w3+fu3j*wscj5(j)+fu3i*wsci5(j)-
534  * fv3j*wscj7(j)-fv3i*wsci7(j)
535  bp3= (wba(i)*wsl3(i)+wbb(i)*wsv3(i)+
536  * w3*a6ai)/wgmc(i)-wbk30(i)
537  s= 1./(wbk11(i)*wbk33(i)-wbk13(i)*wbk31(i))
538  wbs1(i)=(bp1*wbk33(i)-bp3*wbk13(i))*s
539  wbs3(i)=(bp3*wbk11(i)-bp1*wbk31(i))*s
540  wdsd1(i)=wbs1(i)
541  1 CONTINUE
542 C SD1,GL
543  w1= 0.
544  wsd1(1)=0.5*wdsd1(1)
545  s= wgl(na1)*exp(-wdgl(na1))
546 C write(*,*)S,WGL(1),WGL(NA1),WDGL(NA1),EXP(-WDGL(NA1)),WDSD1(1)
547  wgl(1)= s
548  wdgl(1)=wbs2(1)
549  DO 2 i=2,na1
550  j= i-1
551  wgl(i)= s*exp(wdgl(i))
552  wdgl(i)=wbs2(i)
553  w1= w1+(wdsd1(j)+wdsd1(i))*wscj1(j)
554  wsd1(i)=w1/wsaa(i)
555  2 CONTINUE
556 C SD3,SDD3
557  fu3j= 1.+wdgl(na1)*wsaa(na1)
558  wdsd3(na1)=wbs3(na1)+2.*wsd3(na1)*fu3j
559  DO 3 i1=1,na
560  j= na1-i1
561  i= j+1
562  aa= wsaa(j)
563  fu3i= fu3j
564  fu3j= 1.+wdgl(j)*aa
565  wsd3(j)=(wsd3(i)*(wsaa(i)-2.*wsci1(j)*fu3i)-(wbs3(i)+
566  , wbs3(j))*wscj1(j))/(aa+2.*wscj1(j)*fu3j)
567  wdsd3(j)=wbs3(j)+2.*wsd3(j)*fu3j
568  3 CONTINUE
569  wsac1= wsa(na1)*abs(wdsd1(na1)-wsac1)
570  wsac2= wsaa(na1)*abs(wdgl(na1)-wsac2)
571  wsac3= wsa(na1)*abs(wdsd3(na1)-wsac3)
572  wsacc= wsac1+wsac2+wsac3
573  IF(iter.LT.niter.AND.wsacc.GT.acc) go to 200
574 
575  write(*,'(2I5,5F10.7)')iter,niter,wsacc,wsac1,wsac2,wsac3
576  if (iter.EQ.niter) then
577  jj = 7
578  write(*,*)char(jj),
579  > ">>> No convergence in the 3M equilibrium solver >>>"
580  write(*,'(1A)')" You can: (0) Continue anyway"
581  write(*,'(1A)')" (1) Try more iterations"
582  write(*,'(2A)')" (2) Freeze the current ",
583  > "equilibrium and continue"
584  write(*,'(1A)')" (3) Stop the run"
585  write(*,'(A,$)')" Your selection > "
586 !dpc
587  write(*,*) 'Continue anyway'
588  goto 4
589 !cpd
590  read(*,*)jj
591  if (jj .ge. 3) stop
592  if (jj .eq. 0) goto 4
593  iter = 0
594  if (jj .eq. 1) goto 200
595  if (jj .eq. 2) na = 0
596  return
597  endif
598 
599 C SD2D1,GD2L,SD2D3
600  4 wd2sd1(1)=wdsd1(1)
601  fu1i= 0.
602  wd2gl(1)=wdgl(1)*wgl(1)
603  fu2i= -wgl(1)
604  wd2sd3(1)=wdsd3(1)
605  fu3i= 0.
606  DO 5 i=2,na1
607  j= i-1
608  aa= wsaa(i)
609  fu1j= fu1i
610  fu1i= aa*(wdsd1(i)-wsd1(i))
611  wd2sd1(i)=(fu1i-fu1j)/wsci1(j)-wd2sd1(j)
612  fu2j= fu2i
613  fu2i= wgl(i)*(aa*wdgl(i)-1.)
614  wd2gl(i)=(fu2i-fu2j)/wsci1(j)-wd2gl(j)
615  fu3j= fu3i
616  fu3i= aa*(wdsd3(i)-wsd3(i))
617  wd2sd3(i)=(fu3i-fu3j)/wsci1(j)-wd2sd3(j)
618  5 CONTINUE
619  RETURN
620  END subroutine eqab3
621 C======================================================================|
622 
623  subroutine eqgb3(BR00,SA0,GL0,GD30,NA)
624 C----------------------------------------------------------------------|
625 C Parameters used:
626 C Input:
627 C Variables: NA,BR00,SA0,GL0,GD30
628 C Arrays:
629 C Output
630 C Arrays defined here and used elsewhere
631 C WSA = j*SA0/NA current minor radius, a
632 C WSAA = WSA**2 a^2
633 C WSCI1 =
634 C WSCI3 =
635 C WSCI5 =
636 C WSCI7 =
637 C WSCJ1 =
638 C WSCJ3 =
639 C WSCJ5 =
640 C WSCJ7 =
641 C Arrays first defined here and modified in EQAB3
642 C WGL = \lambda
643 C WSD3 = \delta/a^2
644 C WDSD3
645 C----------------------------------------------------------------------|
646  implicit none
647  integer na,na1,i,j
648  double precision
649  1 br00,sa0,gl0,gd30,da,d3,d32,aai,aaj,h,c0,ci,cj,cii,cjj
650 C double precision WBR0,WBR00
651 C double precision WSA(1),WSAA(1),WGL(1),WSD3(1),WDSD3(1)
652 C double precision WSCI1(1),WSCI3(1),WSCI5(1),WSCI7(1)
653 C double precision WSCJ1(1),WSCJ3(1),WSCJ5(1),WSCJ7(1)
654  wbr00=br00
655  wbr0=br00
656  da=sa0/na
657  d3=gd30/sa0**2
658  d32=2.*d3
659  na1=na+1
660  DO 1 i=1,na1
661  j=i-1
662  wsa(i)=da*j
663  wsaa(i)=wsa(i)**2
664  wgl(i)=gl0
665  wsd3(i)=d3
666  wdsd3(i)=d32
667 C---C(J,N),C(I,N)
668  IF(j.EQ.0) go to 1
669  aai=wsaa(i)
670  aaj=wsaa(j)
671  h=(aai-aaj)*0.5
672  wscj1(j)=h*0.5
673  wsci1(j)=wscj1(j)
674  cj=aai+2.*aaj
675  ci=aaj+2.*aai
676  c0=h/6.
677  wscj3(j)=c0*cj
678  wsci3(j)=c0*ci
679 C---CJJ=AJ**2N
680  cjj=aaj*aaj
681 C---CII=AI**2N
682  cii=aai*aai
683  cj=aai*cj+3.*cjj
684  ci=aaj*ci+3.*cii
685  c0=c0*0.5
686  wscj5(j)=c0*cj
687  wsci5(j)=c0*ci
688  c0=h*0.05
689  wscj7(j)=c0*(aai*cj+4.*cjj*aaj)
690  wsci7(j)=c0*(aaj*ci+4.*cii*aai)
691  1 CONTINUE
692  RETURN
693  END subroutine eqgb3
694 
695 
696 C - EQK3.FOR - LEFT HAND SIDE AVERAGER AUGUST 17,1988
697 
698  subroutine eqk3(NA,NT)
699 C----------------------------------------------------------------------|
700 C Parameters used:
701 C Input:
702 C Variables: NA,WBR0
703 C Arrays:
704 C Output (arrays):
705 C
706 C----------------------------------------------------------------------|
707  implicit none
708  integer na,nt,na1,nt1,i,k
709  double precision
710  1 cgp,a,aa,aspa,ee,sc1,sc2,t,den,
711  2 da1,da2,ra,ri,xr,g22a2,r0rd,sf2,sf20,sf21,sf3,sf30,sf31,
712  3 xrxr,sk10,sk11,sk13,sk0,sk01,sk2,sr,sx,sxx,sx1,sxx1,szz,
713  4 s,ss,sy1,c,cc,bm1,bm2,bmi,bn,bn0,bn1,bn2,bn12,sdt,sdt0
714  save cgp
715  data cgp/3.14159265359d0/
716  na1= na+1
717  nt1= nt+1
718  sdt0= 1./nt
719  ri= 1./wbr0
720  DO 1 i=1,na1
721  a= wsa(i)
722  aa= a*a
723  aspa= aa*ri
724  ee= wgl(i)**2
725  sc1= 0.25*(ee-1.)
726  sc2= 0.5*(ee+1.)
727  wbk02(i)=0.
728  wbg332(i)=0.
729  wbg222(i)=0.
730  wsu1(i)=0.
731  wsu2(i)=0.
732  wbk20(i)=0.
733  wbk22(i)=0.
734  wbk10(i)=0.
735  wbk11(i)=0.
736  wbk13(i)=0.
737  wbk30(i)=0.
738  wbk31(i)=0.
739  wbk33(i)=0.
740 C------------------------------------------
741 C KONOVALOV LIKES TO DO SOMETHING BY HIMSELF
742 C------------------------------------------
743  skga(i)=0.
744  skdr(i)=0.
745  sqg22r(i)=0.
746 C------------------------------------------
747  DO 2 k=1,nt1
748  t= sdt0*cgp*(k-1)
749  sdt= sdt0
750  IF(k.EQ.1.OR.k.EQ.nt1) sdt=sdt0*0.5
751 C--------------------------------------
752  c= cos(t)
753  s= sin(t)
754  ss= s*s
755  cc= c*c
756  sx1= -wsd1(i)-wsd3(i)*ss
757  sxx1= sx1*sx1
758  szz= ee*ss
759  sx= c+a*sx1
760  sxx= sx*sx
761  sr= wbr0+a*sx
762  sy1= 2.*c*wsd3(i)
763 C---BN0=N0*DT
764  bn0= (ee*cc+ss)*sdt
765 C---BN1=N1*DT
766  bn1= 2.*sy1*ss*sdt
767 C---BN2=N2*DT
768  bn2= ss*sy1*sy1*sdt
769  bn12= bn1+a*bn2
770  bn= bn0+a*bn12
771  bm1= -(wbs1(i)+wbs3(i)*ss)*c
772  bm2= wbs2(i)*ss
773  bmi= 1./(1.+a*bm1+aa*bm2)
774 C EVEN EQUATIONS
775  den= 1./(1.+a*bm1)
776  sk01= bn1-bn0*bm1
777  sk0= (bn2-bn1*(bm1+a*bm2))*bmi+bn0*bm1*bm1*den
778  sk2= -bn0*ss*den*bmi
779  sf20= cc-szz+sc1
780  sf21= 2.*c*sx1
781  sf2= sf20+a*sf21+aa*sxx1
782  wbk02(i)=wbk02(i)+sk0+wbs2(i)*sk2
783  wbk20(i)=wbk20(i)+bn0*sxx1+sk01*sf21+sk0*sf2
784  wbk22(i)=wbk22(i)+sk2*sf2
785 C ODD EQUATIONS
786  den= 1./(1.+aa*bm2)
787  sk0= bn0*den
788  sk10= bn12*bmi
789  sk11= sk0*c*bmi
790  sk13= sk11*ss
791  sf3= cc-3.*szz
792  sf30= c*sf3
793  sf31= sx1*(sxx+c*sx+sf3)
794  sf3= sf30+a*sf31
795  wbk10(i)=wbk10(i)+sk0*sx1+sk10*sx
796  wbk30(i)=wbk30(i)+sk0*sf31+sk10*sf3
797  wbk11(i)=wbk11(i)+sk11*sx
798  wbk31(i)=wbk31(i)+sk11*sf3
799  wbk13(i)=wbk13(i)+sk13*sx
800  wbk33(i)=wbk33(i)+sk13*sf3
801  xr= sx*ri
802  xrxr= xr*xr
803  r0rd= wbr0/sr
804  wbg222(i)=wbg222(i)+bn*bmi*xrxr*r0rd
805  r0rd=r0rd*sdt
806  wbg332(i)=wbg332(i)+(xrxr+bm2-bm1*xr)*r0rd
807  r0rd= r0rd*sx*c
808  wsu1(i)=wsu1(i)+r0rd*szz
809  wsu2(i)=wsu2(i)+r0rd*sxx
810 C HERE ZAKHAROV'S AUTHOR RIGHTS ARE CANCELED WITHOUT ANY CEREMONY
811 c DRDA=-WDSD1(I)*A+C-WDSD3(I)*A*SS
812 c DZDA=S*WGL(I)*(AA*WDGL(I)+1.)
813 c DRDT=-A*S-2.*AA*WSD3(I)*C*S
814 c DZDT=WGL(I)*A*C
815 C METRIC (SUBSRIPT) TENSOR COMPONENTS
816 c G11=DRDA**2+DZDA**2
817 c G22=DRDT**2+DZDT**2
818 c G12=DRDA*DRDT+DZDA*DZDT
819 c G22/A**2
820  g22a2=ss+4.*a*wsd3(i)*ss*c+(2.*a*wsd3(i)*s*c)**2+(wgl(i)*c)**2
821 C D/A
822  da1=wgl(i)*(cc-a*(wdsd1(i)+wdsd3(i)*ss)*c+
823  + ss*(wdgl(i)*aa+1.)*(1.+2.*a*c*wsd3(i)))
824 C---- EMEQ METRIC COMBINATIONS -------
825  skga(i)=skga(i)+g22a2*sr*sdt/da1
826  sqg22r(i)=sqg22r(i)+sqrt(g22a2)*sr*sdt
827 C D**2/A**2
828  da2=(c-a*(wdsd1(i)+wdsd3(i)*ss))*wgl(i)*c+
829  + wgl(i)*ss*(wdgl(i)*aa+1.)*(1.+2.*a*c*wsd3(i))
830  skdr(i)=skdr(i)+sdt*da2/sr
831  2 CONTINUE
832  wbg332(i)=wbg332(i)+(wsd1(i)+0.5*wsd3(i))*ri
833  wbg222(i)=wbg222(i)+wbk02(i)-(wbk10(i)+
834  * wbs1(i)*wbk11(i)+wbs3(i)*wbk13(i))*ri
835  wbd02(i)= 0.5*wbs2(i)
836  ra= 1.-wsd1(i)*aa*ri
837  wbg02(i)=0.5*wbs2(i)*(ra-0.75*wsd3(i)*aspa)-
838  * (wsd1(i)+0.5*(wsd3(i)+wbs1(i)+0.25*wbs3(i)))*ri
839  wbd12(i)= wbg02(i)-wbg332(i)
840  wbk0(i)= sc2+aa*wbk02(i)
841  wdbk00(i)=ee*wdgl(i)
842  wbd0(i)= 1.+aa*wbd02(i)
843  wbg0(i)= 1.+aa*wbg02(i)
844  wbg33(i)= 1.+aa*wbg332(i)
845  wbg22(i)= sc2+aa*wbg222(i)
846  1 CONTINUE
847  END subroutine eqk3
848 C - EQLVU3.FOR - RIGHT HAND SIDE AVERAGER AUGUST 17,1988
849 
850  subroutine eqlvu3(NA)
851 C----------------------------------------------------------------------|
852 C Parameters used:
853 C Input:
854 C Variables: NA,WBR0
855 C Arrays:
856 C WSAA(*)
857 C WGL(*)
858 C WSD1(*)
859 C WSD3(*)
860 C Output
861 C Arrays defined here and used elsewhere
862 C WSL0
863 C WSL1
864 C WSL2
865 C WSL3
866 C WSL22
867 C WSV0
868 C WSV1
869 C WSV2
870 C WSV3
871 C Arrays modified here, first defined in EQK3, used in EQAB3
872 C WSU0
873 C WSU1
874 C WSU2
875 C WSU3
876 C----------------------------------------------------------------------|
877  implicit none
878 C double precision WBR0,WSAA(1),WGL(1),WSD1(1),WSD3(1)
879 C double precision WSL0(1),WSL1(1),WSL2(1),WSL3(1),WSL22(1)
880 C double precision WSV0(1),WSV1(1),WSV2(1),WSV3(1)
881 C double precision WSU0(1),WSU1(1),WSU2(1),WSU3(1)
882  integer na,na1,i
883  double precision
884  1 cr2,cr3,cr4,cr6,cr8,aa,e,ee,t1,t3,t5,rbr0,
885  2 s,s1,s2,x20,x22,x30,x32,x40,x42,x50,x60,x5x1,x6x1,
886  3 y0,y2,y4,ux0,ux2,ux12,ux30,ut0,ut1,ut2,ut4,ut6,
887  4 uy0,uy1,uy2,uy3,uy4,uz1,uz3
888  data cr2 /5.d-1/,cr3/3.33333333d-1/,cr6/1.66666666d-1/,
889  , cr4/2.5d-1/,cr8/1.25d-1/
890  na1= na+1
891  rbr0=1./wbr0
892  do 1 i=1,na1
893  aa= wsaa(i)
894  e= wgl(i)
895  ee= e*e
896  ux2= cr2*wsd3(i)
897  ux0= -wsd1(i)-ux2
898 C-------------------------------
899  uy4= cr2*ux2*ux2
900  uy0= ux0*ux0+uy4
901  y0= cr2+uy0*aa
902  uy1= 2.*ux0+ux2
903  uy2= 2.*ux0*ux2
904  y2= cr2+uy2*aa
905  uy3= ux2
906  y4= uy4*aa
907  ut0= ux0*y0+cr2*(uy1+ux2*y2)
908  ut1= ux0*uy1+uy0+cr2*(uy2+ux2*(uy3+uy1))
909  t1= 0.75+ut1*aa
910  ut2= ux0*y2+y0*ux2+cr2*(uy3+uy1+ux2*y4)
911  t3= aa*ux0*uy3+cr2*(y2+y4+aa*ux2*uy1)
912  ut4= ux0*y4+cr2*(uy3+ux2*y2)
913  t5= cr2*(y4+aa*ux2*uy3)
914  ut6= cr2*ux2*y4
915  uz1= uy1*(2.*y0+y2)+uy3*(y2+y4)
916  uz3= uy3*2.*y0+uy1*(y4+y2)
917  x5x1= t1*y0+aa*ut0*uy1+cr2*(y2*(t1+t3)+
918  , y4*(t3+t5)+aa*(ut2*uy1+uy3*(ut2+ut4)))
919  x6x1= t1*(2.*ut0+ut2)+t3*(ut2+ut4)+t5*(ut4+ut6)
920  s= ee*cr8
921  x20= cr4*uy1
922  x22= s*(uy1-uy3)*cr2
923  x30= cr6*t1
924  x32= s*(t1-t3)*cr3
925  ux30= cr6*ut1
926  x40= cr8*uz1
927  x42= s*(uz1-uz3)*0.25
928  x50= x5x1*0.1
929  x60= cr2*x6x1*cr6
930 C-------------------------------
931  s1=(1.-ee)*cr4
932  wsl0(i)= e*cr2
933  wsl22(i)= e*ux30
934  wsl2(i)= wsl0(i)*s1+aa*wsl22(i)
935  wsl1(i)= e*x20
936  wsl3(i)= e*(x40-3.*x22)
937  s= 2.*e*rbr0
938  s1= cr2*rbr0
939  s2= aa*s1
940  wsv0(i)= s*(x20+s1*x30)
941  wsv1(i)= s*(x30+s2*x40)
942  wsv2(i)= s*(x40-x22+s1*(x50-x32))
943  wsv3(i)= s*(x50-3.*x32+s2*(x60-3.*x42))
944 
945  ux30= wsu2(i)
946  ux12= wsu1(i)
947  s= ee*rbr0
948  s2= aa*rbr0
949  wsu1(i)= s*(cr2+s2*(ux30*rbr0-2.*x20))
950  wsu0(i)= -rbr0*wsu1(i)
951  wsu2(i)= s*(2.*x20+rbr0*(ux12-ux30))
952  wsu3(i)= s*(ux30-3.*ux12)
953  1 CONTINUE
954  RETURN
955  END subroutine eqlvu3
956 
957  subroutine eqc1(NA)
958 C----------------------------------------------------------------------|
959 C Parameters used:
960 C Input:
961 C Variables: NA,WBR0,WBR00
962 C Arrays: WSA,WSP,WSJP
963 C Output:
964 C Arrays:
965 C WBA Zakharov's function A
966 C WBB Zakharov's function B
967 C WDBA h/a*dA/da
968 C WDBB
969 C----------------------------------------------------------------------|
970  implicit none
971 C double precision WBR0,WBR00
972 C double precision WSA(1),WSP(1),WSJP(1),WBA(1),WBB(1),WDBA(1),WDBB(1)
973 
974  integer na,na1,i
975  double precision h, s, rs, ss
976  na1=na+1
977  wdba(1)=2.*(wsjp(2)-wsjp(1))/wsa(2)**2
978  wba(1)=wsjp(1)
979  wbb(1)=wsp(1)
980  wdbb(1)=2.*(wsp(2)-wsp(1))/wsa(2)**2
981  h=.5/wsa(2)
982  do 1 i=2,na
983  wba(i)=wsjp(i)
984  wbb(i)=wsp(i)
985  wdbb(i)=(wsp(i+1)-wsp(i-1))*h/wsa(i)
986 1 wdba(i)=(wsjp(i+1)-wsjp(i-1))*h/wsa(i)
987  wba(na1)=wsjp(na1)
988  wdba(na1)=(wsjp(na1)-wsjp(na))/wsa(2)/wsa(na1)
989  wbb(na1)=wsp(na1)
990  wdbb(na1)=(wsp(na1)-wsp(na))/wsa(2)/wsa(na1)
991  s=wbr0/wbr00
992  rs=1./s
993  ss=s-rs
994  do 2 i=1,na1
995  wba(i)=wba(i)*rs+wbb(i)*ss
996  wbb(i)=wbb(i)*s
997  wdba(i)=wdba(i)*rs+wdbb(i)*ss
998  wdbb(i)=wdbb(i)*s
999 2 continue
1000  end subroutine eqc1
1001 
1002  subroutine eqppab(NA)
1003 C----------------------------------------------------------------------|
1004 C gm0 = 0.4*cgp
1005 C WSA(I) = a[m], WSAA(I) = a**2
1006 C WBR00 = (RTOR+SHIFT) = R[m]
1007 C WBR0 = R0[m]
1008 C WBJ0 = 0.2*I[MA]
1009 C WBBS0 = Bs[T],
1010 C WBF(I)=0.2*F/R0[MA/m], WBFF(I) = WBF(I)**2
1011 C WSJ(I)=gm0*<j>[MA/m**2], WSJP(I) = <jB>/B
1012 C WSJSL(I) = j(r,gt=cgp), WSJSR(I) = j(r,gt=0)
1013 C WSP(I) = gm0*p[MJ/m**3]
1014 C WGP(I) = gP[VS]/(2*cgp*R0)
1015 C WDSQRQ(I) = aq'/q
1016 C WGPINT, WGPRES, WSLI are defined but not used (not present in commons
1017 C----------------------------------------------------------------------|
1018 C Parameters used:
1019 C Input:
1020 C Variables: NA,WBR0,WBR00,WBJ0,WBBS0
1021 C Arrays:
1022 C WSA(*)
1023 C WSAA(*)
1024 C WGL(*)
1025 C WSD1(*)
1026 C WSD3(*)
1027 C Output
1028 C Variables:
1029 C WGB
1030 C WGB0
1031 C WGBD
1032 C WGBJ
1033 C WGBST
1034 C WGMJ
1035 C WGMJEX
1036 C WSQC
1037 C WGPINT
1038 C WGPRES
1039 C WSLI
1040 C Arrays:
1041 C WSL0
1042 C WSL1
1043 C WSL2
1044 C WSL3
1045 C WSL22
1046 C WSV0
1047 C WSV1
1048 C WSV2
1049 C WSV3
1050 C----------------------------------------------------------------------|
1051  implicit none
1052 C double precision WBR0,WBR00,WBJ0,WBBS0
1053 C double precision WGB,WGB0,WGBD,WGBJ,WGBST,WGMJ,WGMJEX,WSQC
1054 C double precision WSP(1),WSA(1),WSAA(1),WGL(1),WGP(1),WSJ(1),WSJP(1)
1055 C double precision WBA(1),WBB(1),WBF(1),WBFF(1),WSD1(1)
1056 C double precision WDBA(1),WDBB(1),WDGL(1),WBK0(1),WBK02(1),WDBK00(1)
1057 C double precision WDSP(1),WDSJ(1),WSQ(1),WDSQ(1),WDSQRQ(1),WDSJP(1)
1058 C double precision WBD0(1),WBD12(1),WBD02(1),WBG22(1),WBG33(1),WBG332(1)
1059 C double precision WSCI3(1),WSCJ3(1),WSCJ1(1)
1060 C double precision WSL0(1),WSU0(1),WSW0(1),WSV0(1),WGMC(1),WDGMC(1)
1061 C double precision WSJSL(1),WSJSR(1)
1062  integer na,na1,i,j
1063  double precision
1064  1 s,sli,slj,g33i,g33j,cgp,dg33,bj,bji,bjj,rr,rl,
1065  2 dl0,dl0i,dl0j,d12i,d12j,dd12,dv0,ff0,ffi,ffj,v0i,v0j,
1066  3 gb,gbi,gbj,gbs,gbsi,gbsj,gbd,gbdi,gbdj,gpi,gpj,gmji,gmjj,
1067  4 wgpint,wgpres,wsli
1068  save cgp
1069  data cgp/3.14159265359d0/
1070  na1=na+1
1071 C --- GP,SJ,DSJ,BJ,SJL,SJR
1072 
1073  g33i=0.
1074  dg33=2.*(wba(1)*(wbg332(1)-wbd02(1))+wbb(1)*wbd12(1))
1075  wsj(1)=wba(1)
1076  wdsj(1)=wdba(1)+dg33
1077  bji=wsj(1)*wgl(1)
1078  bj=0.
1079 
1080  DO 1 i=2,na1
1081  j=i-1
1082  g33j=g33i
1083  g33i=(wba(i)*(wbg332(i)-wbd02(i))+wbb(i)*wbd12(i))*
1084  , wsaa(i)/wbd0(i)
1085  dg33=(g33i-g33j)/wscj1(j)-dg33
1086  wsj(i)=wba(i)+g33i
1087  wdsj(i)=wdba(i)+dg33
1088  bjj=bji
1089  bji=wsj(i)*wbd0(i)*wgl(i)
1090  bj=bj+(bjj+bji)*wscj1(j)
1091 1 CONTINUE
1092 
1093  s=wbj0/bj
1094 
1095  DO 2 i=1,na1
1096  wba(i) =wba(i)*s
1097  wdba(i) =wdba(i)*s
1098  wbb(i) =wbb(i)*s
1099  wdbb(i) =wdbb(i)*s
1100  wsj(i) =wsj(i)*s
1101  wdsj(i) =wdsj(i)*s
1102  wgmc(i) =wgmc(i)*s
1103  wsw0(i) =wsw0(i)*s
1104 
1105  rr=wbr0-wsd1(i)*wsaa(i)
1106  rl=wbr0/(rr-wsa(i))
1107  rr=wbr0/(rr+wsa(i))
1108  wsjsl(i)=wba(i)*rl+wbb(i)*(1./rl-rl)
1109  wsjsr(i)=wba(i)*rr+wbb(i)*(1./rr-rr)
1110 
1111  wdsp(i)=-wbb(i)*wgmc(i)*wgl(i)
1112 2 CONTINUE
1113 
1114  gpi=wgmc(1)*wgl(1)
1115  wgp(1)=0.
1116  sli=wgp(1)*wsj(1)*wgl(1)*wbd0(1)
1117  wsli=0.
1118 
1119  ffi=(wba(1)-wbb(1))*gpi
1120  wbff(1)=0.
1121  wsp(1)=0.
1122 
1123  gbi=wbb(1)*gpi*wsl0(1)
1124  gb=0.
1125  gbsi=0.
1126  gbs=0.
1127  gbdi=(wba(1)-wbb(1))*wgmc(1)*wsu0(1)
1128  gbd=0.
1129 
1130  DO 3 i=2,na1
1131  j=i-1
1132 
1133  gpj=gpi
1134  gpi=wgmc(i)*wgl(i)
1135  wgp(i)=wgp(j)-(gpj+gpi)*wscj1(j)
1136  slj=sli
1137  sli=wgp(i)*wsj(i)*wgl(i)*wbd0(i)
1138  wsli=wsli+(slj+sli)*wscj1(j)
1139 
1140  ffj=ffi
1141  ffi=(wba(i)-wbb(i))*gpi
1142  wbff(i)=wbff(j)+(ffj+ffi)*wscj1(j)
1143  wsp(i)=wsp(j)+(wdsp(j)+wdsp(i))*wscj1(j)
1144 
1145  gbj=gbi
1146  gbi=wbb(i)*gpi*wsl0(i)
1147  gb=gb+gbj*wscj3(j)+gbi*wsci3(j)
1148 
1149  gbsj=gbsi
1150  gbsi=wsp(i)*gbi
1151  gbs=gbs+gbsj*wscj3(j)+gbsi*wsci3(j)
1152 
1153  gbdj=gbdi
1154  gbdi=(wba(i)-wbb(i))*wgmc(i)*wsu0(i)
1155  gbd=gbd+gbdj*wscj3(j)+gbdi*wsci3(j)
1156 3 CONTINUE
1157  wgpint=-2.*cgp*wbr0*wgp(na1)
1158 c
1159 c inernal iductance
1160 c
1161  wsli=2.*wbr0/(wbr00*wbj0**2)*(wsli-wgp(na1)*wbj0)
1162  wgpres=cgp*wbr00*wsli*wbj0
1163 c
1164 c betta j
1165 c
1166  wgbj=4.*gb/(wbj0**2)
1167  s=wsl0(na1)*wsaa(na1)
1168 c
1169 c betta
1170 c
1171  wgb=2.*gb/(wbbs0**2*s)
1172 c
1173 c betta*
1174 c
1175  wgbst=2.*sqrt(2.*(gbs-wsp(na1)*gb)/s)/wbbs0**2
1176 !DPC ff0 here is not yet defined, moved ff0 calculation up
1177  ff0=(wbbs0*wbr00/wbr0)**2
1178  wgbd=2.*gb/(wbr0**2*(2.*gbd-ff0*wsu0(na1)*wsaa(na1)/wgl(na1)))
1179  wgmj=4.*gbd*(wbr0/wbj0)**2
1180  wsp(1)=-wsp(na1)
1181 !moved up FF0=(WBBS0*WBR00/WBR0)**2
1182  wbff(1)=ff0+2.*wbff(na1)
1183  wbf(1)=sqrt(wbff(1))
1184  wgb0=2.*wsp(1)/wbbs0**2
1185  wsqc=wbbs0*wsaa(na1)*(wgl(na1)**2+1.)*.5/(wbj0*wbr00)
1186 
1187 C --- SJP,SDJP,SP,BFF,DGMC,SQ,SDQ,GMJ,GMJEX
1188  dv0=2.*(wbb(1)*wsv0(1)+wsw0(1)-wgmc(1)*wbk02(1))
1189  v0i=0.
1190  dl0i=(wdgl(1)-2.*wbd02(1))*wgl(1)
1191  dl0=0.
1192  wdgmc(1)=(wdba(1)*wsl0(1)+wba(1)*(wgl(1)*wbd02(1)+.25*dl0i)-
1193  , wgmc(1)*wdbk00(1)+dv0)/wbk0(1)
1194 
1195  ffi=(wba(1)-wbb(1))*wgmc(1)*wgl(1)
1196  wsq(1)=wbf(1)/(wbr0*wgmc(1))
1197  dg33=2.*wbg332(1)*wsq(1)
1198  g33i=0.
1199  wdsq(1)=wsq(1)*(-ffi/wbff(1)-wdgmc(1)/wgmc(1))+dg33
1200  wdsqrq(1)=0.
1201 
1202  dd12=2.*((wba(1)-wbb(1))*wbg22(1)/(wbr0*wsq(1))**2+
1203  , wbb(1)*wbd12(1))
1204  d12i=0.
1205  wsjp(1)=wba(1)
1206  wdsjp(1)=wdba(1)+dd12
1207 
1208  bji=wgl(1)
1209  bj=0.
1210  gmji=ffi/wbf(1)*bji*.5
1211  wgmjex=0.
1212 
1213  DO 4 i=2,na1
1214  j=i-1
1215  wsp(i)=wsp(i)-wsp(na1)
1216 C --- DERIVATIVES
1217  s=1./wscj1(j)
1218  gpi=wgmc(i)*wgl(i)
1219 
1220  wbff(i)=wbff(1)-2.*wbff(i)
1221  wbf(i)=sqrt(wbff(i))
1222 
1223  dl0j=dl0i
1224  dl0i=(wdgl(i)-2.*wbd02(i))*wgl(i)
1225  dl0=dl0+dl0j*wscj3(j)+dl0i*wsci3(j)
1226  v0j=v0i
1227  v0i=(wbb(i)*wsv0(i)+wsw0(i)-wgmc(i)*wbk02(i))*wsaa(i)
1228  dv0=(v0i-v0j)*s-dv0
1229  wdgmc(i)=(wba(i)*(wgl(i)*wbd02(i)+dl0/wsaa(i)**2)+
1230  , wdba(i)*wsl0(i)-wgmc(i)*wdbk00(i)+dv0)/(wgl(i)**2+1.)*2.
1231 
1232  ffi=(wba(i)-wbb(i))*wgl(i)*wgmc(i)
1233  g33j=g33i
1234  wsq(i)=wbf(i)/(wbr0*wgmc(i))
1235  g33i=wsq(i)*wbg332(i)*wsaa(i)
1236  dg33=(g33i-g33j)*s-dg33
1237  wdsq(i)=wsq(i)*(-ffi/wbff(i)-wdgmc(i)/wgmc(i))+dg33
1238  wsq(i)=wsq(i)+g33i
1239  wdsqrq(i)=wsaa(i)*wdsq(i)/wsq(i)
1240 
1241  d12j=d12i
1242  d12i=((wba(i)-wbb(i))*wbg22(i)*wbg33(i)/(wbr0*wsq(i))**2+
1243  , wbb(i)*wbd12(i)/wbg33(i))*wsaa(i)
1244  dd12=(d12i-d12j)*s-dd12
1245  wsjp(i)=wba(i)+d12i
1246  wdsjp(i)=wdba(i)+dd12
1247 
1248  bjj=bji
1249  bji=wbg33(i)*wgl(i)
1250  bj=bj+(bjj+bji)*wscj1(j)
1251  gmjj=gmji
1252  gmji=ffi*bj/(wbf(i)*wsaa(i))
1253  wgmjex=wgmjex+gmjj*wscj3(j)+gmji*wsci3(j)
1254 
1255 4 CONTINUE
1256 
1257  wgmjex=4.*wgmjex*wbbs0/wbj0**2
1258 
1259  END subroutine eqppab
1260 C======================================================================|
1261 
1262  end module emeq_equil
emeq_equil::dealloc_emeq_equil
subroutine, public dealloc_emeq_equil
Definition: emeq_equil.f:93
emeq_equil::alloc_emeq_equil
subroutine, public alloc_emeq_equil(NP_in)
Definition: emeq_equil.f:53
emeq_equil::e3astr
subroutine, public e3astr
Definition: emeq_equil.f:178
emeq_equil
replacement for the EMEQ common block
Definition: emeq_equil.f:8