4 USE descriptor_mod,
ONLY: iam
9 USE nscalingtools,
ONLY: startglobrow, endglobrow
17 INTEGER :: irad_scale = 11
18 INTEGER,
PUBLIC :: niter_max
19 REAL(dp),
ALLOCATABLE :: jdotb(:,:,:)
20 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:) :: bmnc_res
21 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:) :: bmns_res
23 PUBLIC :: init_data, add_perturb
50 CHARACTER(LEN=10) :: date0
51 CHARACTER(LEN=10) :: time0
52 CHARACTER(LEN=10) :: zone0
53 CHARACTER(LEN=256) :: temp
54 CHARACTER(LEN=100) :: command_arg
55 CHARACTER(LEN=256) :: filename
60 CALL getcarg(1, command_arg, numargs)
62 IF (numargs .GT. 0)
THEN
63 IF (len_trim(temp) .NE. 0)
THEN
64 IF (command_arg(1:7) .NE.
'siesta_' .AND.
65 command_arg(1:7) .NE.
'SIESTA_')
THEN
66 temp =
'siesta_' // trim(command_arg)
68 istat = len_trim(temp)
69 IF (temp(istat - 3:istat) .NE.
'.jcf' .AND.
70 temp(istat - 3:istat) .NE.
'.JCF')
THEN
71 temp = trim(temp) //
'.jcf'
81 irad_scale = min(irad_scale,
nsin)
93 &
' is not a valid wout file name.')
101 OPEN (unit=
unit_out, file=filename, iostat=istat)
103 CALL assert_eq(istat, 0,
'ERROR WRITING SIESTA OUTPUT FILE')
104 IF (iam .NE. 0)
RETURN
107 WRITE (*,1001)
'4.0', 100917
113 WRITE (*,1003) trim(temp)
117 CALL date_and_time(date0,time0,zone0)
118 READ (date0(5:6),2000) imon
119 WRITE (
unit_out,1004) months(imon), date0(7:8), date0(1:4),
120 time0(1:2), time0(3:4), time0(5:6)
132 1000
FORMAT(
'output_',a,
'_',i4.4,2(
'X',i3.3),
'.txt')
133 1001
FORMAT(62(
'-'),/,
134 1x,
'SIESTA MHD EQUILIBRIUM CODE v',a,
' (',i6,
')', /,1x,
135 'Scalable Island Equilibrium Solver for Toroidal Applications'
137 1002
FORMAT(
'CASE: ',a)
138 1003
FORMAT(
'OUTPUT FILE (SCREEN DUMP): output_',a,
'.txt')
139 1004
FORMAT(
' DATE = ',a3,
' ',a2,
',',a4,
' ',
' TIME = ',2(a2,
':'),a2,/)
140 1005
FORMAT(i2,
' mres: ',i4,
' HelPert: ',1pe9.2,
' HelPertA: ',1pe9.2)
141 1006
FORMAT(/,
' ngmres_type: ', i4,
' iOrtho: ', i4,
' lColScale: ', l2)
147 SUBROUTINE add_perturb (xc, getwmhd)
157 REAL(dp),
EXTERNAL :: getwmhd
161 INTEGER,
PARAMETER :: ipert_sign=2
162 INTEGER :: js, l, iprint, mres0, nres0, isign1, icount,
164 REAL(dp) :: w0, w1, wmhd, eta_factor_save, p_width_min
165 REAL(dp) :: normal, HelPert0, HelPert0A, HP, rad, &
167 LOGICAL :: lresist_save
175 IF (ntor.EQ.0 .OR. (all(
helpert.EQ.zero)
176 .AND. all(
helperta.EQ.zero)))
RETURN
179 nsmin = max(1,startglobrow - 1)
180 nsmax = min(ns,endglobrow + 1)
182 ALLOCATE (
buv_res(ntheta,nzeta,nsmin:nsmax),
183 jdotb(ntheta,nzeta,nsmin:nsmax), stat=istat)
184 CALL assert(istat.eq.0,
'ALLOCATION ERROR IN add_perturbation')
195 IF (.NOT.lverbose .AND. iprint.EQ.6) cycle
196 WRITE (iprint,
'(/,a,/,2a)')
197 ' Adding helical magnetic field perturbations',
198 ' 10^6 X Del-W mres nres HelPert',
199 ' rad |m*chip+n*phip| iota radial width'
203 #if defined(JDOTB_PERT)
205 jdotb = ksupsijf0(:,:,nsmin:nsmax)*bsubsijf(:,:,nsmin:nsmax)
206 + ksupuijf0(:,:,nsmin:nsmax)*bsubuijf(:,:,nsmin:nsmax)
207 + ksupvijf0(:,:,nsmin:nsmax)*bsubvijf(:,:,nsmin:nsmax)
209 ALLOCATE (bmnc_res(0:mpol,-ntor:ntor, nsmin:nsmax),
210 bmns_res(0:mpol,-ntor:ntor, nsmin:nsmax), stat=istat)
211 CALL assert(istat.eq.0,
'ISTAT != 0 IN ADD_PERT')
224 res_test:
DO icount = 1,
SIZE(
helpert)
225 helpert0 = abs(
helpert(icount)/normal)
226 helpert0a= abs(
helperta(icount)/normal)
227 mres0 = abs(
mres(icount))
228 IF ((helpert0 .EQ. zero .AND.
229 helpert0a .EQ. zero) .OR. mres0.GT.mpol) cycle
233 nres_test:
DO nres0 = -ntor, ntor
235 jsurf_test:
DO WHILE (jstart .LT. ns-1)
237 IF (.NOT.findresonance(mres0, nres0, rad, jstart))
EXIT
242 pert_sign:
DO isign1 = 1, 2
249 scale_test:
DO irscale = 1, irad_scale
250 CALL getrespert(irscale, isign1, mres0, nres0, rad,
251 helpert0, helpert0a, chip0, phip0, p_width
254 IF (wmhd .LT. w1)
THEN
258 p_width_min = p_width
271 IF (isign1 .EQ. 2) hp = -hp
274 IF (.NOT.lverbose .AND. iprint.EQ.6) cycle
275 WRITE (iprint, 100) 10**6*(wmhd-w0), mres0, nres0,
276 hp*normal, rad, abs(mres0*chip0+nres0*nfp*phip0
277 chip0/phip0, p_width_min
279 100
FORMAT(1p,e12.3,2(i8),3x,1pe10.2,0p,f8.2,f11.2,4x,2(f11.2))
283 CALL getrespert(irscale, isign1, mres0, nres0, rad,
284 helpert0, helpert0a, chip0, phip0, p_width
287 IF (abs(wmhd-w1) .GT. 1.e-12_dp)
THEN
288 CALL siesta_error_set_error(siesta_error_general,
292 CALL update_state(.false., zero, zero)
302 IF (.NOT.lverbose .AND. iprint.EQ.6) cycle
307 DEALLOCATE (
buv_res, jdotb, bmnc_res, bmns_res, stat=istat)
311 END SUBROUTINE add_perturb
314 LOGICAL FUNCTION findresonance(mres0, nres0, resrad, jstart)
318 INTEGER,
INTENT(in) :: mres0
319 INTEGER,
INTENT(in) :: nres0
320 INTEGER,
INTENT(inout) :: jstart
321 REAL(dp),
INTENT(out) :: resrad
326 REAL(dp) :: del1, del2, rho1, rho2, delmin
334 IF (mres0*nres0 .EQ. 0)
THEN
335 findresonance = .false.
340 DO js = jstart, ns - 1
341 del1 = mres0*chipf_i(js) + nres0*nfp*phipf_i(js)
342 del1 = del1/max(abs(mres0*chipf_i(js)),abs(nres0*nfp*phipf_i(js
343 IF (delmin .EQ. -1 .OR. abs(del1) .LT. delmin)
THEN
344 IF (delmin .EQ. -1) del2 = del1
349 IF (del1*del2 < zero)
THEN
353 resrad = (rho2*abs(del1) + rho1*abs(del2))/(abs(del1) + abs(del2
361 IF (delmin .LT. 1.e-3_dp)
THEN
363 resrad = hs_i*(jmin-1)
368 findresonance=.false.
371 END FUNCTION findresonance
374 SUBROUTINE getrespert(irscale, isign1, mres0, nres0, rad, &
375 HelPert0, HelPert0A, chip0, phip0, p_width)
385 INTEGER,
INTENT(in) :: irscale, isign1
386 INTEGER,
INTENT(in) :: mres0, nres0
387 REAL(dp),
INTENT(in) :: rad
388 REAL(dp),
INTENT(out) :: chip0, phip0
389 REAL(dp),
INTENT(in) :: HelPert0, HelPert0A
390 REAL(dp),
INTENT(out) :: p_width
395 REAL(dp) :: rho, rhores, pert_prof(ns), HelP_local,
397 REAL(dp),
ALLOCATABLE :: bmn_res(:,:)
398 REAL(dp),
ALLOCATABLE :: resc(:,:,:), ress(:,:,:)
402 help_local = helpert0
403 helpa_local = helpert0a
404 IF (isign1 .EQ. 2)
THEN
405 help_local = -helpert0
406 helpa_local = -helpert0a
410 js = int(rad/hs_i) + 1
411 IF (js.LT.1 .OR. js.GT.ns)
RETURN
414 locrad = real(irscale - 1,dp)/(irad_scale - 1)
415 p_width = 0.1_dp + 0.9_dp*locrad
420 rhores = locrad*rhores
421 pert_prof(js) = one/(one + rhores*rhores)
426 #if defined(JDOTB_PERT)
427 ALLOCATE(resc(0:mpol,-ntor:ntor, nsmin:nsmax),
428 ress(0:mpol,-ntor:ntor, nsmin:nsmax), stat=istat)
430 rho = max(maxval(abs(bmnc_res(mres0,nres0,:))),
431 maxval(abs(bmns_res(mres0,nres0,:))))
435 resc(mres0,nres0,:) = help_local*bmnc_res(mres0,nres0,:)
440 ress(mres0,nres0,:) = helpa_local*bmns_res(mres0,nres0,:)
451 DEALLOCATE (resc, ress, stat=istat)
453 ALLOCATE (bmn_res(0:mpol,-ntor:ntor), stat=istat)
454 CALL assert(istat.eq.0,
'ISTAT != 0 IN GETRESPERT')
457 bmn_res(mres0,nres0) = help_local
461 bmn_res(mres0,nres0) = helpa_local
464 DO js = nsmin + 1, nsmax
470 DEALLOCATE (bmn_res, stat=istat)
473 END SUBROUTINE getrespert
475 SUBROUTINE compute_resonant_epar(xc, epar)
479 REAL(dp),
INTENT(out) :: epar(ns,ntheta,nzeta)
480 REAL(dp),
INTENT(in) :: xc(ns,0:mpol,-ntor:ntor,3)
484 REAL(dp),
ALLOCATABLE :: eparmn(:,:,:)
486 INTEGER :: mres, nres, istat, jstart
493 ALLOCATE (eparmn(ns,0:mpol,-ntor:ntor), stat=istat)
494 CALL assert(istat.eq.0,
'Allocation error in Compute_Resonant_Epar'
501 DO WHILE (jstart .lt. ns-1)
502 IF (.NOT.findresonance(mres, nres, rad, jstart))
EXIT
503 CALL assert(.false.,
'AddResonantE not implemented')
511 END SUBROUTINE compute_resonant_epar
514 END MODULE perturbation