1 SUBROUTINE readin(input_file, iseq_count, ier_flag, lscreen)
7 USE mgrid_mod,
ONLY: nextcur, curlabel, nfper0, read_mgrid
9 USE parallel_include_module,
ONLY: grank, mgrid_file_read_time,
11 USE parallel_vmec_module,
ONLY: runvmec_comm_world
16 INTEGER :: iseq_count, ier_flag
18 CHARACTER(LEN=*) :: input_file
22 INTEGER :: iexit, ipoint, n, iunit, ier_flag_init,
23 & i, ni, m, nsmin, igrid, mj, isgn, ioff, joff,
25 REAL(dp),
DIMENSION(:,:),
POINTER ::
26 & rbcc, rbss, rbcs, rbsc, zbcs, zbsc, zbcc, zbss
27 REAL(dp) :: rtest, ztest, tzc, trc, delta
28 REAL(dp),
ALLOCATABLE :: temp(:)
29 CHARACTER(LEN=100) :: line, line2
30 CHARACTER(LEN=1) :: ch1, ch2
246 CALL second0(treadon)
248 lwrite = (grank .EQ. 0)
249 ier_flag_init = ier_flag
250 ier_flag = norm_term_flag
251 IF (ier_flag_init .EQ. more_iter_flag)
GOTO 1000
256 CALL read_indata(input_file, iunit, ier_flag)
257 IF (ier_flag .NE. norm_term_flag)
RETURN
259 IF (tensi2 .EQ. zero ) tensi2 = tensi
266 CALL heading(input_extension, time_slice,
267 & iseq_count, lmac, lscreen, lwrite)
273 rewind(iunit, iostat=iexit)
275 DO WHILE(iexit .EQ. 0)
276 READ (iunit,
'(a)', iostat=iexit) line
277 IF (iexit .NE. 0)
EXIT
278 iexit = index(line,
'INDATA')
279 iexit = iexit + index(line,
'indata')
280 ipoint = index(line,
'!')
281 IF (ipoint .EQ. 1)
WRITE (nthreed, *) trim(line)
292 CALL read_mgrid (mgrid_file, extcur, nzeta, nfp,
293 & lscreen, ier_flag, comm = runvmec_comm_world)
295 mgrid_file_read_time = mgrid_file_read_time + (tzc - trc)
297 IF (lfreeb .AND. lscreen .AND. lwrite)
THEN
298 WRITE (6,
'(2x,a,1p,e10.2,a)')
'Time to read MGRID file: ',
300 IF (ier_flag .ne. norm_term_flag)
RETURN
301 IF (lwrite)
WRITE (nthreed,20) nr0b, nz0b, np0b, rminb,
302 & rmaxb, zminb, zmaxb, trim(mgrid_file)
303 20
FORMAT(//,
' VACUUM FIELD PARAMETERS:',/,1x,24(
'-'),/,
304 &
' nr-grid nz-grid np-grid rmin rmax zmin',
305 &
' zmax input-file',/,3i9,4f10.3,5x,a)
314 IF (ns_array(1) .eq. 0)
THEN
315 ns_array(1) = min(nsin,nsd)
317 ns_array(multi_ns_grid) = ns_default
320 DO WHILE (ns_array(multi_ns_grid) .ge. nsmin .and.
321 & multi_ns_grid .lt. 100)
322 nsmin = max(nsmin, ns_array(multi_ns_grid))
323 IF (nsmin .le. nsd)
THEN
324 multi_ns_grid = multi_ns_grid + 1
326 ns_array(multi_ns_grid) = nsd
329 print *,
' NS_ARRAY ELEMENTS CANNOT EXCEED ',nsd
330 print *,
' CHANGING NS_ARRAY(',multi_ns_grid,
') to ',
335 multi_ns_grid = multi_ns_grid - 1
337 IF (ftol_array(1) .eq. zero)
THEN
338 ftol_array(1) = 1.e-8_dp
339 IF (multi_ns_grid .eq. 1) ftol_array(1) = ftol
340 DO igrid = 2, multi_ns_grid
341 ftol_array(igrid) = 1.e-8_dp * (1.e8_dp * ftol)**
342 & ( real(igrid-1,dp)/(multi_ns_grid-1) )
359 IF (nvacskip .LE. 0) nvacskip = nfp
361 proc0:
IF (lwrite)
THEN
363 & ns_array(multi_ns_grid),ntheta1,nzeta,mpol,ntor,nfp,
365 & gamma,spres_ped,phiedge,curtor,bcrit,lrfp
367 & gamma,spres_ped,phiedge,curtor,lrfp
369 100
FORMAT(/,
' COMPUTATION PARAMETERS: (u = theta, v = zeta)'/,
371 &
' ns nu nv mu mv',/,
372 & 5i7,//,
' CONFIGURATION PARAMETERS:',/,1x,39(
'-'),/,
373 &
' nfp gamma spres_ped phiedge(wb)'
375 &
' curtor(A) BCrit(T) lRFP',
376 & /,i7,1p,e11.3,2e15.3,2e14.3,l12/)
379 & /,i7,1p,e11.3,2e15.3,e14.3,l12/)
381 WRITE (nthreed,110) ncurr,niter_array(multi_ns_grid),
382 & ns_array(1),nstep,nvacskip,
383 & ftol_array(multi_ns_grid),tcon0,lasym,lforbal,lmove_axis,
384 & lconm1,mfilter_fbdy,nfilter_fbdy,lfull3d1out,
385 & max_main_iterations,lgiveup,fgiveup
386 110
FORMAT(
' RUN CONTROL PARAMETERS:',/,1x,23(
'-'),/,
387 &
' ncurr niter nsin nstep nvacskip ftol tcon0',
388 &
' lasym lforbal lmove_axis lconm1',/,
389 & 4i7,i10,1p,2e10.2,4l9,/,
390 &
' mfilter_fbdy nfilter_fbdy lfull3d1out max_main_iterations',
391 &
' lgiveup fgiveup',/,
392 & 2(6x,i7),l12,10x,i10,l8,e9.1,/)
394 WRITE (nthreed,120) precon_type, prec2d_threshold
395 120
FORMAT(
' PRECONDITIONER CONTROL PARAMETERS:',/,1x,34(
'-'),/,
396 &
' precon_type prec2d_threshold',/,2x,a10,1p,e20.2,/)
398 IF (nextcur .gt. 0)
THEN
399 WRITE(nthreed,
"(' EXTERNAL CURRENTS',/,1x,17('-'))")
401 IF (
ALLOCATED(curlabel))
THEN
402 ni = maxval(len_trim(curlabel(1:nextcur)))
405 WRITE (line,
'(a,i2.2,a)')
"(5a",ni,
")"
406 WRITE (line2,
'(a,i2.2,a)')
"(5(",ni-12,
"x,1p,e12.4))"
408 ni = min(i+4, nextcur)
409 IF (
ALLOCATED(curlabel))
THEN
410 WRITE (nthreed, line, iostat=mj)
411 & (trim(curlabel(n)),n=i,ni)
412 WRITE (nthreed, line2,iostat=mj) (extcur(n), n=i,ni)
418 IF (bloat .ne. one)
THEN
419 WRITE (nthreed,
'(" Profile Bloat Factor: ",1pe11.4)') bloat
420 phiedge = phiedge*bloat
423 IF (pres_scale .ne. one)
THEN
424 WRITE (nthreed,121) pres_scale
426 121
FORMAT(
' Pressure profile factor: ',1pe11.4,
427 &
' (multiplier for pressure)')
430 WRITE(nthreed,131) trim(pmass_type)
432 130
FORMAT(
' MASS PROFILE COEFFICIENTS - newton/m**2',
433 &
' (EXPANSION IN NORMALIZED RADIUS):')
434 131
FORMAT(
' PMASS parameterization type is ''', a,
'''')
435 132
FORMAT(1x,35(
'-'))
437 135
FORMAT(1p,6e12.3)
439 SELECT CASE(trim(pmass_type))
440 CASE (
'Akima_spline',
'cubic_spline')
441 WRITE(nthreed,
"(' am_aux_s is' )")
442 n = nonzerolen(am_aux_s,
SIZE(am_aux_s))
443 WRITE(nthreed,135)(am_aux_s(i),i=1, n)
444 n = nonzerolen(am_aux_f,
SIZE(am_aux_f))
445 WRITE(nthreed,
"(' am_aux_f is' )")
446 WRITE(nthreed,135)(am_aux_f(i),i=1, n)
448 n = nonzerolen(am,
SIZE(am))
449 WRITE(nthreed,135)(am(i-1),i=1,n)
460 WRITE(nthreed,143) trim(piota_type)
461 SELECT CASE(trim(piota_type))
462 CASE (
'Akima_spline',
'cubic_spline')
463 n = nonzerolen(ai_aux_s,
SIZE(ai_aux_s))
464 WRITE(nthreed,
"(' ai_aux_s is' )")
465 WRITE(nthreed,135)(ai_aux_s(i),i=1, n)
466 n = nonzerolen(ai_aux_f,
SIZE(ai_aux_f))
467 WRITE(nthreed,
"(' ai_aux_f is' )")
468 WRITE(nthreed,135)(ai_aux_f(i),i=1, n)
470 n = nonzerolen(ai,
SIZE(ai))
471 WRITE(nthreed,135)(ai(i-1),i=1, n)
476 WRITE(nthreed,146) trim(pcurr_type)
479 SELECT CASE(trim(pcurr_type))
480 CASE (
'Akima_spline_Ip',
'Akima_spline_I',
481 &
'cubic_spline_Ip',
'cubic_spline_I')
482 n = nonzerolen(ac_aux_s,
SIZE(ac_aux_s))
483 WRITE(nthreed,
"(' ac_aux_s is' )")
484 WRITE(nthreed,135)(ac_aux_s(i),i=1, n)
485 n = nonzerolen(ac_aux_f,
SIZE(ac_aux_f))
486 WRITE(nthreed,
"(' ac_aux_f is' )")
487 WRITE(nthreed,135)(ac_aux_f(i),i=1, n)
489 n = nonzerolen(ac,
SIZE(ac))
490 WRITE(nthreed,135)(ac(i-1),i=1, n)
494 140
FORMAT(/
' IOTA PROFILE COEFFICIENTS',
495 &
' (EXPANSION IN NORMALIZED RADIUS):',/,1x,35(
'-'))
496 142
FORMAT(/
' SAFETY-FACTOR (q) PROFILE COEFFICIENTS ai',
497 &
' (EXPANSION IN NORMALIZED RADIUS):',/,1x,35(
'-'))
498 143
FORMAT(
' PIOTA parameterization type is ''', a,
'''')
499 145
FORMAT(/
' TOROIDAL CURRENT DENSITY (*V'') COEFFICIENTS',
500 &
' ac (EXPANSION IN NORMALIZED RADIUS):')
501 146
FORMAT(
' PCURR parameterization type is ''', a,
'''')
502 147
FORMAT(1x,38(
'-'))
505 n = nonzerolen(aphi,
SIZE(aphi))
506 WRITE(nthreed,135)(aphi(i),i=1,n)
507 150
FORMAT(/
' NORMALIZED TOROIDAL FLUX COEFFICIENTS aphi',
508 &
' (EXPANSION IN S):',/,1x,35(
'-'))
510 IF (any(ah .ne. zero))
THEN
512 n = nonzerolen(ah,
SIZE(ah))
513 WRITE(nthreed,135)(ah(i-1),i=1, n)
515 n = nonzerolen(at,
SIZE(at))
516 WRITE(nthreed,135)(at(i-1),i=1, n)
519 160
FORMAT(
' HOT PARTICLE PRESSURE COEFFICIENTS ah',
520 &
' (EXPANSION IN TOROIDAL FLUX):',/,1x,35(
'-'))
521 165
FORMAT(
' HOT PARTICLE TPERP/T|| COEFFICIENTS at',
522 &
' (EXPANSION IN TOROIDAL FLUX):',/,1x,35(
'-'))
527 180
FORMAT(/,
' R-Z FOURIER BOUNDARY COEFFICIENTS AND',
528 &
' MAGNETIC AXIS INITIAL GUESS',/,
529 &
' R = RBC*COS(m*u - n*v) + RBS*SIN(m*u - n*v),',
530 &
' Z = ZBC*COS(m*u - n*v) + ZBS*SIN(m*u-n*v)'/1x,86(
'-'),
531 & /,
' nb mb rbc rbs zbc ',
533 &
' raxis(c) raxis(s) zaxis(c) zaxis(s)')
543 delta = atan((rbs(0,1) - zbc(0,1))/
544 & (abs(rbc(0,1)) + abs(zbs(0,1))))
545 IF (delta .ne. zero)
THEN
548 trc = rbc(n,m)*cos(m*delta) + rbs(n,m)*sin(m*delta)
549 rbs(n,m) = rbs(n,m)*cos(m*delta)
550 & - rbc(n,m)*sin(m*delta)
552 tzc = zbc(n,m)*cos(m*delta) + zbs(n,m)*sin(m*delta)
553 zbs(n,m) = zbs(n,m)*cos(m*delta)
554 & - zbc(n,m)*sin(m*delta)
579 rbcc => rmn_bdy(:,:,rcc)
580 zbsc => zmn_bdy(:,:,zsc)
582 rbss => rmn_bdy(:,:,rss)
583 zbcs => zmn_bdy(:,:,zcs)
587 rbsc => rmn_bdy(:,:,rsc)
588 zbcc => zmn_bdy(:,:,zcc)
590 rbcs => rmn_bdy(:,:,rcs)
591 zbss => zmn_bdy(:,:,zss)
595 rmn_bdy = 0; zmn_bdy = 0
597 ioff = lbound(rbcc,1)
598 joff = lbound(rbcc,2)
603 & (mfilter_fbdy.gt.1 .and. m.gt.mfilter_fbdy)) cycle
606 & (nfilter_fbdy.gt.0 .and. abs(n).gt.nfilter_fbdy)) cycle
610 ELSE IF (n .gt. 0)
THEN
615 rbcc(ni,mj) = rbcc(ni,mj) + rbc(n,m)
616 IF (m .gt. 0) zbsc(ni,mj) = zbsc(ni,mj) + zbs(n,m)
619 IF (m .gt. 0) rbss(ni,mj) = rbss(ni,mj) + isgn*rbc(n,m)
620 zbcs(ni,mj) = zbcs(ni,mj) - isgn*zbs(n,m)
624 IF (m .gt. 0) rbsc(ni,mj) = rbsc(ni,mj) + rbs(n,m)
625 zbcc(ni,mj) = zbcc(ni,mj) + zbc(n,m)
627 rbcs(ni,mj) = rbcs(ni,mj) - isgn*rbs(n,m)
628 IF (m .gt. 0) zbss(ni,mj) = zbss(ni,mj) + isgn*zbc(n,m)
632 IF (ier_flag_init .ne. norm_term_flag) cycle
633 trc = abs(rbc(n,m)) + abs(rbs(n,m))
634 & + abs(zbc(n,m)) + abs(zbs(n,m))
637 IF (trc.eq.zero .and. abs(raxis_cc(n)).eq.zero .and.
638 & abs(zaxis_cs(n)).eq.zero) cycle
639 IF (lwrite)
WRITE (nthreed,195) n, m, rbc(n,m), rbs(n,m),
640 & zbc(n,m), zbs(n,m), raxis_cc(n), raxis_cs(n),
641 & zaxis_cc(n), zaxis_cs(n)
643 IF (trc .eq. zero) cycle
644 IF (lwrite)
WRITE (nthreed,195) n, m, rbc(n,m), rbs(n,m),
649 195
FORMAT(i5,i4,1p,8e12.4)
656 rtest = sum(rbcc(1:ntor1,mj))
657 ztest = sum(zbsc(1:ntor1,mj))
658 lflip=(rtest*ztest .lt. zero)
660 IF (lflip)
CALL flip_theta(rmn_bdy, zmn_bdy)
669 IF (lconm1 .and. (lthreed .or. lasym))
THEN
670 ALLOCATE (temp(
SIZE(rbcc,1)))
674 rbss(:,mj) = p5*(temp(:) + zbcs(:,mj))
675 zbcs(:,mj) = p5*(temp(:) - zbcs(:,mj))
680 rbsc(:,mj) = p5*(temp(:) + zbcc(:,mj))
681 zbcc(:,mj) = p5*(temp(:) - zbcc(:,mj))
683 IF (
ALLOCATED(temp))
DEALLOCATE (temp)
689 precon_type = trim(adjustl(precon_type))
692 ch1 = precon_type(1:1); ch2 = precon_type(2:2)
699 IF (ch2 ==
'g' .or. ch2 ==
'G') itype_precon = 1
703 IF (ch2 ==
'm' .or. ch2 ==
'M') itype_precon = 2
704 IF (len_trim(precon_type) == 6) itype_precon = 3
708 IF (ch2 ==
'f' .or. ch2 ==
'F') itype_precon = 4
717 CALL second0(treadoff)
718 timer(tread) = timer(tread) + (treadoff-treadon)
719 CALL mpi_bcast(lprecond,1,mpi_logical,0,runvmec_comm_world,
721 readin_time = timer(tread)
723 END SUBROUTINE readin
725 INTEGER FUNCTION nonzerolen(array, n)
726 USE vmec_main,
ONLY: dp, zero
728 INTEGER,
INTENT(IN) :: n
729 REAL(dp),
INTENT(IN) :: array(n)
733 IF (array(k) .NE. zero)
EXIT
738 END FUNCTION nonzerolen