-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathadd-ab_flags.f
More file actions
1766 lines (1735 loc) · 65.1 KB
/
add-ab_flags.f
File metadata and controls
1766 lines (1735 loc) · 65.1 KB
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
c add-ab_flags cloned from add-msk-col
c add-msk-col - cloned from SampleCols & add-opt1-cols
c add msk column to an mdex file
c
c vsn 1.0 B80810: initial version
c 1.1 B80813: switched to 32-bit mask image (from 16-bit)
c 1.2 B80814: added processing of bits 21 and 22
c 1.3 B80816: added translation of mask bits to cc_flags;
c changed name to add-cc_flags
c 1.4 B80822: changed initial tagging from DPHO to dpho
c 1.5 B80824: added halo bits tot he processing
c 1.6 B80904: changed name to add-ab_flags; added special
c processing for the "_af" mdex format
c 1.6 B80913: added command-line flags
c 1.6 B80914: added multi-match processing
c 1.6 B80916: fixed cc_flags symbol priority
c 1.6 B80920: separated parent-star-pixel bit clearing by band
c 1.7 B80925: added geometric diffraction spike bits
c 1.8 B80926: major changes to accommodate mrgad not averaging
c w?x and w?y
c 1.8 B80930: added file output to warn if nSrc mismatch
c 1.9 B81005: added alarm/band-aid fix for NaNs in w?mpro &c.
c 1.91 B81005: added counters for NaN combinations
c 1.92 B81101: added option to recompute magnitudes from fluxes
c 1.93 B81103: added many debug checkpoints to trace NaNs
c 1.94 B81105: removed NaN debugging and band-aid NaN fix
c 1.95 B81106: added "pipe" table header conversion to "double" for
c ra, dec, MJDs, elon, & elat
c 1.96 B81120: installed mag upper limit logic
c 1.97 B81207: restored 4-band cc_flags; added w?sat processing and
c w?mcor values
c 1.98 B81210: installed w?sat & w?cov processing
c 1.99 B81214: installed fix for |dec| > 90
c 1.99 B81215: installed fix for |dec_pm| > 90;
c changed w?cov to f8.1 to accomodate > 9999
c 1.99 B81216: defaulted MaxDec to 60 for fixup of mrgad RA Dec;
c changed ab_flags to upper case; copy ra & dec to
c ra_pm and dec_pm when doing MaxDec fix
c 2.0 B81220: removed all attempts to fix |Dec| > 90 problem (now
c fixed in mrgad); defaulted DoMags to false.
c 2.1 B81230: bullet-proofed for incomplete CatWISE/IRSA matching
c 2.2 B90114: added -999.0 clipping for w?snr; set DoMags back to T
c 2.3 B90118: added checking for non-blank characters under "pipes"
c 2.4 B90124: added table header lines identifying version, date
c and time of run, and mask used; removed "cleanup"
c and "equinox" from header pass-through, added "mrgad"
c 2.5 B90127: changed handling of input headers; set DoMags=F
c 2.6 B90204: fixed bug in ERROR_DATA_UNDER_PIPR file name
c
c=======================================================================
c
Integer*4 MaxFld
Parameter (MaxFld = 1000)
c
character*150000 Hfits
Character*5000 Line, HdrLine
Character*500 InFNam, MskNam, OutFNam, NLNam, Cov1Nam, Cov2Nam,
+ ErrNam
Character*50 w1cc_map_str, w2cc_map_str, cc_flags, w1cc_map,
+ w2cc_map, dist_str
Character*25 Field(MaxFld), Field2(MaxFld)
Character*13 w1ab_map_str, w2ab_map_str
Character*11 Vsn, NumStr
Character*9 ab_flags, w1ab_map, w2ab_map
Character*8 CDate, CTime
Character*5 nAWstr
Character*3 Flag, Flag0
Character*1 w3cc, w4cc, w3cc2, w4cc2
Real *8 ra, dec, x8, y8, flux, sigflux, mag, sigmag,
+ w1m0, w2m0, CoefMag, wsnr
Real*4 w1x, w2x, w1y, w2y, dist, dist2, w1mcor, w2mcor,
+ wsat, wcov
Integer*4 IArgC, LNBlnk, FileID1, nHead, MskBitHist(32), msk,
+ NPlanes, NRows, NCols, I, J, N, IStat, NpixPL, k,
+ NPix, status, nSrc, IFa(MaxFld), IFb(MaxFld), NF,
+ w1abmap, w2abmap, w1ccmap, w2ccmap, notZero, nMag,
+ nSrcHdr, nAW, w1ccmap2, w2ccmap2, IOr, nArg, nArgs,
+ nw394, nw395, nw396, nw397, nw398, nw385, wcs,
+ offscl, nNaN, nn11, nn12, nn21, nn22, kBadw3,
+ kBadw4, kBad2w3, kBad2w4, kBadness, i1PSF, j1PSF,
+ iPix, jPix, nPSF, i2PSF, j2PSF, n1Sat, n2Sat,
+ n1Cov, n2Cov, src, src0, nIncomp, IFa2(MaxFld),
+ IFb2(MaxFld), NF2, nUnderPipe
Logical*4 NeedHelp, anynull, SanityChk, GoodXY1, GoodXY2,
+ BitSet, dbg, useWCS, NaNwarn, NaNstat1,
+ NaNpm1, NaNstat2, NaNpm2, doMags, doCov, GotN1,
+ GotN2, ReadAlready, GotEndHdr1
Integer*4 nullval
Integer*4, allocatable :: array1(:,:)
Integer*2 cov1(2048,2048), cov2(2048,2048)
c
Data Vsn/'2.6 B90204'/, nSrc/0/, nHead/0/, SanityChk/.true./,
+ doMags/.false./, useWCS/.true./, NaNwarn/.false./,
+ nn11,nn12,nn21,nn22/4*0/, w1m0,w2m0/2*22.5/, nPSF/2/,
+ NeedHelp/.False./, MskBitHist/32*0/, dbg/.false./,
+ notZero/0/, CoefMag/1.085736205d0/, w1mcor/0.145/,
+ w2mcor/0.177/, kBadw3,kBadw4,kBad2w3,kBad2w4/4*0/,
+ doCov/.true./, GotN1,GotN2/2*.false./, ReadAlready/.false./,
+ nIncomp/0/, nSrcHdr/-9/, nUnderPipe/0/, GotEndHdr1/.false./
c
Common / VDT / CDate, CTime, Vsn
c
namelist / abflagin / doCov, doMags, nPSF, w1m0, w1mcor,
+ w2m0, w2mcor
c
c=======================================================================
c
nArgs = IArgc()
NeedHelp = (nArgs .lt. 6)
1 If (NeedHelp) then
print *,'add-ab_flags vsn ', Vsn
print *
print *,'Usage: add-ab_flags <flags specifications>'
print *
print *,'Where the REQUIRED flags and specifications are:'
print *,
+ ' -i name of a gsa-matched mdex/CatWISE-IRSA-matched file'
print *,' -m name of a CatWISE mask FITS file'
print *,
+ ' -o name of the output CatWISE mdex file with ab_flags'
+ //' columns'
print *,' appended; this file must not already exist'
print *
print *,'The OPTIONAL flags are:'
print *,' -n1 name of a W1 "-n-" coverage image'
print *,' -n2 name of a W2 "-n-" coverage image'
print *,' -m0 turn off magnitude recomputation (default)'
print *,' -m1 turn on magnitude recomputation'
print *,' -n name of an abflagin namelist file'
print *,' -d (enable debug mode)'
Print *
print *,
+ 'If either "-n1" or "-n2" is specified, the other must also be.'
stop
end if
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
call signon('add-ab_flags')
c
NArg = 0
2 NArg = NArg + 1
call GetArg(NArg,Flag)
Flag0 = Flag
call UpCase(Flag)
c ! input CatWISE/IRSA file
If (Flag .eq. '-I') then
call NextNarg(NArg,Nargs)
Call GetArg(NArg,InFNam)
if (Access(InFNam(1:LNBlnk(InFNam)),' ') .ne. 0) then
print *
print *,'ERROR: file not found: ', InFNam(1:LNBlnk(InFNam))
print *
NeedHelp = .True.
Go to 1
end if
c ! Turn debug prints on
else if (Flag .eq. '-D') then
dbg = .true.
print *,'Debug prints enabled'
c ! Turn magnitude recomputation off
else if (Flag .eq. '-M0') then
DoMags = .false.
if (dbg) print *,'magnitude recomputation off'
c ! Turn magnitude recomputation on
else if (Flag .eq. '-M1') then
DoMags = .true.
if (dbg) print *,'magnitude recomputation on'
c
else if (Flag .eq. '-M') then
call NextNarg(NArg,Nargs)
Call GetArg(NArg,MskNam)
if (Index(MskNam,'.fits') .eq. 0)
+ MskNam = MskNam(1:LNBlnk(MskNam))//'.fits'
if (Access(MskNam(1:LNBlnk(MskNam)),' ') .ne. 0) then
print *
print *,'ERROR: File not found: ', MskNam(1:LNBlnk(MskNam))
print *
NeedHelp = .True.
Go to 1
end if
c
else if (Flag .eq. '-N1') then
call NextNarg(NArg,Nargs)
Call GetArg(NArg,Cov1Nam)
if (Index(Cov1Nam,'.fits') .eq. 0)
+ Cov1Nam = Cov1Nam(1:LNBlnk(Cov1Nam))//'.fits'
if (Access(Cov1Nam(1:LNBlnk(Cov1Nam)),' ') .ne. 0) then
print *
print *,'ERROR: File not found: ', Cov1Nam(1:LNBlnk(Cov1Nam))
print *
NeedHelp = .True.
Go to 1
end if
GotN1 = .true.
c
else if (Flag .eq. '-N2') then
call NextNarg(NArg,Nargs)
Call GetArg(NArg,Cov2Nam)
if (Index(Cov2Nam,'.fits') .eq. 0)
+ Cov2Nam = Cov2Nam(1:LNBlnk(Cov2Nam))//'.fits'
if (Access(Cov2Nam(1:LNBlnk(Cov2Nam)),' ') .ne. 0) then
print *
print *,'ERROR: File not found: ', Cov2Nam(1:LNBlnk(Cov2Nam))
print *
NeedHelp = .True.
Go to 1
end if
GotN2 = .true.
else if (Flag .eq. '-O') then
call NextNarg(NArg,Nargs)
Call GetArg(NArg,OutFNam)
if (Index(OutFNam,'.tbl') .eq. 0)
+ OutFNam = OutFNam(1:LNBlnk(OutFNam))//'.tbl'
if (Access(OutFNam(1:LNBlnk(OutFNam)),' ') .eq. 0) then
print *
print *,'ERROR: Output file already exists: ',
+ OutFNam(1:LNBlnk(OutFNam))
print *
NeedHelp = .True.
Go to 1
end if
else if (Flag .eq. '-N') then
call NextNarg(NArg,Nargs)
Call GetArg(NArg,NLNam)
if (Access(NLNam(1:LNBlnk(NLNam)),' ') .ne. 0) then
print *
print *,'ERROR: File not found: ', NLNam(1:LNBlnk(NLNam))
print *
NeedHelp = .True.
Go to 1
end if
open(10, file = NLnam)
read(10, abflagin, end = 3017, err = 3018)
write (6, abflagin)
close(10)
Else
print *,'ERROR: unrecognized command-line specification: '
+ //Flag0
end if
c
If (NArg .lt. NArgs) Go to 2
DoCov = DoCov .and. GotN1 .and. GotN2
if ((GotN1 .or. GotN2) .and. .not.DoCov) then
print *,
+ 'WARNING: "-n1" and/or "-n2" specified, but either DoCov = F'
print *,
+ 'in namelist or only one coverage image specified;'
print *,'coverage image(s) ignored.'
end if
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c ! open the mask file and read in
call readFhead(MskNam,Hfits)
Call GetNAX(MskNam,NCols,NRows,NPlanes,1,FileID1)
if (dbg) print *,'GetNAX returned NCols,NRows,NPlanes,FileID1:',
+ NCols,NRows,NPlanes,FileID1 ! dbg
if (NPlanes .ne. 1) then
print *,'ERROR: this program handles only 2-dimensional images'
stop
end if
c
wcs = -1
call wcsinit(Hfits,WCS)
c if (wcs .le. 0) then
c print *,'ERROR: file has no usable WCS - ',
c + MskNam(1:lnblnk(MskNam))
c useWCS = .false.
c end if
c
allocate(array1(NCols,NRows))
if (.not.allocated(array1)) then
print *,'ERROR: allocation of array1 failed for NAXIS1 = ',NCols,
+ ', NAXIS2 = ',NRows
call exit(64)
end if
c print *,'array1 allocated successfully' ! dbg
NPixpL = NCols*NRows
c
NPix = NPixpL*NPlanes
status = 0
call ftgpvj(FileID1,1,1,NPix,nullval,array1,anynull,status)
if (status .ne. 0) then
print *,'ERROR reading file1; status = ',status
stop
end if
c print *,'ftgpvj returned ',NPix,' pixels for file1' ! dbg
status = 0
c
C The FITS file must always be closed before exiting the program.
C Any unit numbers allocated with FTGIOU must be freed with FTFIOU.
call ftclos(FileID1, status)
call ftfiou(FileID1, status)
if (dbg) then
do 4 j = 1, 2048
do 3 i = 1, 2048
if (array1(i,j) .ne. 0) notZero = notZero + 1
3 continue
4 continue
print *,'No. of nonzero mask pixels: ',notZero,
+ '; fraction = ',float(notZero)/16793604.0
end if
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c ! read in the W1 coverage file
if (.not.DoCov) go to 5
call readFhead(Cov1Nam,Hfits)
Call GetNAX(Cov1Nam,NCols,NRows,NPlanes,1,FileID1)
if (dbg) print *,
+ 'Cov1 - GetNAX returned NCols,NRows,NPlanes,FileID1:',
+ NCols,NRows,NPlanes,FileID1 ! dbg
if (NPlanes .ne. 1) then
print *,'ERROR: this program handles only 2-dimensional images'
stop
end if
c
NPix = NCols*NRows
status = 0
call ftgpvi(FileID1,1,1,NPix,nullval,Cov1,anynull,status)
if (status .ne. 0) then
print *,'ERROR reading Cov1; status = ',status
stop
end if
status = 0
call ftclos(FileID1, status)
call ftfiou(FileID1, status)
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c ! read in the W2 coverage file
call readFhead(Cov2Nam,Hfits)
Call GetNAX(Cov2Nam,NCols,NRows,NPlanes,1,FileID1)
if (dbg) print *,
+ 'Cov2 - GetNAX returned NCols,NRows,NPlanes,FileID1:',
+ NCols,NRows,NPlanes,FileID1 ! dbg
if (NPlanes .ne. 1) then
print *,'ERROR: this program handles only 2-dimensional images'
stop
end if
c
NPix = NCols*NRows
status = 0
call ftgpvi(FileID1,1,1,NPix,nullval,Cov2,anynull,status)
if (status .ne. 0) then
print *,'ERROR reading Cov2; status = ',status
stop
end if
status = 0
call ftclos(FileID1, status)
call ftfiou(FileID1, status)
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c ! Sanity Check
5 open (10, file = InFNam)
6 read (10,'(a)', end = 3000) Line
if (Line(1:1) .eq. '\') go to 6
call GetFlds(Line,Field,IFa,IFb,NF)
HdrLine = Line
rewind(10)
c
if (SanityChk) then
if (NF .ne. 401) then
print *,'ERROR: input is not a two-band wphot mdex file'
print *,' no. of fields =', NF,'; should be 401'
call exit(64)
end if
c ! verify some fields in mdex file
call ChkFld(Field(3), 'ra',3)
call ChkFld(Field(4), 'dec',4)
call ChkFld(Field(5), 'sigra',5)
call ChkFld(Field(6), 'sigdec',6)
call ChkFld(Field(7), 'sigradec',7)
call ChkFld(Field(8), 'w1x',8)
call ChkFld(Field(9), 'w1y',9)
call ChkFld(Field(10),'w2x',10)
call ChkFld(Field(11),'w2y',11)
end if
c
open(20, file = OutFNam)
nw385 = IFb(385) - IFa(385) + 1 ! dist_x
nw394 = IFb(394) - IFa(394) + 1 ! cc_flags
nw395 = IFb(395) - IFa(395) + 1 ! w1cc_map
nw396 = IFb(396) - IFa(396) + 1 ! w1cc_map_str
nw397 = IFb(397) - IFa(397) + 1 ! w2cc_map
nw398 = IFb(398) - IFa(398) + 1 ! w2cc_map_str
if (dbg) then
print *,'dist_x field width:', nw385
print *,'cc_flags field width:', nw394
print *,'w1cc_map field width:', nw395
print *,'w1cc_map_str field width:', nw396
print *,'w2cc_map field width:', nw397
print *,'w2cc_map_str field width:', nw398
end if
c ! filter out unwanted header lines
10 read (10, '(a)', end=1000) Line
if (Line(1:1) .eq. '\') then
if (index(Line,'"Joined"') .gt. 0) GotEndHdr1 = .true.
if (.not.GotEndHdr1) write(20,'(a)') Line(1:lnblnk(Line))
if ((index(Line,'\Nsrc =') .gt. 0) .and. (nSrcHdr .lt. 0)) then
n = index(Line,'=') + 1
read (Line(n:lnblnk(Line)), *, err = 3002) nSrcHdr
if (dbg) print *,'Header \nSrc line: "'
+ //Line(1:lnblnk(Line))//'"'
end if
if (index(Line,'\DATETIME =') .gt. 0) then
n = index(Line,'"') + 1
Line = '\ AllWISE flags retrieved from IRSA using 2.75'
+ //' arcsec radius on '//Line(n:n+18)
write(20,'(a)') Line(1:lnblnk(Line))
end if
if (index(Line,'\ mrgad vsn') .gt. 0) GotEndHdr1 = .true.
go to 10
end if
c
if (Line(1:1) .eq. '|') then
Line = Line(1:Ifb(191))//Line(IFa(385):IFb(385))
+ //Line(IFa(394):IFb(398))//'|'
nHead = nHead + 1
if (nHead .eq. 1) then
write(20,'(a)') '\ add-ab_flags vsn '//Vsn//' run on '
+ //CDate//' at '//CTime
write(20,'(a)') '\ artifact bitmasks from '
+ //MskNam(1:LNBlnk(MskNam))
Line = Line(1:lnblnk(Line))
+ //'n_aw|ab_flags|w1ab_map|w1ab_map_str|w2ab_map|w2ab_map_str|'
write(20,'(a)') Line(1:lnblnk(Line))
HdrLine = Line
call GetFlds(HdrLine,Field2,IFa2,IFb2,NF2) ! for "under-pipe" check
end if
if (nHead .eq. 2) write(20,'(a)') '| char '
+ //'| i | double | double | r | r | r '
+ //'| r | r | r | r | r | r '
+ //'| r | r | r | r | r | r | r '
+ //'| r | r | r | r | r '
+ //'| r | r | r | r | r | r '
+ //'| r | i | i | c | c | r | r | char '
+ //'| r | r | i | r | r | r | r | i '
+ //'| r | r | r | r | i | r '
+ //'| r | i | r | r | i | r '
+ //'| r | i | r | r | i | r '
+ //'| r | i | r | r | i | r '
+ //'| r | i | r | r | i | r '
+ //'| r | i | r | r | i | r '
+ //'| r | i | r | r | i | r '
+ //'| r | i | r | r | i | r '
+ //'| r | i | i | i | r | r | r '
+ //'| r | i | r | double '
+ //'| double | double | i | i | r '
+ //'| r | r | r | i | r | double '
+ //'| double | double | i | i | r '
+ //'| r | r | r | r | r | i | r '
+ //'| r | r | r | r | i | i | i | i '
+ //'| r | r | double | double | double '
+ //'| r | r | r | r | r '
+ //'| r | r | r | r | r '
+ //'| r | r | r | r '
+ //'| r | r | r | r | r '
+ //'| r | c | i | i | real | r '
+ //'| r | r | r | double | r | double '
+ //'| r | r | r | r | r '
+ //'| r | r | r | r | i| i| i| i|c'
+ //'| r | r | r | r | double '
+ //'| char | int | char '
+ //'| int | char | int| char | int '
+ //'| char | int | char |'
if (nHead .eq. 3) write(20,'(a)') Line(1:lnblnk(Line))
+ //' | | | | | |'
if (nHead .eq. 4) write(20,'(a)') Line(1:lnblnk(Line))
+ //'null| null | null | null | null | null |'
go to 10
end if
c
12 nSrc = nSrc + 1
nNaN = index(Line,' NaN')
if (nNaN .gt. 0) NaNwarn = .true.
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c ! check for multi-matches
15 if ((index(Line(IFa(401):IFb(401)),'null') .gt. 0) .or.
+ (index(Line(IFa(401):IFb(401)),'-999') .gt. 0)) then ! not multi
if ((index(Line(IFa(385):IFb(385)),'null') .gt. 0) .or.
+ (index(Line(IFa(385):IFb(385)),'-999') .gt. 0)) then ! no match
nAWstr = ' 0'
dist_str = 'null'
cc_flags = 'null'
w1cc_map = 'null'
w1cc_map_str = 'null'
w2cc_map = 'null'
w2cc_map_str = 'null'
else ! solo match
nAWstr = ' 1'
dist_str = AdjustL(Line(IFa(385):IFb(385)))
cc_flags = AdjustL(Line(Ifa(394):IFb(394)))
w1cc_map = AdjustL(Line(Ifa(395):IFb(395)))
w1cc_map_str = AdjustL(Line(Ifa(396):IFb(396)))
w2cc_map = AdjustL(Line(Ifa(397):IFb(397)))
w2cc_map_str = AdjustL(Line(Ifa(398):IFb(398)))
end if
if (dbg .and. (nSrc .lt. 5)) then
print *,'nSrc:', nSrc
print *,
+ 'dist_str = |'//dist_str(1:lnblnk(dist_str))//'|'
print *,
+ 'cc_flags = |'//cc_flags(1:lnblnk(cc_flags))//'|'
print *,
+ 'w1cc_map = |'//w1cc_map(1:lnblnk(w1cc_map))//'|'
print *,
+ 'w1cc_map_str = |'//w1cc_map_str(1:lnblnk(w1cc_map_str))//'|'
print *,
+ 'w2cc_map = |'//w2cc_map(1:lnblnk(w2cc_map))//'|'
print *,
+ 'w2cc_map_str = |'//w2cc_map_str(1:lnblnk(w2cc_map_str))//'|'
end if
go to 100
end if
c ! process multi-matches
read (Line(IFa(2):IFb(2)) , *, err = 3003) src0
read (Line(IFa(401):IFb(401)), *, err = 3003) nAW
write (nAWstr,'(i5)') nAW ! note that w1cc_map and w2cc_map
read (Line(IFa(395):Ifb(395)), *, err = 3004) w1ccmap ! should never
read (Line(IFa(397):Ifb(397)), *, err = 3005) w2ccmap ! be null for
read (Line(IFa(385):Ifb(385)), *, err = 3006) dist ! matched srcs.
n = IFa(394) + lnblnk(Line(IFa(394):Ifb(394))) - 2 ! monitor W3
w3cc = Line(n:n) ! & W4 cc_flags
w4cc = Line(n+1:n+1)
kBadw3 = kBadness(w3cc)
kBadw4 = kBadness(w4cc)
c
do 20 n = 2, nAW
read (10, '(a)', end=3006) Line
read (Line(IFa(2):IFb(2)), *, err = 3003) src
if (src .ne. src0) then ! Inconsistent Group; missing
ReadAlready = .true. ! at least one member; assume
nIncomp = nIncomp + 1 ! members WOLD BE contiguous
go to 25 ! in input stream
end if
read (Line(IFa(395):Ifb(395)), *, err = 3004) w1ccmap2
read (Line(IFa(397):Ifb(397)), *, err = 3005) w2ccmap2
w1ccmap = IOr(w1ccmap, w1ccmap2)
w2ccmap = IOr(w2ccmap, w2ccmap2)
read (Line(IFa(385):Ifb(385)), *, err = 3006) dist2
if (dist2 .gt. dist) dist = dist2
kBad2w3 = 0
kBad2w4 = 0
k = IFa(394) + lnblnk(Line(IFa(394):Ifb(394))) - 2 ! monitor W3
w3cc2 = Line(k:k) ! & W4 cc_flags
w4cc2 = Line(k+1:k+1)
kBad2w3 = kBadness(w3cc2)
kBad2w4 = kBadness(w4cc2)
if (kBad2w3 .gt. kBadw3) then
w3cc = w3cc2
kBadw3 = kBad2w3
end if
if (kBad2w4 .gt. kBadw4) then
w4cc = w4cc2
kBadw4 = kBad2w4
end if
20 continue
25 write (dist_str,'(f8.5)') dist
30 if (lnblnk(dist_str) .lt. nw385) then
dist_str = ' '//dist_str
go to 30
end if
c
c| cc_flags | w1cc_map | w1cc_map_str | w2cc_map | w2cc_map_str |
c1234567890123456123456789012312345678901234567890123456789012312345678901234567890
cc________1______c________1___c________1_________2c________1___c________1_________2
c
c O H - P D - - o h - p d
c 11 10 09 08 07 06 05 04 03 02 01 00
c
c flag priority: D,P,H,O,d,p,h,o for cc_flags
c flag priority: D,d,P,p,H,h,O,o for w?cc_map_str
c
40 if (w1ccmap .eq. 0) then
cc_flags = '.0'
w1cc_map = '0'
w1cc_map_str = '.null'
else
if (dbg) then
print *,'-------------------------'
print *,'nSrc: ', nSrc,' ',trim(Line(IFa(1):IFb(1))),
+ ', w1ccmap: ', w1ccmap
end if
cc_flags = '.'
write(w1cc_map,'(i9)') w1ccmap
if (BitSet(w1ccmap,7)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'D'
else if (BitSet(w1ccmap,8)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'P'
else if (BitSet(w1ccmap,10)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'H'
else if (BitSet(w1ccmap,11)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'O'
else if (BitSet(w1ccmap,0)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'d'
else if (BitSet(w1ccmap,1)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'p'
else if (BitSet(w1ccmap,3)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'h'
else if (BitSet(w1ccmap,4)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'o'
end if
if (dbg) print *,'cc_flags: ',cc_flags(1:lnblnk(cc_flags))
w1cc_map_str = '.'
if (BitSet(w1ccmap,7))
+ w1cc_map_str = w1cc_map_str(1:lnblnk(w1cc_map_str))//'D'
if (BitSet(w1ccmap,0) .and. .not.BitSet(w1ccmap,7))
+ w1cc_map_str = w1cc_map_str(1:lnblnk(w1cc_map_str))//'d'
if (BitSet(w1ccmap,8))
+ w1cc_map_str = w1cc_map_str(1:lnblnk(w1cc_map_str))//'P'
if (BitSet(w1ccmap,1) .and. .not.BitSet(w1ccmap,8))
+ w1cc_map_str = w1cc_map_str(1:lnblnk(w1cc_map_str))//'p'
if (BitSet(w1ccmap,10))
+ w1cc_map_str = w1cc_map_str(1:lnblnk(w1cc_map_str))//'H'
if (BitSet(w1ccmap,3) .and. .not.BitSet(w1ccmap,10))
+ w1cc_map_str = w1cc_map_str(1:lnblnk(w1cc_map_str))//'h'
if (BitSet(w1ccmap,11))
+ w1cc_map_str = w1cc_map_str(1:lnblnk(w1cc_map_str))//'O'
if (BitSet(w1ccmap,4) .and. .not.BitSet(w1ccmap,11))
+ w1cc_map_str = w1cc_map_str(1:lnblnk(w1cc_map_str))//'o'
if (dbg) print *,'w1cc_map_str: ',w1cc_map_str(1:lnblnk(w1cc_map_str))
end if
c
if (w2ccmap .eq. 0) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'0'
w2cc_map = '0'
w2cc_map_str = '.null'
else
if (dbg) then
print *,'-------------------------'
print *,'nSrc: ', nSrc,' ',trim(Line(IFa(1):IFb(1))),
+ ', w2ccmap: ', w2ccmap
end if
write(w2cc_map,'(i9)') w2ccmap
if (BitSet(w2ccmap,7)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'D'
else if (BitSet(w2ccmap,8)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'P'
else if (BitSet(w2ccmap,10)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'H'
else if (BitSet(w2ccmap,11)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'O'
else if (BitSet(w2ccmap,0)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'d'
else if (BitSet(w2ccmap,1)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'p'
else if (BitSet(w2ccmap,3)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'h'
else if (BitSet(w2ccmap,4)) then
cc_flags = cc_flags(1:lnblnk(cc_flags))//'o'
end if
if (dbg) print *,'cc_flags: ',cc_flags(1:lnblnk(cc_flags))
w2cc_map_str = '.'
if (BitSet(w2ccmap,7))
+ w2cc_map_str = w2cc_map_str(1:lnblnk(w2cc_map_str))//'D'
if (BitSet(w2ccmap,0) .and. .not.BitSet(w2ccmap,7))
+ w2cc_map_str = w2cc_map_str(1:lnblnk(w2cc_map_str))//'d'
if (BitSet(w2ccmap,8))
+ w2cc_map_str = w2cc_map_str(1:lnblnk(w2cc_map_str))//'P'
if (BitSet(w2ccmap,1) .and. .not.BitSet(w2ccmap,8))
+ w2cc_map_str = w2cc_map_str(1:lnblnk(w2cc_map_str))//'p'
if (BitSet(w2ccmap,10))
+ w2cc_map_str = w2cc_map_str(1:lnblnk(w2cc_map_str))//'H'
if (BitSet(w2ccmap,3) .and. .not.BitSet(w2ccmap,10))
+ w2cc_map_str = w2cc_map_str(1:lnblnk(w2cc_map_str))//'h'
if (BitSet(w2ccmap,11))
+ w2cc_map_str = w2cc_map_str(1:lnblnk(w2cc_map_str))//'O'
if (BitSet(w2ccmap,4) .and. .not.BitSet(w2ccmap,11))
+ w2cc_map_str = w2cc_map_str(1:lnblnk(w2cc_map_str))//'o'
if (dbg) print *,'w2cc_map_str: ',w2cc_map_str(1:lnblnk(w2cc_map_str))
end if
c
cc_flags = cc_flags(1:lnblnk(cc_flags))//w3cc//w4cc
cc_flags(1:1) = ' ' ! remove sentinel character
w1cc_map_str(1:1) = ' '
w2cc_map_str(1:1) = ' '
c
100 if (lnblnk(cc_flags) .lt. nw394) then
cc_flags = ' '//cc_flags
go to 100
end if
110 if (lnblnk(w1cc_map_str) .lt. nw396) then
w1cc_map_str = ' '//w1cc_map_str
go to 110
end if
120 if (lnblnk(w2cc_map_str) .lt. nw398) then
w2cc_map_str = ' '//w2cc_map_str
go to 120
end if
130 if (lnblnk(w1cc_map) .lt. nw395) then
w1cc_map = ' '//w1cc_map
go to 130
end if
140 if (lnblnk(w2cc_map) .lt. nw397) then
w2cc_map = ' '//w2cc_map
go to 140
end if
150 if (lnblnk(dist_str) .lt. nw385) then
dist_str = ' '//dist_str
go to 150
end if
c
Line(IFa(385):IFb(385)) = dist_str
Line(Ifa(394):IFb(394)) = cc_flags
Line(Ifa(395):IFb(395)) = w1cc_map
Line(Ifa(396):IFb(396)) = w1cc_map_str
Line(Ifa(397):IFb(397)) = w2cc_map
Line(Ifa(398):IFb(398)) = w2cc_map_str
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c ! recompute mags if requested
if (doMags) then
c
if (index(Line(IFa(22):IFb(23)),'null') .eq. 0) then ! w1flux
nMag = 1
read (Line(Ifa(22):Ifb(22)), *, err=333) flux
read (Line(Ifa(23):Ifb(23)), *, err=333) sigflux
wsnr = flux/sigflux
if (wsnr .gt. 9999.0) wsnr = 9999.0
if (wsnr .lt. -999.0) wsnr = -999.0
write (Line(IFA(20):IFB(20)), '(F7.1)') wsnr
if ((flux .gt. 0.0d0) .and. (wsnr .ge. 2.0d0)) then
mag = w1m0 - 2.5*dlog10(flux)
sigmag = CoefMag*sigflux/flux
write (Line(Ifa(26):Ifb(26)),'(f7.3)') mag
write (Line(Ifa(27):Ifb(27)),'(f10.3)') sigmag
else
if (flux .gt. 0.0d0) then
mag = w1m0 - 2.5*dlog10(flux+2.0d0*sigflux)
else
mag = w1m0 - 2.5*dlog10(2.0d0*sigflux)
end if
write (Line(Ifa(26):Ifb(26)),'(f7.3)') mag
Line(Ifa(27):Ifb(27)) = ' 9.99'
end if
if (index(Line(IFa(26):IFb(26)),'NaN') .ne. 0) then
print *, 'ERROR: NaN produced for w1mpro on row no.', nSrc
print *,' w1flux, w1sigflux:', flux, sigflux
NaNwarn = .true.
nn11 = nn11 + 1
end if
end if
c
if (index(Line(IFa(24):IFb(25)),'null') .eq. 0) then ! w2flux
nMag = 2
read (Line(Ifa(24):Ifb(24)), *, err=333) flux
read (Line(Ifa(25):Ifb(25)), *, err=333) sigflux
wsnr = flux/sigflux
if (wsnr .gt. 9999.0) wsnr = 9999.0
if (wsnr .lt. -999.0) wsnr = -999.0
write (Line(IFA(21):IFB(21)), '(F7.1)') wsnr
if ((flux .gt. 0.0d0) .and. (wsnr .ge. 2.0d0)) then
mag = w2m0 - 2.5*dlog10(flux)
sigmag = CoefMag*sigflux/flux
write (Line(Ifa(29):Ifb(29)),'(f7.3)') mag
write (Line(Ifa(30):Ifb(30)),'(f10.3)') sigmag
else
if (flux .gt. 0.0d0) then
mag = w2m0 - 2.5*dlog10(flux+2.0d0*sigflux)
else
mag = w2m0 - 2.5*dlog10(2.0d0*sigflux)
end if
write (Line(Ifa(29):Ifb(29)),'(f7.3)') mag
Line(Ifa(30):Ifb(30)) = ' 9.99'
end if
if (index(Line(IFa(29):IFb(29)),'NaN') .ne. 0) then
print *, 'ERROR: NaN produced for w2mpro on row no.', nSrc
print *,' w2flux, w2sigflux:', flux, sigflux
NaNwarn = .true.
nn21 = nn21 + 1
end if
end if
c
if (index(Line(IFa(152):IFb(153)),'null') .eq. 0) then ! w1flux_pm
nMag = 3
read (Line(Ifa(152):Ifb(152)), *, err=333) flux
read (Line(Ifa(153):Ifb(153)), *, err=333) sigflux
wsnr = flux/sigflux
if (wsnr .gt. 9999.0) wsnr = 9999.0
if (wsnr .lt. -999.0) wsnr = -999.0
write (Line(IFA(150):IFB(150)), '(f9.1)') wsnr
if ((flux .gt. 0.0d0) .and. (wsnr .ge. 2.0d0)) then
mag = w1m0 - 2.5*dlog10(flux)
sigmag = CoefMag*sigflux/flux
write (Line(Ifa(156):Ifb(156)),'(f10.3)') mag
write (Line(Ifa(157):Ifb(157)),'(f13.3)') sigmag
else
if (flux .gt. 0.0d0) then
mag = w1m0 - 2.5*dlog10(flux+2.0d0*sigflux)
else
mag = w1m0 - 2.5*dlog10(2.0d0*sigflux)
end if
write (Line(Ifa(156):Ifb(156)),'(f10.3)') mag
Line(Ifa(157):Ifb(157)) = ' 9.99'
end if
if (index(Line(IFa(156):IFb(156)),'NaN') .ne. 0) then
print *,'ERROR: NaN produced for w1mpro_pm on row no.', nSrc
print *,' w1flux_pm, w1sigflux_pm:', flux, sigflux
NaNwarn = .true.
nn12 = nn12 + 1
end if
end if
c
if (index(Line(IFa(154):IFb(155)),'null') .eq. 0) then ! w2flux_pm
nMag = 4
read (Line(Ifa(154):Ifb(154)), *, err=333) flux
read (Line(Ifa(155):Ifb(155)), *, err=333) sigflux
wsnr = flux/sigflux
if (wsnr .gt. 9999.0) wsnr = 9999.0
if (wsnr .lt. -999.0) wsnr = -999.0
write (Line(IFA(151):IFB(151)), '(F9.1)') wsnr
if ((flux .gt. 0.0d0) .and. (wsnr .ge. 2.0d0)) then
mag = w2m0 - 2.5*dlog10(flux)
sigmag = CoefMag*sigflux/flux
write (Line(Ifa(159):Ifb(159)),'(f10.3)') mag
write (Line(Ifa(160):Ifb(160)),'(f13.3)') sigmag
else
if (flux .gt. 0.0d0) then
mag = w2m0 - 2.5*dlog10(flux+2.0d0*sigflux)
else
mag = w2m0 - 2.5*dlog10(2.0d0*sigflux)
end if
write (Line(Ifa(159):Ifb(159)),'(f10.3)') mag
Line(Ifa(160):Ifb(160)) = ' 9.99'
end if
if (index(Line(IFa(159):IFb(159)),'NaN') .ne. 0) then
print *,'ERROR: NaN produced for w2mpro_pm on row no.', nSrc
print *,' w2flux_pm, w2sigflux_pm:', flux, sigflux
NaNwarn = .true.
nn22 = nn22 + 1
end if
end if
c
end if
c
go to 400
c
333 print *,'ERROR: read error on flux information group #',
+ nMag,', data row #',nSrc
print *,
+ ' magnitude recomputation terminated for this source'
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c ! process CatWISE mask bits
400 if (useWCS) then
read(Line(IFa(3):Ifb(3)), *, err = 3008) ra
read(Line(IFa(4):Ifb(4)), *, err = 3009) dec
offscl = -1
call wcs2pix(wcs, ra, dec, x8, y8, offscl)
if (offscl .ne. 0) then
w1x = 1.
w1y = 1.
else
w1x = x8
w1y = y8
endif
c iPix = NInt(w1x) ! don't round off; RA & Dec of
c jPix = NInt(w1y) ! a pixel are at its center
iPix = w1x
jPix = w1y
write(Line(IFa(8):IFb(8)), '(f9.3)') w1x
write(Line(IFa(9):IFb(9)), '(f9.3)') w1y
write(Line(IFa(10):IFb(10)),'(f9.3)') w1x
write(Line(IFa(11):IFb(11)),'(f9.3)') w1y
if (dbg .and. (nSrc .le. 10)) then
print *,'-------------------------'
print *,'nSrc: ', nSrc,' ',trim(Line(IFa(1):IFb(1))),
+ ', ra & dec: ', ra, dec
print *,' w1x & w1y: ', w1x, w1y
end if
go to 500
end if
c
GoodXY1 = index(Line(IFA(8):IFB(8)), 'null') .eq. 0
GoodXY2 = index(Line(IFA(10):IFB(10)),'null') .eq. 0
if (GoodXY1 .and. GoodXY2) then
k = 8
read (Line(IFa(k):IFb(k)), *, err = 3001) w1x
k = 10
read (Line(IFa(k):IFb(k)), *, err = 3001) w2x
iPix = (w1x+w2x)/2.0
else if (GoodXY1) then
k = 8
read (Line(IFa(k):IFb(k)), *, err = 3001) w1x
iPix = w1x
else if (GoodXY2) then
k = 10
read (Line(IFa(k):IFb(k)), *, err = 3001) w2x
iPix = w2x
else
ab_flags = ' null'
w1ab_map = ' null'
w1ab_map_str = ' null'
w2ab_map = ' null'
w2ab_map_str = ' null'
go to 700
end if
c
GoodXY1 = index(Line(IFA(9):IFB(9)), 'null') .eq. 0
GoodXY2 = index(Line(IFA(11):IFB(11)),'null') .eq. 0
if (GoodXY1 .and. GoodXY2) then
k = 9
read (Line(IFa(k):IFb(k)), *, err = 3001) w1y
k = 11
read (Line(IFa(k):IFb(k)), *, err = 3001) w2y
jPix = (w1y+w2y)/2.0
else if (GoodXY1) then
k = 9
read (Line(IFa(k):IFb(k)), *, err = 3001) w1y
jPix = w1y
else if (GoodXY2) then
k = 11
read (Line(IFa(k):IFb(k)), *, err = 3001) w2y
jPix = w2y
else
ab_flags = ' null'
w1ab_map = ' null'
w1ab_map_str = ' null'
w2ab_map = ' null'
w2ab_map_str = ' null'
go to 700
end if
c
500 msk = array1(iPix,jPix)
if (dbg .and. (msk .gt. 0)) then
if (.not.BitSet(msk,6) .and. .not.BitSet(msk,21) .and.
+ .not.BitSet(msk,9) .and. .not.BitSet(msk,10) .and.
+ .not.BitSet(msk,22)) then
print *,'-------------------------'
print *,'nSrc: ', nSrc,' ',trim(Line(IFa(1):IFb(1))),
+ ', msk: ', msk
end if
end if
call DoHist(msk,MskBitHist)
c
if (BitSet(msk,21)) then ! don't flag the source
call ClrBit(msk,0) ! that caused the
call ClrBit(msk,1) ! diffraction spikes
call ClrBit(msk,7)
call ClrBit(msk,21)
call ClrBit(msk,23)
call ClrBit(msk,27)
call ClrBit(msk,29)
end if
if (BitSet(msk,22)) then
call ClrBit(msk,2)
call ClrBit(msk,3)
call ClrBit(msk,8)
call ClrBit(msk,22)
call ClrBit(msk,24)
call ClrBit(msk,28)
call ClrBit(msk,30)
end if
c
w1abmap = 0 ! O H - P D - - o h - p d
w2abmap = 0 ! 11 10 09 08 07 06 05 04 03 02 01 00
c
if (BitSet(msk,17) .or. BitSet(msk,18))
+ call SetBit(w1abmap,1) ! 'p' for W1
if (BitSet(msk,19) .or. BitSet(msk,20))
+ call SetBit(w2abmap,1) ! 'p' for W2
if (BitSet(msk,13) .or. BitSet(msk,14))
+ call SetBit(w1abmap,1) ! 'p' for W1
if (BitSet(msk,15) .or. BitSet(msk,16))
+ call SetBit(w2abmap,1) ! 'p' for W2
if (BitSet(msk,25) .or. BitSet(msk,26))
+ call SetBit(w1abmap,4) ! 'o' for W1
if (BitSet(msk,11) .or. BitSet(msk,12))
+ call SetBit(w2abmap,4) ! 'o' for W2
if (BitSet(msk, 0) .or. BitSet(msk, 1) .or. BitSet(msk, 29)
+ .or. BitSet(msk, 7) .or. BitSet(msk, 27))
+ call SetBit(w1abmap,0) ! 'd' for W1
if (BitSet(msk, 2) .or. BitSet(msk, 3) .or. BitSet(msk, 30)
+ .or. BitSet(msk, 8) .or. BitSet(msk, 28))
+ call SetBit(w2abmap,0) ! 'd' for W2
if (BitSet(msk, 23))
+ call SetBit(w1abmap,3) ! 'h' for W1
if (BitSet(msk, 24))
+ call SetBit(w2abmap,3) ! 'h' for W2
c
c flag priority: D,d,P,p,H,h,O,o
c
if (w1abmap .eq. 0) then
ab_flags = '.0 '
w1ab_map = ' 0'
w1ab_map_str = '.null '
else
if (dbg) print *,'nSrc: ', nSrc,' ',trim(Line(IFa(1):IFb(1))),
+ ', w1abmap: ', w1abmap
ab_flags = '. '
write(w1ab_map,'(i9)') w1abmap
if (BitSet(w1abmap,7)) then
ab_flags = ab_flags(1:lnblnk(ab_flags))//'D'
else if (BitSet(w1abmap,8)) then
ab_flags = ab_flags(1:lnblnk(ab_flags))//'P'
else if (BitSet(w1abmap,10)) then
ab_flags = ab_flags(1:lnblnk(ab_flags))//'H'
else if (BitSet(w1abmap,11)) then
ab_flags = ab_flags(1:lnblnk(ab_flags))//'O'
else if (BitSet(w1abmap,0)) then