FTN77,S
$FILES 3,3
      PROGRAM ANTPC(3,99)
***********************************************************************
*
*       PROGRAM TO READ ANALYSE DATA COLLECTED BY TNFPC.
*         INFORMATION WRITTEN ON BY SUBROUTINE WINFO IN PROGRAM TNFPC
*
*        This program writes data to final files
*
*            THIS IS COPY OF NWNTS AS ON 1984/11/05
*               ie can use other calibrators and get old nos from file
*                UPDATED <910118.0903>
*
*       DOES NOT PRESERVE DVM READINGS.
*
***********************************************************************
*
      DOUBLE PRECISION DATE, STOBS
*
      LOGICAL  AOK, PROB, RPTFQ, HKEEP, ENDFIL, RPRINT, DUALBM,
     *         FLXPNT, WVTB, CALSRS, WRTFNL, DONESX, DECDAT, DATF10,
     *         FSTIM
* 
      INTEGER  RESFIL(3), IDCB(144), RESLT(21), SORSNM(4), IPRAM(5),
     *         ITIME(5), TMDCB(144), TMFILE(3), IBLNK(1), NAMFIN(3) 
* 
      REAL     RDFLX(11), ERFLX(11), PEAK(3), PEAKER(3),
     *         DECOF(3), DECER(3) 
* 
      COMMON / FLUXS / RDFLX, ERFLX, WAVE, SORSNM 
* 
      DATA  JREC / 0 /, JS / 0 /, LU / 1 /, 
     *      TMFILE / 2HTM, 2HFI, 2HLE / 
* 
*       In/output unit is unit 1
           CALL FSYSU (LU,LU) 
* 
         CALL EXEC (11,ITIME,IYEAR) 
           WRITE (LU,1010) ITIME(5),IYEAR 
 1010      FORMAT(//" PROGRAM ANTPC. Day ",I3,":",I5/ 
     *            1X, 28"="/, 
     *              3X,"Program to analyse pointing data."/)
* 
           WRITE (LU,1000)
1000      FORMAT(/" Program assumes that no more than 150 sets of "/
     *              " readings were taken."/
     *              " Check at end of TNFPC tracking program that"/ 
     *              " this number is not exceeded."/
     *              "If it is, increase all  200-dimensional arrays"/ 
     *              " in this program and recompile.")
* 
  105      WRITE (LU,1100)
 1100      FORMAT(/3X," Name of results file?_")
             READ (LU,1200) (RESFIL(L), L=1,3)
 1200        FORMAT(3A2)
* 
           WRITE (LU,1300) (RESFIL(L), L=1,3) 
 1300      FORMAT(/3X," Reading results from disc file ",3A2, 
     *            /42X," OK (T)?_") 
* 
              READ (LU,1500) AOK
 1500         FORMAT(L6)
        IF ( .NOT. AOK )  GO TO 105 
* 
C              PRINT '(4X,"Date on file in decimal form (T or F) ? _")' 
C              READ *, DECDAT 
               DECDAT = .TRUE.
C              PRINT '(4X,"Data after D313/84 (T or F) ? _")' 
C              READ *, DATF10 
               DATF10 = .TRUE.
* 
C              PRINT '(4X,"Data after D325/84 (T or F) ? _")' 
C              READ *, FSTIM
               FSTIM = .TRUE. 
* 
C              WRITE (LU,1350)
C1350          FORMAT(/3X," Necessary to print out raw data ? Yes(T)_") 
C                READ (LU,1500) RPRINT
                 RPRINT = .FALSE. 
* 
C            WRITE (LU,1360)
C1360        FORMAT(/3X" Must NET FLUXES etc. be written out? Yes(T)_") 
C              READ (LU,1500) FLXPNT
               FLXPNT = .FALSE. 
* 
         WRITE (LU,1370)
 1370    FORMAT(/3X," Writing results to final file ? Yes(T)_") 
           READ (LU,1500) WRTFNL
* 
*     Open file, RESFIL.
* 
  110    CALL OPEN (IDCB, IERR, RESFIL) 
            IF ( IERR .GE. 0)  GO TO 130
* 
*        Error checking routine.
* 
  120    CALL CLOSE (IDCB)
            WRITE (LU,1400)  IERR 
 1400       FORMAT(" FMGR ERROR ",I4," in MAIN. CHECK STATUS. OK(T)? _")
                READ (LU,1500)  AOK 
                  IF ( .NOT. AOK)   GO TO 900 
            GO TO 110 
* 
*        Position pointer at begining of file.
* 
  130    CALL POSNT (IDCB, IERR, 1, 1)
             IF ( IERR .LT. 0)  GO TO 120 
* 
*        Read housekeeping, first resetting buffer. 
* 
  135          DO 133 IB = 1,21 
                   RESLT(IB) = 2H   
  133          CONTINUE 
* 
         CALL READF (IDCB, IERR, RESLT, 21, LEN)
* 
       JS = JS + 1
       DONESX = .FALSE. 
* 
           IF (( LEN .LE. 0 ) .OR. ( IERR .EQ. -012))  GO TO 150
              IF (IERR .LT. 0)  GO TO 120 
         GO TO 160
* 
  150      WRITE (LU,1600)
 1600      FORMAT(/," END OF FILE ")
             GO TO 500
* 
*      Assign to program parameters.
* 
  160      CALL CODE
                READ (RESLT,1740) HKEEP 
 1740           FORMAT(L1)
* 
*       If line is line of housekeeping, continue.
               IF ( HKEEP )  GO TO 1750 
* 
*            Otherwise move back and write line of data.
                 CALL POSNT (IDCB, IERR, -1)
                    IF ( IERR .LT. 0 )  GO TO 120 
                    GO TO 1850
* 
*       If line is hkeeping, assign to program parameters 
* 
 1750     If (DECDAT) then
              If (DATF10) then
                  If (FSTIM) then 
                      CALL CODE 
                      READ (RESLT,1706) HKEEP, NSORS, IBLNK, SORSNM,
     *                                  DATE, KOUNT, CALSRS, SIDTIM 
                  else
                      CALL CODE 
                      READ (RESLT,1705) HKEEP, NSORS, IBLNK, SORSNM,
     *                                  DATE, KOUNT, CALSRS 
                  endif 
              else
                   CALL CODE
                   READ (RESLT, 1700) HKEEP, NSORS,IBLNK,SORSNM,
     *                                DATE, KOUNT, CALSRS 
              endif 
 1700      FORMAT(L1, I3, A2,4A2, F8.4, I2, L2) 
 1705      FORMAT(L1, I3, A2, 4A2, F10.6, I2, L2) 
 1706      FORMAT(L1, I3, A2, 4A2, F10.6, I2, L2, F12.6)
* 
                 IYEAR = IDINT(DATE)  
                 IDAY = IDINT ((DATE - DBLE(IYEAR)) * 365.) 
                 IYEAR = IYEAR + 1900 
* 
           else 
             CALL CODE
             READ (RESLT, 1710) HKEEP, NSORS, 
     *                        IBLNK, SORSNM, IDAY, IYEAR, KOUNT, CALSRS 
 1710        FORMAT(L1, I3, A2, 4A2, I4, I4, I2, L2)
                 DATE = FLOAT(IYEAR-1900) + FLOAT(IDAY)/365.
* 
           endif
* 
C             CALL NAMNO (NSORS, SORSNM, RAHR, *900,k)
* 
*       Find Hour Angle 
* 
                      DATSG = SNGL(DATE)
                      STOBS = DBLE(SIDTIM)
          CALL HANGL (SORSNM, NSORS, DATE, HA, DATF10, FSTIM, 
     *                STOBS, *900,k)
* 
*         Want greater time resolution for 1144-379 (ie. >1 day)
* 
           If (.NOT. DECDAT) then 
                If (NSORS .EQ. 123 .OR. NSORS .EQ. 916) then
                       PRINT '("Source",I4,2X,4A2,". Enter UT of",
     *                          " observation in hh, mm : _")', 
     *                                       NSORS, SORSNM
                       READ *, IHR, IMIN
                       DTIM = (FLOAT(IHR) + FLOAT(IMIN)/60.)/(24.*365.) 
                       DATE = DATE + DTIM 
                endif 
           endif
* 
*     1144-37 and 6 calibrators/controls common to both lists 
* 
              IF (NSORS .EQ. 916) NSORS = 123 
              IF (NSORS .EQ. 944) NSORS = 144 
              IF (NSORS .EQ. 909) NSORS = 108 
              IF (NSORS .EQ. 915) NSORS = 119 
              IF (NSORS .EQ. 921) NSORS = 127 
              IF (NSORS .EQ. 929) NSORS = 132 
              IF (NSORS .EQ. 937) NSORS = 141 
* 
*       Write headers.
          IF ( .NOT. RPRINT )  GO TO 190
              WRITE (LU,1900) 
 1900         FORMAT(//," No.",1X,"Source",3X,"Day:",6X,"No. freq.")
* 
       WRITE (LU,1800) NSORS, IBLNK, SORSNM, IDAY, IYEAR, KOUNT 
 1800  FORMAT(I3, A2, 4A2, I4, I5, I4)
* 
*        Write headers. 
 1850      IF ( .NOT. RPRINT )  GO TO 190 
           WRITE (LU,2100)
 2100      FORMAT(/6X,"Freq",1X,"Place",7X,"Flux",19X,"DVM",/ 
     *             6X," cm.",9X   , 5X,"Jy  ",12X,"Tube     Source",
     *            /6X,4"=",1X,5"=",1X,17"=",5X,"=======   =======") 
* 
*     Start of loop for each freq.
* 
  190   DO 300  I = 1,KOUNT 
* 
*            Start point for repeat of freq.
* 
*      READ DATA
* 
  170     DO 310 J = 1,11 
* 
             DUALBM = ((WAVE .EQ. 3.4) .OR. (WAVE .EQ. 6.0))
*              Only 5 readings taken at 13 and 18 single beam.
             IF (( J .LT. 6 ) .OR. DUALBM )  GO TO 320
* 
                   DO 325  JJ = 6,11
                       RDFLX(JJ) = 0.0
                       ERFLX(JJ) = 0.0
  325              CONTINUE 
               GO TO 330
* 
  320       CALL RESRD (WAVE, RDFLX(J),ERFLX(J),IDCB, RESLT,
     *                  LU, PROB, RPTFQ, HKEEP, ENDFIL, RESFIL, RPRINT, 
     *                  WVTB, IFQ, JS, TUBE)  
* 
                 IF ( PROB )  GO TO 900 
                 IF (ENDFIL)  GO TO 500 
* 
*               If line was housekeeping, set read flux values to zero. 
                 IF ( HKEEP )   then
                    RDFLX(J) = 0.0
                    ERFLX(J) = 0.0
                    WAVE = 0.0
*                       If reading not complete, ignore 
                     IF ((DUALBM .AND. (J .GT. 1 .AND. J .LT. 11)) .OR. 
     *                   (.NOT. DUALBM .AND. (J .GT. 1 .AND. J .LT. 5)))
     *                JS = JS - 1 
                      GO TO 135 
                  else
                 endif
* 
  310       CONTINUE
* 
              IF ( DONESX .AND. WAVE .EQ. 6.0)  JS = JS + 1 
* 
*        Calculate net fluxes by subtracting off FN values. 
* 
  330       CALL NETTS (LU, FLXPNT, PEAK, PEAKER, AVPK, DECOF, JS,
     *                  RPTFQ, RAOFF, RAFER, PKMEAS, PKMSER, DECER, 
     *                  DATSG, TUBE, IFQ) 
* 
             If ( WAVE .EQ. 6.0 )  DONESX = .TRUE.
* 
*          Write info onto temporary file for later manipulation. 
* 
               CALL TEMPF (LU, NSORS, SORSNM, WAVE, CALSRS, RPTFQ,
     *                     PEAK, PEAKER, AVPK, DECOF, DATSG,
     *                     PROB, TMDCB, TMFILE, JREC, JS, PKMEAS, 
     *                     PKMSER, DECER, RAOFF, RAFER, HA, TUBE) 
                     IF (PROB)  GO TO 900 
* 
            IF (RPTFQ)  GO TO 170 
* 
  300  CONTINUE 
* 
*        Return for next source.
         GO TO 135
* 
  500      CALL CLOSE (IDCB)
* 
*        Correct for RA offset and calbrate.
* 
       CALL CALIB (LU, TMDCB, TMFILE, WRTFNL, RESFIL, 
     *             *900,k)
* 
*        Write data to mag tape 
* 
               PRINT '(//4X,"Write data on disc file ",(3A2),"to ", 
     *                      "mag tape ? Yes(T) _")', (RESFIL(L),L=1,3)
                     READ *, AOK
                         IF (AOK)  CALL MAGTP (RESFIL, *900,k)
* 
  900   CONTINUE
* 
C          WRITE (LU,9100)
C9100      FORMAT(" Want to purge TEMPORARY FILE ? _")
C            READ (LU,1500) AOK 
C              IF ( .NOT. AOK )  GO TO 910
* 
C          CALL PURGE (TMDCB, IERR, TMFILE, 0, 22)
C              IF ( IERR .LT. 0 ) GO TO 120 
* 
             WRITE (LU,9999)
 9999        FORMAT(///"*** END ANTPC ***"///)
      END 
* 
*      SUBROUTINE TO READ RESULTS OFF DISC FILE RESFIL. 
* 
*********************************************************************** 
* 
      SUBROUTINE RESRD (WAVE, SORSJY, ERROR, IDCB, RESLT, LU, 
     *                  PROB, RPTFQ, HKEEP, ENDFIL, RESFIL, RPRINT, 
     *                  WVTB, IFQ, JS, TUBE)  
* 
        LOGICAL PROB, AOK, RPTFQ, HKEEP, ENDFIL, RPRINT, WVTB, CALCTB 
* 
        INTEGER IDCB(1), RESLT(1), RESFIL(1), KPLACE(1) 
* 
          PROB = .FALSE.
          ENDFIL = .FALSE.
          CALCTB = .FALSE.
* 
*    Read off file, first resetting buffer. 
* 
  100         DO  99  IB = 1,21 
                RESLT(IB) = 2H  
   99         CONTINUE
* 
         CALL READF (IDCB, IERR, RESLT, 21, LEN)
             IF (( LEN .LE. 0 ) .OR. (IERR .EQ. -012))  GO TO 110 
              IF (IERR .LT. 0)  GO TO 120 
          GO TO 200 
* 
*        Error routine. 
  110            WRITE (LU,1100)
 1100            FORMAT(" END OF FILE RESFIL ") 
                   ENDFIL = .TRUE.
                   GO TO 900
* 
  120      CALL CLOSE (IDCB)
           WRITE (LU,1200) IERR 
 1200      FORMAT(" FMGR ERROR ",I4," in RESRD. CHECK STATUS. OK(T)? _")
* 
                 READ (LU,1300) AOK 
 1300            FORMAT(L6) 
* 
                   IF ( AOK )  GO TO 180
* 
                        PROB = .TRUE. 
                         GO TO 900
* 
*           File error corrected, so....
  180         CALL OPEN (IDCB, IERR, RESFIL)
                  IF (IERR .LT. 0)  GO TO 120 
* 
*        Assign read values to parameters.
* 
  200       CALL CODE 
                 READ (RESLT,1430)  HKEEP 
 1430            FORMAT(L1) 
* 
*      If hkeeping, move back and leave subroutine
* 
              IF (HKEEP)  GO TO 250 
* 
         CALL CODE
      READ (RESLT,1440) HKEEP, WVTB 
 1440 FORMAT(2L1) 
* 
            IF ( .NOT. HKEEP )  GO TO 300 
* 
  250          CALL POSNT (IDCB, IERR, -1)
                  IF ( IERR .LT. 0 )  GO TO 120 
                  GO TO 900 
* 
*        Read off wavelength and tube DVM values. 
* 
  300        IF ( .NOT. WVTB )  GO TO 310 
C              CALL POSNT (IDCB, IERR, -1)
C                    IF (IERR .LT. 0)  GO TO 120
C              CALL READF (IDCB, IERR, RESLT, 21, LEN)
             CALL CODE
                  READ (RESLT, 1450)  HKEEP, WVTB, WAVE, AVTUB
 1450             FORMAT(L1, L1, F6.1, G12.6) 
* 
               CALCTB = .TRUE.
* 
         GO TO 100
* 
  310    CALL CODE
         READ (RESLT, 1400) HKEEP, WVTB, KPLACE, SORSJY, ERROR, 
     *                      RPTFQ, AVDVM
 1400    FORMAT(L1, L1, 1A2, G12.6, G12.6, L2, G12.6) 
* 
              If (CALCTB) then
               If (WAVE .EQ. 3.4) IFQ = 1 
               If (WAVE .EQ. 6.0) IFQ = 2 
               If (WAVE .EQ. 13.) IFQ = 3 
               If (WAVE .EQ. 18.) IFQ = 4 
               TUBE = AVTUB - AVDVM 
              endif 
* 
C        GO TO 100
* 
               IF ( .NOT. RPRINT )  GO TO 900 
           WRITE (LU,1500) WAVE, KPLACE, SORSJY, ERROR, AVTUB, AVDVM
 1500      FORMAT(4X,F6.1,2X, 1A2,2X, G10.4,"+-", 3G10.4) 
* 
  900 RETURN
       END
* 
*        SUBROUTINE TO COMPUTE NET FLUXES AND OTHER PEARLS OF WISDOM. 
**********************************************************************
* 
         SUBROUTINE NETTS (LU, FLXPNT, PEAK, PEAKER, AVPK, DECOF, JS, 
     *                     RPTFQ, RAOFF, RAFER, PKAV, PKAVER, DECER,
     *                     DATE, TUBE, IFQ) 
* 
           LOGICAL DUALBM, FLXPNT, RPTFQ
* 
           INTEGER TPLACE(7), SORSNM(4) 
* 
           REAL    AVFN(3), AVFNER(3), NETFLX(11), NETER(11), 
     *             RDFLX(11), ERFLX(11), PEAK(3), PEAKER(3), DECOF(3),
     *             HPAV(2), HPAVER(2), HPRAT(2),
     *             DECER(3), TUFLUX(4), ERFNDV(4), DE(4), JYDVM,
     *             NETDVM, NTERDV 
* 
           COMMON / FLUXS / RDFLX, ERFLX, WAVE, SORSNM
* 
           DATA  TPLACE / 2HON, 2HHP, 2HHP, 2HON, 2HHP, 2HHP, 2HCE /, 
     *           TUFLUX / 175.,   24.,   35.,  140.0 /
* 
             If (DATE .LT. 83.6466) then
                 TUFLUX(1) = 167. 
                 TUFLUX(2) = 22.8 
                 TUFLUX(3) = 32.8 
             endif
* 
             DUALBM = (( WAVE .EQ. 3.4) .OR. (WAVE .EQ. 6.0)) 
* 
*            Average of first two FN's. 
               AVFN(1) = (RDFLX(1) + RDFLX(5))/2. 
               ERROR = 0.5 * SQRT( (ERFLX(1)*ERFLX(1))
     *                            + (ERFLX(5)*ERFLX(5)) ) 
               AVFNER(1) = ERROR
* 
*            Net flux of first peak.
               NETFLX(1) = RDFLX(3) - AVFN(1) 
               ERROR = SQRT( (ERFLX(3)*ERFLX(3))
     *                      + (AVFNER(1)*AVFNER(1)) ) 
               NETER(1) = ERROR 
* 
*            Net HP fluxes. 
               NETFLX(2) = RDFLX(2) - AVFN(1) 
               ERROR = SQRT( (ERFLX(2)*ERFLX(2))
     *                      + (AVFNER(1)*AVFNER(1)) ) 
               NETER(2) = ERROR 
* 
               NETFLX(3) = RDFLX(4) - AVFN(1) 
               ERROR = SQRT( (ERFLX(4)*ERFLX(4))
     *                      + (AVFNER(1)*AVFNER(1)) ) 
               NETER(3) = ERROR 
* 
*            Ave. of centre FN's. 
               AVFN(2) = (RDFLX(5) + RDFLX(7))/2. 
               ERROR = 0.5 * SQRT( (ERFLX(5)*ERFLX(5))
     *                            + (ERFLX(7)*ERFLX(7)) ) 
               AVFNER(2) = ERROR
* 
*            Net centre flux
               NETFLX(7) = RDFLX(6) - AVFN(2) 
               ERROR = SQRT( (ERFLX(6)*ERFLX(6))
     *                      + (AVFNER(2)*AVFNER(2)) ) 
               NETER(7) = ERROR 
* 
           IF ( DUALBM )  GO TO 430 
                     NETFLX(7) = 0.0
                     NETER(7) = 0.0 
* 
* 
*            Ave. of last FN's. 
  430          AVFN(3) = (RDFLX(7) + RDFLX(11))/2.
               ERROR = 0.5 * SQRT( (ERFLX(7)*ERFLX(7))
     *                            + (ERFLX(11)*ERFLX(11)) ) 
               AVFNER(3) = ERROR
* 
*            Net peak flux. 
               NETFLX(4) = RDFLX(9) - AVFN(3) 
               ERROR = SQRT( (ERFLX(9)*ERFLX(9))
     *                      + (AVFNER(3)*AVFNER(3)) ) 
               NETER(4) = ERROR 
* 
*            Net HP's.
               NETFLX(5) = RDFLX(8) - AVFN(3) 
               ERROR = SQRT( (ERFLX(8)*ERFLX(8))
     *                      + (AVFNER(3)*AVFNER(3)) ) 
               NETER(5) = ERROR 
* 
               NETFLX(6) = RDFLX(10) - AVFN(3)
               ERROR = SQRT( (ERFLX(10)*ERFLX(10))
     *                      + (AVFNER(3)*AVFNER(3)) ) 
               NETER(6) = ERROR 
* 
           IF ( .NOT. DUALBM)  GO TO 400
*            Find average of all FN's.
               TAVFN = (RDFLX(1) + RDFLX(5) + RDFLX(7)
     *                   + RDFLX(11))/4.
               SMFNSQ = RDFLX(1)*RDFLX(1) + RDFLX(5)*
     *                 RDFLX(5) + RDFLX(7)*RDFLX(7)
     *                 + RDFLX(11)*RDFLX(11)
               AVSQFN = SMFNSQ/4.
*
                 GO TO 420
*
*            Averages for single beam systems-
  400          TAVFN = AVFN(1)
               SMFNSQ = RDFLX(1)*RDFLX(1) + RDFLX(5)
     *                        *RDFLX(5)
               AVSQFN = SMFNSQ/2.
*
  420        DEVFN = SQRT ( AVSQFN - TAVFN*TAVFN )
*
*        Calculate the correct errors from the old
*
            JYDVM = TUFLUX(IFQ) / ABS(TUBE)
* 
            DO 425 IE = 1,4 
                   If (IE .EQ. 1) IEE = 1 
                   If (IE .EQ. 2) IEE = 5 
                   If (IE .EQ. 3) IEE = 7 
                   If (IE .EQ. 4) IEE = 11
            ERFNDV(IE) = ERFLX(IEE) / (JYDVM * SQRT(3.))
            CALL MEDIN (ERFNDV(IE), IE, EMED, DE) 
  425       CONTINUE
* 
            ERTUBE = (DE(2) + DE(3)) / 2. 
* 
            DO 426 IN = 1,7 
* 
               If (.NOT. DUALBM .AND. IN .GE. 4)  GO TO 426 
               NETDVM = NETFLX(IN) / JYDVM
               NTERDV = NETER(IN) / JYDVM 
* 
               CDIFF = (NTERDV * NTERDV) +  
     *                   (ERFNDV(1) * ERFNDV(1) 
     *                              * (NETDVM / TUBE) 
     *                              * (NETDVM / TUBE)) -
     *                      (3. * ERTUBE * ERTUBE / 2.) 
               If (CDIFF .GT. 0.) then
               CNETER = SQRT(CDIFF) 
               NETER(IN) = CNETER * JYDVM 
               endif
  426       CONTINUE
* 
*        Find average of peak fluxes. 
           PKAV = ( NETFLX(1) - NETFLX(4))/2. 
           PKAVER = 0.5 * SQRT( (NETER(1)*NETER(1)) 
     *                         + (NETER(4)*NETER(4)) )
* 
              IF ( .NOT. DUALBM)  PKAV = PKAV * 2.
* 
*         Std. dev. of FN fluxes expressed as % of peak flux. 
           PCNTDV = ( DEVFN/PKAV ) * 100.0
* 
*            HP averages. 
               DO   440  L1= 1,2
                         L2 = 3*L1 - 1
                         L3 = L2 + 1
                         L4 = L1 * L1 
                   HPAV(L1) = ( NETFLX(L2) + NETFLX(L3))/2. 
                   ERROR = 0.5 * SQRT( (NETER(L2)*NETER(L2))
     *                                  + (NETER(L3)*NETER(L3)) ) 
                   HPAVER(L1) = ERROR 
* 
*                Ratios.
                     HPRAT(L1) = HPAV(L1) / NETFLX(L4)
                       IF ( .NOT. DUALBM )  HPRAT(2) = 0.0
  440          CONTINUE 
* 
                     CERAT = NETFLX(7) / PKAV 
* 
         IF ( .NOT. FLXPNT )  GO TO 500 
* 
*        Write out calculated values. 
* 
             WRITE (LU,2000)  (SORSNM(L), L=1,4), WAVE
 2000        FORMAT(//" Source ",4A2," at",F6.1,
     *              "cm.",/ 
     *               48"=") 
* 
               WRITE (LU,2100)  (TPLACE(L), L = 1,7)
 2100          FORMAT(/," Place:",4X,A2, 6(8X,A2))
* 
                 WRITE (LU,2200)  (NETFLX(L), L = 1,7)
 2200            FORMAT("Net flux:", 7G10.4)
* 
                   WRITE (LU,2300)  (NETER(L), L = 1,7) 
 2300              FORMAT( 1X,"+-error",1X, 7G10.4) 
* 
                     WRITE (LU,2400) TAVFN, DEVFN, PCNTDV 
 2400                FORMAT(/" Average of FN values - std dev given"
     *                      " as % of PEAK FLUX ie. shows drift ",/ 
     *                        8X, G10.4,"+-",F10.4," ie. ", F10.4,"%")
* 
               WRITE (LU,2500) PKAV, PKAVER, HPAV(1), HPAVER(1),
     *                         HPRAT(1), HPAV(2), HPAVER(2), HPRAT(2),
     *                         CERAT
 2500          FORMAT(/3X," Average peak flux (Jy): ",G10.4,"+-",G10.4/ 
     *                 6X," Half power averages: ", G10.4,"+-",G10.4, 
     *                    " Ratio HP/Peak:", F10.4,/
     *                 28X, G10.4,"+-",G10.4,15X,F10.4,/
     *                 6X," Ratio CE/Peak : ",F10.4)
* 
  500    CALL GSOFF (LU, NETFLX, WAVE, PEAK, PEAKER, DECOF, AVPK, JS, 
     *               RPTFQ, RAOFF, RAFER, NETER, FLXPNT, DECER, 
     *               AVPKER, DATE, IFQ) 
* 
*    Set limit on PKAV. Assume values>100 are erroneous.
* 
             IF ( PKAV .LT. 0.0 ) PKAV = 0.0
             IF ( ABS(PKAV) .GE. 100. )  PKAV = 99.99 
             IF ( ABS(PKAVER) .GE. 10.)  PKAVER = 9.999 
* 
         RETURN 
         END
* 
********************************************************************* 
* 
*         SUBROUTINE TO CALCULATE OFFSET IN DEC USING 3 POINTS
*            ON ASSUMED GAUSSIAN. 
* 
********************************************************************* 
* 
       SUBROUTINE GSOFF (LU, NETFLX, WAVE, PEAK, PEAKER, DECOF, AVPK, 
     *                   JS, RPTFQ, RAOFF, RAFER, NETER, FLXPNT, DECER, 
     *                   AVPKER, DATE, IFQ) 
* 
         LOGICAL DUALBM, RPTFQ, FLXPNT
* 
         REAL    SIGDEC(4), PEAK(3), DECOF(3), NETFLX(1), NETER(1), 
     *           CNPK(3), BMSEP(4), PEAKER(3), SIGASS(4), DCFLIM(4),
     *           NETFD, DECER(3)
* 
         DATA    SIGDEC /  .044,   .082,    .162,    .247   /,
     *           SIGASS /  .048,   .082,    .162,    .247   /,
     *           DCFLIM /  .055,   .140,    .185,    .250   /,
     *           SLPFAC / 10.9 /,   FRCSFC / 0.2028 /,
     *           BMSEP /  .213,   .242,   0.0,  0.0  /, 
     *           FAC1  / 0.693147181  / 
* 
            If (DATE .GT. 85.595 .AND. IFQ .EQ. 3) then  ! Feed off axis
                   SIGDEC(3) = .158 
                   DCFLIM(3) = .180 
            endif 
* 
           DUALBM = (( IFQ .EQ. 1 ) .OR. ( IFQ .EQ. 2 ))
* 
         DO 110  LL = 1,2 
                 L1 = LL*LL 
* 
              IF ( (LL .EQ. 2) .AND. ( .NOT. DUALBM ) )  GO TO 110
* 
*           Check for read fluxes being wrong sign due to low flux
*                or error in FN readings
* 
               RAT1 = NETFLX(L1)/NETFLX(L1+1) 
               RAT2 = NETFLX(L1)/NETFLX(L1+2) 
* 
           IF ((RAT1 .GE. 0.0) .AND. (RAT2 .GE. 0.0))  GO TO 120
* 
              IF ( .NOT. FLXPNT )  GO TO 100
                 WRITE (LU,1000)
 1000            FORMAT(" Incorrect sign. Check. PEAK set to 0.0 ") 
  100            PEAK(1) = 0.0
                 PEAK(2) = 0.0
                 PEAKER(1) = NETER(1) 
                 PEAKER(2) = NETER(4) 
                 DECOF(1) = 0.0 
                 DECOF(2) = 0.0 
                 DECER(1) = 0.0 
                 DECER(2) = 0.0 
C                LL = LL + 1
               GO TO 115
* 
  120        V = (NETFLX(L1) + NETFLX(L1+2)) / NETFLX(L1+1) 
* 
*         Require SQRT(V+1) - 1 > 0 ie V>0
* 
               IF ( V .LE. 0.0 )  GO TO 100 
* 
            X = SIGDEC(IFQ)*SIGDEC(IFQ)/(SIGASS(IFQ)*2.*FAC1) 
            Y = (EXP(2.*FAC1*SIGASS(IFQ)*SIGASS(IFQ)/ 
     *                            (SIGDEC(IFQ)*SIGDEC(IFQ))))/4.
            Z = (EXP(FAC1*SIGASS(IFQ)*SIGASS(IFQ)/
     *                            (SIGDEC(IFQ)*SIGDEC(IFQ))))/2.
            SW = SQRT(Y+V)
            OFFST = X * ALOG(SW - Z)
            OFSER = 0.5 * (X/SW) / (SW - Z) 
* 
                   V1 = NETFLX(L1) + NETFLX(L1+2) 
                   V1ER = SQRT( (NETER(L1)*NETER(L1)) + (NETER(L1+2)
     *                                                 *NETER(L1+2)) )
                  VFRCER =  ((V1ER*V1ER)/(V1*V1)) + ((NETER(L1+1)*
     *                       NETER(L1+1))/(NETFLX(L1+1)*NETFLX(L1+1)))
                  VER = V * SQRT(VFRCER)
C               WRITE (LU,1010) V, VER
C1010           FORMAT(" V ", F10.4,"+-", F10.6)
* 
             OFSER = OFSER * VER
* 
             DOS = OFFST / SIGDEC(IFQ)
* 
*            DOS limited by overflow error. 
* 
                  IF ( DOS .GT. 11. )  DOS = 11.
* 
           PEAK(LL) = NETFLX(L1) * EXP(FAC1 * DOS * DOS)
* 
               EFDD = EXP(FAC1*DOS*DOS) 
               NETFD = NETFLX(L1) * EFDD * 2.* FAC1 * DOS / SIGDEC(IFQ) 
               PERSQ = (EFDD*EFDD*NETER(L1)*NETER(L1)) +
     *                          (NETFD*NETFD*OFSER*OFSER) 
             PEAKER(LL) = SQRT(PERSQ) 
* 
           DECOF(LL) = OFFST
           DECER(LL) = OFSER
* 
*             Check that dec offset is less than limit for applicability
*               of gaussian fit in RA, else set flux to zero
* 
              If (ABS(DECOF(LL)) .GT. DCFLIM(IFQ)) GO TO 100  
* 
                 GO TO (111,112)  LL
  111                    IF ( PEAK(LL) .LT. 0.0 ) GO TO 100 
                         GO TO 110
  112                    IF ( PEAK(LL) .GT. 0.0 )  GO TO 100
* 
  110    CONTINUE 
* 
*             Single beam.....
  115              IF ( DUALBM )  GO TO 130 
                       AVPK = PEAK(1) 
                       AVPKER = PEAKER(1) 
                       DECOF(2) = 0.0 
                       DECER(2) = 0.0 
                   GO TO 140
* 
  130      AVPK = ( PEAK(1) - PEAK(2) )/2.
           AVPKER = 0.5 * SQRT( (PEAKER(1)*PEAKER(1)) + (PEAKER(2)
     *                                            *PEAKER(2)) ) 
* 
  140      IF ( .NOT. FLXPNT )  GO TO 190 
* 
         WRITE (LU,1100) (DECOF(J), DECER(J), J=1,2), 
     *                   ( PEAK(J),PEAKER(J), J=1,2), 
     *                   AVPK, AVPKER 
 1100    FORMAT(" DEC Offsets :  Beam A ", F10.6,"+-", F10.6,2X,
     *          " Beam B ", F10.6,"+-",F10.6,/
     *          " Peak fluxes ",  7X, 2(3X, F10.4,"+-", F8.4)/
     *          " AVERAGE : ", F10.4,"+-", F10.6) 
* 
  190      IF ( IFQ .NE. 2 )  GO TO 900 
*          Calculate offset of CENTRE.
* 
C            IF ( .NOT. RPTFQ )  GO TO 200
* 
*        If 2 readings at 6cm, choose which is to be used.
* 
C          CALL TWOSX (LU, RAOFF, NETFLX(7), AVPK)
C            GO TO 900
* 
  200    SLOPE = -1.0 * SLPFAC * AVPK 
* 
C               ERSLOP = 10. * AVPKER 
         RAOFF = NETFLX(7) / SLOPE
* 
                   FRCNET = NETER(7) / NETFLX(7)
                   FRCPK = AVPKER / AVPK
                 SQOFER = FRCNET*FRCNET + FRCPK*FRCPK + FRCSFC*FRCSFC 
                 RAFER = ABS(RAOFF) * SQRT(SQOFER)
* 
            IF ( AVPK .NE. 0.0 )  GO TO 900 
                   RAOFF = 0.0
                   RAFER = 0.0
* 
  900 RETURN
       END
* 
********************************************************************* 
* 
*        SUBROUTINE TO WRITE INFO ONTO TEMPORARY FILE 
* 
********************************************************************* 
* 
       SUBROUTINE TEMPF (LU, NSORS, SORSNM, WAVE, CALSRS, RPTFQ,
     *                   PEAK, PEAKER, AVPK, DECOF, DATE, 
     *                   PROB, TMDCB, TMFILE, JREC, JS, PKMEAS, 
     *                   PKMSER, DECER, RAOFF, RAFER, HA, TUBE) 
* 
         LOGICAL  CALSRS, RPTFQ, AOK, PROB
* 
         INTEGER  SORSNM(4), TMDCB(144), TMFILE(3), TMBUF(82),
     *            JTMSIZ(4), LBUF(400)
* 
         REAL     PEAK(3), PEAKER(3), DECOF(3), DECER(3)
* 
         DATA  JTMSIZ  / 0,24,0,0 /, JTMTYP / 3 /, JTMSEC / 0 /,
     *         JTMCR    / 52 /, JDCBS / 0 /, JNUM / 1 / 
* 
          CALL LGBUF (LBUF,400) 
* 
*       Set limit on DECOF for writing onto file. 
* 
             DO 300  LL = 1,2 
                 IF ( DECOF(LL) .GE. 10. )  DECOF(LL) = 9.99999 
  300        CONTINUE 
* 
*      Open temp file.
* 
  400    JREC = JREC + 1
  405      CALL OPEN (TMDCB, IERR, TMFILE, 0, JTMSEC, JTMCR, JDCBS) 
               IF ( IERR .EQ. -006 )  GO TO 410 
               IF ( IERR .GE. 0 )  GO TO 420
* 
*        Error routine. 
* 
  415    CALL CLOSE (TMDCB) 
             WRITE (LU, 4100) IERR
 4100        FORMAT(/" FMGR ERROR ",I4," in TEMPF. Check. OK?_")
               READ (LU,4000) AOK 
 4000          FORMAT(L6) 
* 
           IF ( AOK )  GO TO 405
               PROB = .TRUE.
                 GO TO 900
* 
*        Create file if doesn't exist.
  410      CALL CRETS (TMDCB, IERR, JNUM, TMFILE, JTMSIZ, JTMTYP, 
     *                 JTMSEC, JTMCR, JDCBS)
             IF (IERR .LT. 0)  GO TO 415
C          GO TO 405
* 
*          Position at start of file. 
  420    CALL POSNT (TMDCB, IERR, JREC, 1)
               IF (IERR .LT. 0)  GO TO 415
* 
*        Info to buffer.
* 
       CALL CODE
            WRITE (TMBUF, 4200)  NSORS, (SORSNM(L), L=1,4), WAVE, 
     *                           CALSRS, RPTFQ, (PEAK(L), L=1,2), 
     *                           (PEAKER(L), L=1,2), AVPK,
     *                           (DECOF(L), L=1,2), (DECER(L), L=1,2),
     *                           DATE, JS,
     *                           PKMEAS, PKMSER, RAOFF, RAFER, HA,
     *                           TUBE 
 4200  FORMAT(I3, 4A2, F6.1, 2L1, 5G12.6, 4F7.5, F10.6,I3, F5.2,F5.3, 
     *        2F7.5, F7.3, G12.6) 
* 
*          Write onto file. 
       CALL WRITF (TMDCB, IERR, TMBUF, 82)
             IF (IERR .LT. 0)  GO TO 415
* 
* 
       CALL CLOSE (TMDCB) 
* 
  900  RETURN 
         END
* 
********************************************************************* 
* 
*        SUBROUTINE TO CALCULATE SECOND OF TWO RA OFFSET VALUES 
*          IF TWO MEASUREMENTS WERE MADE AT 6 cm. ON THE SAME SOURCE. 
* 
********************************************************************* 
* 
       SUBROUTINE TWOSX (LU, RAOFF, NETFLX, AVPK) 
* 
         LOGICAL  AOK 
* 
         REAL     NETFLX
* 
           SLOPE = -10.0 * AVPK 
           OFST2 = NETFLX / SLOPE 
* 
  110        WRITE (LU,1100) RAOFF, OFST2 
 1100        FORMAT(" 6cm readings repeated."/
     *              " The 2 values of RA offset calculated are :",
     *              2(3X,F9.5),/" Which is wanted ? "/
     *              " Enter   1 if first value wanted,  2 for second"/
     *              " Enter   0 to get average    _") 
* 
         READ (LU,*)  JCHC
           JCHC = JCHC + 1
* 
         GO TO (111, 222, 333) JCHC 
* 
  111      PROVOF = ( RAOFF + OFST2 )/2.
             GO TO 400
  222      PROVOF = RAOFF 
             GO TO 400
  333      PROVOF = OFST2 
* 
  400  WRITE (LU,4100) PROVOF 
 4100  FORMAT(" Offset is : ", F9.5, " OK (T) ? _") 
* 
         READ (LU,4200) AOK 
 4200    FORMAT(L6) 
* 
           IF ( .NOT. AOK )  GO TO 110
* 
       RAOFF  = PROVOF
* 
      RETURN
       END
* 
********************************************************************* 
* 
*      SUBROUTINE TO CORRECT FOR RA OFFSET AND CALIBRATE MEASURED FLUXES
*               Reads from TMFILE.
* 
********************************************************************* 
* 
       SUBROUTINE CALIB (LU, TMDCB, TMFILE, WRTFNL, 
     *                   RESFIL, k,*) 
* 
         PARAMETER (NSCS = 400) 
* 
         CHARACTER    CALNAM(20)*8
* 
         LOGICAL AOK, CALFND, CALSRS, RPTFQ, DUALBM, WRTFNL, NEWCAL,
     *           SPCRAT 
* 
         INTEGER TMDCB(144), TMFILE(3), TMBUF(82), SORSNM(4),NUMCL(20), 
     *            IBLNK(1), NAMFIN(3), NWCLNM(4), RESFIL(3), LBUF(400), 
     *            NDECF(4)
* 
         REAL    PEAK(3), PEAKER(3), DECOF(3), SIGMA(4),
     *           BMSEP(4), CNPK(2), CALFLX(4,20), CALER(4,20), FRAT(4), 
     *           FRATER(4), CNPKER(2), OFLIMT(4), 
     *           AVDECF(4,2), 
     *           DCDEV(4,2), WVLEN(4),
     *           DECER(3),
     *           RAFARY(NSCS), RARYR(NSCS), DCFARY(4,NSCS,2),
     *           DCARYR(4,NSCS,2), DCF(NSCS), DCFR(NSCS)
C
        EMA DCFARY,DCARYR
*
        DATA   WVLEN /  3.4,    6.0,    13.0,   18.0 /,
     *         SIGMA / .048,   .082,   .162,    .247 /,
     *         BMSEP / .216,   .254,  0.0,      .250 /,
     *         OFLIMT/ .072,   .123,   .123,    .250 /,
     *         FAC1  / 0.693147181  /,
     *         CALFND / .TRUE.  /,
     *         NDUMY / 0 /,  NCAL / 0 /,
     *         IBLNK / 2H   /
*
               CALL LGBUF(LBUF,400)
*
              NWCALN = 0
              NRAF = 0
            DO 501 NZ = 1,150
                   RAFARY(NZ) = 0.
                   RARYR(NZ) = 0. 
                   DCF(NZ) = 0. 
                   DCFR(NZ) = 0.
            DO 500  NFQ = 1,4 
                  NDECF(NFQ) = 0
                  DO 505 K = 1,2
                     DCFARY(NFQ,NZ,K) = 0.
                     DCARYR(NFQ,NZ,K) = 0.
  505             CONTINUE
  500       CONTINUE
  501       CONTINUE
* 
*     Open temp file. 
* 
  100    CALL OPEN (TMDCB, IERR, TMFILE)
             IF (IERR .GE. 0 )  GO TO 120 
* 
*        Error routine. 
* 
  110       CALL CLOSE (TMDCB)
               WRITE (LU,1100) IERR 
 1100          FORMAT(" FMGR ERROR ", I4," in CALIB. OK (T) ? _") 
                   READ (LU,1200) AOK 
 1200              FORMAT(L6) 
* 
               IF ( .NOT. AOK )  GO TO 800
             GO TO 100
* 
*          Position at start of file. 
* 
  120        CALL POSNT (TMDCB, IERR, 1, 1) 
                  IF (IERR .LT. 0)  GO TO 110 
* 
*        Write headers. 
                   IF ( .NOT. CALFND )  GO TO 125 
* 
                        WRITE (LU,1300) 
 1300                   FORMAT(//" CALIBRATORS AND CONTROLS ")
               WRITE (LU,1360)
 1360          FORMAT(" No.",2X," Source ",2X," Wave ",3X," Flux (Jy)", 
     *                  6X," Dec offsets ",2X," RA off ",2X,"HA") 
         GO TO 130
* 
  125          WRITE (LU,1330)
 1330          FORMAT(//,3X,"FINAL VALUES",/
     *             1X," No.",2X,"Source",4X,"Wave",4X,"Meas. Flux", 
     *                4X,"Corr. Flux",3X,"Dec offsets",5X,"RA off", 
     *                2X,"HA",
     *               /1X,4"-",2X,6"-",4X,4"-",4X,10"-",4X,10"-",3X, 
     *                11"-",5X, 6"-",2X,2"-") 
* 
*            Read data. 
  130      CALL READF (TMDCB, IERR, TMBUF, 82, LEN) 
               IF ((IERR .EQ. -012) .OR. (LEN .LE. 0))  GO TO 300 
               IF (IERR .LT. 0)  GO TO 110
* 
*            Assign values from buffer to program parameters. 
* 
        CALL CODE 
           READ (TMBUF, 1400)  NSORS, SORSNM, WAVE, CALSRS, RPTFQ,
     *                         (PEAK(L), L=1,2), (PEAKER(L), L=1,2),
     *                         AVPK, (DECOF(L), L=1,2), 
     *                         (DECER(L), L=1,2), DATE,JS,
     *                         PKMEAS, PKMSER, RAOFF, RAFER, HA,
     *                         TUBE 
 1400  FORMAT(I3, 4A2, F6.1, 2L1, 5G12.6, 4F7.5, F10.6,I3, F5.2,F5.3, 
     *        2F7.5, F7.3, G12.6) 
* 
*        Introduce IFQ for easier manipulation
* 
            IF ( WAVE .EQ. 3.4 )  IFQ = 1 
            IF ( WAVE .EQ. 6.0 )  IFQ = 2 
            IF ( WAVE .EQ. 13.0)  IFQ = 3 
            IF ( WAVE .EQ. 18.0)  IFQ = 4 
         DUALBM = ( (IFQ .EQ. 1) .OR. (IFQ .EQ. 2) )
* 
         If ( CALFND ) then 
           IF ( .NOT. CALSRS .AND. (NSORS .NE. NWCALN))  GO TO 130
* 
*           Put calibrator values into array for later manipulation.
            If ( NSORS .NE. NDUMY .OR. (IFQ .EQ. KDNFQ)) then 
* 
              NCAL = NCAL + 1 
              NDUMY = NSORS 
              KDNFQ = IFQ 
            endif 
         endif
* 
*        Correct for RA offset using value obtained from 6 cm.
* 
         CALL RACOR (RAOFF, RAFER, PKMEAS, PKMSER, IFQ, DUALBM, 
     *               PEAK, PEAKER, FLXAV, FLXER, WVLEN, DATE) 
* 
*      Multiply by calibration factor.
        IF ( CALFND )  GO TO 250
* 
           FINFLX = FRAT(IFQ) * FLXAV 
           FINER = FRAT(IFQ) * FLXER
* 
          SQFR1 = (FLXER*FLXER) / (FLXAV*FLXAV) 
          SQFR2 = (FRATER(IFQ)*FRATER(IFQ)) / (FRAT(IFQ)*FRAT(IFQ)) 
* 
           FINER = FINFLX * SQRT( SQFR1 + SQFR2 ) 
* 
*          Set limit to FINER as written onto file in F5.3
             IF (FINER .GT. 99.99)  FINER = 99.99999
* 
*          If measurement repeated choose which is wanted.
  250      IF (RPTFQ)  WRITE (LU,2222)
 2222                  FORMAT(" TEST. RPTFQ TRUE AFTER AVERAG CALC ") 
* 
*        If searching for CALIBRATORS don't write onto file.
* 
         IF ( CALFND )  GO TO 290 
* 
             IF ( .NOT. WRTFNL )  GO TO 260 
*          Write results onto file
* 
           CALL WFINL (NSORS, SORSNM, WAVE, DATE, FINFLX, FINER,
     *                 DECOF, DECER, RAOFF, HA, RESFIL, *900,k) 
* 
C        CALL WFIN2 (LU, NSORS, SORSNM, WAVE, DATE, 
C    *               FINFLX, FINER, PKMEAS, PKMSER, DECOF, RAOFF, 
C    *               NAMFIN)
* 
  260    WRITE (LU,3000)  NSORS, IBLNK, SORSNM, WAVE, PKMEAS, PKMSER, 
     *                    FINFLX, FINER, (DECOF(L), L=1,2),RAOFF, 
     *                    HA
 3000    FORMAT(I4, A2, 4A2, F6.1, 2(2X,F5.2,"+-",F5.3),F7.3,2(2X,F7.3),
     *            2X,F4.1)
* 
*        Set limit at 2 Jy for use of values to calc ave offset 
* 
       If   (FINFLX .GE. 2.0)  then 
* 
         If (RAOFF .NE. 0.0 .AND. WAVE .EQ. 6.0) then 
              NRAF = NRAF + 1 
              If (NRAF .GT. NSCS) then
                 PRINT '(//"Too many sources. Change parameter NSCS ",
     *                    "in subroutine CALIB and recompile")' 
                   k = 1
                   GO TO 900
              endif 
              RAFARY(NRAF) = RAOFF
              RARYR(NRAF) = RAFER 
         endif
* 
         If (DECOF(1) .NE. 0.0 .OR. DECOF(2) .NE. 0.0) then 
                  NDECF(IFQ) = NDECF(IFQ) + 1 
              If (NDECF(IFQ) .GT. NSCS) then
                 PRINT '(//"Too many sources. Change parameter NSCS ",
     *                    "in subroutine CALIB and recompile")' 
                   k = 1
                   GO TO 900
              endif 
              DO  510  K = 1,2
                  DCFARY(IFQ,NDECF(IFQ),K) = DECOF(K) 
                  DCARYR(IFQ,NDECF(IFQ),K) = DECER(K) 
  510         CONTINUE
         endif
       endif
* 
           GO TO 130
* 
  290  CALFLX(IFQ,NCAL) = FLXAV 
       CALER(IFQ, NCAL) = FLXER 
       NUMCL(NCAL) = NSORS
* 
          Assign 295 to label 
             NSTMT = 295
  295     WRITE (CALNAM(NCAL),'(4A2)', iostat=ierr, err=820)
     *                    (SORSNM(LS), LS=1,4)
         WRITE (LU, 3100) NSORS, IBLNK, SORSNM, WAVE, FLXAV, FLXER, 
     *                    (DECOF(L), L=1,2), RAOFF, HA
 3100    FORMAT(I4, A2, 4A2, F6.1, 3X, F5.2,"+-", F6.3, 3(2X,F7.3), 
     *          2X,F4.1)
* 
*            Return to read next line.
       GO TO 130
* 
  300      IF ( .NOT. CALFND )  GO TO 800 
* 
*               If there are no calibrators in list, choose another 
* 
               If (NCAL .EQ. 0) then
                   PRINT '(/"No calibrators in list."/"Want to specify",
     *                     " another source as a calibrator? Yes(T) _")'
                    READ *, NEWCAL
* 
                       If (NEWCAL) then 
                           PRINT '("Name of new calibrator ? _")' 
                            READ '(4A2)', (NWCLNM(LN), LN=1,4)
                                CALL NAMNO (NWCALN, NWCLNM, RAHR,*900,k)
                                   NDUMY = 0
                       else 
                         PRINT '("Want to specify cal ratio? Yes (T)_")'
                            READ *, SPCRAT
                              If (SPCRAT) then
  310                            DO 315  LF = 1,4 
                                    PRINT '("Enter cal ratio and error",
     *                                       " for",(F5.1)," cm : _")', 
     *                                          WVLEN(LF) 
                                    READ *, FRAT(LF), FRATER(LF)
  315                            CONTINUE 
* 
                                    DO 320 LF = 1,4 
                                       PRINT '(/"Cal ratio for",(F5.1), 
     *                                         " cm is,"(F7.4)," +-", 
     *                                         (F7.4))',
     *                                         WVLEN(LF), FRAT(LF), 
     *                                         FRATER(LF) 
  320                               CONTINUE
* 
                                      PRINT '(12X,"OK (T)? _")' 
                                        READ *, AOK 
                                        IF (.NOT. AOK )  GO TO 310
                              else
                                  DO 330 LN=1,4 
                                     FRAT(LN) = 1.0 
                                     FRATER(LN) = 0.0 
  330                             CONTINUE
                              endif 
                                     CALFND = .FALSE. 
                       endif
                         GO TO 120
               else 
* 
*                  Choose CALIRATION RATIO. 
* 
                       CALL FCALB (LU, NCAL, CALFLX, CALER, FRAT, 
     *                             FRATER, NUMCL, CALNAM, NEWCAL) 
* 
                          CALFND = .FALSE.
* 
*                  Return to find final corrected and calibrated fluxes.
                          GO TO 120 
               endif
* 
  800    CONTINUE 
* 
*          Calculate averages,std dev 
* 
            CALL WMEAN (RAFARY, RARYR, NRAF, AVRAOF, RADEV, AV, ER, DEV)
* 
            DO 520  NFQ = 1,4 
               DO 525 K = 1,2 
                  DO 526 NZ = 1,NDECF(NFQ)
                     DCF(NZ) = DCFARY(NFQ,NZ,K) 
                     DCFR(NZ) = DCARYR(NFQ,NZ,K)
  526             CONTINUE
               CALL WMEAN (DCF, DCFR, NDECF(NFQ), AVDECF(NFQ,K),
     *                     DCDEV(NFQ,K), AV, ER, DEV) 
  525         CONTINUE
  520      CONTINUE 
* 
           WRITE (LU,5100) AVRAOF, RADEV, NRAF, 
     *                     (WVLEN(NFQ), (AVDECF(NFQ,K), 
     *                     DCDEV(NFQ,K), K=1,2), NDECF(NFQ), NFQ = 1,4) 
 5100      FORMAT(//15X,"Ave RA offset :", 4X, F7.3,"+-",F8.3,3X,I4/
     *              15X,"Ave DEC offsets at :"/ 
     *               4(20X,F6.1,"cm : ",2(3X, F7.3,"+-", F8.3),3X,I4/)) 
* 
            GO TO 900 
* 
*         Error routine 
  820         CALL ERMSG (ierr, NSTMT, *900,k)
                    GO TO label 
* 
  900    RETURN k 
         END
* 
********************************************************************* 
* 
*           Subroutine to correct for RA offset 
* 
********************************************************************* 
* 
         SUBROUTINE RACOR (RAOFF, RAFER, PKMEAS, PKMSER, IFQ, DUALBM, 
     *                     PEAK, PEAKER, FLXAV, FLXER, WVLEN, DATE) 
* 
         LOGICAL DUALBM, OK6
* 
         REAL SIGMA(4), BMSEP(4), CNPK(2), CNPKER(2), PEAK(3),
     *        PEAKER(3), WVLEN(4), SEPER(4), RAFLIM(4)
* 
*                       3.4     6.0     13.0     18.0 
        DATA   SIGMA / .046,   .081,   .162,     .247 /,
     *         BMSEP / .213,   .242,  0.0,      0.0   /,
     *         RAFLIM/ .054,   .093,   .175,     .250 /,
     *         SEPER / .002,   .006,  0.0,      0.0   /,
     *         FAC1  / 0.693147181  / 
* 
               If (DATE .GT. 86.945) then    ! New 3.6 cm feeds 
                   BMSEP(1) = .195
                   SEPER(1) = 0.0 
               endif
* 
           If (DATE .GT. 85.595 .AND. IFQ .EQ. 3) then   !Feed off axis 
                    SIGMA(3) = .166 
                    RAFLIM(3) = .170
             endif
* 
             If (IFQ .EQ. 2) then 
                If (PEAK(1) .EQ. 0.0 .OR. PEAK(2) .EQ. 0.0) then
                 OK6 = .FALSE.
                else
                 OK6 = .TRUE. 
                endif 
             endif
* 
             DUMMY = RAOFF / SIGMA(IFQ) 
* 
*         Limit DUMMY for EXP 
                IF ( ABS(DUMMY) .GT. 11. ) DUMMY = 11.
         A = EXP( -FAC1 * DUMMY * DUMMY)
* 
                     DYER = RAFER / SIGMA(IFQ)
                     EMDF2 = 2.* FAC1 * DUMMY * A 
                     AER = ABS (EMDF2 * DYER) 
C                PRINT '(" A ", F10.4,"+-", F10.6)', A, AER 
* 
* 
*          Single beam..... 
        If ( .NOT. DUALBM )  then 
                        CNPK(1) = PEAK(1) / A 
                        CNPK(2) = 0.0 
                           CNRSQ = (AER*AER)/(A*A) + (PEAKER(1)*
     *                               PEAKER(1))/(PEAK(1)*PEAK(1)) 
                        CNPKER(1) = CNPK(1) * SQRT(CNRSQ) 
                        CNPKER(2) = 0.0 
                        FLXAV = CNPK(1) 
                        FLXER = CNPKER(1) 
* 
         else 
               DMINUS = (RAOFF - BMSEP(IFQ) - SEPER(IFQ)) / SIGMA(IFQ)
* 
                   IF ( ABS(DMINUS) .GT. 11. ) DMINUS = 11. 
         B = EXP( -FAC1 * DMINUS * DMINUS)
* 
                     EMFD2 = 2.* FAC1 * DMINUS * B
                     BER = ABS(EMFD2 * DYER)
* 
              EXPD = (RAOFF - SEPER(IFQ)) / SIGMA(IFQ)
         D = EXP( -FAC1 * EXPD * EXPD)
              DER = ABS(2. * FAC1 * EXPD * D * DYER)
* 
C               PRINT '(" B, D ",2(F10.4,"+-", F10.6))', B, BER, D, DER 
* 
             DPLUS = ( RAOFF + BMSEP(IFQ)) / SIGMA(IFQ) 
* 
                      IF ( ABS(DPLUS) .GT. 11. ) DPLUS = 11.
         C = EXP( -FAC1 * DPLUS * DPLUS)
                     EMFD2 = 2.* FAC1 * DPLUS * C 
                     CER = ABS(EMFD2 * DYER)
* 
                  EXPE = (RAOFF + SEPER(IFQ)) / SIGMA(IFQ)
         E = EXP( -FAC1 * EXPE * EXPE)
             EER = ABS(2. * FAC1 * EXPE * E * DYER) 
C               PRINT '(" C, E ",2(F10.4,"+-", F10.6))', C, CER, E, EER 
* 
                   DMB = D - B
                   DMBER = SQRT( (DER*DER) + (BER*BER) )
                   CME = C - E
                   CMEER = SQRT( (CER*CER) + (EER*EER) )
                       DBRTSQ = (DMBER*DMBER) / (DMB*DMB) 
                       CERTSQ = (CMEER*CMEER) / (CME*CME) 
* 
       CNPK(1) = PEAK(1) / DMB
                   PKRAT = PEAKER(1)/PEAK(1)
                   CNPKER(1) = CNPK(1) * SQRT( DBRTSQ + (PKRAT*PKRAT))
* 
       CNPK(2) = PEAK(2) / CME
                   PKRAT = PEAKER(2)/PEAK(2)
                   CNPKER(2) = CNPK(2) * SQRT( CERTSQ + (PKRAT*PKRAT))
* 
       FLXAV = ( CNPK(1) + CNPK(2) ) / 2. 
       FLXER = 0.5 * SQRT( (CNPKER(1)*CNPKER(1)) + (CNPKER(2)*
     *                                         CNPKER(2)) ) 
* 
         endif
* 
            If (.NOT. OK6) then 
               DUM =  RAFLIM(2) / (2. * SIGMA(IFQ))   
               A2 = EXP ( -FAC1 * DUM * DUM)
               FLXER = FLXER / A2 
            endif 
* 
         If (ABS(RAOFF) .GT. RAFLIM(IFQ)) then
             FLXAV = 0.0
             FLXER = 0.0
         endif
* 
         IF ( FLXAV .LT. 0.0 )  FLXAV = 0.0 
* 
        RETURN
        END 
* 
********************************************************************* 
* 
*       SUBROUTINE TO WRITE FINAL DATA ONTO FILE
* 
********************************************************************* 
* 
       SUBROUTINE WFINL (NSORS, SORSNM, WAVE, DATE, FINFLX, FINER,
     *                   DECOF, DECER, RAOFF, HA, RESFIL, k,*)  
* 
      CHARACTER NAMFIN*14 
* 
      INTEGER   SORSNM(4), RESFIL(3), FRSFIL(3) 
* 
      REAL      DECOF(2), FDECOF(2), DECER(2) 
* 
      DATA   NAMFIN(1:1) / 'C' /, NAMFIN(7:14) / '::27:4:4' / 
* 
          NWAVE = IFIX(WAVE*10.)
          IF (WAVE .EQ. 13.) NWAVE = 13 
          IF (WAVE .EQ. 18.) NWAVE = 18 
          NSORS = NSORS + 100 
          WRITE (NAMFIN(2:3),'(I2)',iostat=ierr, err=820) NWAVE 
          WRITE (NAMFIN(4:6),'(I3)',iostat=ierr,err=820) NSORS
* 
C         PRINT '("Name ",(A14))', NAMFIN 
* 
*       File operations 
* 
*    Open file,if it exists 
* 
       Assign 110 to label
         NSTMT=110
  110  OPEN (10, file=NAMFIN, iostat=ierr, err=810, status='OLD') 
* 
*         Find EOF
       Assign 120 to label
         NSTMT=120
  120  DO WHILE (ierr .ne. -1)
          READ (10, 1100, iostat=ierr, err=820, end=150) FDATE,FLUX,
     *              ERROR, (FDECOF(L),L=1,2), FRAOFF, FHA,
     *              (FRSFIL(L),L=1,3) 
        end do
        Assign 150 to label 
           NSTMT=150
  150   BACKSPACE(10, iostat=ierr, err=820) 
* 
*          Write info to file 
* 
         Assign 160 to label
          NSTMT=160 
  160    WRITE (10, 1100, iostat=ierr, err=820)  DATE, FINFLX, FINER, 
     *                         (DECOF(L),L=1,2), RAOFF, HA, 
     *                         (RESFIL(L),L=1,3)
* 
*         Write EOF 
          Assign 170 to label 
              NSTMT=170 
  170     ENDFILE (10, iostat=ierr, err=820)
* 
              Assign 180 to label 
                   NSTMT=180
  180         CLOSE (10, iostat=ierr, err=820)
* 
       GO TO 900
* 
*      Error handling 
* 
*         Check if error due to file not existing.
* 
  810     If (ierr .eq. 462) then 
              PRINT '("Creating file ",(A14))', NAMFIN
           OPEN (10, file=NAMFIN, iostat=ierr, err=820, status='NEW') 
* 
               Assign 815 to label
                  NSTMT = 815 
  815          WRITE (10, '(4A2)',iostat=ierr, err=820)(SORSNM(L),L=1,4)
           GO TO 120
         endif
* 
*     General error checking. 
* 
  820     Call ERMSG (ierr, NSTMT, *900, k) 
          GO TO label 
* 
 1100    FORMAT (F10.6, F5.2, F5.3, 4F6.4, 3A2) 
* 
  900   RETURN k
      END 
* 
********************************************************************* 
* 
*        Subroutine to display errors 
* 
********************************************************************* 
* 
       SUBROUTINE ERMSG (ierr, NSTMT, k,*)
* 
       LOGICAL AOK
* 
         PRINT '("Error ",(I4)," at statement",(I4))', ierr, NSTMT
         PRINT '("OK(T) or end program(F) ? _")'
            READ '(L1)', aok
            If (.NOT. aok)  k=1 
       RETURN k 
       END
* 
********************************************************************* 
* 
*        SUBROUTINE TO CALCULATE RATIO OF MEASURED FLUX TO KNOWN
*              CALIBRATOR OR CONTROL FLUX.
* 
********************************************************************* 
* 
       SUBROUTINE FCALB (LU, NCAL, CALFLX, CALER, FRAT, FRATER, NUMCL,
     *                   CALNAM, NEWCAL)
*
         CHARACTER*8  CALNAM(20)
*
         LOGICAL AOK, NEWCAL
*
         INTEGER NWANT(20), NUMCL(20)
*
         REAL    CALRAT(20), CRATER(20), CALFLX(4,20), CALER(4,20),
     *           TRUFLX(4,6), FRAT(4), FRATER(4), VAL(20), ERR(20)
*
*              WAVE        3.4       6.0        13.0        18.0
*  0915-11
         DATA  TRUFLX /    8.1,     13.5,      27.3,         0.0,
*  1648+05
     *                     6.36,    11.33,     27.9,         0.0,
*  2203-18
     *                     3.38,    4.02,      5.48,         0.0,
*  0500+01
     *                     1.50,    1.95,      2.32,         0.0,
*  1421-49
     *                     3.60,    5.54,      7.51,         0.0,
*  Specified during running
     *                     0.0,     0.0,       0.0,          0.0   /
*
           WRITE (LU,1100)  NCAL, (NUMCL(L), L=1,NCAL)
 1100      FORMAT(/I4," Calibrators used.  Source nos :: ",20I4)
*
*        Loop for each freq.
*
         DO 120  IFQ = 1,4
*
                 IF ( IFQ .EQ. 1 )  WAVE = 3.4
                 IF ( IFQ .EQ. 2 )  WAVE = 6.0
                 IF ( IFQ .EQ. 3 )  WAVE = 13.0
                 IF ( IFQ .EQ. 4 )  WAVE = 18.0
*
             WRITE (LU,1200) WAVE 
 1200        FORMAT(/" Calibrators at ", F5.1," cm."/ 
     *               1X, 24"*"/ 
     *             /"   No.   Meas. flux     True flux      Ratio "/
     *              "  ====   ==========     =========      ===== ")
* 
*              Loop for each source.
             DO 130  MCAL = 1,NCAL
* 
                       IF (CALNAM(MCAL) .EQ. '0915-11 ')  KCAL = 1
                       IF (CALNAM(MCAL) .EQ. '1648+05 ')  KCAL = 2
                       IF (CALNAM(MCAL) .EQ. '2203-18 ')  KCAL = 3
                       IF (CALNAM(MCAL) .EQ. '0500+01 ')  KCAL = 4
                       IF (CALNAM(MCAL) .EQ. '1421-49 ')  KCAL = 5
                       IF (NEWCAL)                        KCAL = 6
* 
                  If ( CALFLX(IFQ,MCAL) .EQ. 0.0 )  then
                       CALRAT(MCAL) = 0.0 
                       CRATER(MCAL) = 0.0 
                     GO TO 125
                  endif 
* 
                If (NEWCAL .AND. MCAL .LT. 2)  then 
                          PRINT '("Enter TRUE FLUX for ",(A8)," at ", 
     *                             (F5.1)," cm. :: _")', CALNAM(MCAL),
     *                                                   WAVE 
                             READ *, TRUFLX(IFQ,KCAL) 
                endif 
* 
*         Calibration ratio for this source.
* 
  115  CALRAT(MCAL) = TRUFLX(IFQ,KCAL) / CALFLX(IFQ,MCAL) 
       CRATER(MCAL) = CALRAT(MCAL) * CALER(IFQ,MCAL) / CALFLX(IFQ,MCAL) 
* 
  125        WRITE (LU,1400) NUMCL(MCAL), CALFLX(IFQ,MCAL), 
     *                    TRUFLX(IFQ,KCAL), CALRAT(MCAL), CRATER(MCAL)
 1400        FORMAT(1X,I4, 3(3X,F10.4),"+-", F6.4)
* 
  130        CONTINUE 
* 
*            Set NWANT to zero
            GO TO 210 
* 
  150    WRITE (LU,1500)
 1500    FORMAT(/" Which are wanted ? "/
     *         18X," Enter 11 if there are no readings at this freq.",
     *          /8X, " Enter 1 for first, 2 for second, etc, ", / 
     *          10X," More than one no. entered causes those to be" 
     *              " averaged :  _") 
       READ *, NWANT
C      READ '(20I2)' , NWANT
          IF ( NWANT(1) .NE. 11 )  GO TO 155
                 FRAT(IFQ) = 1.0
              GO TO 120 
* 
*        Calc average 
  155  DO 160  L = 1,20 
             NUM = NWANT(L) 
               IF ( NUM .EQ. 0 )  GO TO 170 
             VAL(L) = CALRAT(NUM) 
             ERR(L) = CRATER(NUM) 
  160  CONTINUE 
* 
  170    CALL WMEAN (VAL, ERR, (L-1), FRAT(IFQ), FRATER(IFQ), AV,ER,DV) 
* 
         WRITE (LU,1700) WAVE, FRAT(IFQ), FRATER(IFQ) 
 1700    FORMAT(" Calibrator ratio for ",F5.1," cm. measurements" 
     *          " is :", F10.4,"+-",F6.4/ 
     *          8X, " OK(T) ? _") 
* 
             READ (LU,2000) AOK 
 2000        FORMAT(L6) 
* 
       IF ( AOK )  GO TO 120
  210        DO 180 LL = 1,20 
                 NWANT(LL) = 0
  180        CONTINUE 
           GO TO 150
* 
  120  CONTINUE 
* 
         RETURN 
         END
* 
********************************************************************* 
* 
*          SUBROUTINE TO WRITE FINAL DATA ONTO DISC FILE. 
* 
********************************************************************* 
* 
         SUBROUTINE WFIN2 (LU, NSORS, SORSNM, WAVE, DATE, 
     *                     FINFLX, FINER, PKMEAS, PKMSER, DECOF,
     *                     RAOFF, NAMFIN) 
* 
         LOGICAL AOK
* 
         INTEGER SORSNM(4), NAMFIN(3), FNDCB(144), FNBUF(42), IBLNK(1)
* 
         REAL    DECOF(3) 
* 
         DATA  IBLNK / 2H   / 
* 
*        Open file
  100   CALL OPEN (FNDCB, IERR, NAMFIN) 
           IF ( IERR .GE. 0 )  GO TO 120
           IF ( IERR .EQ. -006)  GO TO 777
* 
*            Error routine. 
  110       CALL CLOSE (FNDCB)
               WRITE (LU,1100) IERR 
 1100          FORMAT(" FMGR ERROR ", I4," in WFINL. Check. OK(T)?_") 
                   READ (LU,1200) AOK 
 1200              FORMAT(L2) 
* 
              IF ( .NOT. AOK )  GO TO 900 
              GO TO 100 
* 
*       CREATE FILE 
  777          CALL CREAT (FNDCB, IERR, NAMFIN, 128, 3, 0, 15)
                    IF ( IERR .LT. 0 )  GO TO 110 
* 
*          Position at EOF. 
* 
  120        CALL POSNT (FNDCB, IERR, 32766)
                 IF ( IERR .EQ. -012 )  GO TO 140 
                 IF ( IERR .LT. 0 )  GO TO 110
* 
*        Data to buffer.
* 
  140      CALL CODE
              WRITE (FNBUF, 1400) NSORS, IBLNK,(SORSNM(L), L=1,4), WAVE,
     *                            DATE, FINFLX, FINER, PKMEAS, PKMSER,
     *                            (DECOF(J), J=1,2), RAOFF
 1400       FORMAT(I4, A2, 4A2, F6.1, F12.6, 2G12.6, F5.2, F5.3, 3F6.4) 
* 
*           Write onto file.
* 
               CALL WRITF (FNDCB, IERR, FNBUF, 42)
                     IF ( IERR .LT. 0 )  GO TO 110
* 
*              Write EOF. 
* 
                  CALL WRITF (FNDCB, IERR, FNBUF, -1) 
                         IF ( IERR .LT. 0 )  GO TO 110
* 
*            Close file.
* 
             CALL CLOSE (FNDCB) 
* 
  900    RETURN 
         END
* 
* 
********************************************************************* 
* 
*    WRITE DISC DATA FILE NF... TO MAG TAPE 
* 
********************************************************************* 
* 
       SUBROUTINE MAGTP (RESFIL, k,*) 
* 
       CHARACTER DATFIL*10, DBUFR*42, TPBUF*29
* 
       LOGICAL   AOK, EOLF
* 
       INTEGER   ITIME(5), RESFIL(3)
* 
       DATA      DATFIL (7:10) / '::20' /, NCR / 20 /, MTOUT / 8 /
* 
             WRITE (DATFIL(9:10),'(I2)', iostat=ierr, err=820) NCR
             WRITE (DATFIL(1:6),'(3A2)', iostat=ierr, err=820)
     *                                    (RESFIL(L), L=1,3)
*         Get time
* 
                 CALL EXEC (11, ITIME, IYEAR) 
* 
              GO TO 101 
* 
  100        PRINT '("Data on cartridge no :? _")'
             READ *, NCR
             WRITE (DATFIL(9:10),'(I2)', iostat=ierr, err=820) NCR
* 
             PRINT '("Magtape LU no : _")'
             READ *, MTOUT
* 
 101   PRINT '("Writing data from ",(A10)," to magtape",(I3),/,4X,
     *         " at",(2I5,4I3))',DATFIL,MTOUT,IYEAR,(ITIME(L),L=5,1,-1) 
* 
        PRINT '(9X,"OK (T)? _")'
        READ *, AOK 
        IF ( .NOT. AOK )  GO TO 100 
* 
* 
*      Start tape with EOF mark to be compatible with rest
* 
             PRINT '("Is this a new tape, Yes(T)? _")'
             READ *, AOK
             If ( AOK ) then
                      ENDFILE (MTOUT, iostat=ierr, err=820) 
                      GO TO 121 
                 else 
                      PRINT '("Tape at end of last file ? Yes(T) _")' 
                      READ *, EOLF
             endif
* 
*      If just written into previous file ie. at end of 2X EOF, write data
* 
               If (EOLF)  GO TO 115 
* 
*      Data on tape terminated by double EOF mark. Find it
* 
  111  Do while (ierr .ne. -1)
                 Assign 110 to label
                   nstmt=110
  110            READ (MTOUT, iostat=ierr, err=820, end=112) TPBUF
       end do 
  112              Assign 113 to label
                    nstmt=113 
  113              READ (MTOUT, iostat=ierr, err=820, end=115) TPBUF
       GO TO 111
* 
  115          Assign 116 to label
                NSTMT=116 
  116          BACKSPACE (MTOUT, iostat=ierr, err=820)
* 
*      Write name, date as 1st record 
* 
  121   Assign 120 to label 
           NSTMT = 120
  120    WRITE (MTOUT, 1100, iostat=ierr, err=820) DATFIL, IYEAR, 
     *                                             (ITIME(L), L=5,1,-1) 
* 
*         Dump file to tape, line by line 
* 
        Assign 130 to label 
             NSTMT = 130
  130      OPEN (20, file=DATFIL, iostat=ierr, err=820, status='OLD') 
* 
         Do while (ierr .ne. -1)
              Assign 132 to label 
                NSTMT=132 
  132          READ (20, 1300, iostat=ierr, err=820, end=140) DBUFR 
* 
                Assign 135 to label 
                   NSTMT=135
  135          WRITE (MTOUT, 1300, iostat=ierr, err=820) DBUFR
          end do
* 
*           Write double EOF
* 
  140           Assign 145 to label 
                   NSTMT=145
  145        ENDFILE (MTOUT, iostat=ierr, err=820)
             ENDFILE (MTOUT, iostat=ierr, err=820)
* 
       GO TO 850
* 
  820    CALL ERMSG (ierr, NSTMT, *900,k) 
         GO TO label
* 
 1100    FORMAT (A10, I4, 5I3)
 1300    FORMAT (A42) 
* 
  850    CONTINUE 
          PRINT '("Rewind tape ? Yes(T) _")'
                READ *, AOK 
                If (AOK) then 
                        Assign 855 to label 
                          nstmt=855 
  855                   REWIND (MTOUT, iostat=ierr, err=820)
                endif 
* 
  900  RETURN k 
       END
* 
*********************************************************************** 
* 
*       SUBROUTINE TO GET CORRECT SOURCE NO. FOR FINAL FILE NAME
*         ie. return to old nos.
* 
*********************************************************************** 
* 
         SUBROUTINE NAMNO (NSORS, SORSNM, RAHR,k,*) 
* 
         CHARACTER NAMSRS*8, SNAMFL*10, FILSNM*8, IBLNK*2 
* 
         LOGICAL   TWOTRY 
* 
         INTEGER   SORSNM(4)
* 
          SNAMFL = 'PNTING::16' 
          TWOTRY = .FALSE.
* 
  100      NSRSOL = NSORS 
* 
             Assign 110 to label
                NSTMT = 110 
  110        WRITE (NAMSRS, '(4A2)', iostat=ierr, err=820)
     *               (SORSNM(LS), LS=1,4) 
* 
C           PRINT '("Change no of source ",(A8)," obs no ",(I4))',
C    *                     NAMSRS, NSRSOL 
* 
               Assign 120 to label
                 NSTMT = 120
  120     OPEN (30, file=SNAMFL, iostat=ierr, err=820, status='OLD')
* 
           Assign 130 to label
             NSTMT = 130
        Do while (NAMSRS .NE. FILSNM) 
  130   If (SNAMFL .NE. 'PNTING::16') then
           READ (30, 1300, iostat=ierr, err=820, end=200) FILSNO, IBLNK,
     *                              FILSNM, RH, RM, RS
* 
            RAHR = RH + (RM + RS/60.)/60.        ! RA in decimal hours
* 
         else 
           READ (30,1400, iostat=ierr, err=820, end=200) FILSNO, IBLNK, 
     *                                              FILSNM, SRA, SDEC 
            RAHR = SRA/15.
         endif
* 
        end do
* 
C           PRINT '("Source ",(A8)," Observing no and file no resp",
C    *                 (2I6))', FILSNM, NSRSOL, FILSNO
* 
              NSORS = FILSNO
            GO TO 900 
* 
  200    If (.NOT. TWOTRY) then 
C          PRINT '("No match for source ",(A8)," observing no. ",(I4),
C    *             " in file ",(A10))', 
C    *                NAMSRS, NSRSOL, SNAMFL
             SNAMFL = 'POINTC::16'
             TWOTRY = .TRUE.
           GO TO 100
         else 
              PRINT '("No match for either file. Enter RA of source ",
     *                 (A8),"in hh, mm _")', NAMSRS 
              READ *, IHS, IMS
              RAHR = FLOAT(IHS) + (FLOAT(IMS))/60.
           GO TO 900
         endif
* 
*          Error routine
  820         CALL ERMSG (ierr, NSTMT, *900,k)
                 GO TO label
* 
  900        CLOSE (30) 
* 
            RETURN k
* 
 1300   FORMAT (I3,A2,A8, F2.0, F3.0, F5.1) 
 1400   FORMAT (I3, A2, A8, 2F9.3)
        END 
* 
* 
**********************************************************************
* 
*        SUBROUTINE TO CONVERT BETWEEN SAST, UT AND ST
* 
**********************************************************************
* 
       SUBROUTINE TMCON(LU,ITYPE,TIMIN,HMIN,ID,IYR,KTPOUT,TIMOUT,HMOUT) 
* 
       DOUBLE PRECISION TIMIN, TIMOUT, BASEOF, RATE, OFFJD, SLONG,
     *                  FAC1, FAC2, FAC3, FAC4, RJDNT0, SASTIM, UTIME,
     *                  STIME, GAST, DUMY, DECML, CHHMM 
* 
        LOGICAL AOK, LEAPYR 
* 
        INTEGER HMIN, HMOUT 
* 
C       REAL    JDNT0 
* 
        DATA  BASEOF / 84 806.882 /, LSTDAY / 0 /,
     *        RATE / 0.002 737 909 3 /, 
     *        OFFJD / 7051.0 /,      SLONG / -1.8456667 /,
     *        FAC1  / 6.67170278/,      FAC2  / 0.0657098232 /, 
     *        FAC3  / 2433282.5 /,      FAC4  / 2445334.5 / 
* 
*        ITYPE is form of time coming IN.   SAST ..... 1
*                                             UT ..... 2
*                                             ST ..... 3
*        KTPOUT is reqd form of OUTPUT time.
*                                      Values as for ITYPE
* 
*        HMIN,HMOUT are set to 1 for input of form hh.mm, 
*                              0           in decimal hrs.
* 
*     Calculations done in decimal form 
* 
           IF ( HMIN .EQ. 1 )    TIMIN = DECML(TIMIN) 
* 
*       Find Julian day no.; modified 
* 
                  LEAPYR = ( MOD(FLOAT(IYR)/4.,1.) .EQ. 0.0 ) 
                     If (LEAPYR) then 
                               LPY = 1
                             else 
                               LPY = 2
                     endif
* 
                      If (ID .LE. (61-LPY)) then
                          MNB = 0 
                         else 
                          MNB = 1 
                      endif 
* 
               NTYR = IYR - 1960
* 
                 RJDNT0 = FLOAT( 367*NTYR - IFIX(7.* FLOAT(IYR+MNB)/4.) 
     *                         + ID + LPY*MNB + 30) + OFFJD 
* 
*        Do calcs depending on input time type. 
* 
           GO TO (100,200,300)  ITYPE 
* 
  100   SASTIM = TIMIN
              CALL LIM24(SASTIM, TMULT) 
* 
         UTIME = SASTIM - 2.
               CALL LIM24(UTIME, TMULT) 
* 
              GO TO 210 
* 
*          UT input 
* 
  200      UTIME = TIMIN
                    CALL LIM24(UTIME, TMULT)
* 
           SASTIM = UTIME + 2.
                     CALL LIM24(SASTIM, TMULT)
* 
  210      GAST = FAC1 + FAC2*RJDNT0 + (1.+RATE)*UTIME  
           STIME = GAST - SLONG 
                   CALL LIM24(STIME, TMULT) 
* 
       GO TO 400
* 
*          If ST entered
* 
  300    STIME = TIMIN
                  CALL LIM24(STIME, TMULT)
* 
*         Green app. ST 
* 
  310          GAST = STIME + SLONG 
* 
C      WRITE (LU,3200) SLONG, GAST
C3200  FORMAT(" Gren app ST from longitude ",F10.7,"is :",F10.4)
* 
          DUMY = FAC2*RJDNT0  
        UTIME = (GAST - FAC1 - DUMY) / (1. + RATE)
* 
             If (ABS(UTIME) .GT. 24.) then
                  CALL LIM24 (UTIME, TMULT) 
                  STIME = STIME + 24.*TMULT 
                    GO TO 310 
             endif
                CALL LIM24(STIME, TMULT)
* 
  330   SASTIM = UTIME + 2. 
                   CALL LIM24(SASTIM, TMULT)
* 
  400       CONTINUE
* 
*         Write 
* 
C         WRITE (LU,4100) SASTIM, UTIME, STIME
C4100     FORMAT(/5X,"SAST", 6X,"UT", 6X,"ST",/ 
C    *               3(3X,F6.2))
* 
* 
*        Output time form as reqd. Set by value of KTPOUT 
* 
             GO TO (511,522,533) KTPOUT 
* 
  511      TIMOUT = SASTIM
               GO TO 900
* 
  522      TIMOUT = UTIME 
               GO TO 900
* 
  533      TIMOUT = STIME 
* 
  900  IF ( HMOUT .EQ. 0 )  GO TO 910 
* 
*        Else convert to hh.mm
* 
          ATIM = TIMOUT 
          TIMOUT = CHHMM(ATIM)
* 
  910  RETURN 
       END
* 
**********************************************************************
* 
*         SUBROUTINE TO GET  0 < TIME < 24
* 
**********************************************************************
* 
       SUBROUTINE LIM24 (XTIME, TMULT)
* 
       DOUBLE PRECISION XTIME 
* 
          TMULT = 0.0 
* 
  110      IF ( XTIME .LT. 24. )  GO TO 120 
                XTIME = XTIME - 24. 
                TMULT = TMULT + 1.
                GO TO 110 
* 
  120      IF ( XTIME .GT. 0.0 )  GO TO 130 
                XTIME = XTIME + 24. 
                TMULT = TMULT + 1.
                GO TO 120 
  130   RETURN
       END
* 
* 
*********************************************************************** 
* 
*            Subroutine to calculate HOUR ANGLE from decimal date 
* 
*********************************************************************** 
* 
         SUBROUTINE HANGL (SORSNM, NSORS, DATE, HA, DATF10, FSTIM,
     *                     STOBS, k,*)
* 
           DOUBLE PRECISION UTOBS, STOBS, UTCALC, DATE, DATNW 
* 
           LOGICAL    DATF10, FSTIM 
* 
           INTEGER    SORSNM(4) 
* 
*           Date into year,day,UT 
* 
             IYEAR = IDINT(DATE)  
             DAY = 365. * (DATE - FLOAT(IYEAR)) 
             IDAY = IFIX(DAY) 
             IYEAR = IYEAR + 1900 
             UTOBS = 24. * (DAY - FLOAT(IDAY))
* 
C       If (DATF10)  then 
* 
*           Convert UT to ST
* 
C            CALL TMCON (LU, 2, UTOBS, 0, IDAY, IYEAR,
C    *                       3, STOBS, 0) 
C              CALL LIM24 (STOBS,TMULT) 
* 
C        else 
* 
         If (.NOT. FSTIM) then
             PRINT '("Enter ST of obs of ",I4,2X,4A2," _")', NSORS, 
     *                                             (SORSNM(L),L=1,4)
             READ *, SHR, SMIN
* 
               STOBS = SHR + SMIN/60. 
          endif 
* 
*             Calculate UT from ST to get better precision for UT 
* 
                 CALL TMCON (LU, 3, STOBS, 0, IDAY, IYEAR,
     *                           2, UTCALC, 0)
                    CALL LIM24 (UTCALC,TMULT) 
* 
               If (UTOBS .LT. 2. .AND. UTCALC .GT. 22.) IDAY = IDAY - 1 
* 
*            New decimal year from this UT
          DATNW = (UTCALC/24. + DBLE(IDAY))/365. + DBLE(IYEAR-1900) 
          DATE = DATNW
* 
C        endif
* 
*         Source coords 
* 
             CALL NAMNO (NSORS, SORSNM, RAHR, *900,k) 
* 
* 
          HA = STOBS - RAHR 
            IF ( HA .GT. 12. ) HA = HA - 24.
            IF ( HA .LT. -12.) HA = HA + 24.
* 
           HA = HA * 15.                       ! HA in degrees
* 
C          PRINT '("UTOBS, UTCALC, ST, RA, HA",5(2X,F8.6),/,
C    *              "Filed date , New date",2(2X,F12.6))',
C    *                  UTOBS, UTCALC, STOBS, RAHR, HA, DATE, DATNW 
  900    RETURN k 
         END
********************************************************************* 
* 
*      FUNCTION TO GIVE DECIMAL TIME
* 
********************************************************************* 
* 
       FUNCTION DECML (ATIM)
* 
        DOUBLE PRECISION ATIM, DECML
* 
           NTIM = IDINT(ATIM) 
           ATIM = NTIM + ((ATIM - NTIM) * 100./60.) 
           DECML = ATIM 
* 
       RETURN 
       END
* 
********************************************************************* 
* 
*             FUNCTION TO CONVERT FROM DECIMAL TO HH.MM 
* 
********************************************************************* 
* 
       FUNCTION CHHMM (ATIM)
* 
        DOUBLE PRECISION ATIM, CHHMM
* 
         NTIM = IDINT(ATIM) 
         ATIM = NTIM + ((ATIM - NTIM) * 60./100.) 
         CHHMM = ATIM 
* 
       RETURN 
       END
* 
*********************************************************************** 
* 
*           Subroutine to calculate the weighted mean from arrays of
*             value and error. Also gives arithmetic mean.
*
***********************************************************************
*
         SUBROUTINE WMEAN (VALUE, ERROR, NDATA, WTMEAN, WTERR, AVERAG,
     *                     ARERR, DEVN)
*
           PARAMETER (MXRD=150)
*
           REAL   VALUE(MXRD), ERROR(MXRD), WEIGHT(MXRD), WTVAL(MXRD)
*
        IF (NDATA .EQ. 0) GO TO 900
*
*        Set sum values to zero
*
             SUMWT = 0.
             SUMWV = 0.
*
                  SUMVAL = 0.
                  SUMSQS = 0.
                  SMSQER = 0.
*
         DO 200  J = 1,NDATA
*
              WEIGHT(J) = 1./(ERROR(J)*ERROR(J))
              WTVAL(J)  = WEIGHT(J) * VALUE(J)
*
              SUMWV = SUMWV + WTVAL(J)
              SUMWT = SUMWT + WEIGHT(J)
*
                  SUMVAL = SUMVAL + VALUE(J)
                  SUMSQS = SUMSQS + (VALUE(J) * VALUE(J))
                  SMSQER = SMSQER + (ERROR(J) * ERROR(J))
*
  200    CONTINUE
*
             WTMEAN = SUMWV/SUMWT
             WTERR  = SQRT(1./SUMWT)
* 
                  RNDATA = NDATA
                  AVERAG = SUMVAL / RNDATA
                  AVGSQD = AVERAG * AVERAG
                  AVSMSQ = SUMSQS / RNDATA
                  DIFF = AVSMSQ - AVGSQD
                     If (DIFF .GE. 0.) then 
                  DEVN = SQRT (DIFF)
                     else 
                  DEVN = 0.0
                     endif
                  DEVN = DEVN / SQRT(RNDATA)
                  ARERR = (SQRT(SMSQER)) / RNDATA 
* 
  900    RETURN 
         END
* 
********************************************************************* 
* 
      BLOCK DATA
* 
         INTEGER SORSNM 
* 
         COMMON / FLUXS / RDFLX(11), ERFLX(11), WAVE, 
     *                    SORSNM(4) 
* 
      END 
* 
********************************************************************* 
* 
*          SUBROUTINE TO ARRANGE VALUES IN ASCENDING ORDER AND
*               CALCULATE THE MEDIAN. 
*                                 Uses REAL variables. 
* 
********************************************************************* 
* 
         SUBROUTINE MEDIN (RDVM, ICOUNT, MEDVM, DIGVM)
* 
          REAL    DIGVM(11), MEDVM  
* 
           DO 110  J = 1,ICOUNT 
* 
*           Rise thru' array until position of value is found.
               IF ( RDVM .LE. DIGVM(J) )  GO TO 130 
* 
  110    CONTINUE 
* 
           DIGVM(ICOUNT) = RDVM 
                 GO TO 900
* 
*          Shift all subsequent values up the array.
  130        DO 140  K = ICOUNT,J,-1
                     DIGVM(K+1) = DIGVM(K)
  140        CONTINUE 
* 
           DIGVM(J) = RDVM
* 
  900      CONTINUE 
* 
           MIDL = ICOUNT/2
           MEDVM = DIGVM(MIDL)
* 
         RETURN 
      END 
      END$
