-
Notifications
You must be signed in to change notification settings - Fork 3
/
xspech.h
793 lines (559 loc) · 35.2 KB
/
xspech.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!title (main) ! Main program.
!latex \briefly{main program}
!latex \calledby{}
!latex \calls{\link{global}:readin,
!latex \link{global}:wrtend,
!latex \link{preset},
!latex \link{packxi},
!latex \link{volume},
!latex \link{pc00aa},
!latex \link{newton},
!latex \link{dforce},
!latex \link{hesian},
!latex \link{ra00aa},
!latex \link{bnorml},
!latex \link{sc00aa},
!latex \link{jo00aa} and
!latex \link{pp00aa}}
!latex \tableofcontents
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
! uppercase generally indicates macros;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
program xspech
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
use constants, only : zero, one, pi2
use numerical, only : vsmall, logtolerance
use fileunits, only : ounit, zunit, lunit
use inputlist, only : Wmacros, Wxspech, ext, &
Nfp, Igeometry, Nvol, Lrad, &
tflux, pflux, phiedge, pressure, pscale, helicity, Ladiabatic, adiabatic, gamma, &
Rbc, Zbs, Rbs, Zbc, &
Lconstraint, &
Lfreebound, mfreeits, gBntol, gBnbld, &
Lfindzero, &
odetol, nPpts, nPtrj, &
LHevalues, LHevectors, LHmatrix, Lperturbed, Lcheck, &
Lzerovac
use cputiming, only : Txspech
use allglobal, only : readin, wrtend, ncpu, myid, cpus, &
Mvol, &
YESstellsym, NOTstellsym, &
Iquad, &
mn, im, in, &
Ntz, &
LGdof, NGdof, &
iRbc, iZbs, iRbs, iZbc, &
BBe, IIo, BBo, IIe, &
Btemn, Bzemn, Btomn, Bzomn, &
vvolume, &
Lcoordinatesingularity, Lplasmaregion, Lvacuumregion, &
dtflux, dpflux, &
ImagneticOK, &
ForceErr, &
efmn, ofmn, &
iBns, iBnc, iVns, iVnc, &
Ate, Aze, Ato, Azo, & ! only required for debugging; 09 Mar 17;
nfreeboundaryiterations, &
beltramierror
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
LOCALS
LOGICAL :: LComputeDerivatives, LContinueFreeboundaryIterations, exist, LupdateBn
! INTEGER :: nfreeboundaryiterations, imn, lmn, lNfp, lim, lin, ii ! 09 Mar 17;
INTEGER :: imn, lmn, lNfp, lim, lin, ii, ideriv
INTEGER :: vvol, llmodnp, ifail, wflag, iflag, vflag
REAL :: rflag, lastcpu, bnserr, lRwc, lRws, lZwc, lZws, lItor, lGpol, lgBc, lgBs
REAL, allocatable :: position(:), gradient(:)
CHARACTER :: pack
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
MPISTART ! It might be easy to read the source if this macro was removed; SRH: 27 Feb 18;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
if( myid.eq.0 ) then
write(ounit,'("xspech : ", 10x ," : ")')
! COMPILATION ! do not delete; this line is replaced (see Makefile) with a write statement identifying date, time, compilation flags, etc.;
endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
cpus = GETTIME ! set initial time; 04 Dec 14;
cpuo = cpus
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
if( myid.eq.0 ) then
write(ounit,'("xspech : ", 10x ," : ")')
write(ounit,'("xspech : ",f10.2," : begin execution ; ncpu=",i3," ; calling global:readin ;")') cpus-cpus, ncpu
endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!latex \subsection{reading input, allocating global variables}
!latex \begin{enumerate}
!latex \item The input namelists and geometry are read in via a call to \link{global}\verb+:readin+.
!latex A full description of the required input is given in \link{global}.
!latex \item Most internal variables, global memory etc., are allocated in \link{preset}.
!latex \end{enumerate}
WCALL( xspech, readin ) ! sets Rscale; 03 Nov 16;
WCALL( xspech, preset )
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!latex \subsection{preparing output file \type{.ext.sp.its}}
!latex \begin{enumerate}
!latex \item The file \verb+.ext.sp.its+ is created.
!latex This file contains the interface geometry at each iteration, which is useful for constructing movies illustrating the convergence.
!latex The format of this file is:
!latex \begin{verbatim} open(zunit,file="."//trim(ext)//".sp.its",status="unknown",form="unformatted")
!latex write(zunit) mn, Mvol, Nfp ! integer;
!latex write(zunit) im(1:mn) ! integer;
!latex write(zunit) in(1:mn) ! integer;
!latex ! begin do loop over iterations;
!latex write(zunit) wflag, iflag, Energy, rflag ! written in global.h:wrtend; perhaps redundant;
!latex write(zunit) iRbc(1:mn,0:Mvol) ! real;
!latex write(zunit) iZbs(1:mn,0:Mvol) ! real;
!latex write(zunit) iRbs(1:mn,0:Mvol) ! real;
!latex write(zunit) iZbc(1:mn,0:Mvol) ! real;
!latex ! end do loop over iterations; \end{verbatim}
!latex \end{enumerate}
if( myid.eq.0 ) then ! prepare ``convergence evolution'' output; save restart file;
open(zunit, file = "."//trim(ext)//".sp.its", status = "unknown", form = "unformatted" ) ! this file is written to in globals/wrtend; 11 Aug 14;
write(zunit) mn, Mvol, Nfp
write(zunit) im(1:mn)
write(zunit) in(1:mn)
wflag = 0 ; iflag = 0 ; rflag = zero
WCALL( xspech, wrtend, ( wflag, iflag, rflag ) ) ! write restart file etc. ! 14 Jan 13;
endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
FATAL( xspech, NGdof.lt.0, counting error )
SALLOCATE( position, (0:NGdof), zero ) ! position ; NGdof = #geometrical degrees-of-freedom was computed in preset;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
nfreeboundaryiterations = -1
9000 nfreeboundaryiterations = nfreeboundaryiterations + 1 ! this is the free-boundary iteration loop; 08 Jun 16;
if (nfreeboundaryiterations .eq. 0) then ! first iteration only run fixed-boundary
Mvol = Nvol
else
Mvol = Nvol + Lfreebound
endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!latex \subsection{``packing" geometrical degrees-of-freedom into vector}
!latex \begin{enumerate}
!latex \item If \internal{NGdof.gt.0}, where \internal{NGdof}$\equiv $ counts the geometrical degrees-of-freedom, i.e. the $R_{bc}$, $Z_{bs}$, etc.,
!latex then \link{packxi} is called to ``pack" the geometrical degrees-of-freedom into \internal{position(0:NGdof)}.
!latex \end{enumerate}
if( NGdof.gt.0 ) then ! pack geometry into vector; 14 Jan 13;
pack = 'P'
WCALL( xspech, packxi, ( NGdof, position(0:NGdof), Mvol, mn, iRbc(1:mn,0:Mvol), iZbs(1:mn,0:Mvol), iRbs(1:mn,0:Mvol), iZbc(1:mn,0:Mvol), pack ) )
endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!latex \subsection{initialize adiabatic constants}
!latex \begin{enumerate}
!latex \item If \inputvar{Ladiabatic.eq.0}, then the ``adiabatic constants'' in each region, $P_v$, are calculated as
!latex \be P_v \equiv p_v V_v^\gamma, \label{eq:adiabatic}
!latex \ee
!latex where $p_v\equiv$ \inputvar{pressure(vvol)}, the volume $V_v$ of each region is computed by \link{volume},
!latex and the adiabatic index $\gamma\equiv$ \inputvar{gamma}.
!latex \end{enumerate}
do vvol = 1, Mvol
vflag = 0
WCALL( xspech, volume, ( vvol, vflag ) ) ! compute volume;
if( Ladiabatic.eq.0 ) adiabatic(vvol) = pressure(vvol) * vvolume(vvol)**gamma ! initialize adiabatic constants using supplied pressure profile;
enddo ! end of do vvol = 1, Mvol;
if( Mvol.gt.Nvol ) then ; adiabatic(Mvol) = zero ; pressure(Mvol) = zero ! these are never used; 15 May 13;
endif
if( Wxspech .and. myid.eq.0 ) then
cput = GETTIME
write(ounit,'("xspech : ",f10.2," : myid=",i3," ; adiabatic constants = "999es13.5)') cput-cpus, myid, adiabatic(1:Mvol)
endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!latex \subsection{solving force-balance}
!latex \begin{enumerate}
!latex \item If there are geometrical degress of freedom, i.e. if \internal{NGdof.gt.0}, then \bi
!l tex \item[] if \inputvar{Lminimize.eq.1}, call \link{pc00aa} to find minimum of energy functional
!l tex using quasi-Newton, preconditioned conjugate gradient method, \nag{www.nag.co.uk/numeric/FL/manual19/pdf/E04/e04dgf_fl18.pdf}{E04DGF};
!latex \item[] If \inputvar{Lfindzero.gt.0}, call \link{newton} to find extremum of constrained energy functional using a Newton method,
!latex \nag{www.nag.co.uk/numeric/FL/manual19/pdf/C05/c05pdf_fl19.pdf}{C05PDF};
!latex \ei
!latex \end{enumerate}
! if( Lminimize.eq.1 ) then
!
! ifail = 1 ! this is probably not required; 26 Feb 13;
!
! WCALL(xspech,pc00aa,( NGdof, position(1:NGdof), Mvol, mn, ifail ))
!
! endif
if( NGdof.gt.0 ) then
if( Lfindzero.gt.0 ) then
ifail = 1
WCALL( xspech, newton, ( NGdof, position(0:NGdof), ifail ) )
endif
pack = 'U' ! unpack geometrical degrees of freedom; 13 Sep 13;
WCALL( xspech, packxi, ( NGdof, position(0:NGdof), Mvol, mn, iRbc(1:mn,0:Mvol), iZbs(1:mn,0:Mvol), iRbs(1:mn,0:Mvol), iZbc(1:mn,0:Mvol), pack ) )
endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!latex \subsection{post diagnostics}
!latex \begin{enumerate}
!latex \item The pressure is computed from the adiabatic constants from \Eqn{adiabatic}, i.e. $p=P/V^\gamma$.
!latex \item The Beltrami/vacuum fields in each region are re-calculated using \link{dforce}.
!latex \item If \inputvar{Lcheck.eq.5 .or. LHevalues .or. LHevectors .or. Lperturbed.eq.1}, then the ``Hessian'' matrix is examined using \link{hesian}.
!latex \end{enumerate}
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
! set/reset input variables;
! if( Lconstraint.lt.2 ) helicity(1:Nvol) = lABintegral(1:Nvol) ! updated ``input'' quantity;
#ifdef DEBUG
do vvol = 1, Mvol
FATAL( xspech, vvolume(vvol).lt.vsmall, error dividing adiabatic by volume )
enddo
#endif
pressure(1:Mvol) = adiabatic(1:Mvol) / vvolume(1:Mvol)**gamma ! this matches construction of adiabatic above;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
SALLOCATE( gradient, (0:NGdof), zero )
lastcpu = GETTIME
LComputeDerivatives = .false.
! vvol = Mvol ; ideriv = 0 ; ii = 1
! write(ounit,'("xspech : ", 10x ," : sum(Ate(",i3,",",i2,",",i2,")%s) =",99es23.15)') vvol, ideriv, ii, sum(Ate(vvol,ideriv,ii)%s(0:Lrad(vvol)))
WCALL( xspech, dforce, ( NGdof, position(0:NGdof), gradient(0:NGdof), LComputeDerivatives ) ) ! (re-)calculate Beltrami fields;
DALLOCATE(gradient)
#ifdef DEBUG
do vvol = 1, Mvol-1
; FATAL( xspech, BBe(vvol).lt.logtolerance, underflow )
if( Igeometry.eq.3 ) then ! include spectral constraints; 04 Dec 14;
;FATAL( xspech, IIo(vvol).lt.logtolerance, underflow )
endif
if( NOTstellsym ) then
;FATAL( xspech, BBo(vvol).lt.logtolerance, underflow )
if( Igeometry.eq.3 ) then ! include spectral constraints; 04 Dec 14;
FATAL( xspech, IIe(vvol).lt.logtolerance, underflow )
endif
endif
enddo
#endif
if( myid.eq.0 ) then
cput = GETTIME
write(ounit,1000) cput-cpus, nfreeboundaryiterations, ForceErr, cput-lastcpu, "|BB|e", alog10(BBe(1:min(Mvol-1,28)))
if( Igeometry.ge.3 ) then ! include spectral constraints; 04 Dec 14;
write(ounit,1001) "|II|o", alog10(IIo(1:min(Mvol-1,28)))
endif
if( NOTstellsym ) then
write(ounit,1001) "|BB|o", alog10(BBo(1:min(Mvol-1,28)))
if( Igeometry.ge.3 ) then ! include spectral constraints; 04 Dec 14;
write(ounit,1001) "|II|e", alog10(IIe(1:min(Mvol-1,28)))
endif
endif
endif
1000 format("xspech : ",f10.2," : #freeits=",i3," ; ":"|f|="es12.5" ; ":"time=",f10.2,"s ;":" log"a5,:"="28f6.2" ...")
1001 format("xspech : ", 10x ," : ",3x," ; ":" " 12x " ":" ", 10x ," ;":" log"a5,:"="28f6.2" ...")
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
if( Lcheck.eq.5 .or. LHevalues .or. LHevectors .or. LHmatrix .or. Lperturbed.eq.1 ) then ! check construction of Hessian; 01 Jul 14;
if( myid.eq.0 ) then
cput = GETTIME
write(ounit,'("xspech : ", 10x ," : ")')
write(ounit,'("xspech : ",f10.2," : myid=",i3," ; calling hesian ; see .ext.hessian.myid ;")') cput-cpus, myid
endif
WCALL( xspech, hesian, ( NGdof, position(0:NGdof), Mvol, mn, LGdof ) )
endif ! end of if( Lcheck.eq.5 ) ; 01 Jul 14;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
select case( Igeometry ) ! 08 Feb 16;
case( 1 ) ; tflux(1) = dtflux(1) ; pflux(1) = dpflux(1) ! 08 Feb 16;
case( 2:3 ) ; tflux(1) = dtflux(1) ; pflux(1) = zero ! 08 Feb 16;
end select ! 08 Feb 16;
do vvol = 2, Mvol ; tflux(vvol) = tflux(vvol-1) + dtflux(vvol) ! 01 Jul 14;
; pflux(vvol) = pflux(vvol-1) + dpflux(vvol) ! 01 Jul 14;
enddo
tflux(1:Mvol) = tflux(1:Mvol) * pi2 / phiedge ! this is the "inverse" operation defined in preset; 19 Jul 16;
pflux(1:Mvol) = pflux(1:Mvol) * pi2 / phiedge
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!latex \subsection{free-boundary: re-computing normal field}
!latex \begin{enumerate}
!latex \item If \inputvar{Lfreebound.eq.1} and \inputvar{Lfindzero.gt.0} and\inputvar{mfreeits.ne.0},
!latex then the magnetic field at the computational boundary produced by the plasma currents is computed using \link{bnorml}.
!latex \item The ``new'' magnetic field at the computational boundary produced by the plasma currents is updated using a Picard scheme:
!latex \be \verb+Bns+_i^{j} = \lambda \, \verb+Bns+_i^{j-1} + (1-\lambda) \verb+Bns+_i, \label{eq:blending}
!latex \ee
!latex where $j$ labels free-boundary iterations, the ``blending parameter'' is $\lambda\equiv $ \inputvar{gBnbld},
!latex and \verb+Bns+$_i$ is computed by virtual casing.
!latex The subscript ``$i$'' labels Fourier harmonics.
!latex \item If the new (unblended) normal field is {\em not} sufficiently close to the old normal field, as quantified by \inputvar{gBntol},
!latex then the free-boundary iterations continue.
!latex This is quantified by
!latex \be \sum_i | \verb+Bns+_i^{j-1} - \verb+Bns+_i | / N, \label{eq:gBntol}
!latex \ee
!latex where $N$ is the total number of Fourier harmonics.
!latex \item There are several choices that are available:
!latex \begin{enumerate}
!latex \item if \inputvar{mfreeits} $= -2$ : the vacuum magnetic field
!latex (really, the normal component of the field produced by the external currents at the computational boundary)
!latex required to hold the given equlibrium is written to file.
!latex This information is required as input by
!latex \paper{FOCUS}{Caoxiang Zhu, Stuart Hudson et al.}{10.1088/1741-4326/aa8e0a}{Nucl. Fusion}{58}{016008}{2017}
!latex for example. (This option probably needs to revised.)
!latex \item if \inputvar{mfreeits} $= -1$ : after the plasma field is computed by virtual casing,
!latex the vacuum magnetic field is set to exactly balance the plasma field
!latex (again, we are really talking about the normal component at the computational boundary.)
!latex This will ensure that the computational boundary itself if a flux surface of the total magnetic field.
!latex \item if \inputvar{mfreeits} $= 0$ : the plasma field at the computational boundary is not updated; no ``free-boundary'' iterations take place.
!latex \item if \inputvar{mfreeits} $> 0$ : the plasma field at the computational boundary is updated according to the above blending \Eqn{blending},
!latex and the free-boundary iterations will continue until either the tolerance condition is met (see \inputvar{gBntol} and \Eqn{gBntol})
!latex or the maximum number of free-boundary iterations, namely \inputvar{mfreeits}, is reached.
!latex For this case, \inputvar{Lzerovac} is relevant:
!latex if \inputvar{Lzerovac} $= 1$, then the vacuum field is set equal to the normal field at every iteration,
!latex which results in the computational boundary being a flux surface.
!latex (I am not sure if this is identical to setting \inputvar{mfreeits}$= -1$; the logic etc. needs to be revised.)
!latex \end{enumerate}
!latex \end{enumerate}
LContinueFreeboundaryIterations = .false.
; LupdateBn = .false. ! default;
! if( Lfreebound.eq.1 .and. Lfindzero.gt.0 ) then
if( Lfreebound.eq.1) then ! removed Lfindzero check; Loizu Dec 18;
if( mfreeits.gt.0 .and. nfreeboundaryiterations.lt.mfreeits ) LupdateBn = .true.
if( mfreeits.lt.0 ) LupdateBn = .true.
endif
if( LupdateBn ) then
Mvol = Nvol + Lfreebound
lastcpu = GETTIME
WCALL( xspech, bnorml, ( mn, Ntz, efmn(1:mn), ofmn(1:mn) ) ) ! compute normal field etc. on computational boundary;
FATAL( xspech, mn-1.le.0, divide by zero )
if( YESstellsym ) bnserr = sum( abs( iBns(2:mn) - ofmn(2:mn) ) ) / (mn-1)
if( NOTstellsym ) bnserr = sum( abs( iBns(2:mn) - ofmn(2:mn) ) ) / (mn-1) &
+ sum( abs( iBnc(1:mn) - efmn(1:mn) ) ) / (mn )
if( bnserr.gt.gBntol ) then
LContinueFreeboundaryIterations = .true.
select case( mfreeits )
case( -2 ) ! mfreeits = -2 ; shall set plasma normal field at computational boundary ; 24 Nov 16;
inquire( file=trim(ext)//".Vn", exist=exist )
FATAL( xspech, .not.exist, ext.Vn does not exist : cannot set vacuum field)
if( myid.eq.0 ) then ! only myid = 0 reads in the vacuum field; 04 Jan 17;
; iBns(2:mn) = iVns(2:mn) - iBns(2:mn) ! temporary storage of the total field; 07 Dec 16;
if( NOTstellsym) iBnc(1:mn) = iVnc(1:mn) - iBnc(1:mn)
open(lunit, file = trim(ext)//".Vn", status="old" )
read(lunit,*) lmn, lNfp
read(lunit,*) ( lim, lin, lRwc, lRws, lZwc, lZws, imn = 1, lmn ) ! control surface; should confirm agreement; 04 Jan 17;
read(lunit,*) lmn, lNfp, lItor, lGpol
do imn = 1, lmn
read(lunit,*) lim, lin, lgBc, lgBs
do ii = 1, mn
if( lim.eq.im(ii) .and. lin*lNfp.eq.in(ii) ) then
; iVns(ii) = lgBs
if( NOTstellsym ) iVnc(ii) = lgBc
endif
enddo
enddo
close(lunit)
endif ! end of if( myid.eq.0 ) ; 07 Dec 16;
;RlBCAST( iVns(1:mn), mn, 0 ) ! only required for ii > 1 ;
if( NOTstellsym ) then
RlBCAST( iVnc(1:mn), mn, 0 )
endif
; iBns(2:mn) = - iBns(2:mn) - iVns(2:mn) ! updated vacuum field ; 24 Nov 16;
if( NOTstellsym) iBnc(1:mn) = - iBnc(1:mn) - iVnc(1:mn)
if( myid.eq.0 ) then
write(ounit,'("xspech : " 10x " : oBns=[",999(es11.03,","))') iBns(1:mn) ! 17 Jan 17;
write(ounit,'("xspech : " 10x " : nBns=[",999(es11.03,","))') ofmn(1:mn) ! 17 Jan 17;
if( NOTstellsym ) then
write(ounit,'("xspech : " 10x " : oBnc=[",999(es11.03,","))') iBnc(1:mn) ! 17 Jan 17;
write(ounit,'("xspech : " 10x " : nBnc=[",999(es11.03,","))') efmn(1:mn) ! 17 Jan 17;
endif
endif
case( -1 ) ! mfreeits = -1 ; shall set vacuum normal field at computational boundary ; 24 Nov 16;
; iVns(2:mn) = iVns(2:mn) - iBns(2:mn) ! total normal field ; 24 Nov 16;
if( NOTstellsym) iVnc(1:mn) = iVnc(1:mn) - iBnc(1:mn)
; iBns(2:mn) = ofmn(2:mn) ! updated normal field due to plasma; 24 Nov 16;
if( NOTstellsym) iBnc(1:mn) = efmn(1:mn)
; iVns(2:mn) = iVns(2:mn) + iBns(2:mn) ! updated vacuum field ; 24 Nov 16;
if( NOTstellsym) iVnc(1:mn) = iVnc(1:mn) + iBnc(1:mn)
case( 0 ) ! mfreeits = 0 ; 09 Mar 17;
FATAL( xspech, .true., illegal mfreeits logic )
case( 1: ) ! mfreeits > 0 ; 09 Mar 17;
select case( Lzerovac )
case( 0 ) ! Lzerovac = 0 ; 09 Mar 17;
; iBns(2:mn) = gBnbld * iBns(2:mn) + ( one - gBnbld ) * ofmn(2:mn)
if( NOTstellsym) iBnc(1:mn) = gBnbld * iBnc(1:mn) + ( one - gBnbld ) * efmn(1:mn)
case( 1 ) ! Lzerovac = 1 ; 09 Mar 17;
; iBns(2:mn) = ofmn(2:mn) ! no blend; 27 Feb 17;
if( NOTstellsym) iBnc(1:mn) = efmn(1:mn)
; iVns(2:mn) = + iBns(2:mn) ! update vacuum field to cancel plasma field on computational boundary; 27 Feb 17;
if( NOTstellsym) iVnc(1:mn) = + iBnc(1:mn)
case default ! Lzerovac; 09 Mar 17;
FATAL( xspech, .true., invalid Lzerovac )
end select ! end select case( Lzerovac ) ; 27 Feb 17;
end select ! end select case( mfreeits ) ; 27 Feb 17;
else
LContinueFreeboundaryIterations = .false.
endif ! end of if( bnserr.gt.gBntol ) ; 24 Nov 16;
cput = GETTIME
if( myid.eq.0 ) then ; write(ounit,1003)
; ; write(ounit,1004) cput-cpus, nfreeboundaryiterations, mfreeits, gBntol, bnserr, cput-lastcpu
; ; write(ounit,1003)
endif ! end of if( myid.eq.0 ) ; 24 Nov 16;
1003 format("xspech : " 10x " : ")
1004 format("xspech : "f10.2" : nfreeboundaryiterations = "i6" / "i6.5" ; gBntol ="es8.1" ; bnserr =",es12.5," ; bnorml time ="f10.2"s ;")
endif ! end of if( LupdateBn ) ;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!latex \subsection{output files: vector potential}
!latex \begin{enumerate}
!latex \item The vector potential is written to file using \link{ra00aa}.
!latex \end{enumerate}
WCALL( xspech, ra00aa, ('W') ) ! this writes vector potential to file;
if( myid.eq.0 ) then ! write restart file; note that this is inside free-boundary iteration loop; 11 Aug 14;
wflag = 0 ; iflag = 0 ; rflag = zero
WCALL( xspech, wrtend, ( wflag, iflag, rflag ) ) ! write restart file; save initial input;
endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
! if( LContinueFreeboundaryIterations .and. Lfindzero.gt.0 .and. nfreeboundaryiterations.lt.mfreeits ) goto 9000
if( LContinueFreeboundaryIterations .and. nfreeboundaryiterations.lt.mfreeits ) goto 9000 ! removed Lfindzero check; Loizu Dec 18;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
! FREE-BOUNDARY ITERATIONS HAVE FINISHED;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
! SALLOCATE( gradient, (0:NGdof), zero )
!
! lastcpu = GETTIME
!
! LComputeDerivatives = .false.
!
! WCALL( xspech, dforce, ( NGdof, position(0:NGdof), gradient(0:NGdof), LComputeDerivatives ) ) ! (re-)calculate Beltrami fields;
!
! DALLOCATE(gradient)
!
!#ifdef DEBUG
! do vvol = 1, Mvol-1
! ; FATAL( xspech, BBe(vvol).lt.logtolerance, underflow )
! if( Igeometry.eq.3 ) then ! include spectral constraints; 04 Dec 14;
! ;FATAL( xspech, IIo(vvol).lt.logtolerance, underflow )
! endif
! if( NOTstellsym ) then
! ;FATAL( xspech, BBo(vvol).lt.logtolerance, underflow )
! if( Igeometry.eq.3 ) then ! include spectral constraints; 04 Dec 14;
! FATAL( xspech, IIe(vvol).lt.logtolerance, underflow )
! endif
! endif
! enddo
!#endif
!
! if( myid.eq.0 ) then
! cput = GETTIME
! write(ounit,1000) cput-cpus, nfreeboundaryiterations, ForceErr, cput-lastcpu, "|BB|e", alog10(BBe(1:min(Mvol-1,28)))
! if( Igeometry.ge.3 ) then ! include spectral constraints; 04 Dec 14;
! write(ounit,1001) "|II|o", alog10(IIo(1:min(Mvol-1,28)))
! endif
! if( NOTstellsym ) then
! write(ounit,1001) "|BB|o", alog10(BBo(1:min(Mvol-1,28)))
! if( Igeometry.ge.3 ) then ! include spectral constraints; 04 Dec 14;
! write(ounit,1001) "|II|e", alog10(IIe(1:min(Mvol-1,28)))
! endif
! endif
! endif
!
!2000 format("finish : ",f10.2," : finished ",i3," ; ":"|f|="es12.5" ; ":"time=",f10.2,"s ;":" log"a5,:"="28f6.2" ...")
!2001 format("finish : ", 10x ," : ",3x," ; ":" " 12x " ":" ", 10x ," ;":" log"a5,:"="28f6.2" ...")
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
if( myid.eq.0 ) then ! this is just screen diagnostics; 20 Jun 14;
cput = GETTIME
if( nPpts.gt.0 ) then
write(ounit,'("xspech : ", 10x ," :")')
write(ounit,'("xspech : ",f10.2," : myid=",i3," ; Poincare plot ; odetol="es8.1" ; nPpts="i7" ;":" nPtrj="24(i5",")" ...")') &
cput-cpus, myid, odetol, nPpts, nPtrj(1:min(Mvol,24))
endif
if( Lcheck.eq.1 ) then
write(ounit,'("xspech : ", 10x ," :")')
write(ounit,'("xspech : ",f10.2," : myid=",i3," ; calling jo00aa; computing error in field ;")') cput-cpus, myid
endif
endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!latex \subsection{final diagnostics}
!latex \begin{enumerate}
!latex \item \link{sc00aa} is called to compute the covariant components of the magnetic field at the interfaces;
!latex these are related to the singular currents;
!latex \item if \inputvar{Lcheck} $= 1$, \link{jo00aa} is called to compute the error in the Beltrami equation;
!latex \item \link{pp00aa} is called to construct the \Poincare plot;
!latex \end{enumerate}
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
do vvol = 1, Mvol
LREGION(vvol)
if( myid.ne.modulo(vvol-1,ncpu) ) cycle ! the following is in parallel; 20 Jun 14;
if( .not.ImagneticOK(vvol) ) then ; cput = GETTIME ; write(ounit,1002) cput-cpus ; write(ounit,1002) cput-cpus, myid, vvol, ImagneticOK(vvol) ; cycle
endif
WCALL( xspech, sc00aa, ( vvol, Ntz ) ) ! compute covariant field (singular currents);
if( Lcheck.eq.1 ) then
WCALL( xspech, jo00aa, ( vvol, Ntz, Iquad(vvol), mn ) )
endif
if( nPpts .gt.0 ) then
WCALL( xspech, pp00aa, ( vvol ) ) ! Poincare plots in each volume
endif
enddo ! end of do vvol = 1, Mvol; ! end of parallel diagnostics loop; 03 Apr 13;
1002 format("xspech : ",f10.2," :":" myid=",i3," ; vvol=",i3," ; IBeltrami="L2" ; construction of Beltrami field failed ;")
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
do vvol = 1, Mvol ; llmodnp = modulo(vvol-1,ncpu)
#ifdef DEBUG
FATAL( xspech, .not.allocated(Btemn), error )
FATAL( xspech, .not.allocated(Bzemn), error )
FATAL( xspech, .not.allocated(Btomn), error )
FATAL( xspech, .not.allocated(Bzomn), error )
#endif
RlBCAST( Btemn(1:mn,0:1,vvol), mn*2, llmodnp ) ! this is computed in sc00aa; 07 Dec 16;
RlBCAST( Bzemn(1:mn,0:1,vvol), mn*2, llmodnp )
RlBCAST( Btomn(1:mn,0:1,vvol), mn*2, llmodnp )
RlBCAST( Bzomn(1:mn,0:1,vvol), mn*2, llmodnp )
RlBCAST( beltramierror(vvol,1:3), 3, llmodnp ) ! this is computed in jo00aa; 21 Aug 18;
enddo ! end of do vvol = 1, Mvol; 01 Jul 14;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!latex \subsection{restart files}
!latex \begin{enumerate}
!latex \item \link{global}:\type{wrtend} is called to write the restart files.
!latex \end{enumerate}
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
if( myid.eq.0 ) then
wflag = 1 ; iflag = 0 ; rflag = zero
WCALL( xspech, wrtend, ( wflag, iflag, rflag ) ) ! write restart file; save initial input;
close(zunit) ! this file is written to in globals/wrtend; 11 Aug 14;
cput = GETTIME
write(ounit,'("xspech : ", 10x ," :")')
write(ounit,'("xspech : ",f10.2," : myid=",i3," : time="f8.2"m = "f6.2"h = "f5.2"d ;")') cput-cpus, myid, (cput-cpus) / (/ 60, 60*60, 24*60*60 /)
endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
9999 continue
WCALL( xspech, ending )
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
stop
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
end program xspech
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!latex \subsection{subroutine : ending}
!latex \begin{enumerate}
!latex \item Closes output files, writes screen summary.
!latex \end{enumerate}
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
subroutine ending
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
use constants, only : zero
use fileunits, only : ounit
use inputlist, only : Wmacros, Wxspech, ext, Ltiming
use cputiming
use allglobal, only : myid, cpus, mn
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
LOCALS
REAL :: Ttotal, dcpu, ecpu
CHARACTER :: date*8, time*10
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
cpui = GETTIME ; cpuo = cpui ! see macro expansion for begin; 11 Aug 14;
#ifdef DEBUG
if( Wxspech ) write(ounit,'("ending : ",f10.2," : myid=",i3," ; start ;")') cpui-cpus, myid
#endif
cput = GETTIME
! SUMTIME ! this is expanded by Makefile, and then again by macros; do not remove;
cput = GETTIME ; dcpu = cput-cpus
if( Ltiming .and. myid.eq.0 ) then
Ttotal = zero
! PRTTIME ! this is expanded by Makefile, and then again by macros; do not remove;
write(ounit,'("ending : ",f10.2," : time spent in wrtend =",f10.2," ;")') dcpu, Twrtend ; Ttotal = Ttotal + Twrtend
write(ounit,'("ending : ",f10.2," : time spent in readin =",f10.2," ;")') dcpu, Treadin ; Ttotal = Ttotal + Treadin
ecpu = Ttotal-dcpu ! error in actual cpu time and calculated cpu time; 7 Mar 13;
write(ounit,'("ending : ",f10.2," : Ttotal =",f10.2," s = "f8.2" m = "f6.2" h ; Timing Error = ",f10.2,"s = ",f10.2,"%")') &
dcpu, Ttotal / (/ 1, 60, 3600 /), ecpu, 100*ecpu/dcpu
endif ! end of if( Ltiming .and. myid.eq.0 ) then; 01 Jul 14;
if( myid.eq.0 ) then
call date_and_time(date,time)
write(ounit,'("ending : ", 10x ," : ")')
write(ounit,1000) dcpu, myid, dcpu / (/ 1, 60, 60*60, 24*60*60 /), date(1:4), date(5:6), date(7:8), time(1:2), time(3:4), time(5:6), ext
write(ounit,'("ending : ", 10x ," : ")')
if( myid.eq.0 ) then
call hdfint ! 18 Jul 14;
endif
endif ! end of if( myid.eq.0 ) ; 14 Jan 15;
MPIFINALIZE
1000 format("ending : ",f10.2," : myid=",i3," ; completion ; time=",f10.2,"s = "f8.2"m = "f6.2"h = "f5.2"d ; date= "&
a4"/"a2"/"a2" ; time= "a2":"a2":"a2" ; ext = "a60)
stop
end subroutine ending
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!