FTN77 
$FILES 1,2
      Program MAPR3 (1,10)
C     ============= 
C 
C     Make a map of the sun/moon/etc using the PC-Steer control system. 
C 
C     USAGE:    XQ MAPR3 <lu> <object> <suppress_plot?> 
C 
C               where:   <lu>             is the LU of the log terminal 
C                        <object>         is the object to be mapped
C                        <suppress_plot?> is <>0 to suppress the plotter
C 
C 
C     MODIFICATIONS :  1992-04-29 JQ
C                    * Now assumes a certain PC-Steer Feed System 
C                      arrangement in SysName, NDcode, DScode, TND
C                      Tilted etc. and performs automatic setup 
C                      accordingly. 
C                    * Checks consistency of Feed System definition 
C                      against PC-Steer and also compares clocks. 
C                    * Will run to System console by default but can be 
C                      set to the Esprit by typing XQ,SUN6,7
C 
C                      1992-06-25 JQ
C                    * Now suggests more intelligent defaults for scan
C                      spacing and rate for each Feed system. 
C                    * The map need no longer be square and the centre
C                      may also be offset in both Ecliptic Long and Lat.
C                    * PC-Steer Errors now cause the program to suspend 
C                      and ask for user intervention. 
C 
C    EXTENDED VERSION  1993-02-16 JQ   SUN7 -> MAPR1
C                    * Extended to allow mapping of objects other than
C                      the Sun, in particular the Moon. 
C 
C                      1993-03-18 JQ   MAPR1 -> MAPR2 
C                    * Now saves and retrieves all suitable parameters
C                      in/from parameter file MAPPRM::16, which is
C                      created but not maintained internally. 
C 
C                      1994-03-16 JQ   MAPR2 -> MAPR3 
C                    * Now estimates ecliptic longitude tracking errors 
C                      and outputs real track number including rough
C                      average longitude built in.
C 
C                      1994-03-28 JQ
C                    * Now scans in opposite direction in ecliptic coords 
C                      during first six months of year to reduce HA scanning
C                      effects
C 
C                      1995-06-30 JQ
C                    * Added a second IPRAM which suppresses use of the 
C                      HP9872S plotter when non zero
C 
C                      1997-03-07 JQ
C                    * Now allows differing scan & bin spacings 
C 
      Implicit  None
C 
C   Remember to update DATA statements below if MaxSys or Other is changed
C 
      Integer   Direction, MaxSys, MaxBin, Sun, Moon, Other, NCHAR
C 
      Parameter ( MaxSys=10 , MaxBin=1001, Sun=1, Moon=2, Other=3 ) 
C 
      Logical   Tilted(MaxSys), FilExist, FAULT, Found, LEAPYR, 
     *          Plot
C 
      Character SysName(MaxSys)*14, Reply*166, junk*1, FILNAM*10, 
     *          SysonPC*14, Names(Other)*15, FILPRE(Other)*2, 
     *          Coordtype*10, Long*5, Lat*5, Answer, Today*9, 
     *          Tomorrow*9, Type*15 
C 
      Integer*2 IPRAM(5), ITIME(5), IYEAR, UThrs, UTmin, UTsec
C 
      Integer   Ad, Am, As, UNcmd(2), Lbuf(500), AntennaLU, FLukeLU,
     *          PlotLU, FeedSys, Nscan, iscan, Count(MaxBin), Nsys, 
     *          NDcode(MaxSys), DScode(MaxSys), I, Nfile, iSrate, 
     *          ISWIT, SWIN, Istat, Nbins, Object, DyTm, YrTm 
C 
      Real      Atoday, Atomorrow, ScanRate, SpacLat, A0obj, A1obj, 
     *          DeltaT, Tscale, Zscale, Out(MaxBin), DVMmax, Base,
     *          FOCUS, TILT, Rdvm, TND(MaxSys), FracDay, SizeLong,
     *          SizeLat, A0off, B0off, SpacDef(MaxSys), a(6), b(6), 
     *          p(5), B0obj, Prlax, Old_DeltaT, ScanLong(MaxBin), 
     *          SpacLong
C 
      Double Precision MJDPC, MJDHP 
C 
      Data      AntennaLU /39/, FlukeLU /35/, PlotLU /31/ 
C 
C            Feed System names : must agree with those in PC-STEER
C 
      Data      SysName / '2.5cm Single  ',  '3.6cm Beam A  ',
     *                    '3.6cm Beam B  ',  '4.5cm Single  ',
     *                    '6.0cm Beam A  ',  '6.0cm Beam B  ',
     *                    '6.0cm Single  ',  '13.0cm Single ',
     *                    '18.0cm Single ',  '18.0cm Untilt ' / 
C 
C            Noise Diode select codes to be output to DIO channel 3 
C 
      Data      NDcode  / 2, 1, 1, 2, 3, 3, 2, 4, 5, 5 /
C 
C            DVM select codes to be output to DIO channal 5 
C 
      Data      DScode  / 6, 1, 1, 6, 3, 3, 3, 5, 7, 7 /
C 
C            Equivalent Temperature of Noise Diodes in Kelvin 
C 
      Data      TND     / 158.0,  10.0,  10.0,  82.0,  1.54,
     *                     1.54,  16.0,  3.96,  33.0,  33.0 / 
C 
C            Required Sub-Reflector position i.e. Tilted ? Y/N
C 
      Data      Tilted  / .FALSE., .FALSE., .FALSE., .FALSE., 
     *                    .FALSE., .FALSE., .FALSE., .FALSE., 
     *                    .FALSE.,  .TRUE.  / 
C 
C            Suitable scan spacing ( approximately FWHM/2 of beam ) 
C 
      Data      SpacDef / 0.030, 0.050, 0.050, 0.060, 0.080,
     *                    0.080, 0.080, 0.160, 0.250, 0.250 / 
C 
C            Total Number of usable Feed Systems
C 
      Data      Nsys    / 10 /
C 
C            Default output filename ( 000 replaced by Day Number ) 
C 
      Data      FILNAM   / '  000A::15' /,
     *          FILPRE   / 'SU', 'MN', 'OB' /,
     *          Today    / 'd000,0000' /, 
     *          Tomorrow / 'd000,0000' /
C 
C            Object Descriptors ( Another object replaced later ) 
C 
      Data      Names / 'The Sun        ',
     *                  'The Moon       ',
     *                  'Another object ' / 
C 
      Call LGBUF (Lbuf,250) 
C 
C     recall system console LU from parameter list. Defaults to 1 
C 
      Call RMPAR (IPRAM)
      Call FSYSU (IPRAM(1),IPRAM(1))
      Object = IPRAM(2) 
      Plot = (IPRAM(3).EQ.0)
      Print '(a)',CHAR(12)
C 
C     check PC vs HP clock and obtain Day Number and Year 
C 
      Call EXEC (11,ITIME,IYEAR)
      Read (AntennaLU,'(27x,f13.6)') MJDPC
      Call UNLT 
      Call MJDAY (IYEAR,ITIME,MJDHP)
C 
      If ( DABS(MJDPC-MJDHP) .gt. 0.58D-5 ) then
C 
C        clocks differ by more than half a second - warn user 
C 
         Print *,' ' 
         Print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 
         Print *,'!!! PC-STEER & HP1000 Clocks do not agree !!!' 
         Print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 
         Print *,' ' 
      endif 
C 
      Print '(" General Purpose Mapping Program MAPR3 Day ",I3.3, 
     *        ", ",I4)', ITIME(5), IYEAR
      Print *, '==================================================='
      Print *, '              Version <970701.1101>'
      If (Plot) then
         Print *, 'The HP9872S Plotter MUST be connected to the HPIB' 
         Print *, '    of the CONTROL computer during this program.'
      else
         Print *,'Use of HP9872S Plotter disabled for this run' 
      endif 
      Print *, ' '
      Print *, 'Ephemeris Ref : "The Astronomical Almanac" Sec. C,D&K'
      Print *, ' '
C 
10101 If (Object .gt. 0 .and. Object .le. Other) go to 20202
         Object = Sun 
         Print *, ' ' 
         Print *, 'Available Targets for mapping :' 
         Print *, ' ' 
         DO I=1,Other/2*2,2 
             Print *,I,' : ',Names(I),'     ',I+1,' : ',Names(I+1)
         enddo
         If (Other/2*2 .lt. Other)
     *       Print *,Other,' : ',Names(Other) 
         Print *, ' ' 
         Print *,'Enter Object (1-',Other,' / =',Object,'): _'
         Read *, Object 
      GO TO 10101 
C 
20202 If (Object .eq. Other) then 
         Print *, ' ' 
         Print *, 'Enter Object Name ( 15 chars ): _' 
         Read '(a15)', Names(Object)
      endif 
C 
C     set up default file name using Day Number and Object Prefix 
C 
      Write (FILNAM(1:2),'(A2)') FILPRE(Object) 
      Write (FILNAM(3:5),'(I3.3)') ITIME(5) 
C 
      DyTm = ITIME(5) + 1 
      YrTm = IYEAR
      If ( (.not. LEAPYR(IYEAR) .and. DyTm .gt. 365)
     *                          .OR. (DyTm.GT.366) ) then 
         DyTm = 1 
         YrTm = IYEAR + 1 
      endif 
      Write (Today(2:4),'(I3.3)') ITIME(5)
      Write (Today(6:9),'(I4.4)') IYEAR 
      Write (Tomorrow(2:4),'(I3.3)') DyTm 
      Write (Tomorrow(6:9),'(I4.4)') YrTm 
C 
      Print *, ' '
      Print *, 'Mapping ', Names(Object)
C 
C     Default direction is as before
C 
      Direction = 1 
C 
      If (Object .eq. Sun ) then
C 
         Coordtype = 'Ecliptic' 
         Long      = 'Long.'
         Lat       = 'Lat. '
C 
         Call Set_Type ( Type, 'SUN', Coordtype, Long ) 
C 
         Print *, ' ' 
         Call Get_Parm ( Type, ITIME(5), IYEAR, 1, Atoday, Found )
         If ( Found ) then
            Print *, 'Eclip. Long. of Sun at 0h DT today ', 
     *           Today, ' = _'
            Print '(F9.5)', Atoday
         else 
            Print *, 'Enter Eclip. Long. of Sun for 0h DT today ',
     *           Today, ' (dd mm ss): _'
            Read *, Ad, Am, As
            Atoday = Ad + Am / 60.0 + As / 3600.0 
            Call Put_Parm ( Type, ITIME(5), IYEAR, 1, Atoday )
         endif
C 
         Print *, ' ' 
         Call Get_Parm ( Type, DyTm, YrTm, 1, Atomorrow, Found )
         If ( Found ) then
            Print *, 'Eclip. Long. of Sun at 0h DT tomorrow ',
     *           Tomorrow, ' = _' 
            Print '(F9.5)', Atomorrow 
         else 
            Print *, 'Enter Eclip. Long. of Sun for 0h DT tomorrow ', 
     *           Tomorrow, ': _'
            Read *, Ad, Am, As
            Atomorrow = Ad + Am / 60.0 + As / 3600.0
            Call Put_Parm ( Type, DyTm, YrTm, 1, Atomorrow )
         endif
C 
         B0obj = 0.0
         Prlax = 0.0
C 
C        If DOY < 174 ie. before 23 June, reverse scanning direction
C 
         Direction = -1 
C 
      elseif (Object .eq. Moon) then
C 
         Coordtype = 'Apparent' 
         Long      = 'RA '
         Lat       = 'DEC'
C 
         Call Set_Type ( Type, 'MOON', CoordType, Long )
C 
         Print *, ' ' 
         Call Get_Parm ( Type, ITIME(5), IYEAR, 6, a, Found ) 
         If ( Found ) then
            Print *, 'RA  poly. coef. of Moon today ',Today,' are :'
            Write (*,'(5(2X,F11.7,/),2X,F11.7)') (a(I),I=1,6) 
         else 
            Print *, 'Enter RA  poly. coef. of Moon for today ',
     *           Today, ' : _'
            Read *, (a(I),I=1,6)
            Call Put_Parm ( Type, ITIME(5), IYEAR, 6, a ) 
         endif
C 
         Call Set_Type ( Type, 'MOON', CoordType, Lat ) 
C 
         Print *, ' ' 
         Call Get_Parm ( Type, ITIME(5), IYEAR, 6, b, Found ) 
         If ( Found ) then
            Print *, 'DEC poly. coef. of Moon today ',Today,' are :'
            Write (*,'(5(2X,F11.7,/),2X,F11.7)') (b(I),I=1,6) 
         else 
            Print *, 'Enter DEC poly. coef. of Moon for today ',
     *           Today, ' : _'
            Read *, (b(I),I=1,6)
            Call Put_Parm ( Type, ITIME(5), IYEAR, 6, b ) 
         endif
C 
         Call Set_Type ( Type, 'MOON', CoordType, 'PRLX' )
C 
         Print *, ' ' 
         Call Get_Parm ( Type, ITIME(5), IYEAR, 5, p, Found ) 
         If ( Found ) then
            Print *, 'Parallax poly. coef. for today ',Today,' are :' 
            Write (*,'(4(4X,F10.8,/),4X,F10.8)') (p(I),I=1,5) 
         else 
            Print *, 'Enter Parallax poly. coef. for today ', 
     *           Today, ' : _'
            Read *, (p(I),I=1,5)
            Call Put_Parm ( Type, ITIME(5), IYEAR, 5, p ) 
         endif
C 
      else
C 
30303    Print *, ' ' 
         Print *, 'Enter Co-ord Type : (G)alactic, (B)1950, (M)ean ', 
     *            'or (J)2000: _' 
         Read *, Answer 
         Call Upcase(Answer)
         If (INDEX('GBMJ',Answer) .eq. 0) go to 30303 
C 
         If (Answer .eq. 'G') then
            Coordtype = 'Galactic'
            Long      = 'Long.' 
            Lat       = 'Lat. ' 
         elseif (Answer .eq. 'B') then
            Coordtype = 'B1950' 
            Long      = 'RA ' 
            Lat       = 'DEC' 
         elseif (Answer .eq. 'M') then
            Coordtype = 'Mean'
            Long      = 'RA ' 
            Lat       = 'DEC' 
         elseif (Answer .eq. 'J') then
            Coordtype = 'J2000' 
            Long      = 'RA ' 
            Lat       = 'DEC' 
         endif
C 
         Print *, ' ' 
         Print *, 'Enter ',  CoordType(1:NCHAR(Coordtype)), 
     *            ' ',            Long(1:NCHAR(Long)),
     *            ' of ',Names(Object)(1:NCHAR(Names(Object))),'_'
         If (Long .eq. 'RA ') then
            Print *, '(hh mm ss): _'
         else 
            Print *, '(dd mm ss): _'
         endif
         Read *, Ad, Am, As 
         A0obj = Ad + Am / 60.0 + As / 3600.0 
         If (Long .eq. 'RA ') A0obj = 15.0 * A0obj
C 
         Print *, ' ' 
         Print *, 'Enter ',  Coordtype(1:NCHAR(Coordtype)), 
     *            ' ',             Lat(1:NCHAR(Lat)), 
     *            ' of ',Names(Object)(1:NCHAR(Names(Object))), 
     *            ' (dd mm ss): _'
         Read *, Ad, Am, As 
         B0obj = Ad + Am / 60.0 + As / 3600.0 
C 
         Prlax = 0.0
C 
      endif 
C 
      Print *, ' '
      Call Get_Parm ( 'TAI-UTC', 0, IYEAR, 1, Old_DeltaT, Found ) 
      If (Found) then 
         DeltaT = Old_DeltaT
         Print '(1X,"Enter the difference TAI-UTC from page K9 ", 
     *           "( / = ",I2," ): _")', INT(DeltaT) 
         Read *, DeltaT 
         Found = (Old_DeltaT .eq. DeltaT) 
      else
         DeltaT = 0.0 
         Print *, 'Enter the difference TAI-UTC from page K9: _'
         Read *, DeltaT 
      endif 
      If (.not. Found)
     *   Call Put_Parm ( 'TAI-UTC', ITIME(5), IYEAR, 1, DeltaT )
C 
C     Allow for 32.184s offset between TAI and TDT of ephemerides 
C 
      DeltaT = DeltaT + 32.184
C 
40404 FeedSys = 1 
      Print *, ' '
      Print *, 'PC-Steer Feed systems as implemented here :'
      Print *, ' '
      DO I=1,Nsys/2*2,2 
          Print *,I,' : ',SysName(I),'     ',I+1,' : ',SysName(I+1) 
      enddo 
      If (Nsys/2*2 .lt. Nsys) 
     *    Print *,Nsys,' : ',SysName(Nsys)
      Print *, ' '
      Print *, 'Enter Feed System (1-8 / =',FeedSys,'): _'
      Read *, FeedSys 
      If (FeedSys .le. 0 .or. FeedSys .gt. Nsys) go to 40404
C 
C     check subreflector position 
C 
      FAULT = .FALSE. 
      CALL RHYP (FOCUS, TILT) 
      If (.not. Tilted(FeedSys)) then 
         If ( FOCUS .lt. 7.4 ) then 
            PRINT *,' ' 
            PRINT *,'Hyperbola focus = ',FOCUS,' not 8.0' 
            FAULT = .TRUE.
         endif
         If ( TILT .gt. -0.5 ) then 
            PRINT *,' ' 
            PRINT *,'Hyperbola tilt = ',TILT,' not -0.66' 
           FAULT = .TRUE. 
         endif
      else
         If ( FOCUS .gt. 5.1 ) then 
            PRINT *,' ' 
            PRINT *,'Hyperbola focus = ',FOCUS,' not 5.0' 
            FAULT = .TRUE.
         endif
         If ( TILT .lt.  9.7 ) then 
            PRINT *,' ' 
            PRINT *,'Hyperbola tilt = ',TILT,' not 9.80'
            FAULT = .TRUE.
         endif
      endif 
      If (FAULT) goto 40404 
C 
      SizeLong = 1.0
      Print *, ' '
      Print '(1X,"Enter Map Size in ",A," ",A," (deg / = ",F4.2,
     *           "): _")', Coordtype(1:NCHAR(Coordtype)), 
     *                          Long(1:NCHAR(Long)), SizeLong 
      Read *, SizeLong
C 
      SizeLat = SizeLong
      Print *, ' '
      Print '(1X,"Enter Map Size in ",A," ",A," (deg / = ",F4.2,
     *           "): _")', Coordtype(1:NCHAR(Coordtype)), 
     *                           Lat(1:NCHAR(Lat)), SizeLat 
      Read *, SizeLat 
C 
      A0off = 0.0 
      Print *, ' '
      Print '(1X,"Enter Map Centre ",A," offset  (deg / = ",F4.2, 
     *           "): _")', Long(1:NCHAR(Long)), A0off 
      Read *, A0off 
C 
      B0off = 0.0 
      Print *, ' '
      Print '(1X,"Enter Map Centre ",A," offset  (deg / = ",F4.2, 
     *           "): _")', Lat(1:NCHAR(Lat)), B0off 
      Read *, B0off 
C 
      SpacLong = SpacDef(FeedSys) 
      Print *, ' '
      Print '(1X,"Enter Scan Spacing (deg / = ",F4.3,"): _")',SpacLong
      Read *, SpacLong
C 
      SpacLat=SpacLong
      Print *, ' '
      Print '(1X,"Enter Bin Spacing (deg /= ",F4.3,"): _")',SpacLat 
      Read *, SpacLat 
C 
      ScanRate = SpacLat * 1.2
      Print *, ' '
      Print '(1X,"Enter Scan Rate (deg/sec / = ",F4.3,
     *           "): _")',ScanRate
      Read *, ScanRate
C 
      Zscale = 7.5
      Print *, ' '
      Print *, 'Enter central height (cm / =',Zscale,'): _' 
      Read *, Zscale
C 
C     convert to Plot-Units ( normalised to sun after calibration ) 
C 
      Zscale = Zscale * 400.0 
C 
C     set Fluke to required sample rate 
C 
      iSrate = 1
      Print *, ' '
      Print *,'Enter DVM sampling rate ( 0=S, 1=M, 2=F, / =',iSrate,
     *        '): _'
      Read *, iSrate
      Write (FlukeLU,'(a,i1)') '*S',iSrate
C 
C     select DVM inputs according to FeedSys
C 
      Call EXEC (2, 113B, DScode(FeedSys)*11B, 1, 5, 5) 
      Call EXEC (12,0,2,0,-1) 
C 
C     kill the TEMPS program if it is running 
C 
      Call MESSS (10HOF,TEMPS,1,10) 
C 
      Nscan = (Nint(SizeLong / SpacLong) / 2) * 2 + 1 
      SizeLong = SpacLong * (Nscan - 1) 
      Nbins = (Nint(SizeLat / SpacLat) / 2) * 2 + 1 
      If (Nbins .gt. MaxBin) then 
         Print *,' ' 
         Print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
         Print *,'!!! Nbins >',MaxBin,' : Latitude size corrected !!'
         Print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
         Print *,' ' 
         Nbins=MaxBin 
      endif 
      SizeLat = SpacLat * (Nbins - 1) 
C 
C     calculate co-ordinates of required object 
C 
      Read (AntennaLU,'(27x,f13.6)') MJDPC
      Call UNLT 
      FracDay = Dmod(MJDPC,1.0D0) + DeltaT/86400.0
      If (Object .eq. Sun) then 
C 
         If (Atomorrow .lt. Atoday) Atomorrow = Atomorrow + 360.0 
         A1obj = (Atomorrow - Atoday) / 86400.0 
         A0obj = Atoday + FracDay * A1obj * 86400.0 
C 
      elseif (Object .eq. Moon) then
C 
         A0obj = ((((a(6)+FracDay+a(5))*FracDay+a(4))*FracDay+a(3)) 
     *            *FracDay+a(2))*FracDay+a(1) 
         B0obj = ((((b(6)*FracDay+b(5))*FracDay+b(4))*FracDay+b(3)) 
     *            *FracDay+b(2))*FracDay+b(1) 
         Prlax =  (((p(5)*FracDay+p(4))*FracDay+p(3))*FracDay+p(2)) 
     *            *FracDay+p(1) 
C 
      endif 
C 
C     drive to calculated position of object ( Note : no Centre offset )
C 
      Write (AntennaLU,'(a)') 'Begin FindObj' 
      Write (AntennaLU,'(a)') 'Clear_Fifo'
      Write (AntennaLU,'(a)') 'Duration 900'
      Write (AntennaLU,'(a," ",a)') 'Coordinate_System', Coordtype
      Write (AntennaLU,'(a)') 'Pointing On' 
      Write (AntennaLU,'(a,i2)') 'Feed_System ', FeedSys
      Write (AntennaLU,'(a)') 'A0 0.0'
      Write (AntennaLU,'(a)') 'A1 0.0'
      Write (AntennaLU,'(a)') 'B0 0.0'
      Write (AntennaLU,'(a)') 'B1 0.0'
      Write (AntennaLU,'(a,f7.3)') 'Psi ', A0obj
      Write (AntennaLU,'(a,f7.3)') 'Phi ',-B0obj
      Write (AntennaLU,'(a)') 'Theta 0.0' 
      Write (AntennaLU,'(a,f7.3)') 'Parallax ', Prlax 
      Write (AntennaLU,'(a)') 'End' 
      Call UNLT 
C 
      Print *, ' '
      Print *, 'Driving to ',Names(Object)(1:NCHAR(Names(Object))), 
     *         ' .....' 
C 
C     check PC Selected Feed System Name against internal name
C 
      Read (AntennaLU,'(40x,a14)') SysonPC
      Call UNLT 
      If (SysonPC .ne. SysName(FeedSys)) then 
C 
C        PC-STEER Feed name differs from internal - warn user 
C 
         Print *,' ' 
         Print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
         Print *,'!!! PC-STEER Feed System list inconsistent !!!'
         Print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
         Print *,' ' 
      endif 
C 
C        dither till PC-STEER is 'On Source' or an error occurs 
C 
1        Call EXEC (12,0,2,0,-3)
         Read (AntennaLU,'(a)') Reply 
         Call UNLT
         If (Reply(1:1) .eq. '1') goto 1
C 
      If (Reply(1:1) .ne. '0') then 
         read (Reply,'(I1)') Istat
         call ERROR (Istat,FAULT) 
         If (FAULT) goto 999
         goto 1 
      endif 
C 
      Print *, ' '
      Print *, 'Tracking Source, Adjust radiometer attenuation' 
      Print *, 'Press RETURN when ready: _' 
      Read (*,'(a)') junk 
C 
C     ensure that all noise diodes are turned off 
C 
      Call EXEC(2,113B,0,1,3,3) 
      Call EXEC(12,0,2,0,-2)
C 
C     sample signal peak on the source for scaling purposes 
C 
      Read (FlukeLU,*) DVMmax 
C 
      Do iscan = -(Nscan/2) , +(Nscan/2)
C 
C       recalculate co-ordinates of Object for present time 
C 
        Read (AntennaLU,'(27x,f13.6)') MJDPC
        Call UNLT 
        FracDay = Dmod(MJDPC,1.0D0) + DeltaT/86400.0
        UThrs = INT(  Dmod(MJDPC,1.0D0)*24.0  ) 
        UTmin = INT( (Dmod(MJDPC,1.0D0)*24.0-UThrs)*60.0 )
        UTsec = INT(((Dmod(MJDPC,1.0D0)*24.0-UThrs)*60.0-UTmin)*60.0) 
        If (Object .eq. Sun) then 
C 
           A0obj = Atoday + FracDay * A1obj * 86400.0 
C 
        elseif (Object .eq. Moon) then
C 
           A0obj = ((((a(6)+FracDay+a(5))*FracDay+a(4))*FracDay+a(3)) 
     *              *FracDay+a(2))*FracDay+a(1) 
           B0obj = ((((b(6)*FracDay+b(5))*FracDay+b(4))*FracDay+b(3)) 
     *              *FracDay+b(2))*FracDay+b(1) 
           Prlax =  (((p(5)*FracDay+p(4))*FracDay+p(3))*FracDay+p(2)) 
     *              *FracDay+p(1) 
C 
        endif 
C 
C       drive to start of next scan ( Map centre offset from Object ) 
C 
        Write (AntennaLU,'(a,i2.2)') 'Begin Start', iscan+Nscan/2+1 
        Write (AntennaLU,'(a)') 'Duration 400'
        Write (AntennaLU,'(a,f7.3)') 'Psi ', A0obj + A0off
        Write (AntennaLU,'(a,f7.3)') 'Phi ',-B0obj - B0off
        Write (AntennaLU,'(a,f7.3)') 'Parallax ', Prlax 
        Write (AntennaLU,'(a,f7.3)') 'A0 ', SpacLong*iscan
        Write (AntennaLU,'(a,f7.3)') 'B0 ', -Direction*SizeLat / 2.0
        Write (antennaLU,'(a)') 'B1 0.0'
        Write (AntennaLU,'(a)') 'End' 
        Call UNLT 
        Print *,' ' 
        Print *,'Driving to scan',iscan+Nscan/2+1,'_' 
        Print '(3X,I2,2(":",I2.2)," ... _" )', UThrs, UTmin, UTsec
C 
        If (iscan .ne. (-Nscan/2)) then 
C 
C          plot and write to disc the previous scan 
C 
           CALL OutScan (iscan-1, Nscan, Nbins, Tscale, Zscale, 
     *                   SpacLong, Out, ScanLong, Count, Plot)
        else
           If (Plot) then 
C 
C             set up the plotter
C 
              Write (PlotLU,'("in;")')
              Write (PlotLU,'("sp 1;pu;pa 4000,10000")')
              Write (PlotLU,'("si 1.5,1.5;lb ",a,a)') 
     *              Names(Object)(1:NCHAR(Names(Object))), Char(3)
              Write (PlotLU,'(";pu;pa 400,2500")')
              Write (PlotLU,'("si 0.2,0.3;lb Using ",4a)')
     *              SysName(FeedSys), Char(10), Char(13), Char(3) 
           Write (PlotLU,'("lb Map Size ",F3.1," by ",F3.1," deg",3a)') 
     *              SizeLong, SizeLat, Char(10), Char(13), Char(3)
              Write (PlotLU,'("lb Scan Spacing ",F5.3," deg",3a)')
     *              SpacLong, Char(10), Char(13), Char(3) 
              Write (PlotLU,'("lb Bin Spacing ",F5.3," deg",3a)') 
     *              SpacLat, Char(10), Char(13), Char(3)
              Write (PlotLU,'("lb Scan Rate ",F5.3," deg/s",3a)') 
     *              ScanRate, Char(10), Char(13), Char(3) 
              If (UThrs .lt. 10) then 
            Write (PlotLU,'("lb ",I1,2(":",I2.2)," Day ",I3.3,I5,3a)')
     *               UThrs, UTmin, UTsec, ITIME(5), IYEAR, Char(10),
     *               Char(13), Char(3)
              else
            Write (PlotLU,'("lb ",I2,2(":",I2.2)," Day ",I3.3,I5,3a)')
     *               UThrs, UTmin, UTsec, ITIME(5), IYEAR, Char(10),
     *               Char(13), Char(3)
              endif 
           endif
C 
C          set up the standard output disc file 
C 
           FilExist=.TRUE.
           Nfile=0
           DO WHILE (FilExist)
               FILNAM(6:6)=CHAR(ICHAR('A')+Nfile) 
               Inquire (FILE=FILNAM, EXIST=FilExist)
               Nfile = Nfile+1
           enddo
           Print *,' '
           Print *,' '
           Print *,'Creating results file ',FILNAM,' ...' 
           Open (FILE=FILNAM, UNIT=20, STATUS='NEW')
           If (Plot) then 
              Write (PlotLU,'("si 0.2,0.3;lb Output File ",4a)')
     *              FILNAM, Char(10), Char(13), Char(3) 
              Write (PlotLU,'("sp 0;")')
           endif
        endif 
C 
C       wait here till drive to start of scan is complete 
C 
2          Call EXEC (12,0,2,0,-2)
           Read (AntennaLU,'(a)') Reply 
           Call UNLT
           If (Reply(1:1) .eq. '1') goto 2
C 
        If (Reply(1:1) .ne. '0') then 
           read (Reply,'(I1)') Istat
           call ERROR (Istat,FAULT) 
           If (FAULT) goto 999
           goto 2 
        endif 
C 
C 
        If (iscan .eq. (-Nscan/2)) then 
C 
C          calibrate on the first scan
C 
           If (Object .eq. Sun) then
              Print *, ' '
              Print *, 'Ready to Calibrate : Remove 20dB of attenuation'
              Print *, 'Press RETURN when ready: _' 
              Read (*,'(a)') junk 
           endif
C 
           Print *,' '
           CALL CALIBRATE (TND(FeedSys), NDcode(FeedSys), Tscale, FAULT)
           Print *, ' ' 
           If (FAULT) goto 999
C 
           If (Object .eq. Sun) then
              Print *, 'Done ! Put back 20dB attenuation' 
              Print *, 'Press RETURN when ready: _' 
              Read (*,'(a)') junk 
              Print *, ' '
              Tscale = Tscale*100.0 
           endif
C 
C          sample signal at map corner for base level 
C 
           Base = 0.0 
           DO I=1,20
               Read (FlukeLU,*) Rdvm
           Base = Base + Rdvm 
           enddo
           Base = Base / 20.0 
C 
C          Normalise central height by signal on sun
C 
           Zscale = Zscale / (DVMmax - Base)
C 
        endif 
C 
C       reselect Fluke DVM inputs 
C 
        Call EXEC (2, 113B, DScode(FeedSys)*11B, 1, 5, 5) 
        Call EXEC (12,0,2,0,-1) 
C 
C       duration extended by 5 Seconds to ensure Scan is completed
C 
        Write (AntennaLU,'(a,i2.2)') 'Begin ObjScn', iscan+Nscan/2+1
        Write (AntennaLU,'(a,f8.3)') 'Duration ', SizeLat/ScanRate+5
        Write (AntennaLU,'(a,f7.3)') 'B1 ', Direction*ScanRate
        Write (AntennaLU,'(a)') 'End' 
        Call UNLT 
        Print *,'Doing scan', iscan+Nscan/2+1,' of',Nscan,' _'
        CALL DoScan (Nbins, SizeLat, Base, Out, ScanLong, 
     *               Count, FAULT)
        If (FAULT) goto 999 
      enddo 
C 
      CALL OutScan (iscan-1, Nscan, Nbins, Tscale, Zscale,
     *              SpacLong, Out, ScanLong, Count) 
C 
      Close (20,err=900)
C 
  900 Print *, ' '
      Print *, '|======================|' 
      Print *, '|  Map scans finished  |' 
      Print *, '|======================|' 
      If (Plot) then
C 
C        Chart Advance on the Plotter 
C 
         Write (PlotLU,'("af;")') 
      endif 
C 
C     Send telescope to Zenith if bit 11 is set 
C 
      If (ISSW(11) .lt. 0) then 
         ISWIT = SWIN(0)
         ISWIT = ISWIT .AND. 173777B
         CALL SWOUT(ISWIT)
         Write (AntennaLU,'(a)') 'Begin ZENITH' 
         Write (AntennaLU,'(a)') 'Clear_Fifo' 
         Write (AntennaLU,'(a)') 'Duration 9999'
         Write (AntennaLU,'(a)') 'Coordinate_System HORIZON'
         Write (AntennaLU,'(a)') 'Pointing Off' 
         Write (AntennaLU,'(a)') 'Feed_System 0'
         Write (AntennaLU,'(a)') 'A0 0.0' 
         Write (AntennaLU,'(a)') 'A1 0.0' 
         Write (AntennaLU,'(a)') 'B0 90.0'
         Write (AntennaLU,'(a)') 'B1 0.0' 
         Write (AntennaLU,'(a)') 'Psi 0.0'
         Write (AntennaLU,'(a)') 'Phi 0.0'
         Write (AntennaLU,'(a)') 'Theta 0.0'
         Write (AntennaLU,'(a)') 'Parallax 0.0' 
         Write (AntennaLU,'(a)') 'End'
         Call UNLT
      endif 
C 
  999 Print *,'Bye !' 
      Stop
      End 
C 
C 
      SUBROUTINE CALIBRATE (TND, NDcode, Tscale, FAULT) 
C     ====================
C     To obtain the Kelvins / Volt calibration factor 
C     Modified to take out linear slope in sampling 
C 
C     inputs :    TND     real        Cal Noise Diode in Kelvins
C                 NDcode  integer     Cal Noise Diode select code in DIO
C 
C     outputs:    Tscale  real        Kelvins / Volt calibration factor 
C                 FAULT   logical     .TRUE. if unable to calibrate 
C 
      Implicit None 
C 
      Logical          FAULT
C 
      Integer          I, J, NVALID(4), NDcode, FlukeLU, Retries, Number
C 
      Real             TND, Tscale, DTScale 
C 
      Double Precision DVM(4), DVMSQD(4), ADVM(4), DMS(4), RMS(4),
     *                 RDVM, RMS_SUM, DMS_SUM, cal_COUNT
C 
      Data             FlukeLU / 35 / 
C 
      Retries = 0 
      Number = 150
      FAULT = .FALSE. 
      Print *,' System calibration :' 
      Print *,' ' 
C 
C         restart point after error  - allow for up to 10 failures
C 
   10     Retries = Retries + 1 
          DMS_SUM = 0D0 
          RMS_SUM = 0D0 
C 
C         sample with noise cal off (J = 1,4) or on (J = 2,3) 
C 
          DO J = 1, 4 
              DVM(J)  = 0D0 
              DVMSQD(J) = 0D0 
              NVALID(J)  = 0
              If (J .eq. 2) then
C 
C                 turn on the noise diode 
C 
                  Call EXEC (2,113B,NDcode,1,3,3) 
              endif 
              If (J .eq. 1 .or. J .eq. 4) then
C 
C                 turn off the noise diode
C 
                  Call EXEC (2,113B,0,1,3,3)
              endif 
              If (J .ne. 3) then
C 
C                 let the system settle 
C 
                  Call EXEC (12,0,2,0,-2) 
              endif 
C 
C             sample 0.1 sec filtered 'direct' radiometer output
C 
              DO I = 1,Number 
                  read (FlukeLU,*) RDVM 
C 
C                 Store DVM and DVM SQUARED readings
C 
                  NVALID(J) = NVALID(J) + 1 
                  DVM(J) = DVM(J) + RDVM
                  DVMSQD(J) = DVMSQD(J) + RDVM * RDVM 
              enddo 
C 
C             sampling done : calculate average DVM reading and rms error 
C 
              If (NVALID(J) .gt. 2) then
C 
C                 average DVM reading 
C 
                  ADVM(J) = DVM(J) / NVALID(J)
                  DMS(J) = (DVMSQD(J) - DVM(J)*DVM(J)/NVALID(J))
     *                      / (NVALID(J) - 1.0) 
                  DMS_SUM = DMS_SUM + DMS(J)
C 
C                 rms error in reading
C 
                  RMS(J) = DSQRT(DMS(J))
                  RMS_SUM = RMS_SUM + RMS(J)
C 
                  If (J .eq. 1 .or. J .eq. 4) then
                      Print *,'   Noise cal. off : _' 
                  endif 
                  If (J .eq. 2 .or. J .eq. 3) then
                      Print *,'   Noise cal.  on : _' 
                  endif 
                  Print '(F12.6,'' +-'',F12.6,'' V'')',ADVM(J),RMS(J) 
              endif 
          enddo 
C 
C         error trap
C 
          If (Retries .ge. 10) then 
              FAULT = .TRUE.
              Print *,' ' 
              Print *,' Abandoning calibration' 
              RETURN
          endif 
C 
C         check noise tube > 4 * average rms noise
C 
          cal_COUNT = ((ADVM(2) + ADVM(3)) - (ADVM(1) + ADVM(4)))/2 
          If (DABS(cal_COUNT) .lt. RMS_SUM) then
              Print *,' ' 
              Print *,' Noise cal not seen - trying again'
              Print *,' ' 
              go to 10
          endif 
C 
C     calculate Kelvins per DVM unit conversion factor
C 
      Tscale = TND / cal_COUNT
C 
C     fractional error in Tscale is half of the square root of the sum of 
C     the squares of the fractional errors in the four DVM averages 
C 
      DTscale = Tscale * DSQRT (DMS_SUM / (2*cal_COUNT)**2) 
      DTscale = ABS (DTscale) 
      Print *,' ' 
      Print '(1x,"Calibration gives",F12.6," +-",F8.6,
     *           " Kelvins / Volt")', Tscale, DTscale 
C 
      Return
      End 
C 
C 
      SUBROUTINE DOSCAN (Nbins, SizeLat, Base, Out, ScanLong, 
     *                   Count, FAULT)  
C     ================= 
C     To sample the DVM and PC-Steer co-ordinates during a scan 
C     and then bin and average them as required.
C 
C     inputs :    Nbins       integer     No of output bins required
C                 Sizelat     real        Eclip. Lat. length of scan
C                 Base        real        DVM output reading base level 
C 
C     outputs:    Out(*)      real        Binned averaged DVM output
C                 ScanLong(*) real        Binned averaged scan ScanLongitude
C                 Count(*)    integer     No of values averaged in bin
C                 FAULT       logical     Unrecoverable PC-Steer Error
C 
      Implicit None 
C 
      Integer  Nbins, Count(*), AntennaLU, FlukeLU, i,
     *         Istat, Nempty
C 
      Real     X, DVM, Y, Out(*), ScanLong(*), SizeLat, Base
C 
      Logical  FAULT
C 
      DATA     AntennaLU /39/, FlukeLU /35/ 
C 
      Do i = 1 , Nbins
        Count(i) = 0
        Out(i) = 0.0
        ScanLong(i) = 0.0 
      enddo 
C 
1       Read (AntennaLU,'(i1,151x,2f8.3)') Istat, X, Y
        Call UNLT 
        If (Istat .ge. 2) goto 2
        Read (FlukeLU,*) DVM
        if (X.gt.180.0) X = X - 360.0 
        i = Nint ((SizeLat/2.0 + Y) / SizeLat * Nbins)
        If (i .ge. 1 .and. i .le. Nbins) then 
           Count(i) = Count(i) + 1
           Out(i) = Out(i) + DVM-Base 
           ScanLong(i) = ScanLong(i) + X
        endif 
      goto 1
C 
2     If (Istat .ne. 2) then
         call ERROR (Istat,FAULT) 
         If (FAULT) goto 999
      endif 
C 
      Nempty=0
      Do i = 1 , Nbins
        If (Count(i) .ne. 0) then 
           Out(i) = Out(i) / Count(i) 
           ScanLong(i) = ScanLong(i) / Count(i) 
        else
           Nempty=Nempty+1
        endif 
      enddo 
      Print *,'( ',Nempty,' empty bins )' 
C 
  999 Return
      End 
C 
C 
      SUBROUTINE OutScan (iscan, Nscan, Nbins, Tscale, Zscale,
     *                    Spacing, Out, ScanLong, Count, Plot)
C     =================== 
C     To output the binned data for each scan both into the output
C     file and to the plotter in the form of a continuous line. 
C 
C     inputs :    iscan       integer     Present scan number 
C                 Nscan       integer     Total no of scans in map
C                 Nbins       integer     No of output bins in scan 
C                 Tscale      real        Kelvin/Volt calibration factor
C                 Zscale      real        Plot-Unit/Volt scale factor 
C                 Spacing     real        Scan spacing in longitude 
C                 Out(*)      real        Binned averaged DVM output
C                 ScanLong(*) real        Binned averaged scan ScanLongitude
C                 Count(*)    integer     No of values in bin 
C                 Plot        logical     Output to Plotter ? 
C 
C     outputs :   None except to file and plotter 
C 
      Implicit None 
C 
      Integer  iscan, Nscan, Nbins, Count(*), PlotLU, i 
C 
      Logical  Plot 
C 
      Real     Tscale, Zscale, Spacing, Out(*), ScanLong(*) 
C 
      DATA     PlotLU /31/
C 
      If (Plot) then
         Write (PlotLU,'("pu;")') 
         If (iscan .lt. 0) Write (PlotLU,'("sp 2;")') 
         If (iscan .eq. 0) Write (PlotLU,'("sp 3;")') 
         If (iscan .gt. 0) Write (PlotLU,'("sp 4;")') 
      endif 
C 
      Do i = 1, Nbins 
        If (Count(i) .gt. 0) then 
           If (Plot) then 
              Write (PlotLU,'("pa ",i6,",",i6,";pd;")') 
     *           Nint(FLOAT(i)/Nbins*9000+3000),
     *           Nint(Out(i)*Zscale) +
     *           FLOAT(Iscan+Nscan/2)/Nscan*6000 + 500
           endif
           Write (20,'(F8.2,'','',I5,'','',F12.5,'','',I5)')
     *          -ScanLong(i)/Spacing, i-Nbins/2-1,
     *          Out(i)*Tscale, Count(i) 
        else
           Write (20,'(I3,''.00,'',I3,'','',12X,'','',I5)') 
     *          -iscan, i-Nbins/2-1, Count(i) 
        endif 
      enddo 
      If (Plot) then
         Write (PlotLU,'("pu;")') 
         Write (PlotLU,'("sp 0;")') 
      endif 
C 
      Return
      End 
C 
C 
      Subroutine UNLT 
C     =============== 
C     To reset the PC-Steer HPIB card to UNLISTEN mode
C 
C     inputs :    None
C 
C     outputs:    None except on HPIB 
C 
      Integer Cmd(2)
C 
      DATA Cmd /2, 2H_?/
C 
      Call CMDR (30,Cmd,0)
C 
      Return
      End 
C 
C 
      Subroutine MJDAY (IYR, ITIME, MJD)
C     ================
C     To calculate the Modified Julian Day, as defined in the 
C     Astronomical Almanac, for a given year, day & time
C 
C     inputs :    IYR      Integer*2         Year ( from 1983 - ) 
C                 ITIME(5) Integer*2         Day No, Hr, Min, Sec, 100th
C 
C     outputs:    MJD      Double Precision  Modified Julian Day Number 
C 
C 
      Implicit None 
C 
      Integer*2         IYR, ITIME(5) 
C 
      Integer           JYR 
C 
      Double Precision  MJD 
C 
      Logical           LEAPYR
C 
      MJD = ( ( (DBLE(ITIME(1))/100D0 + DBLE(ITIME(2)) ) / 60.0D0 
     *                                  + DBLE(ITIME(3)) ) / 60.0D0 
     *                                    + DBLE(ITIME(4)) ) / 24.0D0 
     *                                      + DBLE(ITIME(5))
C 
      If (IYR .lt. 1983)  then
C 
        Print '(I4," - year not catered for in MJD routine MJDAY")', 
     *    IYR 
        MJD = 0.0D0 
C 
      else
C 
        MJD = MJD + 45334.0D0      ! 1983 offset
        DO 100 JYR = 1984, IYR
          MJD = MJD + 365.0D0 
          If (LEAPYR(JYR-1))  MJD = MJD + 1.0D0 
  100   CONTINUE
C 
      endif 
C 
      Return
      End 
C 
C 
      Logical Function LEAPYR (IYEAR) 
C     ======================= 
C     To determine if any given year is a leap year 
C 
C     inputs :    IYEAR   Integer  Year to be assessed
C 
C     outputs:    LEAPYR  Logical  .TRUE. if IYEAR is a leap year 
C 
      Implicit None 
C 
      Integer IYEAR 
C 
      LEAPYR = (((MOD(IYEAR,4) .eq. 0) .AND. (MOD(IYEAR,100) .ne. 0)) 
     *  .OR. (MOD(IYEAR,400) .eq. 0)) 
C 
      Return
      End 
C 
C 
      Subroutine ERROR (Istat,FAULT)
C     ================
C     To report PC-Steer Errors and request User intervention 
C 
C     inputs :    Istat   Integer  PC-Steer Error Number
C 
C     outputs:    FAULT   Logical  Condition Fatal so terminate 
C 
      Implicit None 
C 
      Integer  Istat
C 
      Character CHAR
  
      Logical  FAULT
C 
      Print *,' '
      Print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 
      Print *,'!!! PC-Steer Error',Istat,' has been detected !!' 
      Print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 
      Print *,' '
C 
      FAULT = .FALSE. 
1     Print *,'Should Map be terminated  ( <cr> = ',FAULT,' ) ? _'
      Read (*,'(a)') CHAR 
      Call Upcase(CHAR) 
      If (CHAR .eq. 'Y' .or. CHAR .eq. 'T' ) then 
         FAULT =.TRUE.
         goto 999 
      endif 
      If (CHAR .ne. 'N' .and. CHAR .ne. 'F' .and. CHAR .ne. '/' .and. 
     *    CHAR .ne. ' ' ) goto 1
C 
  999 Return
      End 
C 
C 
      Subroutine Set_Type ( Type, Object, CoordSys, Coord ) 
C     =================== 
C     To construct the appropriate Parameter Identifier from the
C     various descriptive strings to an underscore separated list.
C 
C     inputs :    Object     Character*(*)    Map Object Identifier 
C                 CoordSys   Character*(*)    Coordinate System 
C                 Coord      Character*(*)    Coordinate Identifier 
C 
C     outputs:    Type       Character*(*)    Parameter Identifier
C 
      Implicit None 
C 
      Character  Type*(*), Object*(*), CoordSys*(*), Coord*(*)
C 
      Integer    Nchar, L1, L2, L3
C 
C     evaluate lengths of three components, maximums 4_5_4
C 
      L1 = MIN( 4, Nchar( Object ) )
      L2 = MIN( 5, Nchar( CoordSys ) )
      L3 = MIN( 4, Nchar( Coord ) ) 
      If (Coord(L3:L3) .eq. '.') L3 = L3 - 1
C 
C     construct Parameter Identify as Object_CoordSys_Coord 
C 
      Type = Object(1:L1) //'_'// CoordSys(1:L2) //'_'// Coord(1:L3)
      Call Upcase ( Type )
C 
      Return
      End 
C 
C 
      Subroutine Open_File ( Filename, LU, Position, IERR ) 
C     ====================
C     To open any sequential file if not already opened and to return 
C     the LU number associated with the file.  New LU's are allocated 
C     sequentially from 7 upwards as necessary.  Note that the file 
C     will be created if it does not exist already. 
C 
C     inputs :    Filename   Character*(*)    Name of file to be opened 
C                 Position   Character*(*)    'START' or 'END' of file
C 
C     outputs:    LU         Integer          Logical Unit attatched
C                 IERR       Integer          I/O IOSTAT error number 
C 
      Implicit None 
C 
      Character Filename*(*), Position*(*), Junk
C 
      Integer   LU, Default, IERR 
C 
      Logical   Exists, OpenNow, Used 
C 
      Data      Default / 65 /
C 
C     first we must determine if the required file is open already
C 
      Call Upcase ( Position )
      Inquire ( FILE=Filename, IOSTAT=IERR, ERR=995,
     *          EXIST=Exists, OPENED=OpenNow, NUMBER=LU ) 
C 
      If (.not. OpenNow) then 
C 
C        if file is not yet opened, determine free LU to use for file 
C 
         Used = .TRUE.
         LU = Default - 1 
         Do While (Used)
            LU = LU + 1 
            Inquire ( UNIT=LU, IOSTAT=IERR, ERR=995, OPENED=Used )
         enddo
C 
C        and open required file FILENAME, creating if necessary 
C 
         Open ( FILE=Filename, UNIT=LU, IOSTAT=IERR, ERR=995, 
     *          STATUS='UNKNOWN', ACCESS='SEQUENTIAL' ) 
C 
C        if file did not exist before then write EOF marker 
C 
         If (.not. Exists) Endfile ( UNIT=LU, IOSTAT=IERR, ERR=995 )
C 
      endif 
C 
C     if Position='START' then simply rewind file and Exit
C 
      If (Position .eq. 'START') then 
         Rewind ( UNIT=LU, IOSTAT=IERR, ERR=995 ) 
         goto 999 
      endif 
C 
C     if Position='END' then rewind, locate EOF 
C 
      If (Position .eq. 'END') then 
         Rewind (UNIT=LU, IOSTAT=IERR, ERR=995 )
         Do While (.TRUE.)
            Read ( LU, *, IOSTAT=IERR, ERR=995, END=999 ) JUNK
         enddo
      endif 
C 
  990 Print '("Unknown Position ",A," in Open_File")', Position
      goto 999
C 
  995 Print '("I/O Error ",I4," on opening file ",A," to ",A)',
     *     IERR, Filename, Position 
C 
  999 Return
      End 
C 
C 
      Subroutine Get_Parm ( Type, DayNo, Year, Nvals, Values, Found ) 
C     =================== 
C     To get the values associated with Type on the day no. DayNo of
C     year Year from the MAPPRM::16 parameter file. 
C     Note: A DayNo of 0 will retrieve the most recent value. 
C 
C     inputs :    Type     Chsracter*(*)    Keyword identifying Parameter 
C                 DayNo    Integer          Julian Day Number and ... 
C                 Year     Integer          Year the values are required for
C                 Nvals    Integer          No of associated values 
C 
C     outputs:    Values   Real(*)          Values retrieved from file
C                 Found    Logical          If vslues were available
C 
      Implicit None 
C 
      Character Type*(*), FType*15
C 
      Integer   DayNo, Year, Nvals, FileLU, IERR, FDyNo, FYear, I 
C 
      Logical   Found 
C 
      Real      Values(*), SValue(20), Epoch, FEpch 
C 
C     first Open the MAPPRM::16 parameter storage file at the START 
C 
      Epoch = 0.0 
      Found = .FALSE. 
      Call Upcase ( Type )
      Call Open_File ( 'MAPPRM::16', FileLU, 'START', IERR )
      If (IERR .gt. 0) goto 995 
C 
C     read next line in file till Parameter is found or EOF detected
C 
      Do While (.NOT. Found .OR. DayNo .eq. 0)
         Read ( FileLU, *, IOSTAT=IERR, ERR=995, END=900 ) FType
C 
C        if Type matches that read from file, backspace and read values 
C 
         If (Index(FType,Type) .eq. 1) then 
            Backspace ( UNIT=FileLU, IOSTAT=IERR, ERR=995 ) 
            Read ( FileLU, *, IOSTAT=IERR, ERR=995, END=900 ) 
     *          FType, FDyNo, FYear, (Values(I),I=1,Nvals)
C 
C           if searching for latest value, compare Epoch's etc
C 
            If (DayNo .eq. 0) then
               FEpch = FYear + (FDyNo-1)/366.0
               If (FEpch .ge. Epoch .or. .not. Found) then
                  Found = .TRUE.
                  Epoch = FEpch 
                  Do I = 1, Nvals 
                     SValue(I) = Values(I)
                  enddo 
               endif
C 
C           otherwise both Day No. and Year should match
C 
            elseif (DayNo .eq. FDyNo .and. Year .eq. Fyear ) then 
               Found = .TRUE. 
            endif 
C 
         endif
      enddo 
C 
C     at EOF must copy back latest values if DAYNO = 0 and then Exit
C 
  900 If (Found .and. DayNo .eq. 0) then
        Do I = 1, Nvals 
           Values(I) = SValue(I)
        enddo 
      endif 
      goto 999
C 
  995 Print '("I/O Error ",I4," on MAPPRM::16 in Get_Parm")', IERR 
C 
  999 Return
      End 
C 
C 
      Subroutine Put_Parm ( Type, DayNo, Year, Nvals, Values )
C     =================== 
C     To put the values associated with Type on the day no. DayNo of
C     year Year into the MAPPRM::16 parameter file.  A total of Nvals 
C     values are saved from the Values array. 
C 
C     inputs :    Type     Chsracter*(*)    Keyword identifying Parameter 
C                 DayNo    Integer          Julian Day Number and ... 
C                 Year     Integer          Year the values are valid on
C                 Nvals    Integer          No of associated values 
C                 Values   Real(*)          Values to be saved on file
C 
C     outputs:    None except to file 
C 
      Implicit None 
C 
      Character Type*(*)
C 
      Integer   DayNo, Year, Nvals, FileLU, IERR, I, Nchar
C 
      Real      Values(*) 
C 
C     first Open the MAPPRM::16 parameter storage file at the END 
C 
      Call Upcase ( Type )
      Call Open_File ( 'MAPPRM::16', FileLU, 'END', IERR )
      If (IERR .gt. 0) goto 995 
C 
C     write suitable line into file followed by an EOF marker and Exit
C 
      Write ( FileLU, 100, IOSTAT=IERR, ERR=995 ) 
     *     Type(1:Nchar(Type)), DayNo, Year, (Values(I),I=1,Nvals)
  100 Format ( "'", A, "'", I4.3, I5, 20F13.8 ) 
      Endfile ( UNIT=FileLU, IOSTAT=IERR, ERR=995 ) 
      goto 999
C 
  995 Print '("I/O Error ",I4," on MAPPRM::16 in Put_Parm")', IERR 
C 
  999 Return
      End 
C 
C=======================
      INCLUDE &RHYP::16 
      INCLUDE &NCHAR::16
      INCLUDE UPCASE::16
C=======================