* Last processed by NICE on 12-Jun-2000 15:53:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      SUBROUTINE FIELD (LINE,ELEVATION,DEVIATION,ERROR)
C ----------------------------------------------------------------------
C       Support routine for command
C
C       FIELD   /SMEAR [time]
C       FIELD   /DECL delta /HOUR hour_angle /SIZE field_size /SKEW azimuth
C               /NOISE rms rmin /NMB scan wobbler_throw channels
C ----------------------------------------------------------------------
      CHARACTER*(*) LINE
      LOGICAL ERROR
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
      INCLUDE 'inc:format.inc'
*
      REAL*8 LATI, DECL, INDEX1, INDEX2, HOUR, COSW, HIGHDEC
      REAL*8 COS_EL_SIN_AZ, COS_EL_COS_AZ, WHEN, SINDEC, COSDEC, LOWDEC
      REAL*8 SIN_EL, AZ, EL, EL1, EL2, SINL, COSL, COSH, SINH, AZ1, AZ2
      REAL*8 EL3, EL4, AZ3, AZ4, A1, A2, A3, DURA, SINW, OBSER
      REAL*8 HOURM, LST, AZN, ELN, SKEW_OPT, USEROFFRA, USEROFFDE
      REAL*8 DEL
*
      REAL*4 ELEVATION(0:180,0:180)      ! elevation
      REAL*4 DEVIATION(0:180,0:180)      ! positional error
      REAL*4 LEV1(7), LEV2(2), DE3, SCANTIM
      REAL*4 OFFDE(251*2), OFFRA(251*2), XPOS, YPOS
      REAL*4 RMS, RMIN, WOB
*
      CHARACTER LIMITS*40, OBSERVATORY*80, INTTIME*80, DECLINATION*80
      CHARACTER DELTA*10, DRAWMARK*40, VELOCITY*80, ELLIPSE*80
      CHARACTER WOBB_MESS*80
      CHARACTER CHAIN*40, NMB_NAME*40
      INTEGER*4 IP_OUTPUT, IP, IADDR_OUTPUT, POINTER
      INTEGER IER, SIC_GETVM, ISIZE
      INTEGER INDEX3, DIM(4), OBSE1, OBSE2, OBSE3, OBSE4, I
      INTEGER N1, DE1, DE2, NSUB, ISCAN, NCHAN_NMB
      LOGICAL SIC_PRESENT, ANYDECLOFF, NMB_OPT
      DATA LEV1 /.0,.1,.2,.3,.4,.5,1./ LEV2 /20,80/
*
      LATI = SITE_LATITUDE
      IF (LATI.EQ.0) THEN
         LATI = (37D0+4/6D1+5.6/36D2)*DEG_TO_RAD
*        CALL MESSAGE (2,1,'FIELD','Assuming Pico Veleta Observatory')
      ENDIF
      SINL = SIN(LATI)
      COSL = COS(LATI)
      OBSER = LATI/DEG_TO_RAD
      OBSE1 = INT(OBSER)
      OBSE2 = INT((OBSER-OBSE1)*60)
      OBSE3 = INT((OBSER-OBSE1-OBSE2/60.)*36D2)
      OBSE4 = NINT((OBSER-OBSE1-OBSE2/60.-OBSE3/36D2)*36D5)
      WRITE (OBSERVATORY,'(A,I3.2,1H:,I2.2,1H:,I2.2,1H.,I3.3,A)')
     $'DRA TEX 1 1.02 "OBS LATITUDE: ',OBSE1,OBSE2,
     $OBSE3,OBSE4,'" 7 /user'
      ERROR = .FALSE.
*
* Option /SMEAR
*
      DURA = 1/4.                        ! default is 20 minutes (left and right)
      IF (SIC_PRESENT(1,1)) THEN
         CALL SIC_R8(LINE,1,1,DURA,.TRUE.,ERROR)
         IF (ERROR) GOTO 99
         DURA = DURA/120.                ! left and right
         IF (DURA.LE.1E-5.OR.DURA.GT.5) THEN
            CALL MESSAGE (2,1,'FIELD','Invalid time')
            ERROR = .TRUE.
            RETURN
         ENDIF
      ENDIF
      IF (SIC_PRESENT(1,0).AND.(SIC_PRESENT(2,0).OR.SIC_PRESENT(3,0)
     $.OR.SIC_PRESENT(4,0).OR.SIC_PRESENT(5,0).OR.SIC_PRESENT(6,0)
     $.OR.SIC_PRESENT(7,0).OR.SIC_PRESENT(8,0))) THEN
         CALL MESSAGE (2,1,'FIELD','Incompatible options')
         ERROR = .TRUE.
         RETURN
      ENDIF
      WRITE (INTTIME,'(A,F5.1,A)')
     $'DRA TEX 1 1.1 "INTEGRATION TIME: ',DURA*120,' MIN" 7 /user'
*
* Option /NMB
*
      NMB_OPT = .FALSE.
      NCHAN_NMB = 37
      IF (SIC_PRESENT(7,0)) THEN
         NMB_NAME = '9999.nmb'
         IF (SIC_PRESENT(7,1)) THEN
            CALL SIC_I4(LINE,7,1,ISCAN,.TRUE.,ERROR)
            CHAIN = ' '
            WRITE(CHAIN(1:4),'(i4.4)') ISCAN
            CALL SIC_PARSEF(CHAIN,NMB_NAME,' ','.nmb')
         ENDIF
         WOB = -1
         IF (SIC_PRESENT(7,2)) THEN
            CALL SIC_R4(LINE,7,2,WOB,.TRUE.,ERROR)
         ELSE
            ERROR = .FALSE.
            CALL SIC_GET_REAL ('WOBBLER_THROW',WOB,ERROR)
            IF (ERROR) THEN
               CALL MESSAGE (2,1,'FIELD','Wobbler throw is not defined')
               ERROR = .TRUE.
               RETURN
            ELSE
               WOBB_MESS = 'Using xx.x arcs for wobbler throw'
               WRITE (WOBB_MESS(7:10),'(F4.1)') WOB
               CALL MESSAGE (2,1,'FIELD',WOBB_MESS)
            ENDIF
         ENDIF
         NMB_OPT = .TRUE.
         IF (SIC_PRESENT(7,3)) THEN
            CALL SIC_I4(LINE,7,3,NCHAN_NMB,.TRUE.,ERROR)
         ENDIF
      ENDIF
*
* Option /SKEW
*
      SKEW_OPT = 0                       ! default is no SKEW_OPT
      IF (SIC_PRESENT(6,1)) THEN
         CALL SIC_R8(LINE,6,1,SKEW_OPT,.TRUE.,ERROR)
         IF (ERROR) GOTO 99
         IF (ABS(SKEW_OPT).GT.50) THEN
            CALL MESSAGE (2,1,'FIELD','Crazy SKEW_OPT')
            ERROR = .TRUE.
            RETURN
         ENDIF
      ENDIF
*
      WRITE (INTTIME,'(A,F5.1,A)')
     $'DRA TEX 1 1.1 "INTEGRATION TIME: ',DURA*120,' MIN" 7 /user'
*
* Case: ON-THE-FLY
*
      IF (.NOT.SIC_PRESENT(1,0)) THEN
         IF (SINT.LT.1E-5) SINT = 4.
         CALL SIC_GET_REAL ('VAZM',VAZM,ERROR)
         IF (ERROR) THEN
            VAZM = 4/60.
            VELOCITY = 'Assuming x.x arcs/sec in azimuth (VAZM)'
            WRITE (VELOCITY(10:12),'(F3.1)') VAZM*60
         ELSE
            IF (VAZM.LT.1E-5) VAZM = 4/60.
            VELOCITY = 'Using x.x arcs/sec setting (VAZM)'
            WRITE (VELOCITY(7:9),'(F3.1)') VAZM*60
         ENDIF
*
         DE1 = 0
         DE2 = 0
         DE3 = 0
         IF (SIC_PRESENT(2,1)) THEN
            CALL SIC_CH(LINE,2,1,DELTA,N1,.TRUE.,ERROR)
            IF (ERROR) GOTO 99
            N1 = INDEX(DELTA,':')
            IF (N1.NE.0) THEN
               DELTA(N1:N1) = ' '
               N1 = INDEX(DELTA,':')
               IF (N1.NE.0) THEN
                  DELTA(N1:N1) = ' '
                  IF (INDEX(DELTA,':').NE.0) GOTO 99
                  READ (DELTA,*,ERR=99,END=99) DE1,DE2,DE3
                  GOTO 98
               ENDIF
               READ (DELTA,*,ERR=99,END=99) DE1,DE2
               GOTO 98
            ENDIF
            READ (DELTA,*,ERR=99,END=99) DE1
         ENDIF
98       N1 = 1
         IF (DE1.LT.0) N1 = -1
         DECL = N1*(ABS(DE1)+DE2/60.+DE3/36D2)
         USEROFFRA = 0
         USEROFFDE = 0
         ANYDECLOFF = SIC_PRESENT(5,2)
         IF (ANYDECLOFF) THEN
            CALL SIC_R8 (LINE,5,2,USEROFFDE,.TRUE.,ERROR)
            IF (ERROR) GOTO 99
         ENDIF
         DECL = DECL+USEROFFDE/36D2
         IF (ABS(DECL).GT.90) GOTO 99
         ERROR = .FALSE.
*
         AZ = 200
         IF (SIC_PRESENT(3,1)) THEN
            CALL SIC_R8(LINE,3,1,AZ,.TRUE.,ERROR)
            IF (ERROR) GOTO 99
            IF (AZ.LE.1E-5.OR.AZ.GT.1000) THEN
               CALL MESSAGE (2,1,'FIELD','Invalid subscan width')
               ERROR = .TRUE.
               RETURN
            ENDIF
         ENDIF
         EL = 100
         IF (SIC_PRESENT(3,2)) THEN
            CALL SIC_R8(LINE,3,2,EL,.TRUE.,ERROR)
            IF (ERROR) GOTO 99
            IF (EL.LE.1E-5.OR.EL.GT.1000) THEN
               CALL MESSAGE (2,1,'FIELD','Invalid scan height')
               ERROR = .TRUE.
               RETURN
            ENDIF
         ENDIF
         NSUB = INT(EL/SINT+1)/2*2+1     ! to avoid surprises
         SCANTIM = (64./3600*AZ*EL/SINT/(VAZM*60))/60.
         AZ = AZ/2
         EL = EL/2
         DEL = EL
*
         HOURM = 0
         IF (SIC_PRESENT(4,1)) THEN
            CALL SIC_R8(LINE,4,1,HOURM,.TRUE.,ERROR)
            IF (ERROR) GOTO 99
            IF (HOURM.LT.-12.OR.HOURM.GT.12) THEN
               CALL MESSAGE (2,1,'FIELD','Invalid Hour Angle')
               ERROR = .TRUE.
               RETURN
            ENDIF
         ENDIF
         COSH = COS(HOURM/12*PI)
*
* option /NOISE
*
         NRECORD = (NINT(2*AZ/(VAZM*60.))*2+1)*NSUB
         ISIZE = 2*NRECORD*NCHAN_NMB+2*NCHAN_NMB**2+5*NCHAN_NMB
         IER = SIC_GETVM (ISIZE,IADDR_OUTPUT)
         IP_OUTPUT = POINTER(IADDR_OUTPUT,MEMORY)
         IP = IP_OUTPUT+NRECORD*NCHAN_NMB
         IF (IER.NE.1) THEN
            CALL MESSAGE (2,1,'FIELD',
     $      ' *** Unable to allocate memory *** ')
            GOTO 99
         ENDIF
         IF (SIC_PRESENT(8,0)) THEN
            CALL SIC_R4(LINE,8,1,RMS,.TRUE.,ERROR)
            IF (ERROR) GOTO 99
            CALL SIC_R4(LINE,8,2,RMIN,.TRUE.,ERROR)
            IF (ERROR) GOTO 99
            IF (RMIN.LT.0.OR.RMIN.GT.1) THEN
               CALL MESSAGE (2,1,'FIELD',
     $         'Allowed range for Mean_Corr: [0...1] ')
               ERROR = .TRUE.
               RETURN
            ENDIF
            CALL CORRELATED_NOISE(RMS,RMIN,NCHAN_NMB,
     $      MEMORY(IP),MEMORY(IP+NCHAN_NMB),   ! AZ+EL Offset
     $      MEMORY(IP+2*NCHAN_NMB),      ! CHAN_GAIN
     $      MEMORY(IP+3*NCHAN_NMB),      ! CHAN_FLAG
     $      MEMORY(IP+4*NCHAN_NMB),      ! NEW_CHAN_NUM
     $      MEMORY(IP+5*NCHAN_NMB),      ! A
     $      MEMORY(IP+5*NCHAN_NMB+NCHAN_NMB**2),     ! R
     $      MEMORY(IP+5*NCHAN_NMB+2*NCHAN_NMB**2),   ! E
     $      MEMORY(IP_OUTPUT),           ! Output
     $      ERROR)
            IF (ERROR) GOTO 99
         ELSEIF  (.NOT.SIC_PRESENT(8,0)) THEN
            CALL CORRELATED_NOISE(RMS,RMIN,NCHAN_NMB,
     $      MEMORY(IP),MEMORY(IP+NCHAN_NMB),   ! AZ+EL Offset
     $      MEMORY(IP+2*NCHAN_NMB),      ! CHAN_GAIN
     $      MEMORY(IP+3*NCHAN_NMB),      ! CHAN_FLAG
     $      MEMORY(IP+4*NCHAN_NMB),      ! NEW_CHAN_NUM
     $      MEMORY(IP+5*NCHAN_NMB),      ! A
     $      MEMORY(IP+5*NCHAN_NMB+NCHAN_NMB**2),     ! R
     $      MEMORY(IP+5*NCHAN_NMB+2*NCHAN_NMB**2),   ! E
     $      MEMORY(IP_OUTPUT),           ! Output
     $      ERROR)
            IF (ERROR) GOTO 99
         ENDIF
*
         DECL = DECL*DEG_TO_RAD
         DO I=1,NSUB
            LST = -(HOURM-SCANTIM*(NSUB/2.-I)/NSUB)*15
*           AZN = AZ+SKEW_OPT*(I-1)          ! starting unskewed...
            AZN = AZ+SKEW_OPT*(I-(NSUB/2+1))   ! starting skewed...
            ELN = -EL+SINT*(I-1)
            CALL AZTORA (AZN,ELN,LST,LATI,DECL,OFFRA(I),OFFDE(I))
         ENDDO
         DO I=1,NSUB
            LST = -(HOURM+SCANTIM*(NSUB/2.-I)/NSUB)*15
*           AZN = -AZ+SKEW_OPT*(NSUB-I)      ! starting unskewed...
            AZN = -AZ+SKEW_OPT*(NSUB/2+1-I)    ! starting skewed...
            ELN = EL-SINT*(I-1)
            CALL AZTORA (AZN,ELN,LST,LATI,DECL,OFFRA(I+NSUB),
     $      OFFDE(I+NSUB))
         ENDDO
         OFFRA(2*NSUB+1) = OFFRA(1)
         OFFDE(2*NSUB+1) = OFFDE(1)
         IF (SIC_PRESENT(5,1)) THEN
            CALL SIC_R8 (LINE,5,1,USEROFFRA,.TRUE.,ERROR)
            IF (ERROR) GOTO 99
         ENDIF
         DO I=1,2*NSUB+1
            OFFRA(I) = OFFRA(I) + USEROFFRA
            OFFDE(I) = OFFDE(I) + USEROFFDE
         ENDDO
         SINDEC = SIN(DECL)
         COSDEC = COS(DECL)
         SIN_EL        =  SINDEC*SINL+COSDEC*COSH*COSL
         EL = ASIN(SIN_EL)/DEG_TO_RAD
         CALL GR4_GIVE ('X',2*NSUB+1,OFFRA)
         CALL GR4_GIVE ('Y',2*NSUB+1,OFFDE)
         IF (.NOT.SIC_PRESENT(5,0)) THEN
            CALL GR_EXEC ('cl pl')
            CALL GR_EXEC ('set box 4 20 2 18')
            CALL GR_EXEC ('lim 0 1 -1 1')
            CALL GR_EXEC ('label "RA offset (arcsec)" /x')
            CALL GR_EXEC (OBSERVATORY)
            WRITE (INTTIME,'(A,i4,A)')
     $      'DRA TEX 1 1.1 "INTEGRATION TIME: ',
     $      NINT(SCANTIM*12+0.5)*5,' MIN" 7 /user'
            CALL GR_EXEC (INTTIME)
            WRITE (INTTIME,'(A,i4,A)')
     $      'DRA TEX 0 1.1 "MEAN ELEVATION: ',
     $      NINT(EL),' DEG" 9 /user'
            CALL GR_EXEC (INTTIME)
            CALL GR_EXEC ('set orie 90')
            CALL GR_EXEC ('dra tex -.1 0 "DEC offset (arcsec)" 5 /user')
            CALL GR_EXEC ('set orie 0')
            CALL GR_EXEC ('lim /var y x')
            CALL GR_EXEC ('lim < > < > /rev x')
            CALL GR_EXEC ('box')
         ENDIF
         CALL GR_EXEC ('set ang rad')    ! to avoid any problem with GREG
         ELLIPSE = 'elli 5.5 5.5 /user  xxxx.xxxx xxxx.xxxx'
         WRITE (ELLIPSE(20:39),'(2F10.4)') USEROFFRA, USEROFFDE
         CALL GR_EXEC (ELLIPSE)
         CALL GR_EXEC ('conn')
         XPOS = OFFRA(NSUB*2)-OFFRA(1)
         YPOS = OFFDE(1)-OFFDE(NSUB*2)
         DRAWMARK = 'set orie '
         WRITE (DRAWMARK(10:15),'(f6.1)') ATAN2(YPOS,XPOS)/DEG_TO_RAD
         CALL GR_EXEC (DRAWMARK)
         XPOS = (OFFRA(1)+OFFRA(NSUB*2))/2
         YPOS = (OFFDE(1)+OFFDE(NSUB*2))/2
         DRAWMARK = 'dra tex -1000.000 -1000.000 "'
     $   //CHAR(92)//'S|" 5 /user'
         WRITE (DRAWMARK(09:17),'(f9.3)') XPOS
         WRITE (DRAWMARK(19:27),'(f9.3)') YPOS
         CALL GR_EXEC ('set exp 3.')
         CALL GR_EXEC ('pen 0 /wei 5')
         CALL GR_EXEC (DRAWMARK)
         CALL GR_EXEC ('pen 0 /wei 0')
         CALL GR_EXEC ('set exp 1.')
         CALL GR_EXEC ('set orie 0')
*
         IF (ANYDECLOFF) THEN
            DECLINATION = 'Declination set to  xx:xx:xx.x'
            N1 = 1
            IF (DECL.LT.0) N1 = -1
            DECL = ABS(DECL/DEG_TO_RAD)
            DE1 = INT(DECL)
            DE2 = INT((DECL-DE1)*60.)
            DE3 = (DECL-DE1-DE2/60.)*3600
            DE1 = DE1*N1
            WRITE (DECLINATION(20:22),'(I3.2)') DE1
            WRITE (DECLINATION(24:25),'(I2.2)') DE2
            WRITE (DECLINATION(27:28),'(I2.2)') INT(DE3)
            WRITE (DECLINATION(30:30),'(I1)') (DE3-INT(DE3))*10
            CALL MESSAGE (2,1,'FIELD',DECLINATION)
         ENDIF
         CALL MESSAGE (2,1,'FIELD',VELOCITY)
         IF (NINT(VAZM*60).LT.4) THEN
            CALL MESSAGE (2,1,'FIELD','VAZM is not standard')
         ENDIF
         WRITE (INTTIME,'(A,I4,A)') 'Estimated integration time',
     $   NINT(SCANTIM*12+0.5)*5,' min'
         CALL MESSAGE (2,1,'FIELD',INTTIME)
         IF (NINT(EL).GT.10) THEN
            WRITE (INTTIME,'(A,I2,A)') 'Mean elevation of ',
     $      NINT(EL),' deg'
         ELSE
            WRITE (INTTIME,'(A,I2,A)') 'Too low in elevation'
         ENDIF
         CALL MESSAGE (2,1,'FIELD',INTTIME)
*
         IF (NMB_OPT)
     $   CALL WRITE_NMB(NMB_NAME,DECL,AZ,DEL,NSUB,WOB,HOURM,
     $   LATI,NCHAN_NMB,MEMORY(IP_OUTPUT))
         CALL FREE_VM(ISIZE,IADDR_OUTPUT)
         RETURN
      ENDIF
*
* Case: SMEAR
*
      LOWDEC  = -88+90+LATI/DEG_TO_RAD
      HIGHDEC = +88+90
*
      DO INDEX1 = 0,180
         HOUR = (INDEX1-90)/90*PI
         COSH = COS(HOUR)
         SINH = SIN(HOUR)
         INDEX3 = -1
         DO INDEX2 = LOWDEC,HIGHDEC,(HIGHDEC-LOWDEC)/180.
            INDEX3 = INDEX3 + 1
            DECL = (INDEX2-90)*DEG_TO_RAD
            SINDEC = SIN(DECL)
            COSDEC = COS(DECL)
            COS_EL_SIN_AZ = -COSDEC*SINH
            COS_EL_COS_AZ =  SINDEC*COSL-COSDEC*COSH*SINL
            IF (COS_EL_SIN_AZ.NE.0.OR.COS_EL_COS_AZ.NE.0) THEN
               AZ = ATAN2(COS_EL_SIN_AZ,COS_EL_COS_AZ)/DEG_TO_RAD
            ENDIF
            SIN_EL        =  SINDEC*SINL+COSDEC*COSH*COSL
            EL = ASIN(SIN_EL)/DEG_TO_RAD
            ELEVATION (INT(INDEX1),INT(INDEX3)) = EL
*
* Observing an hour: Object turns by 15 degrees
*
            WHEN = HOUR-(DURA+0.01/60.)/12.*PI
            SINW = SIN(WHEN)
            COSW = COS(WHEN)
            COS_EL_SIN_AZ = -COSDEC*SINW
            COS_EL_COS_AZ =  SINDEC*COSL-COSDEC*COSW*SINL
            IF (COS_EL_SIN_AZ.NE.0.OR.COS_EL_COS_AZ.NE.0) THEN
               AZ1 = ATAN2(COS_EL_SIN_AZ,COS_EL_COS_AZ)/DEG_TO_RAD
            ENDIF
            SIN_EL        =  SINDEC*SINL+COSDEC*COSW*COSL
            EL1 = ASIN(SIN_EL)/DEG_TO_RAD
*
            WHEN = HOUR-(DURA-0.01/60.)/12.*PI
            SINW = SIN(WHEN)
            COSW = COS(WHEN)
            COS_EL_SIN_AZ = -COSDEC*SINW
            COS_EL_COS_AZ =  SINDEC*COSL-COSDEC*COSW*SINL
            IF (COS_EL_SIN_AZ.NE.0.OR.COS_EL_COS_AZ.NE.0) THEN
               AZ2 = ATAN2(COS_EL_SIN_AZ,COS_EL_COS_AZ)/DEG_TO_RAD
            ENDIF
            SIN_EL        =  SINDEC*SINL+COSDEC*COSW*COSL
            EL2 = ASIN(SIN_EL)/DEG_TO_RAD
*
            WHEN = HOUR+(DURA-0.01/60.)/12.*PI
            SINW = SIN(WHEN)
            COSW = COS(WHEN)
            COS_EL_SIN_AZ = -COSDEC*SINW
            COS_EL_COS_AZ =  SINDEC*COSL-COSDEC*COSW*SINL
            IF (COS_EL_SIN_AZ.NE.0.OR.COS_EL_COS_AZ.NE.0) THEN
               AZ3 = ATAN2(COS_EL_SIN_AZ,COS_EL_COS_AZ)/DEG_TO_RAD
            ENDIF
            SIN_EL        =  SINDEC*SINL+COSDEC*COSW*COSL
            EL3 = ASIN(SIN_EL)/DEG_TO_RAD
*
            WHEN = HOUR+(DURA+0.01/60.)/12.*PI
            SINW = SIN(WHEN)
            COSW = COS(WHEN)
            COS_EL_SIN_AZ = -COSDEC*SINW
            COS_EL_COS_AZ =  SINDEC*COSL-COSDEC*COSW*SINL
            IF (COS_EL_SIN_AZ.NE.0.OR.COS_EL_COS_AZ.NE.0) THEN
               AZ4 = ATAN2(COS_EL_SIN_AZ,COS_EL_COS_AZ)/DEG_TO_RAD
            ENDIF
            SIN_EL        =  SINDEC*SINL+COSDEC*COSW*COSL
            EL4 = ASIN(SIN_EL)/DEG_TO_RAD
*
            IF ((EL1-EL2).NE.0.OR.(AZ1-AZ2).NE.0)
     $      A1 = ATAN2((EL1-EL2),(AZ1-AZ2))
            IF ((EL3-EL4).NE.0.OR.(AZ3-AZ4).NE.0)
     $      A2 = ATAN2((EL3-EL4),(AZ3-AZ4))
            A3 = ((EL3+EL4)-(EL1+EL2))/2*DEG_TO_RAD
            DEVIATION (INT(INDEX1),INT(INDEX3)) = ABS(TAN((A1-A2)+A3))
         ENDDO
      ENDDO
      DIM(1) = 181
      DIM(2) = 181
      DIM(3) = 1
      DIM(4) = 1
      CALL SIC_DEF_REAL ('field2',DEVIATION,2,DIM,.TRUE.,ERROR)
      CALL GR_EXEC ('clea pl')
      CALL GR_EXEC ('rgdata field2 /var')
      CALL GR_EXEC ('lim /rg')
      CALL GR_EXEC ('set box 4 20 2 18')
      CALL GR4_LEVELS (6,LEV1)
      CALL GR_EXEC ('rgmap quiet /grey 0 2')
      CALL GR_EXEC ('rgmap quiet /abs')
      WRITE (LIMITS,'(a,4f8.2)') 'lim',-12.,12.,REAL(LOWDEC-90),
     $REAL(HIGHDEC-90)
      CALL GR_EXEC (LIMITS)
      CALL GR_EXEC ('box')
      CALL GR_EXEC ('label "HOUR ANGLE" /x')
      CALL GR_EXEC ('lim 0 1 -1 1')
      CALL GR_EXEC (OBSERVATORY)
      CALL GR_EXEC (INTTIME)
      CALL GR_EXEC ('set orie 90')
      CALL GR_EXEC ('dra tex -.1 0 "DECLINATION" 5 /user')
      CALL GR_EXEC ('set orie 0')
      CALL SIC_DELVARIABLE ('field2',.FALSE.,ERROR)
      CALL SIC_DEF_REAL ('field1',ELEVATION,2,DIM,.TRUE.,ERROR)
      CALL GR_EXEC ('pen 0 /wei 5')
      CALL GR_EXEC ('rgdata field1 /var')
      CALL GR_EXEC ('lim /rg')
      CALL GR4_LEVELS (2,LEV2)
      CALL GR_EXEC ('rgmap quiet /abs')
      CALL GR_EXEC ('pen 0 /wei 1')
      CALL SIC_DELVARIABLE ('field1',.FALSE.,ERROR)
99    RETURN
      END
*
      SUBROUTINE AZTORA (ALONG,ELAT,HA,OBSLAT,DECUR,DRA,DDEC)
      INCLUDE 'const.inc'
*
      REAL*8   ALONG, ELAT, DECUR
      REAL*8   OBSLAT, DDEC1, DLO, DLA, HA
      REAL*8   SINR, COSR, SINE, COSE, SINPA, COSPA, SINLO
      REAL*8   COSLO, SINLA, COSLA
      REAL*8   SINDD, COSDD, COSDRA, SINDRA
      REAL*8   PA, SINOL, COSOL, SINDEC, COSDEC
      REAL*4   DRA, DDEC
*
      DLO = ALONG / 3600. * DEG_TO_RAD
      DLA = ELAT /3600. * DEG_TO_RAD
      HA = HA*DEG_TO_RAD
*
      SINDEC   = DSIN(DECUR)
      COSDEC   = DCOS(DECUR)
      SINOL  = DSIN(OBSLAT)
      COSOL  = DCOS(OBSLAT)
*
      SINR   = DSIN(HA)
      COSR   = DCOS(HA)
      SINE   = SINDEC*SINOL+COSDEC*COSOL * COSR
      COSE   = DSQRT(1.D0-SINE*SINE)
      SINPA  = COSOL*SINR/COSE
      COSPA  = (SINOL-SINDEC*SINE)/(COSDEC*COSE)
      PA     = DATAN2(SINPA,COSPA)
      SINLO  = DSIN(DLO)
      COSLO  = DCOS(DLO)
      SINLA  = DSIN(DLA)
      COSLA  = DCOS(DLA)
      SINDD  = COSPA*SINLA-SINPA*COSLA*SINLO
      DDEC1  = DASIN(SINDD)
      COSDD  = DCOS(DDEC1)
      DDEC   = DDEC1/SEC_TO_RAD          ! DDEC in arcsec
      COSDRA = COSLA*COSLO/COSDD
      IF(ABS(SINPA).LT.1D-6) THEN
         DRA = -DLO/SEC_TO_RAD           ! DRA  in arcsec
      ELSE
         SINDRA = -(SINLA-COSPA*SINDD)/(SINPA*COSDD)
         DRA = DATAN2(SINDRA,COSDRA)/SEC_TO_RAD! DDEC in arcsec
      ENDIF
      RETURN
      END
*
C ----------------------------------------------------------------------
C
C     WRITE_NMB()
C
C ----------------------------------------------------------------------
C
      SUBROUTINE WRITE_NMB(NMB_NAME,DECL,AZ,EL,NSUB,WOB,HOURM,
     $LATI,NCHAN_NMB,OUTPUT)
*
      INCLUDE 'nmb.inc'
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
      REAL*8   AZ ,EL, DECL, HOURM, LATI
      REAL*4   WOB, OUTPUT(*)
      INTEGER  NSUB, NCHAN_NMB
      CHARACTER*(*) NMB_NAME
*
      INTEGER*4 IP_OUTPUTSTREAM,POINTER,IP,IADDR
      INTEGER*4 OUTPUTSTREAM
      REAL*8       SCANTIM
      REAL*8       SINDEC,COSDEC,COSH,SINL,COSL,HA,SIN_EL
      INTEGER      LENC,NFILE
      INTEGER      I,J,K
      INTEGER      IREC,NREC_SUB
      CHARACTER*40 FILENAME
      INTEGER      NUMBER, ILAST, LEN_RECORD, ISUBSCAN
      REAL*4       SUBSCAN_EPOCH,SUBSCAN_STARTTIME, DX
      REAL*4       SCAN_COORDX, SCAN_COORDY,LST,EL_RECORD
      INTEGER      SIC_GETVM, IER
      INTEGER      IM,IMON(12),ID,N,IDATE
      CHARACTER*40 DATE_OF_REDUC
      CHARACTER*3  MONTH(12)
      CHARACTER*120 STRING
      INTEGER     NREF
      PARAMETER   (NREF=19)
      INTEGER      DATE_REF(NREF)        ! copied from nic_input_nmb
      INTEGER     NCHAN_REF(NREF)        ! "
      DATA   DATE_REF    /19930208, 19940401, 19941201,
     $19950201, 19950308, 19960111, 19960201, 19960228,
     $19961127, 19970123, 19970423, 19980107, 19980125,
     $19981201, 19990205, 19990208.5, 19990211.06, 19991203,
     $20000124/
*
      DATA   NCHAN_REF
     $/ 7,  7, 19,  7, 19,  7, 19, 19,  7, 19,
     $37, 37, 37, 37, 37, 37, 37, 37, 37/
*
      DATA MONTH /'JAN','FEB','MAR','APR','MAY','JUN',
     $'JUL','AUG','SEP','OCT','NOV','DEC'/
      DATA     IMON/31,28,31,30,31,30,31,31,30,31,30,31/
*
*
      NFILE = LENC(NMB_NAME)
      FILENAME = NMB_NAME
      FILENAME(NFILE+1:NFILE+1) = CHAR(0)
      CALL OPEN_NEW_STREAM(FILENAME,OUTPUTSTREAM)
      IP_OUTPUTSTREAM=POINTER(OUTPUTSTREAM,MEMORY)
*
* Header section
* --------------
      DO I=1,500                         ! raz header structure
         HEADER_STRUCT(I) = 0
      ENDDO
*
c      HEADER_PREHEADER(1)=
c     $'>DATE 00-00-00  processed by NIC'//CHAR(13)//CHAR(0)
      HEADER_PREHEADER(2)=
     $'>TELESCOPE Pico Veleta  30m MRT'//CHAR(13)//CHAR(0)
      HEADER_PREHEADER(3)=
     $'>RECEIVER'//CHAR(13)//CHAR(0)
      HEADER_PREHEADER(4)=
     $'>WAVELENGTH'//CHAR(13)//CHAR(0)
      HEADER_PREHEADER(5)=
     $'>FORMAT LST AZ EL SCAN_COORD SIGNAL'//CHAR(13)//CHAR(0)
      HEADER_PREHEADER(6)=
     $'>INTEGRATION'//CHAR(13)//CHAR(0)
      HEADER_PREHEADER(7)=
     $'>CHANNELS    '//CHAR(13)//CHAR(0)
      WRITE(HEADER_PREHEADER(7)(11:13),'(i3)') NCHAN_NMB
      HEADER_PREHEADER(8)=
     $'>DUMPS 1'//CHAR(13)//CHAR(0)
      HEADER_PREHEADER(9)=
     $'>DOUBLE'//CHAR(13)//CHAR(0)
      HEADER_PREHEADER(10)=
     $'>CLOCK LST'//CHAR(13)//CHAR(0)
      HEADER_PREHEADER(11)=
     $'>PROJECT unknown'//CHAR(13)//CHAR(0)
*
      I = INDEX(NMB_NAME,'.')
      IF (I.EQ.0) THEN
         HEADER_SCANNUMBER = 9999
      ELSE
         READ (NMB_NAME(1:I-1),'(i4)') HEADER_SCANNUMBER
      ENDIF
      HEADER_SOURCENAME = 'simulation'//CHAR(0)
      HEADER_FREQUENCY  = 243.
      HEADER_GAIN = 1
      HEADER_BETA(1) = DECL/DEG_TO_RAD
      HEADER_LAMBDA(1) = 0               ! ???
      HEADER_AZIMUTH(2) = -AZ/3600.
      HEADER_AZIMUTH(3) = VAZM
      HEADER_ELEVATION(2) = -EL/3600.
      HEADER_SINT = SINT
      HEADER_SBAS = 1                    ! ???
      HEADER_DAZM = 2*AZ/3600.
      HEADER_CHANNELS  = NCHAN_NMB
      IF (WOB.GT.0) THEN
         HEADER_WOBBLINGAMPLI  = WOB
      ELSE
         HEADER_WOBBLINGAMPLI  = WOBBLER_THROW
      ENDIF
      HEADER_WOBBLINGPERIOD = 0.5
      HEADER_WOBBLINGDEVICE = 1
*
*  reduction date
      CALL SIC_DATE (DATE_OF_REDUC)
      DO I=1,12
         IF (DATE_OF_REDUC(4:6).EQ.MONTH(I)) IM = I
      ENDDO
      WRITE(DATE_OF_REDUC(4:),'(I2.2,A5)') IM, DATE_OF_REDUC(7:11)
      HEADER_PREHEADER(1)=
     $'>DATE '//DATE_OF_REDUC(1:10)//
     $'  processed by NIC'//CHAR(13)//CHAR(0)
*
* Date of acquisition: must correspond to a bolometer configuration
      IDATE = 0
      DO I=NREF,1,-1
         IF (NCHAN_REF(I).EQ.NCHAN_NMB.AND.IDATE.EQ.0) THEN
            IDATE = DATE_REF(I)
         ENDIF
      ENDDO
      IF (IDATE.EQ.0) THEN
         STRING = ' Bolometer configuration not available'
         CALL MESSAGE(2,1,'FIELD',STRING)
         IDATE = DATE_REF(NREF)
      ENDIF
      HEADER_DATE(1) = IDATE/10000
      HEADER_DATE(2) = (IDATE - HEADER_DATE(1)*10000)/100
      HEADER_DATE(3) = IDATE - HEADER_DATE(1)*10000
     $- HEADER_DATE(2)*100
*
      SCANTIM = (64./3600*AZ*2/(VAZM*60)*NSUB)/60.   ! in hour
      N = 365
      IF (MOD(HEADER_DATE(1),4).EQ.0) THEN
         IMON(2)=29
         N = 366
      ENDIF
      ID = HEADER_DATE(3)
      IF (HEADER_DATE(2).GT.1) THEN
         DO I = 1,HEADER_DATE(2)-1
            ID = ID+IMON(I)
         ENDDO
      ENDIF
      HEADER_EPOCH  = HEADER_DATE(1)+ ID/FLOAT(N)
      CALL WRITEHEADERINFORMATION(HEADER_STRUCT,MEMORY(IP_OUTPUTSTREAM))
*
* Subscan section
* ---------------
*
      SUBSCAN_EPOCH = HEADER_DATE(1)+ ID/FLOAT(N)
      ILAST = 0
      DO I=1,NSUB
         SUBSCAN_STARTTIME = (HOURM - SCANTIM*(NSUB/2-I+1)/NSUB)*3600.
     $   - 1/60.
         NUMBER = I
         IF (I.EQ.NSUB) ILAST = 1
         CALL WRITESUB(NUMBER,SUBSCAN_EPOCH,SUBSCAN_STARTTIME,
     $   ILAST, MEMORY(IP_OUTPUTSTREAM))
      ENDDO
*
* Record section
* ---------------
      LEN_RECORD = 5 + NCHAN_NMB         ! lst,az,el,2*scan_coord + n channels
      IER = SIC_GETVM(LEN_RECORD,IADDR)
      IP = POINTER(IADDR,MEMORY)
*
      NREC_SUB = NINT(2*AZ/(VAZM*60.))*2 + 1
      DX = VAZM*60.*0.5
!
      ILAST = 0
      SINDEC = SIN(DECL)
      COSDEC = COS(DECL)
      SINL   = SIN(LATI)
      COSL   = COS(LATI)
      DO ISUBSCAN=1,NSUB
         SCAN_COORDY = (-EL+(ISUBSCAN-1)*SINT)*10.
         SUBSCAN_STARTTIME = (HOURM - SCANTIM*(NSUB/2-ISUBSCAN+1)/NSUB)
     $   *3600.
         DO IREC=1,NREC_SUB
            K = (ISUBSCAN-1)*NREC_SUB+IREC
            IF (MOD(ISUBSCAN,2).EQ.0) THEN
               SCAN_COORDX = (AZ-(IREC-1)*DX)*10.
            ELSE
               SCAN_COORDX = (-AZ+(IREC-1)*DX)*10.
            ENDIF
            LST = SUBSCAN_STARTTIME + 0.5*(IREC-1)
            CALL R4TOR4(LST,MEMORY(IP),1)
            HA  = (LST/3600.)*15
            COSH = COS(HA*DEG_TO_RAD)
            SIN_EL= SINDEC*SINL+COSDEC*COSH*COSL
            EL_RECORD = ASIN(SIN_EL)/DEG_TO_RAD
            CALL R4TOR4(0.,MEMORY(IP+1),1)
            CALL R4TOR4(EL_RECORD,MEMORY(IP+2),1)
            CALL R4TOR4(SCAN_COORDX,MEMORY(IP+3),1)
            CALL R4TOR4(SCAN_COORDY,MEMORY(IP+4),1)
            IF (K.EQ.NRECORD) ILAST = 1
            DO J=1,NCHAN_NMB
               CALL R4TOR4(OUTPUT(K+(J-1)*NRECORD),
     $         MEMORY(IP+J+4),1)
            ENDDO
            CALL WRITERECORD(MEMORY(IP),NCHAN_NMB,ILAST,
     $      MEMORY(IP_OUTPUTSTREAM))
         ENDDO
      ENDDO
*
      CALL FREE_VM(LEN_RECORD,IADDR)
      CALL CLOSE_STREAM(OUTPUTSTREAM)
      END
*
*
*
      SUBROUTINE CORRELATED_NOISE(RMS,RMIN,NCH,CHAN_AZ_OFF,
     $CHAN_EL_OFF,CHAN_GAIN,CHAN_FLAG,
     $NEW_CHAN_NUM,A,R,E,OUTPUT,ERROR)
C ----------------------------------------------------------------------
C     support routine for option NOISE in command FIELD
C     11-FEB-2000 H. Wiesemeyer
C
C     Arguments:
C
C ----------------------------------------------------------------------
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
      INCLUDE 'inc:format.inc'
*
      INTEGER I, IDUM, II, J, K, L, MAXPAIR, N_CHANPAIRS, NCH
      INTEGER NCH_OLD, NEW_CHAN_NUM(*)
      PARAMETER(MAXPAIR = NREC*(NREC-1)/2)
      REAL*4 CHAN_AZ_OFF(*),CHAN_EL_OFF(*),CHAN_GAIN(*)
      REAL*4 A(NCH,NCH), CHAN_DIST(MAXPAIR), DMIN, D2
      REAL*4 E(NRECORD,NCH), FACTOR, GAUSSDEV, MEAN, ONE
      REAL*4 R(NCH,NCH), RMIN, RMS, SUM, ZERO, ONE_OVER_E_WIDTH
      REAL*4 Z, OUTPUT(*)
      PARAMETER (ONE = 1.0, ZERO = 0.0)
      LOGICAL ERROR, CHAN_FLAG(*)
      DATA IDUM /-12/
*
      INTEGER     IDATE,NREF
      PARAMETER   (NREF=19)
      INTEGER     DATE_REF(NREF), NCHAN_REF(NREF)
      CHARACTER   STRING*120, TEXT*80
      REAL*8      DATE
      DATA   DATE_REF    /19930208, 19940401, 19941201,
     $19950201, 19950308, 19960111, 19960201, 19960228,
     $19961127, 19970123, 19970423, 19980107, 19980125,
     $19981201, 19990205, 19990208.5, 19990211.06, 19991203,
     $20000124/
*
      DATA   NCHAN_REF
     $/ 7,  7, 19,  7, 19,  7, 19, 19,  7, 19,
     $37, 37, 37, 37, 37, 37, 37, 37, 37/
*
* Date of acquisition: must correspond to a bolometer configuration
*
      IDATE = 0
      DO I=NREF,1,-1
         IF (NCHAN_REF(I).EQ.NCH.AND.IDATE.EQ.0) THEN
            IDATE = DATE_REF(I)
         ENDIF
      ENDDO
      IF (IDATE.EQ.0) THEN
         STRING = ' Bolometer configuration not available'
         CALL MESSAGE(2,1,'FIELD',STRING)
         IDATE = DATE_REF(NREF)
      ENDIF
      DATE = DBLE(IDATE)
*
      CALL INIT_OFFSET_GAIN(CHAN_AZ_OFF,CHAN_EL_OFF,CHAN_GAIN,
     $NCH,TEXT,DATE)
*
      ERROR = .FALSE.
      IF (RMS.EQ.ZERO) THEN
         DO L = 1, NRECORD
            DO I = 1, NCH
               CALL R4TOR4(0,OUTPUT(L+(I-1)*NRECORD),1)
            ENDDO
         ENDDO
         RETURN
      ENDIF
c
c     compute vectors of uncorrelated noise
c
      N_CHANPAIRS = NCH*(NCH-1)/2
      DO I = 1, NCH
         SUM = ZERO
         MEAN = ZERO
         DO J = 1, NRECORD
            E(J,I) = GAUSSDEV(IDUM)
            SUM = SUM+E(J,I)*E(J,I)
            MEAN = MEAN+E(J,I)
         ENDDO
         MEAN = MEAN/REAL(NRECORD)
         DO J = 1, NRECORD
            E(J,I) = (E(J,I)-MEAN)/SQRT(SUM)
         ENDDO
      ENDDO
c
c     compute distances of channel pairs
c
      DO I = 1, NCH
         CHAN_FLAG(I) = .FALSE.
         IF (CHAN_GAIN(I).EQ.ZERO) CHAN_FLAG(I) = .TRUE.
      ENDDO
      K = 0
      DMIN = 1.E06
      DO I = 1, NCH-1
         IF (.NOT.CHAN_FLAG(I)) THEN
            DO J = I+1, NCH
               IF (.NOT.CHAN_FLAG(J)) THEN
                  D2 = (CHAN_AZ_OFF(I)-CHAN_AZ_OFF(J))**2
     $            +(CHAN_EL_OFF(I)-CHAN_EL_OFF(J))**2
                  K = K+1
                  CHAN_DIST(K) = SQRT(D2)
                  DMIN = MIN(DMIN,CHAN_DIST(K))
               ENDIF
            ENDDO
         ENDIF
      ENDDO
c
      L = 0
      DO I = 1, NCH
         IF (.NOT.CHAN_FLAG(I)) THEN
            L = L+1
            NEW_CHAN_NUM(I) = L
         ENDIF
      ENDDO
      NCH_OLD = NCH
      NCH = L
c
c     compute correlation coefficients as a function of channel distance
c     according to exponential decay (decay constant determined by
c     observer)
c
      RMIN = ABS(RMIN)
      IF (RMIN.EQ.ZERO.OR.RMIN.EQ.ONE) THEN
         DO I = 1, NCH-1
            DO J = I+1, NCH
               R(I,J) = RMIN
            ENDDO
         ENDDO
      ELSE
         K = 0
         ONE_OVER_E_WIDTH = -DMIN/LOG(RMIN)
         DO I = 1, NCH-1
            DO J = I+1, NCH
               K = K+1
               R(I,J) = EXP(-CHAN_DIST(K)/ONE_OVER_E_WIDTH)
            ENDDO
         ENDDO
      ENDIF
c
c     compute matrix elements for linear transformation of uncorrelated
c     noise vectors to correlated channel signals
c
      A(1,1) = ONE
      DO J = 2, NCH
         A(J,1) = R(1,J)
      ENDDO
      IF (RMIN.NE.ONE) THEN
         A(2,2) = SQRT(ONE-A(2,1)*A(2,1))
         DO I = 2, NCH-1
            DO J = I+1, NCH
               A(J,I) = R(I,J)
               DO K = 1, I-1
                  A(J,I) = A(J,I)-A(I,K)*A(J,K)
               ENDDO
               A(J,I) = A(J,I)/A(I,I)
            ENDDO
            SUM = ONE
            DO K = 1, I
               SUM = SUM-A(I+1,K)*A(I+1,K)
            ENDDO
            A(I+1,I+1) = SQRT(SUM)
         ENDDO
      ELSE
         DO I = 1, NCH
            A(I,1) = ONE
            DO J = 2, I
               A(I,J) = ZERO
            ENDDO
         ENDDO
      ENDIF
c
c     compute output correlated signal vectors
c
      FACTOR = RMS*SQRT(REAL(NRECORD))
      DO I = 1, NCH_OLD
         IF (.NOT.CHAN_FLAG(I)) THEN
            II = NEW_CHAN_NUM(I)
            DO L = 1, NRECORD
               SUM = ZERO
               DO J = 1, II
                  SUM = SUM+A(II,J)*E(L,J)
               ENDDO
               Z = INT(SUM*FACTOR/CHAN_GAIN(I))
               CALL R4TOR4(Z,OUTPUT(L+(I-1)*NRECORD),1)
            ENDDO
         ELSE
            DO L = 1, NRECORD
               Z = 0.
               CALL R4TOR4(Z,OUTPUT(L+(I-1)*NRECORD),1)
            ENDDO
         ENDIF
      ENDDO
      NCH = NCH_OLD
      RETURN
      END
*
*
*
      FUNCTION GAUSSDEV(IDUM)
      INTEGER ISET, IDUM
      REAL*4 GAUSSDEV, V1, V2, R, RAN, FAC, GSET
      SAVE GSET, ISET
      DATA ISET/0/
      IF (ISET.EQ.0) THEN
1        V1=2.*RAN(IDUM)-1.
         V2=2.*RAN(IDUM)-1.
         R=V1**2+V2**2
         IF(R.GE.1.)GO TO 1
         FAC=SQRT(-2.*LOG(R)/R)
         GSET=V1*FAC
         GAUSSDEV=V2*FAC
         ISET=1
      ELSE
         GAUSSDEV=GSET
         ISET=0
      ENDIF
      RETURN
      END
!
      FUNCTION RAN(IDUM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER IDUM,IA,IM,IQ,IR,NTAB,NDIV
      REAL RAN,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     $NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2E-7,RNMX=1.-EPS)
      INTEGER J,K,IV(NTAB),IY
      SAVE IV,IY
      DATA IV /NTAB*0/, IY /0/
      IF (IDUM.LE.0.OR.IY.EQ.0) THEN
         IDUM=MAX(-IDUM,1)
         DO 11 J=NTAB+8,1,-1
            K=IDUM/IQ
            IDUM=IA*(IDUM-K*IQ)-IR*K
            IF (IDUM.LT.0) IDUM=IDUM+IM
            IF (J.LE.NTAB) IV(J)=IDUM
11       CONTINUE
         IY=IV(1)
      ENDIF
      K=IDUM/IQ
      IDUM=IA*(IDUM-K*IQ)-IR*K
      IF (IDUM.LT.0) IDUM=IDUM+IM
      J=1+IY/NDIV
      IY=IV(J)
      IV(J)=IDUM
      RAN=MIN(AM*IY,RNMX)
      RETURN
      END
*
