* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
*
      SUBROUTINE INPUT_FITS(FICH,NFILE,ERROR,DATA_TYPE)
*
      INCLUDE 'nmb.inc'
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
! Parameters
      CHARACTER*(*) FICH
      INTEGER*4 NFILE
      LOGICAL ERROR, DATA_TYPE
! Local variables
! pointers
      INTEGER*4 POINTER
      INTEGER*4 IP
      INTEGER*4 DESC(9)
      INTEGER*4 IP_SKYUT, IP_SKYTAU
*
! real*8
      REAL*8    UT1, UT2
      REAL*8    DX,OBSLONG
*
! others ...
      CHARACTER FILENAME*80, STRING*400, MONTH(12)*3, TEXT*80
      CHARACTER SOFT_VERSION*8
      CHARACTER SAVE_STRING*80, C*1, D*1
      LOGICAL   STATE, COMPARE
*
      INTEGER*4 I, SL1,SL2,SL3,SL4,SB1,SB2,SB3,SB4
      INTEGER*4 SIC_FRELUN,SIC_GETVM,IER,SIC_OPEN,SIC_GETLUN
      INTEGER*4 IDATE, J, IMONTH1, ISCAN, LUN, FOUND_SCAN
      INTEGER   NPARAM_SUBSCAN
      LOGICAL   OFF
      LOGICAL   RECORD_2PHASES, FOUND
*
      INTEGER*4 LEN_RECORD,LENGTH
      INTEGER   OLD_FLAG
      INTEGER   PNT
      REAL*4    OPACITY_GRADIENT, UT3, UT4
      INTEGER   MAX_SEG
      PARAMETER (MAX_SEG=100)
      LOGICAL   JY_STATE
      INTEGER   NC
*
      DATA MONTH /'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG',
     $'SEP','OCT','NOV','DEC'/
      DATA C /':'/
      DATA D /'.'/
*
* Close previous file
*
      IF (INPUTSTREAM.NE.0) THEN
         CALL CLOSE_STREAM(INPUTSTREAM)
         INPUTSTREAM = 0
      ENDIF
*
      IF (ILUN.GT.0) THEN
         CLOSE(UNIT=ILUN)
         IER = SIC_FRELUN(ILUN)
         ILUN = 0
      ENDIF
*
* Open new file
      FILENAME=FICH
      INQUIRE(FILE=FILENAME,EXIST=FOUND)
      IF (.NOT.FOUND) THEN
         OLD_FILENAME = ' '
         CALL MESSAGE(6,3,'INPUT_FITS',
     $   'File not found: '//FILENAME)
         ERROR =.FALSE.
         SCAN_NUMBER = -1
         RETURN
      ENDIF
! Read fits file
      STRING = 'DEFINE FITS F "'//FICH(1:NFILE)//'" read'
      CALL EXEC_ALL(STRING)
! Test if SIC structure created
      CALL SIC_DESCRIPTOR('F',DESC,FOUND)
      IF (.NOT.FOUND) THEN
         OLD_FILENAME = ' '
         CALL MESSAGE(6,3,'INPUT_FITS',
     $   'Unable to read fits file '//FILENAME)
         ERROR =.FALSE.
         SCAN_NUMBER = -1
         RETURN
      ENDIF
      OLD_FILENAME = FICH
      SCAN_NUMBER=-1
      OLD_FLAG = 0
*
      IF (IADDR_SUBFLAG.NE.0) THEN
         IP = POINTER(IADDR_SUBFLAG,MEMORY)
         CALL FREE_SUBSCAN_MEMORY(MEMORY(IPSUBREC),MEMORY(IP))
         LENGTH = NSUBSCAN*NCHAN*2       ! address array (nchan*nsubscan)
         CALL FREE_VM(LENGTH,IADDR_SUBFLAG)
         IADDR_SUBFLAG = 0
      ENDIF
      IF (IADDRRECORD.NE.0) THEN
         OLD_FLAG = NCHAN
         CALL R4TOR4(MEMORY(IP_FLAG),RECEIVER_FLAG(1),NCHAN)
         IF (DATA_SHMID.NE.0)
     $   CALL REMOVE_MEMORY(DATA_SHMID,IADDRRECORD)
         DATA_SHMID=0
         IADDRRECORD = 0
      ENDIF
      IF (IADDRSPIKE.NE.0) THEN
         LENGTH = NCHAN + 2*LENGTH_SPIKE*NCHAN
         CALL FREE_VM(LENGTH,IADDRSPIKE)
         IADDRSPIKE = 0
      ENDIF
      IF (IADDR_PHASES.NE.0) THEN
         LENGTH = 2*NCHAN*NRECORD
         CALL FREE_VM(LENGTH,IADDR_PHASES)
         IADDR_PHASES = 0
      ENDIF
*
C Read Header
*
* Header Section -----------------------------
*
*
      CALL SIC_GET_CHAR('F%SOFTVER',SOFT_VERSION,NC,ERROR)
      IF (ERROR) THEN
         CALL MESSAGE(6,3,'INPUT_FITS',
     $   'Missing SOFTVER parameter in FITS file')
         CALL SIC_DELVARIABLE('F',.FALSE.,ERROR)
         CALL MESSAGE(2,1,'INPUT_FITS','Deleting fits structure F%')
         ERROR =.FALSE.
         SCAN_NUMBER = -1
         RETURN
      ENDIF
      CALL SIC_GET_CHAR('F%PROJECT',PROJECT_NAME,NC,ERROR)
      CALL SIC_GET_CHAR('F%OBJECT',SOURCE_NAME,NC,ERROR)
      CALL SIC_GET_DBLE('F%SLAM',DX,ERROR)
      SOURCE_LAMBDA = DX                 !deg
      CALL SIC_GET_DBLE('F%SBET',DX,ERROR)
      SOURCE_BETA = DX                   ! deg
      CALL SIC_GET_CHAR('F%TELESCOP',TELESCOPE,NC,ERROR)
      CALL SIC_GET_DBLE('F%TELALTIM',DX,ERROR)
      SITE_ALTITUDE = DX                 ! (m)
      CALL SIC_GET_DBLE('F%TELLATID',SITE_LATITUDE,ERROR)
      SITE_LATITUDE = SITE_LATITUDE/180*PI     ! (rad)
      CALL SIC_GET_DBLE('F%TELLONGI',SITE_LONGITUDE,ERROR)
      SITE_LONGITUDE = SITE_LONGITUDE/180*PI   ! (rad)
      CALL SIC_GET_DBLE('F%TELSIZM',DX,ERROR)
      TELESCOPE_DIAMETER = DX            ! (m)
      DATA_FORMAT    = FORMAT_FITS
      CALL SIC_GET_DBLE('F%NSCAN',DX,ERROR)
      SCAN_NUMBER = DX
* Define scan type:
      CALL SIC_GET_DBLE('F%SRP1',DX,ERROR)
      ISRP1 = DX
      CALL SIC_GET_DBLE('F%SRP2',DX,ERROR)
      ISRP2 = DX
      CALL SIC_GET_DBLE('F%SRP1FLAG',DX,ERROR)
      ISR1OPTFLAG = DX
      CALL SIC_GET_DBLE('F%OBSMODE',DX,ERROR)
      SCAN_TYPE = DX
      CALL SIC_GET_DBLE('F%DAZM',DX,ERROR)
      DAZM = DX/3600.
      IF (SCAN_TYPE.EQ.5.AND.DAZM.GT.1E-8)    SCAN_TYPE = 1! ON-OFF
      IF (SOURCE_NAME(1:6).EQ.'SKYPAD')       SCAN_TYPE = 8! SKYPAD
      IF (SOURCE_NAME(1:9).EQ.'RASTERDIP'.AND.SCAN_TYPE.EQ.3)
     $SCAN_TYPE = 9                      ! RASTERDIP
*
      CALL SIC_GET_DBLE('F%EPOCH',DX,ERROR)
      EPOCH = DX
      CALL SIC_GET_DBLE('F%SBAS',DX,ERROR)
      COORDINATE_SYSTEM = DX
      CALL SIC_GET_CHAR('F%DATE-OBS',STRING,NC,ERROR)
      READ(STRING(1:4),'(i4)') IYEAR     ! '2000-02-07T01:41:41'
      READ(STRING(6:7),'(i2)') IMONTH
      READ(STRING(9:10),'(i2)') IDAY
      DATE_OF_OBSERVATION = IYEAR*10000 + IMONTH*100 + IDAY
      IF (SCAN_TYPE.LE.2.AND.(COORDINATE_SYSTEM.LE.4.OR.
     $(COORDINATE_SYSTEM.GE.6.AND.COORDINATE_SYSTEM.LE.9))) THEN
         CALL POSITION (SL1,SL2,SL3,SL4,SB1,SB2,SB3,SB4)
         WRITE (FIELD_CENTER,
     $   '(2(i2.2,a),i2.2,a,i3.3,x,I3.2,2(A,I2.2),A,I2.2)')
     $   SL1,C,SL2,C,SL3,D,SL4,SB1,C,SB2,C,SB3,D,SB4
      ELSE
         FIELD_CENTER = ' '
      ENDIF
c     name_of_reducer
      CALL SIC_GET_DBLE('F%NCHAN',DX,ERROR)
      N_RECEIVERS = DX
      CALL SIC_GET_DBLE('F%NRECAC',DX,ERROR)
      NCHAN       = DX
      NCALPOW = N_RECEIVERS - NCHAN
*
      IS_LST = .TRUE.
      SCAN_COORD_FOUND = .TRUE.
*
      RECORD_2PHASES = .FALSE.
      CALL SIC_DESCRIPTOR('F%DATA%PHASE2',DESC,FOUND)
      IF (FOUND) RECORD_2PHASES = .TRUE.
      VAX_PROCESSED = .TRUE.             ! to separate DRUMBEAT files
*
* Bolometer Section -----------------------
*
      DO I = 1,NCHAN
         RECEIVER_NUMBER(I) = I
      ENDDO
      REFERENCE_RECEIVER = -1
      CALL SIC_GET_DBLE('F%RECCNTRL',DX,ERROR)
      IF (.NOT.ERROR) REFERENCE_RECEIVER = DX
      IF (TELESCOPE(1:3).EQ.'CSO') THEN
         IF (REFERENCE_RECEIVER.EQ.-1) THEN
            CALL MESSAGE (2,1,'FILE','No reference channel defined')
            CALL MESSAGE (2,1,'FILE','Using channel 12 by default')
            REFERENCE_RECEIVER = 12
         ENDIF
      ENDIF
      CALL SIC_GET_DBLE('F%GAIN',DX,ERROR)
      REFERENCE_GAIN     = DX
      CALL SIC_GET_DBLE('F%FREQUENC',DX,ERROR)
      FREQUENCY          = DX
*
* Antenna Section -------------------------
*
      WOBBLER_DEVICE     = 0
      WOBBLER_PERIOD     = 0
      WOBBLER_THROW = -999
      CALL SIC_GET_DBLE('F%WOBTHROW',DX,ERROR)
      IF (.NOT.ERROR) WOBBLER_THROW     = DX
      IF (WOBBLER_THROW.NE.0.AND.WOBBLER_THROW.NE.-999) THEN
         CALL SIC_GET_DBLE('F%WOBDEVCE',DX,ERROR)
         WOBBLER_DEVICE     = DX
         CALL SIC_GET_DBLE('F%WOBTHROW',DX,ERROR)
         WOBBLER_THROW     = DX
         CALL SIC_GET_DBLE('F%WOBCYCLE',DX,ERROR)
         WOBBLER_PERIOD     = DX
      ENDIF
*
* Calibration Section ----------------------
*
*     Gain is initialized by routine INIT_OFFSET_GAIN()
*
* Flag Section ----------------------
*
*     GAIN_ELEV_FLAG      = .FALSE.
*     GAIN_RECEIVER_FLAG  = .FALSE.
*     OPACITY_FLAG        = .FALSE.
*     AMPLI_FLAG          = .FALSE.
*     JANSKY_FLAG         = .FALSE.
      NOISE_REMOVAL_FLAG  = .FALSE.
      DO I=1,NCHAN
         SPIKE_REMOVAL_FLAG(I) = -1
      ENDDO
*
* Total Power Section -----------------------
*
      IF (NCALPOW.NE.0) THEN
         DO I = 1,NCALPOW
            CALPOW_FLAG(I)   = .FALSE.
            CALPOW_GAIN(I)   = 1.
         ENDDO
         CALPOW_NUMBER(1) = 1
         CALPOW_NUMBER(2) = 17
         CALPOW_NUMBER(3) = 32
      ENDIF
*
*
* Parameters used by MAKPLAN
      NBI = 1
      CALL SIC_GET_DBLE('F%OLAM',DX,ERROR)
      OLAM  = DX
      CALL SIC_GET_DBLE('F%VLAM',DX,ERROR)
      VLAM  = DX
      CALL SIC_GET_DBLE('F%OBET',DX,ERROR)
      OBET  = DX
      CALL SIC_GET_DBLE('F%SAZM',DX,ERROR)
      SAZM  = DX
      CALL SIC_GET_DBLE('F%OAZM',DX,ERROR)
      OAZM  = DX/3600
      CALL SIC_GET_DBLE('F%VAZM',DX,ERROR)
      VAZM  = DX/60
      CALL SIC_GET_DBLE('F%OELV',DX,ERROR)
      OELV  = DX/3600
      CALL SIC_GET_DBLE('F%VELV',DX,ERROR)
      VELV  = DX/60
      CALL SIC_GET_DBLE('F%SINT',DX,ERROR)
      SINT  = DX
      CALL SIC_GET_DBLE('F%DAZM',DX,ERROR)
      DAZM  = DX/3600
      CALL SIC_GET_DBLE('F%TEMPERAT',DX,ERROR)
      AMBIENT_TEMPERATURE  = DX
      CALL SIC_GET_DBLE('F%PRESSURE',DX,ERROR)
      AMBIENT_PRESSURE    = DX
      CALL SIC_GET_DBLE('F%HUMIDITY',DX,ERROR)
      HUMIDITY            = DX
      CALL SIC_GET_DBLE('F%REFRACT',DX,ERROR)
      REFRACTION          = DX
      CALL SIC_GET_DBLE('F%NULA',DX,ERROR)
      PNULA = DX
      CALL SIC_GET_DBLE('F%NULE',DX,ERROR)
      PNULE = DX
      CALL SIC_GET_DBLE('F%COLSTAR',DX,ERROR)
      COLSTAR = DX
      IF (SCAN_TYPE.EQ.5.AND.DAZM.GT.1E-8) SCAN_TYPE = 1   ! ONOFF
*
C Read number of subscans
      CALL SIC_DESCRIPTOR('F%SUBS%TIME',DESC,FOUND)
      NSUBSCAN = 0
      IF (FOUND) THEN
         NSUBSCAN = DESC(4)
         WRITE(6,*) 'Number of subscans: ',NSUBSCAN
      ENDIF
      SUBSCAN_NUMBER = NSUBSCAN
      IF (NSUBSCAN.LE.0) THEN
         NRECORD = 0
         RETURN
      ENDIF
C Read number of records
      NRECORD = 0
      CALL SIC_DESCRIPTOR('F%DATA%LST',DESC,FOUND)
      IF (FOUND) THEN
         NRECORD = DESC(4)
      ELSE
         CALL SIC_DESCRIPTOR('F%DATA%UT',DESC,FOUND)
         IF (FOUND) NRECORD = DESC(4)
      ENDIF
      WRITE(6,*) 'Number of records: ',NRECORD
      IF (NRECORD.LE.0) RETURN
*
* Allocate shared memory for :
*    - 3 parameters to describe structure of shared memory :
*        n_receivers, nparam(=5), nparam_subscan(=4)
*    - NPARAM parameters :
*         receiver_az_offset(nchan)
*         receiver_el_offset(nchan)
*         receiver_flag(nchan)
*         receiver_rms(nchan)
*         receiver_gain(nchan)
*    - subscan info : subscan_number,subscan_epoch,subscan_time,nsubrec
*    - LST,AZ,EL,SIGNAL,SCAN_COORD,MDB and EKH support
*
      NRECORD_ALLOCATED = NRECORD
      NPARAM = 5
      NPARAM_SUBSCAN=4
      LEN_RECORD = 5+N_RECEIVERS         ! UT/LST,AZ,EL,2*SCAN_COORD and n chann
      LENGTH = 3+NCHAN*NPARAM + NSUBSCAN*NPARAM_SUBSCAN +
     $NRECORD_ALLOCATED*LEN_RECORD
      IF (SCAN_TYPE.EQ.2 .OR. SCAN_TYPE.EQ.4)
     $LENGTH = LENGTH+ 6*MAX_SEG*NCHAN+2 ! support
      CALL CREATE_SHARED_MEMORY(DATA_SHMID,LENGTH*4)
      CALL ATTACH_MEMORY(DATA_SHMID,IADDRRECORD)
      IP   = POINTER(IADDRRECORD,MEMORY)
      CALL R4TOR4(N_RECEIVERS,MEMORY(IP),1)
      CALL R4TOR4(NPARAM,MEMORY(IP+1),1)
      CALL R4TOR4(NPARAM_SUBSCAN,MEMORY(IP+2),1)
      IPPARAM   = IP +3
      IPSUBSCAN = IPPARAM   + NPARAM*NCHAN
      IPRECORD  = IPSUBSCAN + NPARAM_SUBSCAN*NSUBSCAN
! Allocate memory to store 2 phases
      IF (RECORD_2PHASES) THEN
         IER = SIC_GETVM (2*NRECORD*NCHAN,IADDR_PHASES)
         IF (IER.NE.1) THEN
            CALL MESSAGE (2,1,'FILE','Memory allocation failure')
            RETURN
         ENDIF
         IP_PHASE1 = POINTER(IADDR_PHASES,MEMORY)
         IP_PHASE2 = IP_PHASE1 + NRECORD*NCHAN
      ENDIF
*
*
*
* Update data :
*
*
*   Subscan information
*
*
      CALL SIC_DESCRIPTOR('F%SUBS%NUMBER',DESC,FOUND)
      IP   = POINTER(DESC(2),MEMORY)
      CALL R4TOR4(MEMORY(IP),MEMORY(IPSUBSCAN),NSUBSCAN)
      CALL SIC_DESCRIPTOR('F%SUBS%EPOCH',DESC,FOUND)
      IP   = POINTER(DESC(2),MEMORY)
      CALL R4TOR4(MEMORY(IP),MEMORY(IPSUBSCAN+NSUBSCAN),NSUBSCAN)
      CALL SIC_DESCRIPTOR('F%SUBS%TIME',DESC,FOUND)
      IP   = POINTER(DESC(2),MEMORY)
      CALL R4TOR4(MEMORY(IP),MEMORY(IPSUBSCAN+2*NSUBSCAN),NSUBSCAN)
*
*
*   LST,AZ,EL,SCAN_COORD,SUPPORT,UT
*
      IP_LST    = IPRECORD
      IP_AZ     = IP_LST + NRECORD
      IP_EL     = IP_AZ  + NRECORD
      IP_SIGNAL = IP_EL  + NRECORD
      IP_CALPOW = IP_SIGNAL + NCHAN*NRECORD
      IP_SCAN_COORD = IP_SIGNAL + NRECORD*N_RECEIVERS
      IF (SCAN_TYPE.EQ.2 .OR. SCAN_TYPE.EQ.4) THEN   ! ON-THE-FLY
         IP_MDB_SUPPORT = IP_SCAN_COORD  + NRECORD*2
         IP_EKH_SUPPORT = IP_MDB_SUPPORT + 4*MAX_SEG*NCHAN + 1
         CALL RAZ_MEMORY(MEMORY(IP_MDB_SUPPORT),6*MAX_SEG*NCHAN+2)
      ENDIF
*
      IER = SIC_GETVM (NRECORD,IADDR_UT)
      IP_UT = POINTER (IADDR_UT,MEMORY)
*
      CALL SIC_DESCRIPTOR('F%DATA%LST',DESC,FOUND)
      IF (FOUND) THEN
         IP   = POINTER(DESC(2),MEMORY)
         CALL R4TOR4(MEMORY(IP),MEMORY(IP_LST),NRECORD)
      ENDIF
      CALL SIC_DESCRIPTOR('F%DATA%UT',DESC,FOUND)
      IF (FOUND) THEN
         IP   = POINTER(DESC(2),MEMORY)
         CALL R4TOR4(MEMORY(IP),MEMORY(IP_LST),NRECORD)
         OBSLONG = SITE_LONGITUDE/DEG_TO_RAD   ! in degrees and double precision
         CALL UTST(MEMORY(IP_LST),NRECORD,IDAY,IMONTH,IYEAR,OBSLONG)
      ENDIF
      CALL SIC_DESCRIPTOR('F%DATA%AZ',DESC,FOUND)
      IP   = POINTER(DESC(2),MEMORY)
      CALL R4TOR4(MEMORY(IP),MEMORY(IP_AZ),NRECORD)
      CALL SIC_DESCRIPTOR('F%DATA%EL',DESC,FOUND)
      IP   = POINTER(DESC(2),MEMORY)
      CALL R4TOR4(MEMORY(IP),MEMORY(IP_EL),NRECORD)
      CALL SIC_DESCRIPTOR('F%DATA%SCAN_COO',DESC,FOUND)
      IP   = POINTER(DESC(2),MEMORY)
      CALL R4TOR4(MEMORY(IP),MEMORY(IP_SCAN_COORD),2*NRECORD)
      CALL SIC_DESCRIPTOR('F%DATA%PHASE1',DESC,FOUND)
      IP   = POINTER(DESC(2),MEMORY)
      IF (FOUND) CALL R4TOR4(MEMORY(IP),MEMORY(IP_PHASE1),NCHAN*NRECORD)
      CALL SIC_DESCRIPTOR('F%DATA%PHASE2',DESC,FOUND)
      IP   = POINTER(DESC(2),MEMORY)
      IF (FOUND) CALL R4TOR4(MEMORY(IP),MEMORY(IP_PHASE2),NCHAN*NRECORD)
      CALL SIC_DESCRIPTOR('F%DATA%TOT',DESC,FOUND)
      IP   = POINTER(DESC(2),MEMORY)
      IF (FOUND) CALL R4TOR4(MEMORY(IP),MEMORY(IP_CALPOW),
     $NCALPOW*NRECORD)
* Compute signal (phase1-phase2)
      IF (RECORD_2PHASES) THEN
         CALL DIFF_PHASE1_PHASE2(MEMORY(IP_SIGNAL),MEMORY(IP_PHASE1),
     $   MEMORY(IP_PHASE2),NCHAN*NRECORD)
      ELSE
         CALL SIC_DESCRIPTOR('F%DATA%SIGNAL',DESC,FOUND)
         IP   = POINTER(DESC(2),MEMORY)
         IF (FOUND) CALL R4TOR4(MEMORY(IP),MEMORY(IP_SIGNAL),
     $   NCHAN*NRECORD)
      ENDIF
*
*   Gain-elevation curve and Receivers offsets and gain
*
      IP_AZOFF = IPPARAM
      IP_ELOFF = IP_AZOFF + NCHAN
      IP_FLAG  = IP_ELOFF + NCHAN
      IP_RMS   = IP_FLAG  + NCHAN
      IP_GAIN  = IP_RMS   + NCHAN
      FIELD_ROTATION = 1                 ! default value for Pico Veleta
      IF (TELESCOPE(1:3).EQ.'CSO') THEN
         CALL INIT_OFFSET_GAIN_SHARC(RECEIVER_AZ_OFFSET(1),
     $   RECEIVER_EL_OFFSET(1),RECEIVER_GAIN(1),NCHAN,TEXT)
         GAIN_ELEV_FLAG = .FALSE.        ! no correction applied
         FIELD_ROTATION = 0              ! no rotation
      ELSEIF (TELESCOPE(1:3).EQ.'HHT') THEN
         CALL INIT_OFFSET_GAIN_HHT(RECEIVER_AZ_OFFSET(1),
     $   RECEIVER_EL_OFFSET(1),RECEIVER_GAIN(1),NCHAN,TEXT)
         GAIN_ELEV_FLAG = .FALSE.        ! no correction applied
         FIELD_ROTATION = -1             ! opposite rotation (compared to Pico Veleta)
      ELSE                               ! pico veleta
         CALL GET_SCAN_UT (UT1,MEMORY(IP_LST),1)
         CALL INIT_OFFSET_GAIN(RECEIVER_AZ_OFFSET(1),
     $   RECEIVER_EL_OFFSET(1),RECEIVER_GAIN(1),NCHAN,TEXT,
     $   UT1)
         CALL INIT_GAIN_ELEVATION(GAIN_ELEV)
      ENDIF
*
* Check for file nic.fit
*
      IER = SIC_GETLUN(LUN)
      IF (IER.NE.1) THEN
         WRITE(6,*)' Error getting unit for nic.fit'
         GOTO 99
      ENDIF
      IER = SIC_OPEN(LUN,'nic.fit','old',.FALSE.)
      IF (IER.EQ.0) THEN
         FOUND_SCAN = -1
300      READ (LUN,'(a)',ERR=350,END=400) STRING
         IF (INDEX(STRING,'! Scan:').NE.0) THEN
            SAVE_STRING = STRING
            READ (STRING(9:12),'(i4)',ERR=350,END=350) ISCAN
            I = INDEX(STRING,'Date: ')+6
            READ (STRING(I:I+1),'(i2)',ERR=350,END=350) IDATE
            I = I+3
            DO J=1,12
               IF (STRING(I:I+2).EQ.MONTH(J)) IMONTH1 = J
            ENDDO
            IDATE = IDATE+100*IMONTH1
            I = I+4
            READ (STRING(I:I+3),'(i4)',ERR=350,END=350) J
            IDATE = IDATE+J*10000
            IF (DATE_OF_OBSERVATION.LT.IDATE) GOTO 300
            CALL R4TOR4 (MEMORY(IP_LST),UT4,1)
            CALL STUT (UT4,1,IDAY,IMONTH,IYEAR,
     $      SITE_LONGITUDE/DEG_TO_RAD)
            I = INDEX(STRING,'Time: ')+6
            READ (STRING(I:I+6),'(f7.1)',ERR=350,END=350) UT3
            IF (DATE_OF_OBSERVATION.EQ.IDATE.AND.
     $      UT3-UT4.GT.1) GOTO 300       ! one second enough ?
*
* Get Az, El and Gain
*
310         READ (LUN,'(a)',ERR=350,END=350) STRING
            IF (INDEX(STRING,'! Az:').EQ.0) GOTO 310
            DO J=1,NCHAN
               READ (STRING((J-1)*8+7:J*8+6),'(f8.2)',ERR=350,END=350)
     $         RECEIVER_AZ_OFFSET(J)
            ENDDO
320         READ (LUN,'(a)',ERR=350,END=350) STRING
            IF (INDEX(STRING,'! El:').EQ.0) GOTO 320
            DO J=1,NCHAN
               READ (STRING((J-1)*8+7:J*8+6),'(f8.2)',ERR=350,END=350)
     $         RECEIVER_EL_OFFSET(J)
            ENDDO
330         READ (LUN,'(a)',ERR=350,END=350) STRING
            IF (INDEX(STRING,'! Ga:').EQ.0) GOTO 330
            DO J=1,NCHAN
               READ (STRING((J-1)*8+7:J*8+6),'(f8.2)',ERR=350,END=350)
     $         RECEIVER_GAIN(J)
            ENDDO
            FOUND_SCAN = ISCAN
            TEXT = SAVE_STRING
         ENDIF
         GOTO 300
350      CALL MESSAGE (2,1,'INIT',
     $   'File NIC.FIT is corrupted or incomplete')
400      CLOSE (LUN)
         IF (FOUND_SCAN.NE.-1) THEN
            WRITE (STRING,'(a,i4.4,a)')
     $      'Beam parameter retrieved from scan [',FOUND_SCAN,
     $      '] in NIC.FIT'
            CALL MESSAGE (2,1,'BEAM',STRING)
            TEXT = TEXT(15:)
         ENDIF
      ENDIF
      IER = SIC_FRELUN(LUN)
*
      CALL R4TOR4(RECEIVER_AZ_OFFSET(1),MEMORY(IP_AZOFF),NCHAN)
      CALL R4TOR4(RECEIVER_EL_OFFSET(1),MEMORY(IP_ELOFF),NCHAN)
      CALL R4TOR4(RECEIVER_GAIN(1),MEMORY(IP_GAIN),NCHAN)
      IF (OLD_FLAG.EQ.0 .OR. OLD_FLAG.NE.NCHAN) THEN
         DO I=1,NCHAN
            RECEIVER_FLAG(I) = .FALSE.
         ENDDO
      ENDIF
      CALL R4TOR4(RECEIVER_FLAG(1),MEMORY(IP_FLAG),
     $NCHAN)
      CALL RAZ_MEMORY(MEMORY(IP_RMS),NCHAN)    !receiver_rms
*
* Allocate memory for spike data
*
      LENGTH_SPIKE = 1000
      LENGTH = NCHAN + 2*NCHAN*LENGTH_SPIKE
      IER = SIC_GETVM (LENGTH,IADDRSPIKE)
      IPSPIKE = POINTER(IADDRSPIKE,MEMORY)
      CALL RAZ_MEMORY(MEMORY(IPSPIKE),LENGTH)
*
* Allocate memory for flagged subscans
*
      LENGTH = NCHAN*NSUBSCAN*2
      IER = SIC_GETVM (LENGTH,IADDR_SUBFLAG)
      IF (IER.NE.1) THEN
         CALL MESSAGE (2,1,'FILE','Memory allocation failure')
         RETURN
      ENDIF
      IP = POINTER(IADDR_SUBFLAG,MEMORY)
      CALL RAZ_MEMORY(MEMORY(IP),LENGTH)
! Unit (jansky or counts)
      IF (SCAN_TYPE.NE.3.AND.JY_TO_COUNTS.NE.0.AND.JANSKY_FLAG) THEN
         JY_STATE = .FALSE.
         CALL DO_JY (MEMORY(IP_SIGNAL),JY_STATE,JY_TO_COUNTS)
      ENDIF
*
* Compute number of records per subscan
*
      IPSUBREC = IPSUBSCAN + 3*NSUBSCAN
      CALL RECORDS_IN_SUBSCAN (MEMORY(IPSUBREC),MEMORY(IPSUBSCAN),
     $MEMORY(IPSUBSCAN+2*NSUBSCAN),NSUBSCAN,MEMORY(IP_LST),
     $MEMORY(IP_SCAN_COORD))
      CALL R4TOR4(MEMORY(IP_LST),MEMORY(IP_UT),NRECORD)
      CALL STUT (MEMORY(IP_UT),NRECORD,IDAY,IMONTH,IYEAR,
     $SITE_LONGITUDE/DEG_TO_RAD)
*
* skew
      SKEW = 0
      IF (SCAN_TYPE.EQ.2.AND.TELESCOPE(1:3).NE.'CSO') THEN
         CALL ANY_SKEW (MEMORY(IPSUBREC),MEMORY(IP_SCAN_COORD),
     $   NSUBSCAN)
         IF (OELV.GT.0) SKEW = -SKEW
      ENDIF
      OFF = .FALSE.
      CALL STANDARD_VARIABLES (OFF,.FALSE.)
      IF (SECTION(1).EQ.1) CALL HEADER_VARIABLES(OFF,.FALSE.)
      IF (SECTION(1).EQ.2) CALL HEADER_VARIABLES(OFF,.TRUE.)
      IF (SECTION(2).EQ.1) CALL BOLOMETER_VARIABLES(OFF,.FALSE.)
      IF (SECTION(2).EQ.2) CALL BOLOMETER_VARIABLES(OFF,.TRUE.)
      IF (SECTION(3).EQ.1) CALL ANTENNA_VARIABLES(OFF,.FALSE.)
      IF (SECTION(3).EQ.2) CALL ANTENNA_VARIABLES(OFF,.TRUE.)
      IF (SECTION(4).EQ.1) CALL CALIBRATION_VARIABLES(OFF,.FALSE.)
      IF (SECTION(4).EQ.2) CALL CALIBRATION_VARIABLES(OFF,.TRUE.)
      IF (SECTION(5).EQ.1) CALL FLAG_VARIABLES(OFF,.FALSE.)
      IF (SECTION(5).EQ.2) CALL FLAG_VARIABLES(OFF,.TRUE.)
      IF (SECTION(6).EQ.1) CALL SKYDIP_VARIABLES(OFF,.FALSE.)
      IF (SECTION(6).EQ.2) CALL SKYDIP_VARIABLES(OFF,.TRUE.)
      IF (SECTION(7).EQ.1) CALL FIT_VARIABLES(OFF,.FALSE.)
      IF (SECTION(7).EQ.2) CALL FIT_VARIABLES(OFF,.TRUE.)
      IF (SECTION(8).EQ.1) CALL DATA_VARIABLES(OFF,.FALSE.)
      IF (SECTION(8).EQ.2) CALL DATA_VARIABLES(OFF,.TRUE.)
      IF (SECTION(9).EQ.1) CALL CALPOW_VARIABLES(OFF,.FALSE.)
      IF (SECTION(9).EQ.2) CALL CALPOW_VARIABLES(OFF,.TRUE.)
*
      ANY_BASE = .FALSE.
      STATE = .FALSE.
      IF (AMPLI_FLAG) CALL DO_AM (MEMORY(IP_SIGNAL),STATE,
     $NINT(REFERENCE_GAIN))
      GAIN_REC = .FALSE.
      IF (GAIN_RECEIVER_FLAG) CALL DO_GR (MEMORY(IP_SIGNAL),
     $MEMORY(IP_GAIN))
      IF (COMPARE (MEMORY(IP_GAIN),BEAM_GAIN(1)).AND.
     $GAIN_RECEIVER_FLAG)
     $CALL MESSAGE (2,1,'BEAM',TEXT)
      CALL R4TOR4 (MEMORY(IP_GAIN),BEAM_GAIN(1),NCHAN)
      STATE = .FALSE.
      IF (GAIN_ELEV_FLAG.AND.SCAN_TYPE.NE.3.AND.SCAN_TYPE.NE.8)
     $CALL DO_GE (MEMORY(IP_EL),MEMORY(IP_SIGNAL),GAIN_ELEV,STATE)
      SCAN_OPACITY = 0
      IF (OPACITY_FLAG.AND.SCAN_TYPE.NE.3.AND.SCAN_TYPE.NE.8) THEN
         CALL GET_SCAN_UT (UT1,MEMORY(IP_LST),1)
         CALL GET_SCAN_UT (UT2,MEMORY(IP_LST),NRECORD)
         ATMOSPHERIC_OPACITY = MANUAL_OPACITY
         OPACITY_GRADIENT = 0
         IF (NUMBER_OF_SKYDIPS.GT.0.AND.MANUAL_OPACITY.EQ.0) THEN
            IP_SKYUT = POINTER(SKYDIP_MEM,MEMORY)
            IP_SKYTAU = IP_SKYUT + 2*NUMBER_OF_SKYDIPS     ! double precision
            CALL GET_OPACITY (UT1,MEMORY(IP_SKYUT),
     $      MEMORY(IP_SKYTAU),NUMBER_OF_SKYDIPS,
     $      ATMOSPHERIC_OPACITY,OPACITY_GRADIENT,DATA_TYPE)
         ENDIF
         IF (ATMOSPHERIC_OPACITY.GT.1.AND.DATA_TYPE)
     $   CALL MESSAGE (2,1,'OPACITY','Crazy atmosphere')
         IF (ATMOSPHERIC_OPACITY.LT.0) THEN
            IF (DATA_TYPE) CALL MESSAGE (2,1,'OPACITY',
     $      'no correction applied (Undefined ?)')
         ELSE
            STATE = .FALSE.
            CALL DO_OP (MEMORY(IP_SIGNAL),MEMORY(IP_EL),STATE,
     $      ATMOSPHERIC_OPACITY,OPACITY_GRADIENT,UT1,UT2,DATA_TYPE)
            OLD_OPACITY = ATMOSPHERIC_OPACITY
            OLD_OPACITY_GRADIENT = OPACITY_GRADIENT
            OLD_UT1 = UT1
            OLD_UT2 = UT2
         ENDIF
      ENDIF
*
* Update MDB and EKH support
      IF (SCAN_TYPE.EQ.2 .AND. IADDR_SUPPORT.NE.0) THEN
         CALL R4TOR4(MEMORY(IP_CNV_SUPPORT),PNT,1)
         CALL UPDATE_MDB_EKH_SUPPORT(PNT,
     $   MEMORY(IP_CNV_SUPPORT+1),MEMORY(IP_CNV_SUPPORT+1+PNT),
     $   ERROR)
      ENDIF
*
* Delete FITS structure
      CALL SIC_DELVARIABLE('F',.FALSE.,ERROR)
      CALL MESSAGE(2,1,'INPUT_FITS','Deleting fits structure F%')
      RETURN
*
99    CALL MESSAGE (8,4,'INPUT_FITS','Not enough memory ')
      RETURN
      END
*
***********************************************************************
      SUBROUTINE  DIFF_PHASE1_PHASE2(OUT,PHASE1,PHASE2,N)
*
      INTEGER N
      REAL OUT(1),PHASE1(1),PHASE2(1)
*
      INTEGER I
      DO I=1,N
         OUT(I)=PHASE1(I)-PHASE2(I)
      ENDDO
      END
