V3FIT
datain.f
1  SUBROUTINE datain(extension, iunit_in, iunit, ians)
2 c--
3 c READ and initialize data for BOOTSJ
4 c--
5 C-----------------------------------------------
6 C M o d u l e s
7 C-----------------------------------------------
8  USE parambs
9  USE vmec0
10  USE read_boozer_mod
11  USE bootsj_input
12  IMPLICIT NONE
13 C-----------------------------------------------
14 C D u m m y A r g u m e n t s
15 C-----------------------------------------------
16  INTEGER :: iunit, ians, iunit_in
17  CHARACTER*(*) :: extension
18 C-----------------------------------------------
19 C L o c a l P a r a m e t e r s
20 C-----------------------------------------------
21  REAL(rprec), PARAMETER :: one = 1, zero = 0
22 C-----------------------------------------------
23 C L o c a l V a r i a b l e s
24 C-----------------------------------------------
25  INTEGER :: i, n, m, i1, ntheta_min, nzeta_min,
26  1 mn, ir, idum
27  REAL(rprec) :: temp
28  REAL(rprec), DIMENSION(:), ALLOCATABLE :: work
29  REAL(rprec) :: status, tempe0, tempi0, pres10, pres0, a, a1
30 C-----------------------------------------------
31 !
32 ! Read in boozer |B|, iota, phip, pres AND ALLOCATE needed arrays
33 !
34  CALL read_boozer(extension)
35 
36  jlist = 0
37  READ (iunit_in, *,iostat=idum) jlist
38  IF (idum .gt. 0) stop ' Error reading bootsj input file in DATAIN'
39 
40  CLOSE (iunit_in)
41 
42 
43 !
44 !
45 !- give default values of control data otherwise read from the file 'bootin'
46 !
47 ! nrho is an obsolete variable that is no longer used. it is set to ns-1
48 ! where ns is the number of whole mesh points used in VMEC
49 !
50 ! nrho is an obsolete variable that is no longer used. it is set to ns-1
51 ! where ns is the number of whole mesh points used in VMEC
52 !
53  nrho = 30 !number of rho values to use
54  mbuse = 6 !number of m (poloidal) terms in B field.
55  nbuse = 4 !number of nzetah (toroidal) terms in B field.
56  zeff1 = 1.0_dp !effective ion charge
57  dens0 = 0.3_dp !central electron density in 10**20 m-3
58  teti = 2.0_dp !ratio of Te/Ti for electron to ion
59  !temperature profiles
60  tempres = -one !tempe1(i)=pres(ir)**tempres
61  !IF(tempres.eq.-1.0_dp) THEN
62  !tempe1(i)=SQRT(pres(i))
63  damp = -0.01_dp !superceded by damp_bs
64  damp_bs = -0.01_dp !damping factor for resonances
65  isymm0 = 0 !IF ne.0 THEN force a symmetric-device calculation,
66  tempe1 = one !1 keV DEFAULT temperature
67  tempi1 = zero
68  ate = zero
69  ati = zero
70  ate(0) = -one
71  ati(0) = -one
72 !
73 !- Read control data
74 !
75  CALL read_namelist (iunit, i, 'bootin')
76  IF (i .ne. 0 ) stop 'Error reading bootin NAMELIST'
77 
78 ! If the new damping variable is not READ in, set to to the value this gives
79 ! an unchanged damping factor, namely damp_bs**2 = damp. If, in addition,
80 ! damp is not READ in, THEN set damp_bs = 0.001
81 
82  IF(damp_bs .lt. zero) THEN !in this CASE no damp_bs was READ in
83  IF(damp .gt. zero) THEN
84  damp_bs = sqrt(damp)
85  ELSE
86  damp_bs = 0.001_dp !in this CASE no damp was READ in
87  ENDIF
88  ENDIF
89  teti = abs(teti)
90 
91 !
92 ! CHECK DIMENSION SIZE CONSISTENCY
93 !
94  IF(nboz_b .lt. nbuse) THEN
95  IF (lscreen) WRITE(*,*) 'nbuse > nbos_b, nbuse = nboz_b'
96  nbuse = nboz_b
97  ENDIF
98  IF(mboz_b .lt. mbuse) THEN
99  IF (lscreen) WRITE(*,*) 'mbuse > mbos_b, mbuse = mboz_b'
100  mbuse = mboz_b
101  ENDIF
102 
103  nzeta_min = 2*nbuse + 1
104  ntheta_min = 2*mbuse + 1
105 
106  DO i = 0, 6
107  nzetah = 4*2**i
108  IF(nzetah .gt. nzeta_min) EXIT
109  nzetah = 2*2**i * 3
110  IF(nzetah .gt. nzeta_min) EXIT
111  ENDDO
112 
113  DO i = 0, 6
114  nthetah = 4*2**i
115  IF(nthetah .gt. ntheta_min) EXIT
116  nthetah = 2*2**i * 3
117  IF(nthetah .gt. ntheta_min) EXIT
118  ENDDO
119 
120  IF(lscreen) print *, 'mbuse = ',mbuse,'nbuse = ',
121  1 nbuse,'nthetah = ',nthetah, ' nzetah = ',nzetah
122 
123  90 FORMAT(5e16.8)
124 
125 ! convert bmn's to amnfit's (positive m's ONLY to positive n's only)
126 
127  amnfit = zero
128 
129  lsurf = .false.
130  status = tiny(a1)
131  DO ir = 1, irup
132  DO mn = 1,mnboz_b
133  m = ixm_b(mn)
134  n = ixn_b(mn)/nfp_b
135  IF (m .gt. mboz_b) stop 'boozmn indexing conflict, m'
136  IF (abs(n) .gt. nboz_b) stop 'boozmn indexing conflict, n'
137  IF (n.lt.0) THEN
138  m = -m
139  n = -n
140  END IF
141  IF (m.eq.0 .and. n.eq.0 .and. bmnc_b(mn,ir).gt.status)
142  1 lsurf(ir) = .true.
143  amnfit(ir,m,n) = bmnc_b(mn,ir+1) !!2nd half grid == 1st grid pt. here
144  END DO
145  END DO
146 
147  CALL read_boozer_deallocate
148 
149  zeff1 = max(one,zeff1)
150 c setup s grid. A nonuniform mesh is allowed.
151 
152  psimax = maxval(abs(flux))
153 
154 c we need to keep the sign of psimax
155  IF(flux(irup) .lt. zero) psimax = -psimax
156 
157 c first normalize the mesh to 1
158 c and THEN shift from the full to half mesh with the first point on 1, not 2
159 c also need to calulate the deltas for differencing on the vmec grid
160 c special values at the ENDs will be calculated in bongrid
161 
162  DO ir = 1, irup
163  rhoar(ir) = 0.5_dp*(flux(ir) + flux(ir+1))/psimax
164  d_rho(ir) = (flux(ir+1) - flux(ir))/psimax
165  END DO
166 
167 c WRITE (ijbs, 130) irup, psimax
168 c 130 FORMAT(' Last flux surface used is number ',i2,' with PSIMAX =',1p
169 c 1 e11.4,' WB')
170 c
171 c- Switch sign of q, IF desired.
172 c
173  IF (iotasign .lt. 0) THEN
174  qsafety(:irup) = iotasign*qsafety(:irup)
175  ENDIF
176 
177 c
178 c !Boozer I/g values
179  aiogar(:irup) = aipsi(:irup)/(gpsi(:irup)+1.0e-36_dp)
180 c
181  CALL positiv (pres1, irup, 2) !to be sure that arrays are positive
182  CALL positiv (betar, irup, 2)
183 c evaluate electron and ion temperature profiles on vmec mesh
184  IF (any(ate .ne. zero)) THEN
185  DO ir = 1, irup
186  tempe1(ir) = temp(rhoar(ir), ate)
187  END DO
188  END IF
189  IF (any(ati .ne. zero)) THEN
190  DO ir = 1, irup
191  tempi1(ir) = temp(rhoar(ir), ati)
192  END DO
193  END IF
194 
195  tempe0 = tempe1(1) !central electron temperature in keV
196  tempi0 = tempi1(1) !central ion temperature in keV
197  pres10 = pres1(1) !central total pressure in Pa
198  pres0 = 1.6022e4_dp
199 c If the leading coefficient on either te or ti is negative, the intent is
200 c to assume a maximum density, and then calcuate the profiles using a
201 c power law given by ABS(tempres). this coefficient must originally be
202 c be negative and ca not be greater than one. That is to say, the profiles
203 c are determined by ne(0) and tempres. The ratio of Te to Ti is given by
204 c teti.
205  IF (tempe0.le.zero .or. tempi0.le.zero) tempres = abs(tempres)
206  tempres = min(one,tempres)
207 c
208  IF (tempres .ge. zero) THEN !in that CASE, calculate the temperature
209  teti = teti + 1.e-36_dp
210  a = one + one/(zeff1*teti)
211 c !central electron temperature in keV
212  tempe0 = pres10/(a*pres0*dens0)
213  tempi0 = tempe0/teti !central ion temperature in keV
214  tempe1(:irup) = pres1(:irup)**tempres
215 c !suggested Te and Ti profiles are the same
216  tempi1(:irup) = tempe1(:irup)
217  a = tempe0/tempe1(1)
218  a1 = tempi0/tempi1(1)
219  tempe1(:irup) = tempe1(:irup)*a
220  tempi1(:irup) = tempi1(:irup)*a1
221  ENDIF
222 c
223  CALL positiv (tempe1, irup, 2)
224  CALL positiv (tempi1, irup, 2)
225 
226  dense(:irup) = pres1(:irup)/(pres0*(tempe1(:irup)+tempi1(:irup)/
227  1 zeff1)+1.e-36_dp)
228  ALLOCATE(work(irdim))
229  CALL smooth1 (dense, 1, irup, work, zero)
230  CALL positiv (dense, irup, 2)
231  i1 = irup - 1
232  a = tempe1(irup) + tempi1(irup)/zeff1
233  a1 = tempe1(i1) + tempi1(i1)/zeff1
234  dense(irup) = dense(i1)*a1*betar(irup)/(a*betar(i1)+1.e-36_dp)
235 c the above game is apparently to control the spline derivatives at the
236 c edge.
237  densi(:irup) = dense(:irup)/zeff1
238  dens0 = dense(1) !central electron density
239 
240 
241 c
242 c- Echo input data to output file.
243 c
244  WRITE (ians, bootin)
245 
246  DEALLOCATE (work)
247 
248  END SUBROUTINE datain