-
Notifications
You must be signed in to change notification settings - Fork 86
/
covr.f90
2250 lines (2152 loc) · 63.9 KB
/
covr.f90
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
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
module covm
! provides subroutine covr for NJOY2016
use locale
implicit none
private
public covr
! global variables
! user's input
integer::nin,nout
integer::nplot
integer::icolor
real(kr)::epmin
real(kr)::einc
integer::irelco
integer::ncase
integer::noleg
integer::nstart
integer::ndiv
integer::matype
character(12)::hlibid
character(21)::hdescr
integer::mat,mt,mat1,mt1
integer::mfflg,mf3,mf5,mf35
! other common variables
integer,parameter::ncasemx=100
integer::nscr,nscr1,nscr2
integer::nrow,ncol,itype,nvf,ncf
integer::ixmin,ixmax
integer::izero,ismall
integer::ishade
integer::nlev
integer::nexp
integer::iza,iza1,izap
integer::nwscr,max,nwcf,nwig
real(kr),dimension(:),allocatable::scr
real(kr),dimension(:),allocatable::xlev
real(kr),dimension(:),allocatable::cf
real(kr),dimension(:),allocatable::x,y
real(kr),dimension(:),allocatable::xx,xy
real(kr),dimension(:),allocatable::rsdx,rsdy
integer,dimension(:),allocatable::imtx,imtx1,imatx1
contains
subroutine covr
!--------------------------------------------------------------------
!
! Plot covariance data from errorr or make a condensed library.
!
! In the plot option, covr plots a matrix of correlation
! coefficients and an associated pair of standard deviation
! vectors, i.e.,a covariance matrix. The correlation
! matrix is plotted as a shaded contour plot and the vectors
! are plotted as semi-log plots, one rotated by 90 degrees.
! The log energy grids for the vector plots are identical
! to the grids for the matrix plot. This version plots
! through viewr.
!
! In the library option, covr produces a condensed bcd
! covariance library in the boxer format. This format is
! efficient for matrices of simple blocks.
!
!---input specifications (free format)---------------------------
!
! card 1
! nin input tape unit
! nout output tape unit
! (default=0=none)
! nplot viewr output unit
! (default=0=none)
!
! ---cards 2, 2', 2a, and 3a for nout.le.0 only (plot option)
!
! card 2
! icolor select color or monochrome style
! 0=monochrome (uses cross hatching)
! 1=color background and contours
! 2=color background and contours plus
! card 2' follows.
! (default=0)
! card 2' (only when icolor=2)
! nlev,(tlev(i),i=1,nlev)
! defines the number of correlation matrix
! intervals and their boundaries. Zero is
! assumed as the lower limit of the first
! boundary, but the User must specify unity
! as the upper limit of the last boundary.
! nlev is a positive integer .le. 9.
! default values (when icolor=1) are:
! 6,0.001,0.1,0.2,0.3,0.6,1.0
! card 2a
! epmin lowest energy of interest (default=0.)
! card 3a
! irelco type of covariances present on nin
! 0/1=absolute/relative covariances
! (default=1)
! ncase no. cases to be run (maximum=60)
! (default=1)
! noleg plot legend option
! -1/0/1=legend for first subcase only/
! legend for all plots/no legends
! (default=0)
! nstart sequential figure number
! 0/n=not needed/first figure is figure n.
! (default=1)
! ndiv no. of subdivisions of each of the
! gray shades (default=1)
!
! ---cards 2b, 3b, and 3c for nout gt 0 (library option) only--
!
! card 2b
! matype output library matrix option
! 3/4=covariances/correlations
! (default=3)
! ncase no. cases to be run (maximum=ncasemx=100)
! (default=1)
! card 3b
! hlibid up to 6 characters for identification
! card 3c
! hdescr up to 21 characters of descriptive
! information
!
! ---cards 4 for both options---
!
! card 4
! mat desired mat number
! mt desired mt number
! mat1 desired mat1 number
! mt1 desired mt1 number
! (default for mt, mat1 and mt1 are 0,0,0
! meaning process all mts for this mat
! with mat1=mat)
! (neg. values for mt, mat1, and mt1 mean
! process all mts for this mat, except for
! the mt-numbers -mt, -mat1, and -mt1. in
! general, -n will strip both mt=1 and mt=n.
! -4 will strip mt=1, mt=3, and mt=4, and
! -62, for example, will strip mt=1, mt=62,
! mt=63, ... up to and incl. mt=90.)
! repeat card 4 ncase times
!
! note---if more than one material appears on the input tape,
! the mat numbers must be in ascending order.
!
!--------------------------------------------------------------------
use mainio ! provides nsysi,nsyso,nsyse
use util ! provides timer,openz,error,repoz,mess
use endf ! provides endf routines and variables
! internals
integer,parameter::nlevmx=9
integer::iomit,nfigmx,ncamx,ne,i,nfig,nwlev,ilev,j,n,nb,nw
integer::nwy,ixpn
real(kr)::sec,da,tlow,tx
character(60)::strng
character(6)::hinpid
real(kr),dimension(:),allocatable::xig,yig
real(kr),dimension(nlevmx)::tlev=(/.001e0_kr,.1e0_kr,.2e0_kr,&
.3e0_kr,.6e0_kr,1.0e0_kr,0.0e0_kr,0.0e0_kr,0.0e0_kr/)
integer::imat(ncasemx),imt(ncasemx),imat1(ncasemx),imt1(ncasemx)
real(kr),parameter::rdn=0.999998e0_kr
character(5),dimension(3)::hreas=(/'null ','small','empty'/)
integer,parameter::ntics3=600
!--initialize
call timer(sec)
write(nsyso,'(/&
&'' covr...process covariance data'',38x,f8.1,''s'')') sec
write(nsyse,'(/'' covr...'',61x,f8.1,''s'')') sec
nfigmx=ncasemx
max=1+nfigmx*(nfigmx+1)/2
mf3=3
mf5=5
ncamx=ncasemx
iomit=0
!--read and write out user input
nout=0
nplot=0
read(nsysi,*) nin,nout,nplot
call openz(nin,0)
call openz(nout,1)
if (nout.le.0) then
nscr1=11
nscr2=0
call openz(nscr1,1)
call repoz(-nscr1)
write(nscr1,'(/&
&'' the following plots were omitted or were plotted but '',&
&''empty'',21x/4x,''mat'',5x,''mt'',3x,''mat1'',4x,''mt1'',3x,&
&''reason'',43x/3x,''----'',4x,''---'',3x,''----'',4x,&
&''---'',3x,''------'',43x)')
! input for plot option
icolor=0
read(nsysi,*) icolor
nlev=6
if (icolor.eq.2) then
read(nsysi,*)nlev,(tlev(i),i=1,nlev)
if (nlev.gt.nlevmx) then
write(strng,'('' nlev>'',i2,'' not allowed'')')nlevmx
call error('covr',strng,' ')
endif
do i=2,nlev
if (tlev(i).le.tlev(i-1)) then
call error('covr',&
'tlev array must sequentially increase',&
' ')
endif
enddo
if (tlev(nlev).ne.1) then
write(strng,'(''reset tlev('',i2,'') from '',1pe11.4,&
&'' to 1.0'')')nlev,tlev(nlev)
call mess('covr',strng,' ')
tlev(nlev)=1
endif
endif
epmin=0
read(nsysi,*) epmin
epmin=epmin*rdn
irelco=1
ncase=1
noleg=0
nstart=1
ndiv=1
read(nsysi,*) irelco,ncase,noleg,nstart,ndiv
if (ndiv.eq.0) ndiv=1
if (ncase.gt.ncamx)&
call error('covr','requested too many cases.',' ')
else
! input for library option
nscr1=11
nscr2=12
call openz(-nscr1,1)
call openz(-nscr2,1)
call repoz(-nscr1)
call repoz(-nscr2)
! do not translate abs. covariances to rel. in library option,
! thus allowing stand. devs. and covariances to track the
! errorr output.
irelco=1
matype=3
ncase=1
read(nsysi,*) matype,ncase
if (matype.ne.4) matype=3
if (ncase.le.0) ncase=1
if (ncase.gt.ncamx)&
call error('covr','requested too many cases.',' ')
read(nsysi,*) hinpid
hdescr=' '
read(nsysi,*) hdescr
endif
do i=1,ncase
imt(i)=0
imat1(i)=0
imt1(i)=0
read(nsysi,*) imat(i),imt(i),imat1(i),imt1(i)
enddo
if (nout.eq.0) write(nsyso,'(/&
&'' unit for input covariance tape ....... '',i10/&
&'' unit for output covariance tape ...... '',i10/&
&'' unit for plot output ................. '',i10/&
&'' icolor ............................... '',i10/&
&'' rel. cov. option (0=abs/1=rel) ....... '',i10/&
&'' ncase ................................ '',i10/&
&'' legend option (-1=first/0=all/1=none) '',i10/&
&'' starting fig. no. (0=omit fig. nos.) . '',i10/&
&'' no. subdivisions of each shade ....... '',i10/&
&'' minimum energy to be plotted ......... '',1p,e10.3)')&
nin,nout,nplot,icolor,irelco,ncase,noleg,nstart,ndiv,epmin
if (nout.gt.0) write(nsyso,'(/&
&'' unit for input covariance tape ....... '',i10/&
&'' unit for output covariance tape ...... '',i10/&
&'' output matrix option (3=cov./4=corr.) '',i10/&
&'' no. cases to run ..................... '',i10/&
&'' library identification ............... ''/6x,a6/&
&'' library description .................. ''/6x,a21)')&
nin,nout,matype,ncase,hinpid,hdescr
write(nsyso,'(/12x,''mat'',6x,''mt'',6x,''mat1'',6x,''mt1''/&
&11x,''----'',5x,''---'',6x,''----'',6x,''---''/&
&(11x,i4,5x,i3,6x,i4,6x,i3))')&
(imat(i),imt(i),imat1(i),imt1(i),i=1,ncase)
!--initialize parameters for plots
if (nout.le.0) then
nfig=nstart-1
! set up shade levels
nwlev=nlev*ndiv
allocate(xlev(nwlev))
ilev=0
tlow=0
do i=1,nlev
tx=(tlev(i)-tlow)/ndiv
do j=1,ndiv
ilev=ilev+1
xlev(ilev)=tlow+tx*j
enddo
tlow=tlev(i)
enddo
da=1
da=da/1000
xlev(ilev)=xlev(ilev)+da
nlev=ilev
endif
call repoz(nin)
nscr=10
nsc=0
if (nin.lt.0) nscr=-nscr
call openz(nscr,1)
nwscr=npage*2+50
!--read the first two records to verify that this is a legal tape
!--set the mf35 flag and determine the number of groups for scr
allocate(scr(17))
call tpidio(nin,0,0,scr,nb,nw)
call contio(nin,0,0,scr,nb,nw)
if (mfh.ne.1.or.mth.ne.451) call error('covr','illegal input tape',' ')
mfflg=nint(scr(5))
if (mfflg.eq.-11.or.mfflg.eq.-14) then
mf35=mf3
else if (mfflg.eq.-12) then
mf35=mf5
else
call error('covr','illegal errorr output tape for covr',' ')
endif
call contio(nin,0,0,scr,nb,nw)
nwscr=l1h+7
if (nwscr.lt.17) nwscr=17
deallocate(scr)
allocate(scr(nwscr))
call repoz(nin)
call openz(nplot,1)
if (nplot.ne.0) write(nplot,'(''1 2 .22'',i3,''/'')') icolor
!--loop over cases
do 130 n=1,ncase
mat=imat(n)
! copy this mat from nin to nscr
call repoz(nin)
call repoz(nscr)
call tpidio(nin,0,nscr,scr,nb,nw)
call finds(mat,0,0,nin)
call contio(nin,0,nscr,scr,nb,nw)
iverf=l1h
iza=nint(scr(1))
iza1=iza
call tomend(nin,0,nscr,scr)
allocate(imtx(max))
allocate(imtx1(max))
allocate(imatx1(max))
nexp=1
imtx(1)=imt(n)
imatx1(1)=imat1(n)
imtx1(1)=imt1(n)
if (imat1(n).ne.imat(n).and.imat1(n).gt.0) then
! if needed, also copy mat1 from nin to nscr
call finds(imat1(n),0,0,nin)
call contio(nin,0,nscr,scr,nb,nw)
iza1=nint(scr(1))
call tomend(nin,0,nscr,scr)
endif
call atend(0,nscr)
call repoz(nscr)
!--expand the mt-mt1 list
if (imt(n).le.0) call expndo(nscr)
!--loop over mt-mt1 pairs
do 191 ne=1,nexp
mt=imtx(ne)
mat1=imatx1(ne)
if (mat1.le.0) mat1=mat
mt1=imtx1(ne)
if (mt1.eq.0) mt1=mt
write(nsyso,'(/&
&'' reaction-pairs processed '',i4,2x,i3,3x,i4,2x,i3)')&
mat,mt,mat1,mt1
!--read through input errorr file, calculating relative standard
!--deviations and correlation coefficients
call corr(nscr,mat,mt,mat1,mt1,izap)
if (nout.gt.0) go to 210
if (ismall.gt.0) go to 155
150 continue
iomit=iomit+1
! document null and small matrices on nscr1, then skip plot
write(nscr1,'(4i7,3x,a5,44x)') mat,mt,mat1,mt1,hreas(izero+1)
go to 190
155 continue
nwcf=ixmax*ixmax
nwig=2*(2*(ixmax+1)+ntics3)
allocate(xig(nwig))
nwy=nwig
if (nwig.lt.nwcf) nwy=nwcf
allocate(yig(nwy))
! cf storage is in the middle. move it to the top
do i=1,nwcf
yig(i)=cf(i)
enddo
deallocate(cf)
call truncg(epmin,x,y,xx,xy,ixmin,ixmax)
! cf storage goes into highest part
allocate(cf(nwcf))
do i=1,nwcf
cf(i)=yig(i)
yig(i)=0
enddo
if (izero.eq.0) go to 190
ixpn=ixmax-ixmin+1
ishade=0
call plotit(x(ixmin:),y(ixmin:),xig,yig,ixpn,ixmax,&
rsdx(ixmin:),rsdy(ixmin:),xlev,noleg,ne,izap)
if (ishade.eq.0) izero=2
if (ishade.eq.0) go to 150
if (nfig-nstart+1.lt.nfigmx) go to 190
write(strng,'(''mat='',i4,'' mt='',i3,'' mat1='',i4,'' mt1='',i3,&
&'' nfig='',i3)') mat,mt,mat1,mt1,nfig
call mess('covr','have plotted all that fit.',strng)
go to 310
!--process library option
210 continue
deallocate(y)
deallocate(xy)
if (izero.ne.0) go to 220
write(nsyso,'(/&
&3x,''null covariance matrix. output suppressed.'')')
go to 190
220 continue
if (ne.le.1) then
if (matype.eq.3) write(hlibid,'(a6,''-a-'',i3)') hinpid,ixmax
if (matype.eq.4) write(hlibid,'(a6,''-b-'',i3)') hinpid,ixmax
endif
nrow=ixmax+1
ncol=1
nvf=10
ncf=3
itype=0
! for first sub-case of first case, write the group structure
if (n.eq.1.and.ne.eq.1) call press(0,nout,x,nrow,ncol)
nrow=ixmax
if (mat1.eq.mat.and.mt1.eq.mt) then
itype=1
call press(0,nout,xx,ixmax,1)
itype=2
call press(0,nout,rsdx,ixmax,1)
endif
ncol=nrow
if (mat1.eq.mat.and.mt1.eq.mt) ncol=0
nvf=7
if (matype.eq.3) nvf=10
ncf=6
if (ixmax.le.100) ncf=5
if (ixmax.le.30) ncf=4
itype=matype
call press(0,nout,cf,2,ixmax)
190 continue
if (allocated(x))deallocate(x)
if (allocated(y))deallocate(y)
if (allocated(rsdx)) deallocate(rsdx)
if (allocated(rsdy)) deallocate(rsdy)
if (allocated(xig)) deallocate(xig)
if (allocated(yig)) deallocate(yig)
if (allocated(xx)) deallocate(xx)
if (allocated(xy)) deallocate(xy)
if (allocated(cf)) deallocate(cf)
191 continue
deallocate(imtx)
deallocate(imtx1)
deallocate(imatx1)
130 continue
if (nout.eq.0) go to 310
nscr1=-nscr1
nscr2=-nscr2
go to 330
!--plot the table of contents
310 continue
write(nplot,'(''99/'')')
!--copy the statistics file to output
if (iomit.ne.0) then
write(nsyso,'(/)')
call copyst(nscr1,nsyso)
endif
!--finished
330 continue
if (allocated(scr)) deallocate(scr)
if (allocated(xlev)) deallocate(xlev)
call repoz(nin)
call repoz(nout)
call closz(nin)
call closz(nout)
call closz(nplot)
call closz(nscr)
call closz(nscr1)
call closz(nscr2)
call timer(sec)
write(nsyso,'(69x,f8.1,''s''/&
&1x,7(''**********''),''****'')') sec
return
end subroutine covr
subroutine expndo(nin)
!--------------------------------------------------------------------
! Expand the MT-MT1 list to include all possible combinations
! in ENDF order. Obtain the set of MT-numbers present by
! examining the cross section file on unit nin.
!--------------------------------------------------------------------
use util ! provides error,mess,repoz
use endf ! provides endf routines and variables
! externals
integer::nin
! internals
integer::nlstm,nb,nw,nmt,l,mtn,istr,im,jm
character(60)::strng
integer::mstrip(3)
integer,dimension(:),allocatable::lstm
nlstm=ncasemx
allocate(lstm(nlstm))
call repoz(nin)
call tpidio(nin,0,0,scr,nb,nw)
nmt=0
mstrip(1)=-imtx(1)
mstrip(2)=-imatx1(1)
mstrip(3)=-imtx1(1)
call finds(mat,mf35,0,nin)
130 continue
call listio(nin,0,0,scr,nb,nw)
if (mfh.ne.mf35) go to 160
l=1
mtn=mth
do while (nb.ne.0)
l=l+nw
if (l.gt.nwscr)&
call error('expndo','storage exceeded.',' ')
call moreio(nin,0,0,scr(l),nb,nw)
enddo
! read send card
call contio(nin,0,0,scr,nb,nw)
do 155 istr=1,3
if (mstrip(istr).eq.0) go to 155
if (mtn.eq.1) go to 180
if (mtn.eq.mstrip(istr)) go to 180
if (mtn.eq.3.and.mstrip(istr).eq.4) go to 180
if (mstrip(istr).lt.51.or.mstrip(istr).gt.90) go to 155
if (mtn.gt.mstrip(istr).and.mtn.le.90) go to 180
155 continue
nmt=nmt+1
lstm(nmt)=mtn
go to 130
160 continue
nexp=0
do im=1,nmt
do jm=im,nmt
nexp=nexp+1
if (nexp.gt.max) then
call error('expndo','storage exceeded.',' ')
endif
imtx(nexp)=lstm(im)
imatx1(nexp)=0
imtx1(nexp)=lstm(jm)
enddo
enddo
deallocate(lstm)
return
180 continue
write(strng,'(''mt='',i3,'' stripped from list.'')') mtn
call mess('expndo',strng,' ')
go to 130
end subroutine expndo
subroutine corr(nin,mat,mt,mat1,mt1,izap)
!--------------------------------------------------------------------
! Convert relative covariances to standard deviations
! and correlations.
!--------------------------------------------------------------------
use util ! provides error,repoz
! externals
integer::nin,mat,mt,mat1,mt1,izap
! internals
integer::icall,ixtest,ixp,i,ig,ii,j,nswap
real(kr)::epmn,xcycle
character(60)::strng
real(kr),parameter::zp4=0.4e0_kr
real(kr),parameter::xsize=4.25e0_kr
real(kr),parameter::epmn0=0.499999e-4_kr
real(kr),parameter::zero=0
ismall=0
!--for cross-reaction matrices, check to see if the matrix
!--is null before proceeding with the correlation calculation
icall=0
if (mat.ne.mat1.or.mt.ne.mt1) then
call covard(nin,mat,mt,mat1,mt1,izap,icall)
icall=1
if (izero.eq.0) return
endif
110 continue
call covard(nin,mat1,mt1,mat1,mt1,izap,icall)
ixtest=ixmax
if (icall.eq.0) icall=1
allocate(rsdx(ixmax))
allocate(rsdy(ixmax))
icall=iabs(icall)
go to (130,165,190),icall
130 continue
ixp=ixmax+1
do i=1,ixp
x(i)=y(i)
enddo
if (nout.le.0) then
epmn=epmn0
xcycle=0
do while (xcycle.lt.zp4)
epmn=epmn*2
xcycle=xsize/log10(x(1+ixmax)/epmn)
enddo
! have found epmin value that satifies criterion
if (epmn.gt.epmin) then
epmin=epmn
write(strng,&
'(''epmin reset to'',1p,e12.4,'', xcycle='',1p,e12.4)')&
epmin,xcycle
call mess('corr',strng,' ')
endif
do i=1,ixmax
if (cf(ixmax*(i-1)+i).gt.0) then
rsdy(i)=sqrt(cf(ixmax*(i-1)+i))
else
rsdy(i)=0
endif
rsdx(i)=rsdy(i)
enddo
else
call repoz(-nscr2)
do i=1,ixmax
read(nscr2) ig,(cf(ii),ii=1,ixmax)
rsdy(i)=sqrt(cf(i))
rsdx(i)=rsdy(i)
enddo
endif
if (mat.eq.mat1.and.mt.eq.mt1) go to 190
icall=icall+1
call covard(nin,mat,mt,mat,mt,izap,icall)
if (icall.lt.0) go to 110
165 continue
if (nout.le.0) then
do i=1,ixmax
rsdx(i)=sqrt(cf(ixmax*(i-1)+i))
enddo
else
call repoz(-nscr2)
do i=1,ixmax
read(nscr2) ig,(cf(ii),ii=1,ixmax)
rsdx(i)=sqrt(cf(i))
enddo
endif
icall=icall+1
call covard(nin,mat,mt,mat1,mt1,izap,icall)
if (icall.lt.0) go to 110
190 continue
if (izero.ne.0) then
if (nout.le.0) then
do i=1,ixmax
do j=1,ixmax
if (cf(ixmax*(i-1)+j).ne.zero.and.&
rsdx(i)*rsdy(j).ne.zero) then
cf(ixmax*(i-1)+j)=cf(ixmax*(i-1)+j)/&
(rsdx(i)*rsdy(j))
! test for the presence of any significant (i.e.,
! plotable) values of the correlation coefficient
if (abs(cf(ixmax*(i-1)+j)).ge.xlev(1)) ismall=1
else
cf(ixmax*(i-1)+j)=0
endif
enddo
enddo
else
call repoz(-nscr2)
call repoz(-nscr1)
if (matype.ne.3) then
do i=1,ixmax
read(nscr2) ig,(cf(ii),ii=1,ixmax)
do j=1,ixmax
if (cf(j).ne.zero.and.&
rsdx(i)*rsdy(j).ne.zero) then
cf(j)=cf(j)/(rsdx(i)*rsdy(j))
else
cf(j)=0
endif
enddo
write(nscr1) ig,(cf(ii),ii=1,ixmax)
enddo
call repoz(-nscr2)
else
nswap=nscr1
nscr1=nscr2
nscr2=nswap
endif
endif
endif
! finished
if (ixmax.ne.ixtest)&
call error('corr','group structures do not agree.',' ')
return
end subroutine corr
subroutine covard(nin,mat,mt,mat1,mt1,izap,icall)
!--------------------------------------------------------------------
! Read covariance data in the ENDF-like format output by
! njoy/errorr and transform it into engineering format,
! that is, with full matrices (zeroes included).
!--------------------------------------------------------------------
use mainio ! provides nsysi,nsyso,nsyse
use util ! provides repoz,finds,error,mess
use endf ! provides ENDF routines and variables
! externals
integer::nin,mat,mt,mat1,mt1,izap,icall
! internals
integer::nb,nw,l,i,k,nmt,imt,k1,mat1x,mtx,kgp,kl,lgp1
integer::j1,j2,kj,j,ir,ipflag,ind,kk,n,mf3x,ixp,ixpo
character(60)::strng,strng2
integer::igprnt=0
real(kr),parameter::zero=0
save ixpo,igprnt
if (icall.le.1) ixpo=0
call repoz(nin)
call tpidio(nin,0,0,scr,nb,nw)
call finds(mat,1,451,nin)
if (nout.gt.0) call repoz(-nscr1)
if (nout.gt.0) call repoz(-nscr2)
!--read group structure
call contio(nin,0,0,scr,nb,nw)
call listio(nin,0,0,scr,nb,nw)
ixmax=l1h
ixp=n1h
if (icall.eq.0) go to 120
if (icall.eq.1) go to 130
if (ixpo.eq.0.and.ixp.eq.igprnt) go to 120
if (ixp.le.ixpo) go to 130
deallocate(x)
deallocate(y)
deallocate(cf)
deallocate(xx)
deallocate(xy)
icall=-icall
120 continue
allocate(x(ixp))
allocate(y(ixp))
allocate(xx(ixmax))
allocate(xy(ixmax))
nwcf=ixmax*(ixmax+3)
if (nout.gt.0) nwcf=ixmax*2
allocate(cf(nwcf))
130 continue
l=1
do while (nb.ne.0)
l=l+nw
if (l.gt.nwscr)&
call error('covard','storage exceeded.',' ')
call moreio(nin,0,0,scr(l),nb,nw)
enddo
do i=1,ixp
y(i)=scr(i+6)
enddo
if (igprnt.ne.ixp) write(nsyso,'(/&
&'' group structure''/'' ---------------''/&
&(1x,1p,6e12.4))') (y(i),i=1,ixp)
igprnt=ixp
ixpo=ixp
!--read cross sections for mt and mt1
k=1
call finds(mat,mf35,mt,nin)
200 continue
call listio(nin,0,0,scr,nb,nw)
if (mf35.eq.5) einc=scr(2)
l=1
do while (nb.ne.0)
l=l+nw
if (l.gt.nwscr)&
call error('covard','storage exceeded.',' ')
call moreio(nin,0,0,scr(l),nb,nw)
enddo
if (k.gt.1) go to 240
do i=1,ixmax
xx(i)=scr(i+6)
enddo
if (mth.ne.mt1) go to 260
240 continue
do i=1,ixmax
xy(i)=scr(i+6)
enddo
260 continue
if (k.gt.1) go to 270
if (mat.eq.mat1.and.mt.eq.mt1) go to 270
call finds(mat1,mf35,mt1,nin)
k=k+1
go to 200
270 continue
do i=1,nwcf
cf(i)=0
enddo
!--read covariances.
mf3x=33
if (mfflg.eq.-14) mf3x=40
if (mf35.eq.5) mf3x=35
if (mt.eq.251) mf3x=34
call finds(mat,mf3x,mt,nin)
call contio(nin,0,0,scr,nb,nw)
nmt=n2h
if (nmt.eq.0) go to 450
imt=0
k1=0
! search for desired subsection of this mt
do
imt=imt+1
if (imt.gt.nmt) then
write(strng,'(''did not find file '',i2,'' subsection'')')mf3x
write(strng2,'(''for mt='',i3,'' mat='',i4)') mt1,mat1
call error('covard',strng,strng2)
endif
izap=0
if (mfflg.eq.-14) then
if (mat.eq.mat1.and.mt.eq.mt1) then
izap=nint(c2h)
endif
endif
call contio(nin,0,0,scr,nb,nw)
mat1x=l1h
if (mf3x.eq.34) mat1x=math
if (mat1x.eq.0) mat1x=math
mtx=l2h
if (mf3x.eq.34) mtx=l1h
kgp=n2h
if (mat1x.eq.mat1.and.mtx.eq.mt1) exit
do
call listio(nin,0,0,scr,nb,nw)
kl=n2h
do while (nb.ne.0)
call moreio(nin,0,0,scr,nb,nw)
enddo
if (kl.ge.kgp) exit
enddo
enddo
do
call listio(nin,0,0,scr,nb,nw)
lgp1=l2h
kl=n2h
l=1
do while (nb.ne.0)
l=l+nw
if (l.gt.nwscr)&
call error('covard','storage exceeded.',' ')
call moreio(nin,0,0,scr(l),nb,nw)
enddo
j1=7
j2=l+nw-1
if (nout.le.0) then
kj=0
do j=j1,j2
kj=kj+1
cf(ixmax*(kl-1)+lgp1-1+kj)=scr(j)
enddo
else
nw=j2-j1+1
write(nscr1) kl,lgp1,nw,(scr(j),j=j1,j2)
if (k1.eq.0) k1=kl
endif
if (kl.ge.kgp) exit
enddo
if (nout.eq.0) go to 410
call repoz(-nscr1)
kl=k1
ir=0
!--eliminate spurious covariances in regions of zero xsec,
!--convert absolute covariances to relative,
!--and test for a null covariance matrix
410 continue
ipflag=0
izero=0
do 420 k=1,ixmax
if (nout.eq.0) ind=ixmax*(k-1)
if (nout.eq.0) go to 425
ind=0
do kk=1,ixmax
cf(kk)=0
enddo
if (k.lt.kl) go to 425
if (ir.eq.0) read(nscr1) kl,lgp1,nw,(scr(i),i=1,nw)
ir=ir+1
if (k.lt.kl) go to 425
ir=0
do kk=1,nw
cf(lgp1+kk-1)=scr(kk)
enddo
425 continue
do 428 n=1,ixmax
if (xx(k)*xy(n).ne.zero) go to 430
if (cf(ind+n).eq.zero) go to 428
cf(ind+n)=0
ipflag=ipflag+1
if (ipflag.eq.1) then
write(strng,&
'(''mt='',i3,'' ig='',i3,'' mt1='',i3,'' ig1='',i3)')&
mt,k,mt1,n
call mess('covard',&
'due to zero xsec, covariance zeroed for',strng)
endif
if (ipflag.eq.2) call mess('covard','additional warnings suppressed',' ')
go to 428
430 continue
if (irelco.ne.1) cf(ind+n)=cf(ind+n)/(xx(k)*xy(n))
if (cf(ind+n).ne.zero) izero=1
428 continue
if (nout.ne.0) then
write(nscr2) k,(cf(kk),kk=1,ixmax)
endif
420 continue
450 return
end subroutine covard
subroutine truncg(epmin,x,y,xx,xy,ixmin,ixmax)
!--------------------------------------------------------------------
! Truncate the lower end of the MT and MT1 energy group
! structures to eliminate zero cross section regions from plots.
! Define ethresh to be the lower of either the actual threshold
! or 1 MeV (emin). Then the plots will be truncated at ethresh
! or epmin, whichever is higher.
!--------------------------------------------------------------------
use util ! provides error
! externals
integer::ixmin,ixmax
real(kr)::epmin,x(*),y(*),xx(*),xy(*)
! internals
integer::i,iymin,ielo
real(kr)::rlimx,rlimy,ethrx,ethry,elo
real(kr),parameter::emin=0.9999e6_kr
real(kr),parameter::xslim=1.e-4_kr
real(kr),parameter::zero=0
real(kr),parameter::ten=10.e0_kr
!--define alternative limit as average/xslim
rlimx=0
rlimy=0
ethrx=x(1)
ethry=y(1)
do i=2,ixmax
rlimx=rlimx+xx(i-1)*(x(i)-x(i-1))
rlimy=rlimy+xy(i-1)*(y(i)-y(i-1))
if (xx(i-1).le.zero) ethrx=x(i)
if (xy(i-1).le.zero) ethry=y(i)
enddo
if (rlimx.ne.zero.and.rlimy.ne.zero) then
rlimx=xslim*rlimx/(y(ixmax)-ethry)
rlimy=xslim*rlimy/(x(ixmax)-ethrx)
endif
ixmin=1
do i=1,ixmax
if (x(1+i).gt.epmin) then
if (xx(i).ge.xslim.or.xx(i).ge.rlimx) go to 120
if (x(1+i).gt.emin) go to 120
endif
ixmin=i+1
enddo
call error('truncg','bad data.',' ')
120 continue
iymin=1
do i=1,ixmax
if (y(1+i).gt.epmin) then
if (xy(i).ge.xslim.or.xy(i).ge.rlimy) go to 140
if (y(1+i).gt.emin) go to 140
endif
iymin=i+1
enddo
call error('truncg','bad data.',' ')
140 continue
if (iymin.lt.ixmin) ixmin=iymin
! limit bin width of thermal groupr to one decade
if (ixmin.eq.1) then
if (10*x(1).lt.x(2).and.10*y(1).lt.y(2)) then
elo=log10(x(2)/10)
ielo=nint(elo)