* Last processed by NICE on 12-Jun-2000 15:53:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
* Customized for :  EEEI, HPUX, UNIX, MOTIF, F77
* Customized for : EEEI HPUX UNIX
* Written by r.neri January-1995
* Updated by h.wiesemeyer February-2000
* (with corrected formula for calibration voltage)
*
* Computes opacity for multibeam bolometer channels
* Normally channel 1 is used for opacity evaluation.
*
      SUBROUTINE SOLVE_SKYDIP (LINE,ERROR)
*
*----------------------------------------------------------------------
* NIC Support routine for command
*        SOLVE /CHANNEL /FILE [file_name]
*
* Arguments :
*
*    LINE C*(*) Command line Input
*    ERROR L Logical error flag
*    Output
*
*----------------------------------------------------------------------
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
*
      INTEGER*4 IPSIGNAL,WORK1
      INTEGER*4 IP_SKYUT, IP_SKYTAU, POINTER
      INTEGER*4 IP_DUMUT, IP_DUMTAU
      INTEGER*4 P1, P2, P3
      INTEGER*4 DESC(9)
      INTEGER         PARIS(1), MADRID(1), TOKYO(1), MY_BAND_ELEM
      LOGICAL ERROR, SIC_PRESENT, SKYDIP_NEW, FOUND, FILE_OPT
      LOGICAL ANY_BAND, USE_BAND, USER_DEFINED_BAND
      CHARACTER*(*) LINE
      CHARACTER*70 DRAW_TEXT, DRAW_RELO, DRAW_LINE, FILENAME
      CHARACTER STR*1, MODE*5
      INTEGER CHAN, LEN, LENC, IER, SIC_GETVM, NC
      REAL VAL1, VAL2, USER_YMAX, USER_YMIN, USER_XMAX, USER_XMIN
      REAL TAU, AZI, SIGMA
      REAL*8 CURRENT_UT
      INTEGER ICOUNTS_TO_K,IAZI, DIM(4)
      CHARACTER ITALIC*3,GREEK*2,HAUT*3,BAS*3
      COMMON /BAND/ P1, P2, P3, PARIS, MADRID, TOKYO, MY_BAND_ELEM
      DATA USER_DEFINED_BAND /.TRUE./
      SAVE /BAND/
      SAVE USER_DEFINED_BAND
*------------------------------------------------------------------------
*
* Any file ?
*
      IF (IADDRRECORD.EQ.0) THEN
         CALL MESSAGE(2,1,'SOLVE','No input file opened')
         ERROR = .TRUE.
         RETURN
      ENDIF
*
* Enough subscans ?
*
      IF (NSUBSCAN.LT.3) THEN
         CALL MESSAGE(2,1,'SOLVE','Not enough data points')
         ERROR = .TRUE.
         RETURN
      ENDIF
*
* Option /CHANNEL
*
      CHAN = 1
      IF (SIC_PRESENT(1,0) .AND. SIC_PRESENT(1,1)) THEN
         CALL SIC_I4 (LINE,1,1,CHAN,.TRUE.,ERROR)
         IF (CHAN.GT.NCHAN.OR.CHAN.LT.1) THEN
            CALL MESSAGE (2,1,'SOLVE','Bad channel number ')
            ERROR = .TRUE.
            RETURN
         ENDIF
      ENDIF
*
* Option /BAND
*
      USE_BAND = .FALSE.
      IF (SIC_PRESENT(8,0)) THEN
         IF (SIC_PRESENT(8,1)) THEN
            CALL SIC_KE(LINE,8,1,MODE,NC,.TRUE.,ERROR)
            IF (MODE.NE."MAMBO") THEN
               CALL MESSAGE (2,1,'SOLVE',
     $         'Only Mambo II bandpass available... ')
               ERROR = .TRUE.
               RETURN
            ENDIF
            USE_BAND = .TRUE.
         ELSE
            CALL MESSAGE (2,1,'SOLVE','Using MAMBO II bandpass')
            USE_BAND = .TRUE.
         ENDIF
      ENDIF
      IF (USE_BAND) THEN
         USER_DEFINED_BAND = .FALSE.
         CALL DEFINE_BANDPASS
      ELSE
         IF (.NOT.USER_DEFINED_BAND) THEN
            CALL SIC_DELVARIABLE('BAND_FREQ',.FALSE.,ERROR)
            CALL SIC_DELVARIABLE('BAND_TRANS',.FALSE.,ERROR)
            CALL SIC_DELVARIABLE('BAND_ATMOS',.FALSE.,ERROR)
         ENDIF
      ENDIF
*
* Option /FILE
*
      FILENAME = 'nic.sky'
      FILE_OPT = .FALSE.
      IF (SIC_PRESENT(2,0)) THEN
         FILE_OPT = .TRUE.
         IF (SIC_PRESENT(2,1)) THEN
            CALL SIC_CH (LINE,2,1,FILENAME,NC,.TRUE.,ERROR)
            IF (ERROR) GOTO 99
         ENDIF
!DB:03-99 CALL SIC_LOWER(FILENAME)
         OPEN (UNIT=2,FILE=FILENAME,FORM='formatted',STATUS='old',
     $   IOSTAT=IER)
         IF (IER.NE.0) THEN
            OPEN (UNIT=2,FILE=FILENAME,FORM='formatted',STATUS='new',
     $      IOSTAT=IER)
            WRITE (2,'(a)') '!'
            WRITE (2,'(a)') '!   Universal Time  Opacity   '//
     $      'Water (mm) Azimuth (deg)  Counts/K    Scan'
            WRITE (2,'(a)') '!'
         ELSE
            DO WHILE (.TRUE.)
               READ(2,'(a)',END=10) DRAW_TEXT
            ENDDO
10       ENDIF
         IF (IER.NE.0) THEN
            CALL MESSAGE (2,1,'FIND','Not able to create file')
            RETURN
         ENDIF
      ENDIF
c     CALL MESSAGE (4,1,'SOLVE','Skydip is experimental...')
*
      IPSIGNAL = IP_SIGNAL + NRECORD*(CHAN-1)
*
      DRAW_TEXT = 'dra tex'
      DRAW_RELO = 'dra rel'
      DRAW_LINE = 'dra lin'
      DRAW_RELO(27:30) = '/use'
      DRAW_LINE(27:30) = '/use'
      ITALIC = CHAR(92)//CHAR(92)//'i'
      GREEK = CHAR(92)//'g'
      BAS = CHAR(92)//CHAR(92)//'d'
      HAUT = CHAR(92)//CHAR(92)//'u'
*
      CALL GR_EXEC ('clear plot')
      CALL GR_EXEC1 ('set box 4 20 4 16')
      CALL GR_EXEC1 ('set exp .8')
*
      CALL FIT_OPACITY (MEMORY(IP_EL),MEMORY(IPSIGNAL),CHAN,
     $TAU,WATER_VAPOUR,ANY_BAND,COUNTS_TO_KELVIN,SIGMA,ERROR)
      IF (ERROR) RETURN
      COUNTS_TO_KELVIN = 1/COUNTS_TO_KELVIN
c      CALL SIC_DESCRIPTOR ('SKYDIP_OPACITY',DESC,FOUND)
c      IF (FOUND) THEN
c         CALL SIC_DELVARIABLE('skydip_time',.FALSE.,ERROR)
c         CALL SIC_DELVARIABLE('skydip_opacity',.FALSE.,ERROR)
c      ENDIF
      IF (NUMBER_OF_SKYDIPS.GT.0) THEN
         IP_SKYUT = POINTER(SKYDIP_MEM,MEMORY)
         FOUND = SKYDIP_NEW(MEMORY(IP_SKYUT),MEMORY(IP_LST),
     $   NUMBER_OF_SKYDIPS,CURRENT_UT)
         IF (FOUND) THEN
            IER = SIC_GETVM(3*NUMBER_OF_SKYDIPS,WORK1)
            IF (IER.NE.1) THEN
               CALL MESSAGE (2,1,'SOLVE','Memory allocation failure')
               RETURN
            ENDIF
            IP_DUMUT = POINTER(WORK1,MEMORY)   ! double precision
            IP_DUMTAU = IP_DUMUT + 2*NUMBER_OF_SKYDIPS
            IP_SKYTAU = IP_SKYUT + 2*NUMBER_OF_SKYDIPS     ! double precision
            CALL R8TOR8(MEMORY(IP_SKYUT),MEMORY(IP_DUMUT),
     $      NUMBER_OF_SKYDIPS)
            CALL R4TOR4(MEMORY(IP_SKYTAU),MEMORY(IP_DUMTAU),
     $      NUMBER_OF_SKYDIPS)
            CALL FREE_VM (3*NUMBER_OF_SKYDIPS,SKYDIP_MEM)
            IER = SIC_GETVM(3*(NUMBER_OF_SKYDIPS+1),SKYDIP_MEM)
            IP_SKYUT = POINTER(SKYDIP_MEM,MEMORY)
            IP_SKYTAU = IP_SKYUT + 2*(NUMBER_OF_SKYDIPS+1) ! double precision
            IF (IER.NE.1) THEN
               CALL MESSAGE (2,1,'SOLVE','Memory allocation failure')
               RETURN
            ENDIF
            CALL SKYDIP_NEXT (MEMORY(IP_DUMUT),MEMORY(IP_SKYUT),
     $      MEMORY(IP_LST), MEMORY(IP_DUMTAU), MEMORY(IP_SKYTAU),TAU,
     $      NUMBER_OF_SKYDIPS,NUMBER_OF_SKYDIPS+1)
            CALL FREE_VM (3*NUMBER_OF_SKYDIPS,WORK1)
            NUMBER_OF_SKYDIPS = NUMBER_OF_SKYDIPS+1
         ENDIF
      ELSE
         IER = SIC_GETVM(3,SKYDIP_MEM)
         IF (IER.NE.1) THEN
            CALL MESSAGE (2,1,'SOLVE','Memory allocation failure')
            RETURN
         ENDIF
         IP_SKYUT = POINTER(SKYDIP_MEM,MEMORY)
         IP_SKYTAU = IP_SKYUT + 2        ! double precision
         CALL SKYDIP_FIRST(MEMORY(IP_SKYUT),MEMORY(IP_LST),
     $   MEMORY(IP_SKYTAU),TAU,CURRENT_UT)
         NUMBER_OF_SKYDIPS = 1
      ENDIF
      IF (FILE_OPT) THEN
         STR = ' '
         IF (SIGMA.GT.1.5) STR = '!'
         WRITE (2,'(a,F17.8,3(x,F7.3,4x),2(3x,i6))') STR, CURRENT_UT,
     $   TAU,WATER_VAPOUR,AZI(MEMORY(IP_LST+NRECORD)),
     $   NINT(COUNTS_TO_KELVIN),SCAN_NUMBER
         CLOSE (UNIT=2)
      ENDIF
      IF (FOUND) THEN
         CALL SIC_DESCRIPTOR ('SKYDIP_OPACITY',DESC,FOUND)
         IF (FOUND) THEN
            CALL SIC_DELVARIABLE('skydip_time',.FALSE.,ERROR)
            CALL SIC_DELVARIABLE('skydip_opacity',.FALSE.,ERROR)
         ENDIF
         CALL SIC_DEF_DBLE ('skydip_time',MEMORY(IP_SKYUT),1,
     $   NUMBER_OF_SKYDIPS,.TRUE.,ERROR)
         CALL SIC_DEF_REAL ('skydip_opacity',MEMORY(IP_SKYTAU),1,
     $   NUMBER_OF_SKYDIPS,.TRUE.,ERROR)
         IF (SECTION(6).EQ.1) CALL SKYDIP_VARIABLES(.FALSE.,.FALSE.)
         IF (SECTION(6).EQ.2) CALL SKYDIP_VARIABLES(.FALSE.,.TRUE.)
      ENDIF
      DIM(1)=1
      CALL SIC_DELVARIABLE('ATM_TEMP',.FALSE.,ERROR)
      CALL SIC_DEF_REAL ('ATM_TEMP',TATM,0,DIM,.TRUE.,ERROR)
*
      SCAN_OPACITY = TAU
      CALL SIC_GET_REAL ('USER_YMAX',USER_YMAX,ERROR)
      CALL SIC_GET_REAL ('USER_YMIN',USER_YMIN,ERROR)
      CALL SIC_GET_REAL ('USER_XMAX',USER_XMAX,ERROR)
      CALL SIC_GET_REAL ('USER_XMIN',USER_XMIN,ERROR)
      LEN = (USER_YMAX-USER_YMIN)/2
      VAL1 = (2*USER_YMAX+USER_YMIN)/3
      WRITE (DRAW_TEXT(17:26),'(f10.1)') VAL1
      WRITE (DRAW_RELO(17:26),'(f10.1)') VAL1
      WRITE (DRAW_LINE(17:26),'(f10.1)') VAL1
      VAL2 = (2*USER_XMAX+USER_XMIN)/3
      WRITE (DRAW_TEXT(9:15),'(f7.3)') VAL2
      WRITE (DRAW_RELO(9:15),'(f7.3)') VAL2
      WRITE (DRAW_LINE(9:15),'(f7.3)') VAL2-8*(VAL2-USER_YMIN)/10.
      DRAW_TEXT(28:) = '"'//ITALIC//GREEK//'t'//
     $BAS//'ATM'//HAUT//' =      " 5 /user'
      WRITE (DRAW_TEXT(46:51),'(f6.3)') TAU
      CALL GR_EXEC1 (DRAW_TEXT)
      FITTED_OPACITY = TAU
      WRITE (DRAW_TEXT(9:15),'(f7.4)') .7
      WRITE (DRAW_TEXT(17:24),'(f7.4)') -.152
      DRAW_TEXT(25:) = '"'//ITALIC//'SCAN =      " 6 /box 9'
      WRITE (DRAW_TEXT(35:40),'(i6)') SCAN_NUMBER
      CALL GR_EXEC1 (DRAW_TEXT)
      WRITE (DRAW_TEXT(17:24),'(f7.4)') -.752
      LEN = LENC(SOURCE_NAME)
      DRAW_TEXT(25:) = '"'//ITALIC//'SOURCE = '//
     $SOURCE_NAME(1:LEN)//' " 6 /box 9'
      CALL GR_EXEC1 (DRAW_TEXT)
      WRITE (DRAW_TEXT(17:24),'(f7.4)') -1.352
      DRAW_TEXT(25:) = '"'//ITALIC//'CHANNEL =    " 6 /box 9'
      WRITE (DRAW_TEXT(38:40),'(i3)') CHAN
      CALL GR_EXEC1 (DRAW_TEXT)
      WRITE (DRAW_TEXT(17:24),'(f7.4)') -1.952
      DRAW_TEXT(25:) = '"'//ITALIC//'AZIMUTH =     "     6 /box 9'
      CALL R4TOI4(AZI(MEMORY(IP_AZ)),IAZI,1)
      WRITE (DRAW_TEXT(39:41),'(i3)') IAZI
      CALL GR_EXEC1 (DRAW_TEXT)
      WRITE (DRAW_TEXT(17:24),'(f7.4)') -2.552
      DRAW_TEXT(25:) = '"'//ITALIC//'WATER =      mm" 6 /box 9'
      WRITE (DRAW_TEXT(37:40),'(f4.1)') WATER_VAPOUR
      CALL GR_EXEC1 (DRAW_TEXT)
      WRITE (DRAW_TEXT(17:24),'(f7.4)') -3.152
      DRAW_TEXT(25:) = '"'//ITALIC//'FORWARD_EFF'//
     $' =     " 6 /box 9'
      WRITE (DRAW_TEXT(43:46),'(f4.2)') FORWARD_EFFICIENCY
      CALL GR_EXEC1 (DRAW_TEXT)
      WRITE (DRAW_TEXT(17:24),'(f7.4)') -3.752
      DRAW_TEXT(25:) = '"'//ITALIC//'COUNTS/K'//
     $' =       " 6 /box 9'
      ICOUNTS_TO_K = COUNTS_TO_KELVIN
      WRITE (DRAW_TEXT(40:45),'(i6)') ICOUNTS_TO_K
      CALL GR_EXEC1 (DRAW_TEXT)
      IF (ANY_BAND) THEN
         CALL GR_EXEC1 ('set box 20.5 28 4 7.8')
         CALL GR4_GIVE ('X',MY_BAND_ELEM,MADRID(P1))
         CALL GR4_GIVE ('Y',MY_BAND_ELEM,PARIS(P2))
         CALL GR_EXEC1 ('lim * * -.1 1.1')
         CALL GR_EXEC1 ('box p n')
         CALL GR_EXEC1 ('axi yr /label o ONLY')
         CALL GR_EXEC1 ('conn')
         CALL GR_EXEC1 ('label "Frequency (GHz)" /x')
         CALL GR_EXEC1 ('set orie 270')
         CALL GR_EXEC1 ('lim 0 1 -1 1')
         CALL GR_EXEC1 ('dra tex 1.2 0 "Instr.Bandpass" 5 /user')
         CALL GR_EXEC1 ('set orie 0')
         CALL GR_EXEC1 ('set box 20.5 28 7.9 11.7')
         CALL GR4_GIVE ('X',MY_BAND_ELEM,MADRID(P1))
         CALL GR4_GIVE ('Y',MY_BAND_ELEM,TOKYO(P3))
         CALL GR_EXEC1 ('lim * * -.1 1.1')
         CALL GR_EXEC1 ('box n n')
         CALL GR_EXEC1 ('axi yr /label o ONLY')
         CALL GR_EXEC1 ('curve')
         CALL GR_EXEC1 ('set orie 270')
         CALL GR_EXEC1 ('lim 0 1 -1 1')
         CALL GR_EXEC1 ('dra tex 1.2 0 "Atm.Trans." 5 /user')
         CALL GR_EXEC1 ('set orie 0')
      ENDIF
      CALL GR_EXEC1 ('set exp 1.')
      RETURN
*
99    ERROR = .TRUE.
      RETURN
      END
*
      SUBROUTINE DEFINE_BANDPASS
*
      REAL*4 BANDF(30), BANDT(30), BANDA(30)
      LOGICAL ERROR
      INTEGER DIM(4)
*
      DATA BANDF /178.65,183.33,191.67,198.70,201.56,204.69,211.20,
     $216.93,222.92,230.73,235.94,241.15,246.88,254.17,
     $259.38,265.63,273.70,279.69,284.90,293.23,295.83,
     $301.04,304.95,312.50,320.57,330.21,338.80,348.44,
     $360.94,373.96/
      DATA BANDT / .00,.01,.08,.17,.39,.51,.56,.60,.72,.90,
     $1.00,.99,.95,.93,.87,.81,.77,.70,.58,.51,
     $.45,.37,.26,.20,.17,.12,.05,.02,.01,.00 /
      DATA BANDA / .00,.00,.00,.00,.00,.00,.00,.00,.00,.00,
     $.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,
     $.00,.00,.00,.00,.00,.00,.00,.00,.00,.00/
*
      DIM(1) = 30
      CALL SIC_DELVARIABLE('BAND_FREQ',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('BAND_TRANS',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('BAND_ATMOS',.FALSE.,ERROR)
      CALL SIC_DEF_REAL ('BAND_FREQ',BANDF,1,DIM,.TRUE.,ERROR)
      CALL SIC_DEF_REAL ('BAND_TRANS',BANDT,1,DIM,.TRUE.,ERROR)
      CALL SIC_DEF_REAL ('BAND_ATMOS',BANDA,1,DIM,.TRUE.,ERROR)
      RETURN
      END
      SUBROUTINE FIT_OPACITY (EL,SIGNAL,CHAN,TAU,WATER,FOUND3,
     $CALIB,SIGMA,ERROR)
*----------------------------------------------------------------------
* NIC SOLVE support routine
* Compute air masses
* Modified by Michael Bremer November 1995
*
* Arguments :
*     EL       R*4 elevation in degrees      Input
*     SIGNAL   R*4 channel countrate         Input
*     CHAN     I*4 channel number            Input
*     TAU      R*4 atmospheric opacity       Output
*
*----------------------------------------------------------------------
*
      IMPLICIT NONE
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
      INCLUDE 'inc:format.inc'
*
      REAL SIGNAL(NRECORD), EL(NRECORD)
*
      INTEGER*4 I, J, K, L, MAX_DIP, MAX_K, NPAR, CHAN
      INTEGER*4 POINTER, STATUS
      INTEGER*4 BAND_ELEM, MODE, ITER
      INTEGER*4 ARRAY1(9), ARRAY2(9), ARRAY3(9), ARRAY4(9)
      INTEGER*4 P1, P2, P3, ARRAY5(9), ARRAY6(9)
      INTEGER         MADRID(1),PARIS(1), TOKYO(1), MY_BAND_ELEM
      INTEGER         NPOINT
      PARAMETER (BAND_ELEM = 5, NPOINT = 100)
      PARAMETER (MAX_DIP = 20, NPAR = 3) ! not more than 20 dips in a scan
      REAL SKY(MAX_DIP), ELV(MAX_DIP), WEI(MAX_DIP), WORK(MAX_DIP)
      REAL AIR_MASS, AM
      REAL TAU, MINX, MAXX, MINY, MAXY, AMPL
      REAL XM(100), YM(100)
      REAL CALIB
      REAL SIGMA, T_CHOPPER, WATER
      REAL BAND_FREQ(BAND_ELEM), BAND_TRANS(BAND_ELEM)
      REAL BAND_ATMOS(BAND_ELEM)
      LOGICAL FOUND1, FOUND2, FOUND3, FOUND5
      LOGICAL ERROR, ERROR1, ERROR2, ERROR5
      COMMON /BAND/ P1, P2, P3, PARIS, MADRID, TOKYO, MY_BAND_ELEM
      SAVE /BAND/
*
      I = 1                              ! number of data points
      J = 1                              ! data counter (within subscan)
      K = 1                              ! subscan counter
      L = 1
      FOUND1 = .FALSE.
      DO WHILE (.NOT.FOUND1)
         IF (ABS(SIGNAL(L)-BLANKING).GT.0.5) THEN
            ELV(L) = AIR_MASS(EL(L))
            SKY(L) = SIGNAL(L)
            FOUND1 = .TRUE.
         ENDIF
         L = L+1
      ENDDO
      MINX = ELV(L-1)
      MAXX = MINX
      MINY = ABS(SKY(L-1))
      MAXY = MINY
      DO I = L,NRECORD
         FOUND1 = .FALSE.
         AM = AIR_MASS(EL(I))
         IF (NINT(ELV(K)*100).NE.NINT(AM*100)) THEN  ! to avoid rounding errors
            FOUND1 = .TRUE.
            IF (J.NE.1) THEN             ! blanked subscans ?
               SKY(K) = ABS (SKY(K) / J) ! normalize
               MINX = MIN (MINX,ELV(K))
               MAXX = MAX (MAXX,ELV(K))
               MINY = MIN (MINY,SKY(K))
               MAXY = MAX (MAXY,SKY(K))
               J = 1
               K = K+1
               IF (K.GT.MAX_DIP) THEN
                  CALL MESSAGE (4,1,'SOLVE','Too many DIPs ...')
                  ERROR = .TRUE.
                  RETURN
               ENDIF
            ENDIF
            ELV(K) = AM
            SKY(K) = SIGNAL(I)           ! initialize
         ENDIF
         IF (.NOT.FOUND1.AND.ABS(SIGNAL(I)-BLANKING).GT.0.5) THEN
            SKY(K) = SKY(K) + SIGNAL(I)
            J = J+1
         ENDIF
      ENDDO
      IF (J.NE.1) THEN
         SKY(K) = ABS(SKY(K) / J)        ! normalize
         MINX = MIN (MINX,ELV(K))
         MAXX = MAX (MAXX,ELV(K))
         MINY = MIN (MINY,SKY(K))
         MAXY = MAX (MAXY,SKY(K))
      ELSE
         K = K-1
      ENDIF
      MAX_K = K                          ! = number of subscans !!
      IF (MINY.EQ.MAXY.OR.K.EQ.0) THEN
         CALL MESSAGE (4,1,'SOLVE','Crazy data')
         ERROR = .TRUE.
         RETURN
      ENDIF
*
*
      AMPL = MAXX-MINX
      MINX = MINX-.05*AMPL
      MAXX = MAXX+.05*AMPL
      AMPL = MAXY-MINY
      MINY = MINY-.05*AMPL
      MAXY = MAXY+.05*AMPL
*
      CALL GR_LIMI (4,MINX,MAXX,MINY,MAXY)
      CALL GR_EXEC1 ('box')
      CALL GR_EXEC1 ('label "Air mass" /x')
      CALL GR_EXEC1 ('label "Counts" /y')
      CALL GR_SEGM ('plot',ERROR)
      CALL GR4_MARKER (K,ELV,SKY,BLANKING,10D0)
*
* Look if the bandpass variables exist
*
      CALL SIC_DESCRIPTOR ('BAND_FREQ', ARRAY1, FOUND1)
      CALL SIC_DESCRIPTOR ('BAND_TRANS', ARRAY2, FOUND2)
      CALL SIC_DESCRIPTOR ('BAND_ATMOS', ARRAY5, FOUND5)
*
      T_CHOPPER   = 290.0
      IF (CHOPPER_TEMPERATURE.GT.0) T_CHOPPER = CHOPPER_TEMPERATURE
      CHOPPER_TEMPERATURE = T_CHOPPER
      MODE = 2
      DO I = 1, MAX_K
         WEI(I) = 1.0
      ENDDO
*
* No bandpass variables present: run iteration with a tophat profile
*
      FOUND3 = FOUND1.AND.FOUND2
      IF (.NOT. FOUND3) THEN
*
         DO I = 1, BAND_ELEM
            BAND_FREQ(I)  = 200.0 + (I-1.)/(BAND_ELEM-1.)*60.
            BAND_TRANS(I) = 1.0
         ENDDO
*
         CALL SOLV_SKYD (ELV, SKY, WEI, WORK, MAX_DIP, MAX_K,
     $   FORWARD_EFFICIENCY, T_CHOPPER, TATM,
     $   AMBIENT_TEMPERATURE+273.15, AMBIENT_PRESSURE,
     $   SITE_ALTITUDE/1000.,
     $   BAND_FREQ, BAND_TRANS, BAND_ATMOS, BAND_ELEM,
     $   CALIB, TAU, WATER, SIGMA, MODE, ITER, STATUS)
*
      ELSE
*
* Bandpass variables present:
         IF (FOUND3) THEN
*
* Must have same size and one dimension
            IF ((ARRAY1(3).EQ.1) .AND. (ARRAY2(3).EQ.1) .AND.
     $      (ARRAY1(4).EQ.ARRAY2(4)).AND. (ARRAY1(4).GE.1) .AND.
     $      (ARRAY5(3).EQ.1) .AND. (ARRAY1(4).EQ.ARRAY5(4))) THEN
*
* Must allow to create FMT_R4 incarnations
               CALL SIC_INCARNATE (FMT_R4, ARRAY1, ARRAY3, ERROR1)
               CALL SIC_INCARNATE (FMT_R4, ARRAY2, ARRAY4, ERROR2)
               CALL SIC_INCARNATE (FMT_R4, ARRAY5, ARRAY6, ERROR5)
               IF ((.NOT. ERROR1) .AND. ( .NOT. ERROR2)
     $         .AND. (.NOT. ERROR5)) THEN
*
* OK, this should be safe now
                  P1 = POINTER(ARRAY3(2), MADRID)
                  P2 = POINTER(ARRAY4(2), PARIS)
                  P3 = POINTER(ARRAY6(2), TOKYO)
                  MY_BAND_ELEM = ARRAY3(4)
*
                  CALL SOLV_SKYD (ELV, SKY, WEI, WORK, MAX_DIP, MAX_K,
     $            FORWARD_EFFICIENCY, T_CHOPPER, TATM,
     $            AMBIENT_TEMPERATURE+273.15, AMBIENT_PRESSURE,
     $            SITE_ALTITUDE/1000.,
     $            MADRID(P1), PARIS(P2), TOKYO(P3), MY_BAND_ELEM,
     $            CALIB, TAU, WATER, SIGMA, MODE, ITER, STATUS)
*
               ELSE
                  CALL MESSAGE (7,1,'SOLVE',
     $            'Error in the REAL version of the bandpass variables')
                  ERROR = .TRUE.
                  RETURN
               ENDIF
            ELSE
               CALL MESSAGE (6,1,'SOLVE',
     $         'BAND_FREQ, BAND_TRANS and BAND_ATMOS badly defined')
               ERROR = .TRUE.
               RETURN
            ENDIF
         ELSE
            CALL MESSAGE (5,1,'SOLVE',
     $      'BAND_FREQ, BAND_TRANS or BAND_ATMOS not available')
            ERROR = .TRUE.
            RETURN
         ENDIF
*
      ENDIF
*
c      print*,'Calibration factor,tau,water,sigma ',
c     &        calib,tau,water,sigma
*
* How is the fit? Sigma is in Kelvin.
*
      IF ((SIGMA .GT. 1.5) .AND. (SIGMA .LT. 3))
     $CALL MESSAGE (8,1,'SOLVE','FIT is not good')
      IF (SIGMA .GT. 3.0) CALL MESSAGE (8,1,'SOLVE','FIT is bad')
      IF (STATUS .EQ. 1) CALL MESSAGE (9,1,'SOLVE','No data to fit')
      IF (STATUS .EQ. 2) CALL MESSAGE (10,1,'SOLVE',
     $'Skydip counts must not have average zero.')
      IF (STATUS .EQ. 3) CALL MESSAGE (11,1,'SOLVE',
     $'Bandpass frequencies must lie between 1 and 1000 GHz')
      IF (STATUS .EQ. 4) CALL MESSAGE (11,1,'SOLVE',
     $'BAND_TRANS must be >= 0, with at least one element > 0')
      IF (STATUS .EQ. 11) CALL MESSAGE (12,1,'SOLVE',
     $'Fit did not converge, exit by maximum iterations')
      IF (STATUS .EQ. 10) CALL MESSAGE (13,1,'SOLVE',
     $'Fit was not terminated')
*
* Plot only if result is ok
*
      IF (STATUS .EQ. 0) THEN
*
         DO I = 1, NPOINT
            XM(I) = ELV(1)+((I-1.)/(NPOINT-1.))*(ELV(MAX_K)-ELV(1))
            YM(I) = (T_CHOPPER - FORWARD_EFFICIENCY * (TATM *
     $      (1. - EXP(-TAU*XM(I))) + 2.7*EXP(-TAU*XM(I)))
     $      - (AMBIENT_TEMPERATURE+273.15) * (1.0 - FORWARD_EFFICIENCY))
     $      /CALIB
         ENDDO
*     Compute baseline value for edges
*
         CALL GR4_CONNECT (NPOINT,XM,YM,0.,-1.)
         CALL GR_OUT
*
      ENDIF
      RETURN
      END
*
      FUNCTION AIR_MASS (ELEVATION)
*----------------------------------------------------------------------
* NIC	SOLVE support routine
*	Compute air masses
*
* Arguments :
*	ELEV		R*4	elevation in degrees		Input
*	AIR_MASS	R*4	number of airmasses		Output
*
*----------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'const.inc'
      REAL*4  ELEVATION, AIR_MASS
      REAL*8  HZ, H0, R0, R1
      REAL*8 ELV, EPS, GAM
c      PARAMETER (H0 = 5.5, R0 = 180.0)  ! atmospheric scale height
*						! earth radius (both in km)
      PARAMETER (H0 = 5.5, R0 = 6370.0)  ! atmospheric scale height
*						! earth radius (both in km)
      ELV = ELEVATION*PI/180
      R1  = R0+H0
      EPS = ASIN(R0/R1*COS(ELV))
      GAM = PI/2-ELV-EPS
      HZ  = R0*R0 + R1*R1 - 2*R0*R1*COS(GAM)
      HZ  = SQRT(MAX(HZ,H0*H0))          ! avoid less than one air mass
*
      AIR_MASS = HZ/H0
      RETURN
      END
*
      SUBROUTINE SKYDIP_NEXT (UT_OLD,UT_NEW,LST,TAU_OLD,TAU_NEW,
     $TAU,N_OLD,N_NEW)
*----------------------------------------------------------------------
* NIC	SOLVE support routine
*	Compute air masses
*
* Arguments :
*	TAU		R*4	atm.opacity		Input
*	LST		R*4	local sidereal time	Input
*       N_OLD           R*4     number of OLD skydips   Input
*       TAU_OLD	        R*4	list of opacities	Input
*	UT_OLD    	R*8	list of times 		Input
*       N_NEW           R*4     number of NEW skydips   Input
*       TAU_NEW	        R*4	list of opacities	Output
*	UT_NEW    	R*8	list of times 		Output
*
*----------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
      INTEGER  J, L, N_NEW, N_OLD
      REAL*4 TAU, TAU_OLD(N_OLD), TAU_NEW(N_NEW)
      REAL*4 LST(NRECORD), UT_DUMMY
      REAL*8 UT_OLD(N_OLD), UT_NEW(N_NEW), AM
*
      UT_DUMMY = LST(1)
      CALL STUT (UT_DUMMY,1,
     $IDAY,IMONTH,IYEAR,SITE_LONGITUDE/DEG_TO_RAD)
      AM = DBLE(DATE_OF_OBSERVATION)+UT_DUMMY/86400D0
      IF (AM.GT.UT_OLD(N_OLD)) THEN
         CALL R8TOR8 (UT_OLD,UT_NEW,N_OLD)
         CALL R4TOR4 (TAU_OLD,TAU_NEW,N_OLD)
         UT_NEW(N_NEW) = AM
         TAU_NEW(N_NEW) = TAU
      ELSE
         J = 1
         L = 1
         DO WHILE (L.LE.N_OLD)
            IF (AM.LT.UT_OLD(L)) THEN
               UT_NEW(J) = AM
               TAU_NEW(J) = TAU
               AM = 50000000             ! nic shouldnt survive 5000 a.d.
            ELSE
               UT_NEW(J) = UT_OLD(L)
               TAU_NEW(J) = TAU_OLD(L)
               L = L+1
            ENDIF
            J = J+1
         ENDDO
      ENDIF
      RETURN
      END
*
      SUBROUTINE SKYDIP_FIRST (UT_NEW,LST,TAU_NEW,TAU,AM)
*----------------------------------------------------------------------
* NIC	SOLVE support routine
*	Compute air masses
*
* Arguments :
*	TAU		R*4	atm.opacity		Input
*	LST		R*4	local sidereal time	Input
*       TAU_NEW	        R*4	list of opacities	Output
*	UT_NEW    	R*8	list of times 		Output
*
*----------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
* Local
      REAL*4 TAU, TAU_NEW(1), LST(NRECORD), UT_DUMMY
      REAL*8 UT_NEW(1),AM
*
      UT_DUMMY = LST(1)
      CALL STUT (UT_DUMMY,1,
     $IDAY,IMONTH,IYEAR,SITE_LONGITUDE/DEG_TO_RAD)
      UT_NEW(1) = DBLE(DATE_OF_OBSERVATION)+UT_DUMMY/86400D0
      AM = UT_NEW(1)
      TAU_NEW(1) = TAU
      RETURN
      END
*
      FUNCTION SKYDIP_NEW (UT_NEW,LST,N,AM)
      IMPLICIT NONE
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
* Local
      INTEGER J, N
      LOGICAL SKYDIP_NEW
      REAL*4 LST(NRECORD), UT_DUMMY
      REAL*8 UT_NEW(N),AM
*
      SKYDIP_NEW = .TRUE.
      UT_DUMMY = LST(1)
      CALL STUT (UT_DUMMY,1,
     $IDAY,IMONTH,IYEAR,SITE_LONGITUDE/DEG_TO_RAD)
      AM = DBLE(DATE_OF_OBSERVATION)+UT_DUMMY/86400D0
      DO J=1,N
         IF (ABS(AM-UT_NEW(J)).LT.1E-5) SKYDIP_NEW = .FALSE.
      ENDDO
      RETURN
      END
*
      FUNCTION AZI (ADRESS)
      REAL AZI, ADRESS(1)
*
      AZI = NINT(ADRESS(1))
      RETURN
      END
*---------------------------------------------------------------------------
*
      SUBROUTINE SOLV_SKYD (X, Y, WW, WORK, NDIM, NVAL,
     $FEFF, TCHOP, TATM, TAMB, PAMB, ALTI,
     $BAND_FREQ, BAND_TRANS, BAND_ATMOS, BAND_ELEM,
     $CT, TAU, WATER, SIGMA, MODE, ITER, STATUS)
*
* Solve a skydip in a conservative way
* take forward efficiency and chopper temperature for granted
* Author: Michael Bremer, IRAM 21-Nov-1995
*
* Input: x(nval), y(nval), nval: vectors of airmass and signal, and
*                                their number of elements
*                   feff, tchop: forward efficiency, chopper temperature
*              tamb, pamb, alti: ambiant temperature and pressure,
*                                obseratory altitude
*      band_freq, band_trans, band_elem: bandpass points, transmissions and
*                                number of elements
*                          mode: 1 for GOLDEN, 2 for BRENT
* Output:       ct, tau, status: calibration proportionality constant, tau,
*                                coded return message that is 0 for sucess
*--------------------------------------------------------------------------
*
      IMPLICIT NONE
      INTEGER STATUS, NVAL, BAND_ELEM, I, MODE, ITMAX, ITER, NDIM
      REAL    X(NDIM), Y(NDIM), WW(NDIM), WORK(NDIM),
     $TCHOP, TATM, FEFF, TAMB, PAMB, ALTI, CT, TAU,
     $BAND_FREQ(BAND_ELEM), BAND_TRANS(BAND_ELEM), BAND_ATMOS(BAND_ELEM)
      REAL WATER_MAX, WATER, FREQ_MIN, FREQ_MAX, TRANS_MAX, TRANS_MIN,
     $SUM, CENT_Y, EPSI, SIGMA,
     $X0,X1,X2,X3, F0,F1,F2,F3, AX,BX,CX, R,CGOLD,XMIN, TOL, GOLDEN
      REAL ZEPS,TOL1,TOL2, P,Q,R1, U,V,W, A,B,D,E, FX,FU,FV,FW, XM,ETEMP
      PARAMETER (R=.61803399,CGOLD=.38196602)
      PARAMETER (ITMAX=20,ZEPS=1.0E-10)
*
      TOL  = 0.001
      EPSI = 0.001
      STATUS = 10
*
* Set calibration factors on something stupid but not lethal
      TAU   = 2.0
      CT    = 1.0
*
*
* !!! 1. error condition: Check if there is anything to work with:
*
      IF (NVAL .LE. 0) THEN
         STATUS = 1
         RETURN
      ENDIF
*
* Estimate center for given skydip points:
      SUM = 0.0
      DO I = 1, NVAL
         SUM = SUM + Y(I)
      ENDDO
      CENT_Y = SUM / NVAL
*
* !!! 2. error condition: This center must not be zero:
*
      IF (ABS(CENT_Y) .LT. EPSI) THEN
         STATUS = 2
         RETURN
      ENDIF
*
* Check if the bandpass function is useful
*
      FREQ_MIN = BAND_FREQ(1)
      FREQ_MAX = BAND_FREQ(1)
      TRANS_MIN = BAND_TRANS(1)
      TRANS_MAX = BAND_TRANS(1)
      DO I = 1, BAND_ELEM
         FREQ_MIN = MIN(FREQ_MIN, BAND_FREQ(I))
         FREQ_MAX = MIN(FREQ_MAX, BAND_FREQ(I))
         TRANS_MIN = MIN(TRANS_MIN, BAND_TRANS(I))
         TRANS_MAX = MAX(TRANS_MAX, BAND_TRANS(I))
      ENDDO
      IF ((FREQ_MIN .LT. 1) .OR. (FREQ_MAX .GT. 1000)) THEN
         STATUS = 3
         RETURN
      ENDIF
      IF ((TRANS_MIN .LT. 0) .OR. (TRANS_MAX .LT. 1E-6)) THEN
         STATUS = 4
         RETURN
      ENDIF
*
* Set up atmosphere:
*
      CALL B_ATM_ATMOSP (TAMB, PAMB, ALTI, WATER_MAX)
*
      AX = EPSI
      BX = 0.2
      CX = 1.0
*
* MINIMIZE:
*
*  mode=1: use GOLDEN from Numres <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      IF (MODE .EQ. 1) THEN
*
         X0=AX
         X3=CX
         IF(ABS(CX-BX).GT.ABS(BX-AX))THEN
            X1=BX
            X2=BX+CGOLD*(CX-BX)
            X2 = MIN(MAX(X2, EPSI), 5.0)
         ELSE
            X2=BX
            X1=BX-CGOLD*(BX-AX)
            X1 = MIN(MAX(X1, EPSI), 5.0)
         ENDIF
         CALL TRYVAL (X, Y, WW, WORK, NDIM, NVAL, CENT_Y, X1,
     $   FEFF, TCHOP, TATM, TAMB, WATER_MAX, WATER,
     $   BAND_FREQ, BAND_TRANS, BAND_ATMOS, BAND_ELEM,
     $   CT, F1, STATUS)
         CALL TRYVAL (X, Y, WW, WORK, NDIM, NVAL, CENT_Y, X2,
     $   FEFF, TCHOP, TATM, TAMB, WATER_MAX, WATER,
     $   BAND_FREQ, BAND_TRANS, BAND_ATMOS, BAND_ELEM,
     $   CT, F2, STATUS)
         ITER = 0
1        CONTINUE
         ITER = ITER + 1
         IF(ABS(X3-X0).GT.TOL*(ABS(X1)+ABS(X2)))THEN
            IF(F2.LT.F1)THEN
               X0=X1
               X1=X2
               X2=R*X1+CGOLD*X3
               X2 = MIN(MAX(X2, EPSI), 5.0)
               F0=F1
               F1=F2
               CALL TRYVAL (X, Y, WW, WORK, NDIM, NVAL, CENT_Y, X2,
     $         FEFF, TCHOP, TATM, TAMB, WATER_MAX, WATER,
     $         BAND_FREQ, BAND_TRANS, BAND_ATMOS, BAND_ELEM,
     $         CT, F2, STATUS)
            ELSE
               X3=X2
               X2=X1
               X1=R*X2+CGOLD*X0
               X1 = MIN(MAX(X1, EPSI), 5.0)
               F3=F2
               F2=F1
               CALL TRYVAL (X, Y, WW, WORK, NDIM, NVAL, CENT_Y, X1,
     $         FEFF, TCHOP, TATM, TAMB, WATER_MAX, WATER,
     $         BAND_FREQ, BAND_TRANS, BAND_ATMOS, BAND_ELEM,
     $         CT, F1, STATUS)
            ENDIF
            IF (ITER .LE. ITMAX) GOTO 1
         ENDIF
         IF(F1.LT.F2)THEN
            GOLDEN=F1
            XMIN=X1
         ELSE
            GOLDEN=F2
            XMIN=X2
         ENDIF
*
*  finished  GOLDEN from Numres <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
*
      ELSE
*
*  mode<>1: use BRENT from Numres : +++++++++++++++++++++++++++++++
         A=MIN(AX,CX)
         B=MAX(AX,CX)
         V=BX
         W=V
         X1=V
         X1 = MIN(MAX(X1, EPSI), 5.0)
         E=0.
         CALL TRYVAL (X, Y, WW, WORK, NDIM, NVAL, CENT_Y, X1, FEFF,
     $   TCHOP, TATM, TAMB, WATER_MAX, WATER,
     $   BAND_FREQ, BAND_TRANS, BAND_ATMOS, BAND_ELEM,
     $   CT, FX, STATUS)
         FV=FX
         FW=FX
         DO 11 ITER=1,ITMAX
            XM=0.5*(A+B)
            TOL1=TOL*ABS(X1)+ZEPS
            TOL2=2.*TOL1
            IF(ABS(X1-XM).LE.(TOL2-.5*(B-A))) GOTO 103
            IF(ABS(E).GT.TOL1) THEN
               R1=(X1-W)*(FX-FV)
               Q=(X1-V)*(FX-FW)
               P=(X1-V)*Q-(X1-W)*R1
               Q=2.*(Q-R1)
               IF(Q.GT.0.) P=-P
               Q=ABS(Q)
               ETEMP=E
               E=D
               IF(ABS(P).GE.ABS(.5*Q*ETEMP).OR.P.LE.Q*(A-X1).OR.
     $         P.GE.Q*(B-X1)) GOTO 101
               D=P/Q
               U=X1+D
               IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=SIGN(TOL1,XM-X1)
               GOTO 102
            ENDIF
101         IF(X1.GE.XM) THEN
               E=A-X1
            ELSE
               E=B-X1
            ENDIF
            D=CGOLD*E
102         IF(ABS(D).GE.TOL1) THEN
               U=X1+D
            ELSE
               U=X1+SIGN(TOL1,D)
            ENDIF
            U = MIN(MAX(U, EPSI), 5.0)
            CALL TRYVAL (X, Y, WW, WORK, NDIM, NVAL, CENT_Y, U, FEFF,
     $      TCHOP, TATM, TAMB, WATER_MAX, WATER,
     $      BAND_FREQ, BAND_TRANS, BAND_ATMOS, BAND_ELEM,
     $      CT, FU, STATUS)
            IF(FU.LE.FX) THEN
               IF(U.GE.X1) THEN
                  A=X1
               ELSE
                  B=X1
               ENDIF
               V=W
               FV=FW
               W=X1
               FW=FX
               X1=U
               FX=FU
            ELSE
               IF(U.LT.X1) THEN
                  A=U
               ELSE
                  B=U
               ENDIF
               IF(FU.LE.FW .OR. W.EQ.X1) THEN
                  V=W
                  FV=FW
                  W=U
                  FW=FU
               ELSE
                  IF(FU.LE.FV .OR. V.EQ.X1 .OR. V.EQ.W) THEN
                     V=U
                     FV=FU
                  ENDIF
               ENDIF
            ENDIF
11       CONTINUE
103      CONTINUE
         XMIN=X1
*  finished  BRENT from Numres : ++++++++++++++++++++++++++++++++++
*
      ENDIF
*
*
*
      CALL TRYVAL (X, Y, WW, WORK, NDIM, NVAL, CENT_Y, XMIN, FEFF,
     $TCHOP, TATM, TAMB, WATER_MAX, WATER,
     $BAND_FREQ, BAND_TRANS, BAND_ATMOS, BAND_ELEM,
     $CT, SIGMA, STATUS)
*
      TAU  = XMIN
      SIGMA = SQRT(SIGMA)
      STATUS = 0
      IF (ITER .EQ. ITMAX) STATUS = 11
      END
*
* *********************************************************************
*
      SUBROUTINE TRYVAL (X, Y, W, Y_TRY, NDIM, NVAL, CENT_Y, TAU, FEFF,
     $TCHOP, TATM, TAMB, WATER_MAX, WATER,
     $BAND_FREQ, BAND_TRANS, BAND_ATMOS, BAND_ELEM,
     $CT, SIGMA, STATUS)
      IMPLICIT NONE
      INTEGER STATUS, NDIM, NVAL, BAND_ELEM, IER, I
      REAL    X(NDIM), Y(NDIM), W (NDIM), Y_TRY(NDIM),
     $TCHOP, FEFF, CT, TAU,
     $BAND_FREQ(BAND_ELEM), BAND_TRANS(BAND_ELEM), BAND_ATMOS(BAND_ELEM)
      REAL WATER, WATER_MAX, TEMI, TATM,
     $TAMB, TEMI_B,
     $SUM1, WSUM, CENT_Y, SIGMA, DELTA
*
      CALL B_TAUT_TO_WP (WATER_MAX, BAND_FREQ, BAND_TRANS, BAND_ATMOS,
     $BAND_ELEM,TEMI,TATM, TAU, WATER, IER)
*
* Calculate emission for all given airmasses. Accelerate by using
* exponential formula
*
*
      WSUM = 0.0
      SUM1 = 0.0
*
      DO I = 1, NVAL
         TEMI_B = TATM*(1.0 - EXP(-TAU*X(I)))
     $   + 2.7 *EXP(-TAU*X(I))
         Y_TRY(I) = TCHOP - FEFF*TEMI_B - TAMB*(1.0 - FEFF)
         SUM1 = SUM1 + W(I) * Y_TRY(I)
         WSUM = WSUM + W(I)
      ENDDO
      SUM1 = SUM1 / WSUM
*
* NOTE: it is possible to calculate sigma in one go by separated sums,
*       and to calibrate and combine them later:
*
*       sum (yi-yti)^2 = sum(yi^2)-2*sum(yi*yti)+sum(yti^2)
*
*       But in this case, the terms will differ
*       considerably in order of magnitude, causing round-off errors even
*       in double precision if the conversion factor counts->temperature is
*       several orders of magnitude, which is usually the case.
*       Therefore, it is safer to work with an additional vector buffer.
*
* Divide the centers of the curves by each other. This gives the
* proportionality factor
*
      CT  = SUM1 / CENT_Y
*
* Estimate the sigma
*
      SIGMA = 0.0
      DO I = 1, NVAL
         DELTA = (Y_TRY(I) - Y(I)*CT)
         SIGMA = SIGMA + W(I) * DELTA*DELTA
      ENDDO
*
      SIGMA = SIGMA / WSUM
*
      END
*
***************************************************************************
*
      SUBROUTINE B_ATM_ATMOSP(T0,P0,H0,WATER_MAX)
C----------------------------------------------------------------------
C Compute an atmospheric model, interpolated between standard atmospheres
C of winter and summer (subroutines ase45 and asj 45), to fit with temperature
C t0 (k) and pressure p0 (mbar) at altitude h0 (km).
C 15 layers are used.
C The transmission of the model atmosphere can then be computed by calling
C entry point transm.
C
C input 	t0	R	temperature	(K)
C		po	R	pressure	(mbar)
C		h0	R	altitude	(km)
C----------------------------------------------------------------------
      INTEGER MP
      PARAMETER (MP=15)
      REAL*4  T(80), H(80), P(80), R(80), RR(80), RS(80),
     $PE, TE, DE, H0, T0, P0, PJ, TJ, DJ, APE, ATE, APJ, ATJ, P1,
     $T1, D1, HEIGHT, R1, WATER, AIRMASS, TAUW, TAUOX, TAUT, TEMI,
     $FREQ, TAG, TOX, TAGU, TOXI, TATM, CT_1, WATER_MAX, T_CELSIUS,A,B
      INTEGER NP, J, IER
*
      SAVE T,H,P,R
*
* average summer and winter model atmospheres according to
* given values of temperature and pression
      CALL ASE45(PE,TE,DE,H0)
      CALL ASJ45(PJ,TJ,DJ,H0)
      APE = (P0-PJ)/(PE-PJ)
      APJ = (P0-PE)/(PJ-PE)
      ATE = (T0-TJ)/(TE-TJ)
      ATJ = (T0-TE)/(TJ-TE)
*
* set layers
      DO J=1,6
         H(J) = .5E5
      ENDDO
      DO J=7,12
         H(J) = 2.E5
      ENDDO
      DO J=13,15
         H(J) = 15.E5
      ENDDO
      NP = 15
*
* Set t,p, and r (H2O for 1mm precipitable content) profiles
      HEIGHT = H0
      P1 = P0
      T1 = T0
      R1 = .5
      CT_1 = 0.0
      WATER_MAX = 0.0
      DO J = 1, NP
         HEIGHT = HEIGHT + H(J)/100000.  ! in km.
         P(J) = P1
         T(J) = T1
         R(J) = R1
         CALL ASE45(PE,TE,DE,HEIGHT)
         CALL ASJ45(PJ,TJ,DJ,HEIGHT)
         P1 = APE*PE+APJ*PJ
         T1 = ATE*TE+ATJ*TJ
         D1 = DE*(1+(P1-PE)/PE-(T1-TE)/TE)
         R1 = .5*EXP(-.5*(HEIGHT-H0))
         IF(HEIGHT.GT.15.) R1 = R1 + D1*2E-6
         P(J) = (P(J) + P1)/2.
         T(J) = (T(J) + T1)/2.
         R(J) = (R(J) + R1)/2.
         CT_1 = CT_1 + R(J)*H(J)*1E-5
         T_CELSIUS = T(J) - 273.155
*
* Fleagle, Atmospheric Physics, Academic Press 1980:
* 50% of clouds at T=-8 degree Celsius consist of supercooled water drops
*
         IF (T_CELSIUS .GE. -8.0) THEN
*
* coefficients for saturated vapor pressure over water (Kyle,
* Atmospheric Transmission, Emission and Scattering, Pergamon Press 1991:
            A = 7.5
            B = 237.3
         ELSE
*
* coefficients for saturated vapor pressure over ice:
            A = 9.5
            B = 265.5
         ENDIF
c         rs(j) = 0.0
         RS(J) = 1324/T(J) * 10 ** (A*T_CELSIUS/(B + T_CELSIUS))
         IF (HEIGHT .LE. 30)
     $   WATER_MAX = WATER_MAX + RS(J)*H(J)*1E-5
      ENDDO
c      print*,'Atmospheric water vapor capacity: '
c      print*,'Zenith saturation up to 30 km would be ',
c     &       water_max,' mm vapor'
*
* Ensure that the 1mm water profile is 1mm indeed:
*
      DO J = 1, NP
         R(J) = R(J) / CT_1
      ENDDO
*
      RETURN
*
*
      ENTRY ATM_TRANSM(WATER,AIRMASS,FREQ,TEMI,TATM,TAUOX,TAUW,TAUT
     $,IER)
C----------------------------------------------------------------------
C Compute atmospheric emission and  absorption.
C
C Input:
C 	water 	R	H2O precipitable content(mm)
C	airmass R	Number of air masses
C	freq 	R	Frequency		(GHz)
C
C Output:
C	temi	R	atmosph emission	(K)
C	tatm   	R	mean temperature	(K)
C	tauox	R	Oxygen optical depth  AT ZENITH	(nepers)
C       tauw   	R	Water  optical depth  AT ZENITH	(nepers)
C       taut   	R	Total  optical depth  AT ZENITH	(nepers)
C	IER	I	Error code
C----------------------------------------------------------------------
      NP = 15
      DO J = 1, NP
         RR(J) = R(J) * WATER
      ENDDO
      IER = 0
      CALL KVATM(NP,P,T,RR,H,TAUW,TAUOX,FREQ,TEMI,TATM,TAG,TAGU,TOX,
     $TOXI,0,0,TAUT,AIRMASS,IER)
      TAUOX = TAUOX / AIRMASS            ! RL 14 MAR 86
      TAUW = TAUW / AIRMASS              !
      TAUT = TAUT / AIRMASS              !
      RETURN
      END
*
*-----------------------------------------------------------------------
*
*
***********************************************************************
*
      SUBROUTINE B_TAUT_TO_WP ( WATER_MAX, BAND_F, BAND_T, BAND_A,
     $BAND_E, TEMI, TATM, TAUT, WATER, IER)
*
* Broad band version, using an optional bandpass for the bolometer
* Basic idea: Use the form of the sky dip curve to calibrate the bolometer
* This routine assumes the total atmospheric tau is known from a previous
* gauss fit, and derives the corresponding
*
* Assumes zenith elevation, i.e. one airmass
* Calculates the amount of total precipitable water from ...
* Input:  real*4 water_max amount of precipitable water in mm
*                         that would cause total saturation.
*                 taut    total atmospheric opacity (nepers)
*
* Output:
*         real*4  temi    Sky emissivity
*                 tatm    effective atmospheric temperature
*                 water   total amount of precipitable water (mm)
*         integer ier     error code
*
* !!! BEFORE this routine is executed, set up the atmosphere !!!
* !!! with a call to B_ATM_ATMOSP.                             !!!
*
* Author: 03-Nov-1995  Michael Bremer, IRAM
*----------------------------------------------
*
      IMPLICIT NONE
      INTEGER IER,  BAND_E
      REAL*4 AIRMASS, WATER_MAX,
     $TEMI, TATM, TAUT, WATER,  TAUOX, TAUW,
     $W_1, W_2, W_3, TAUT_M1, TAUT_M2,
     $TAUT_M3, EPSI,
     $TAUW_M1, TAUW_M2, A,B,C,
     $BAND_F(BAND_E), BAND_T(BAND_E), BAND_A(BAND_E)
*
      IER   = 0
      EPSI  = 0.001
*
* Error code ier=4 will be set if the water would be too low,
*            ier=5 if it is too high, i.e. higher as saturation amount
*            ier=6 if the precision does not match epsi
* (code 1-3 are already used in KVATM)
*
      AIRMASS = 1.0
*
      W_1 = 0.001
*
      CALL B_ATM_TRANSM (W_1, AIRMASS, BAND_F,BAND_T,BAND_A,BAND_E,
     $TEMI, TATM, TAUOX, TAUW_M1, TAUT_M1, IER)
*
*
      IF (TAUT_M1 .GT. TAUT) THEN
         PRINT*, TAUT_M1, TAUT
         PRINT*,'Error: given taut would cause water < 0'
         WATER = 0.0
         IER   = 4
         RETURN
      ENDIF
*
      W_2 = 7
*
      CALL B_ATM_TRANSM(W_2, AIRMASS, BAND_F,BAND_T,BAND_A,BAND_E,
     $TEMI, TATM, TAUOX, TAUW_M2, TAUT_M2, IER)
*
* water and taut depend linearly on each other (in good approximation).
* construct a linear relationship and derive the precise amount of water
* for the current atmospheric conditions and the desired frequency
*
      W_3 = (TAUT * (W_2 - W_1) + (TAUT_M2*W_1-TAUT_M1*W_2)) /
     $(TAUT_M2 - TAUT_M1)
*
* Precision not high enough, digresses by some percent.
* Recalculate the exact model for the derived water
*
      CALL B_ATM_TRANSM(W_3, AIRMASS, BAND_F,BAND_T,BAND_A,BAND_E,
     $TEMI, TATM, TAUOX, TAUW, TAUT_M3, IER)
*
*
* make a parabolic interpolation with the third point.
*
      CALL POLI2(TAUT_M1,TAUT_M2,TAUT_M3,W_1,W_2,W_3,A,B,C)
      WATER = A + TAUT * (B + C*TAUT)
      CALL B_ATM_TRANSM(WATER, AIRMASS, BAND_F,BAND_T,BAND_A,BAND_E,
     $TEMI, TATM, TAUOX, TAUW, TAUT_M3, IER)
*
      IF (WATER .GT. WATER_MAX) IER = 5
*
      IF (ABS(TAUT - TAUT_M3) .GT. EPSI) IER = 6
*
      RETURN
*
      END
*
***********************************************************************
*
      SUBROUTINE B_ATM_TRANSM (WATER, AIRMASS,
     $BAND_FREQ,BAND_TRANS,BAND_ATMTR,BAND_ELEM,
     $TEMI_B, TATM_B, TAUOX_B, TAUW_B, TAUT_B, IER)
*
* calculates the atmospheric parameters for a broad band instrument
* Version 2.0 Michael Bremer, IRAM 02-Nov-1995
* Version 3.0 MB 20-Mar-1997 include a real BAND_ATMTR(BAND_ELEM) variable
*                for atmospheric transmission. This variable will be
*                normalized to 1 at band_freq(1). Value is exp(-taut)
*
* Uses trapezoid integration to convolve the atmospheric parameters with
* the bolometer bandpass and to integrate the bandpass itself. Returns the
* convolved parameters normalized with the integral of the bandpass.
*
* Can use arbitrarily sized spectral bins, but the frequency supports must
* be sorted in ascending order
*
* temi, tatm, tauox, tauw, taut, path, di    will be averaged over the band
* nhb determines the number of spectral bins. band_elem atmospherical models
* will be calculated for the specified spectral range.
*
* NOTE: band_elem must at least have one element
*       otherwise an error condition (ier = 100) is set
*
*
      IMPLICIT NONE
      INTEGER IER, I, K
      REAL*4 TEMI_B, TATM_B, TAUOX_B, TAUW_B, TAUT_B
      REAL*4  WATER, AIRMASS, F_M, TEMI,TATM,TAUOX,TAUW,TAUT
      REAL*4 YA(3), YB(3), ZA1, ZA2, ZB1, ZB2, A_INTEG(3),
     $NORM_BIN1, NORM_BIN2, DELT
*
      INTEGER BAND_ELEM
      REAL*4    BAND_FREQ(BAND_ELEM), BAND_TRANS(BAND_ELEM),
     $BAND_ATMTR(BAND_ELEM), REF_TRANS
*
* Check if valid bandpass is given
      IF (BAND_ELEM .LT. 1) THEN
         IER = 100
         RETURN
      ENDIF
*
* intitialize integration array
      DO I = 1, 3
         A_INTEG(I) = 0.0
      ENDDO
      NORM_BIN1 = 0.0
      NORM_BIN2 = 0.0
*
* Memo: who is who
* a_integ(1) =  temi_b
* a_integ(2) =  tatm_b
* a_integ(3) =  tauox_b
*
* arbitrary frequency spacing: do convolution, trapezoidal integration
* Integral over an interval is 0.5*(y(freq_a)+y(freq_b))*(freq_b-freq_a)
* calculate first border, index a
*
* First point
      F_M = BAND_FREQ(1)
      CALL ATM_TRANSM(WATER, AIRMASS, F_M,
     $TEMI, TATM, TAUOX, TAUW, TAUT, IER)
      YA(1) = TEMI  * BAND_TRANS(1)
      YA(2) = TATM  * BAND_TRANS(1)
      YA(3) = TAUOX * TATM * BAND_TRANS(1)
      ZA1    = BAND_TRANS(1)
      ZA2    = TATM * BAND_TRANS(1)
*
      BAND_ATMTR(1) = 1.0
      REF_TRANS     = EXP(-TAUT * AIRMASS)
*
* Only one bandpass element given: routine collapses into narrowband case
*
      IF (BAND_ELEM .EQ. 1) THEN
         TEMI_B  = TEMI
         TATM_B  = TATM
         TAUOX_B = TAUOX
         TAUW_B  = TAUW
         TAUT_B  = TAUT
         RETURN
      ENDIF
*
* More than one frequency for bandpass given: calculate integration
*
      DO I = 2, BAND_ELEM
         F_M = BAND_FREQ(I)
         CALL ATM_TRANSM(WATER, AIRMASS, F_M,
     $   TEMI, TATM, TAUOX, TAUW, TAUT, IER)
*
* calculate second border (index b)
         YB(1) = TEMI  * BAND_TRANS(I)
         YB(2) = TATM  * BAND_TRANS(I)
         YB(3) = TAUOX * TATM * BAND_TRANS(I)
         ZB1    = BAND_TRANS(I)
         ZB2    = TATM * BAND_TRANS(I)
*
* integrate over atmosph. Parameters convolved with bandpass, and  bandpass
*
         DELT  = 0.5 * (BAND_FREQ(I) - BAND_FREQ(I-1))
         DO K = 1, 3
            A_INTEG(K) = A_INTEG(K) + (YB(K)+YA(K)) * DELT
         ENDDO
         NORM_BIN1 = NORM_BIN1 + (ZA1+ZB1) * DELT
         NORM_BIN2 = NORM_BIN2 + (ZA2+ZB2) * DELT
*
* change ya,za into yb,zb for the next frequency step
*
         DO K = 1, 3
            YA(K) = YB(K)
         ENDDO
         ZA1 = ZB1
         ZA2 = ZB2
         BAND_ATMTR(I) = EXP(-TAUT*AIRMASS) / REF_TRANS
      ENDDO
      DO I=1,BAND_ELEM
         BAND_ATMTR(I) = BAND_ATMTR(I)*REF_TRANS
      ENDDO
*
* normalize values with transmission function integral, and for oxygen tau
* with transmission function times tatm
*
      A_INTEG(1) = A_INTEG(1) / NORM_BIN1
      A_INTEG(2) = A_INTEG(2) / NORM_BIN1
      A_INTEG(3) = A_INTEG(3) / NORM_BIN2
      TEMI_B  = A_INTEG(1)
      TATM_B  = A_INTEG(2)
      TAUOX_B = A_INTEG(3)
*
* Calculate broad band total tau from broad band tsky and broad band tatm
      TAUT_B = - LOG( 1.0 - TEMI_B/TATM_B) / AIRMASS
*
* calculate wet tau from total tau - taut oxygen
* make sure that the wet tau becomes not negative (we don't want a masing
* atmosphere, don't we?)
*
      TAUW_B = MAX(TAUT_B - TAUOX_B, 0.0)
*
      END
*
**************************************************************************
*
      SUBROUTINE KVATM(NP,P,T,RHO,H,AGU,OXI,V,TEMI,TATM,TAG,TAGU,TOX,
     $TOXI,ILAG,ILOX,KVAT,AMA,IER)
cc     $TOXI,ILAG,ILOX,KVAT,AMA,path,di,IER)
C----------------------------------------------------------------------
C	opacidad de la atmosfera a la frecuencia 'v' debida al
C	vapor de agua y al oxigeno.
C
C	np .... numero de capas
C	h  .... espesor de las capas     (cm)
C	p  .... presion (milibares)
C	t  .... temperatura (k)
C	rho ... cantidad de vapor de agua (gr/m**3)
C	temi .. emisividad total de la atmosfera     (k)
C	kvat .. opacidad total (nepers)
C	tatm .. temperatura media de la atmosfera   (k)
C	agu ... opacidad debida al vapor de agua (nepers)
C	oxi ... opacidad debida al oxigeno       (   "  )
C	tag ... emisividad debida al vapor de agua (atmosfera sin oxigeno)
C	tox ... idem para el oxigeno (sin vapor de agua)
C	tagu .. temperatura media de la capa de vapor de agua (k)
C	toxi .. temperatura media de la capa de oxigeno (k)
C
C J.Cernicharo
C----------------------------------------------------------------------
      INTEGER NP,ILAG,ILOX,IER
      REAL*4 AGU,OXI,V,TEMI,TATM,TAG,TAGU,TOX,TOXI,KVAT,AMA
      REAL*4 H(*),P(*),T(*),RHO(*)
c      REAL*4 H(np),P(np),T(np),RHO(np)
*
      INTEGER J
      REAL*4 R,PR,TEM,DH,OX,AG,KV
      REAL*4 KH2O, KO2
*
*
      TEMI=0.
      KV=0.
      TAG=0.
      TOX=0.
      AGU=0.
      OXI=0.
      DO 1 J=1,NP
         R=RHO(J)
         PR=P(J)
         TEM=T(J)
         DH=H(J)
*
* partial pressure of water vapor. Rspec = Rgas/M_H2O = 8314/18.02 = 461.4
* Conversion from pascal->mbar 1e-2, g->kg 1e-3:
c         pr_ag = 4.614E-03  * tem * r
*
         AG=KH2O(R,TEM,PR,V,ILAG)*DH*AMA
         OX=KO2(TEM,PR,V,ILOX)*DH*AMA
         TAG=TAG+TEM*EXP(-AGU)*(1.-EXP(-AG))
         AGU=AGU+AG
         TOX=TOX+TEM*EXP(-OXI)*(1.-EXP(-OX))
         OXI=OXI+OX
         TEMI=TEMI+TEM*EXP(-KV)*(1.-EXP(-AG-OX))
1     KV=AGU+OXI
      KVAT=KV
      IF ( KV.LE.1.E-10 ) THEN
         IER = 1
      ELSEIF ( OXI.LE.1.E-20) THEN
         IER = 2
      ELSEIF ( AGU.LE.1.E-20 ) THEN
         IER = 3
      ELSE
         TATM=TEMI/(1.-EXP(-KV))
         TAGU=TAG/(1.-EXP(-AGU))
         TOXI=TOX/(1.-EXP(-OXI))
         IER = 0
      ENDIF
      END
* Last processed by NICE on
* Customized for :  EEEI        HPUX        UNIX
      SUBROUTINE POLI2(X1,X2,X3,Y1,Y2,Y3,A,B,C)
C----------------------------------------------------------------------
C	ESTA SUBRUTINA CALCULA LOS COEFICIENTES A,B,C DEL POLINOMIO DE
C	SEGUNDO GRADO A+BX+CX**2, QUE PASA POR LOS PUNTOS (X1,Y1),
C	(X2,Y2),(X3,Y3)
C J.Cernicharo
C----------------------------------------------------------------------
      REAL*4 X1,X2,X3,Y1,Y2,Y3,A,B,C
*
      C=(Y3-Y2)*(X2-X1)-(Y2-Y1)*(X3-X2)
      B=(X2-X1)*(X3*X3-X2*X2)-(X2*X2-X1*X1)*(X3-X2)
      C=C/B
      B=(Y2-Y1)-C*(X2*X2-X1*X1)
      B=B/(X2-X1)
      A=Y1-C*X1*X1-B*X1
      END
*******************************************************************
*
* Last processed by NICE on
* Customized for :  EEEI        HPUX        UNIX
      SUBROUTINE ASE45(P,T,D,HA)
C J CERNICHARO
C
C	ATMOSFERA U.S. 1962 MES DE ENERO 45 GRADOS DE LATITUD NORTE
C	HA ES LA ALTURA A LA QUE SE QUIEREN CALCULAR LA PRESION,
C	(MILIBARES),LA TEMPERATURA (K) Y LA DENSIDAD (GR/M**3)
C	HA DEBE ESTAR EN KM
C
C	P === PRESION
C	T === TEMPERATURA
C	D === DENSIDAD
C
      REAL*4 P,T,D,HA
      REAL*4 PR(91),TE(91),DEN(91)
      INTEGER I0,I1,I2,J
      REAL*4 X1,X2,X3,Y1,Y2,Y3,A,B,C
      DATA TE/
     $272.150,268.650,265.154,261.654,255.668,249.674,243.686,237.698,
     $231.710,225.728,219.746,219.159,218.661,218.163,217.665,217.167,
     $216.670,216.172,215.675,215.178,215.150,215.150,215.150,215.150,
     $215.150,215.150,215.150,215.150,215.852,216.644,217.437,218.230,
     $219.022,221.723,224.789,227.855,230.921,233.987,237.049,240.112,
     $243.175,246.235,249.294,252.354,255.414,258.470,261.527,264.580,
     $265.650,265.650,265.650,265.650,265.650,264.526,262.560,260.594,
     $258.628,256.664,254.698,252.736,250.772,248.810,246.848,244.886,
     $242.926,240.932,238.874,236.818,234.762,232.706,230.650,228.596,
     $226.543,224.491,222.439,220.387,218.336,216.286,214.236,212.187,
     $210.150,210.150,210.150,210.150,210.150,210.150,210.150,210.150,
     $210.150,210.150,210.150/
      DATA PR/
     $.101800E04,.897340E03,.789746E03,.693761E03,.608132E03,.531322E03,
     $.462749E03,.401636E03,.347333E03,.299257E03,.256837E03,.219907E03,
     $.188249E03,.161092E03,.137803E03,.117839E03,.100747E03,.861032E02,
     $.735727E02,.628431E02,.536674E02,.458314E02,.391458E02,.334355E02,
     $.285581E02,.243961E02,.208407E02,.178034E02,.152146E02,.130097E02,
     $.111307E02,.952848E01,.816277E01,.700062E01,.601741E01,.518290E01,
     $.447305E01,.386792E01,.335147E01,.290933E01,.253005E01,.220437E01,
     $.192389E01,.168188E01,.147270E01,.129175E01,.113478E01,.998513,
     $.879632,.774975,.682770,.601612,.530101,.467019,.411144,.361607,
     $.317729,.278938,.244604,.214337,.187596,.164042,.143293,.125033,
     $.108996,.949100E-01,.825474E-01,.717187E-01,.622342E-01,
     $.539367E-01,.466861E-01,.403641E-01,.348523E-01,.3006E-01,
     $.258872E-01,.223E-01,.191221E-01,.164020E-01,.140484E-01,
     $.120146E-01,.102596E-01,.875578E-02,.747237E-02,.637812E-02,
     $.544411E-02,.464687E-02,.396639E-02,.338610E-02,.289071E-02,
     $.246780E-02,.210710E-02/
      DATA DEN/
     $.130098E04,.116211E04,.103651E04,.922974E03,.828225E03,
     $.741093E03,.661404E03,.588557E03,.522172E03,.461856E03,
     $.407161E03,.349556E03,.299915E03,.257234E03,.220550E03,
     $.189031E03,.161984E03,.138758E03,.118838E03,.101741E03,
     $.868975E02,.742096E02,.633844E02,.541483E02,.462409E02,
     $.395018E02,.337449E02,.288270E02,.245552E02,.209198E02,
     $.178331E02,.152106E02,.129834E02,.109993E02,.932552E01,
     $.792416E01,.674806E01,.575870E01,.492532E01,.422101E01,
     $.362450E01,.311870E01,.268847E01,.232179E01,.200867E01,
     $.174103E01,.151159E01,.131472E01,.115353E01,.101629E01,
     $.895370,.788942,.695164,.615041,.545511,.483404,.427977,
     $.378600,.334562,.295439,.260605,.229681,.202224,.177869,
     $.156306,.137232,.120385,.105501,.923507E-01,.807449E-01,
     $.705135E-01,.615126E-01,.535945E-01,.466437E-01,.405427E-01,
     $.351939E-01,.305104E-01,.264185E-01,.228440E-01,.197255E-01,
     $.170075E-01,.145146E-01,.123870E-01,.105731E-01,.902476E-02,
     $.770318E-02,.657513E-02,.561318E-02,.479197E-02,.409090E-02,
     $.349296E-02/
C
C	INTERPOLACION PARABOLICA ENTRE TRES PUNTOS
C
C	CALL POLI2
C
      I1=HA+1
      I0=I1-1
      I2=I1+1
      IF(I0.GT.0)GO TO 10
      I1=I1+1
      I2=I2+1
      I0=I0+1
10    X1=I0-1
      X2=I1-1
      X3=I2-1
      DO 3 J=1,3
         GO TO (5,6,7)J
5        Y1=PR(I0)
         Y2=PR(I1)
         Y3=PR(I2)
         CALL POLI2(X1,X2,X3,Y1,Y2,Y3,A,B,C)
         P=A+B*HA+C*HA**2
         GO TO 3
6        Y1=TE(I0)
         Y2=TE(I1)
         Y3=TE(I2)
         CALL POLI2(X1,X2,X3,Y1,Y2,Y3,A,B,C)
         T=A+B*HA+C*HA*HA
         GO TO 3
7        Y1=DEN(I0)
         Y2=DEN(I1)
         Y3=DEN(I2)
         CALL POLI2(X1,X2,X3,Y1,Y2,Y3,A,B,C)
         D=A+B*HA+C*HA**2
3     CONTINUE
      RETURN
      END
* Last processed by NICE on
* Customized for :  EEEI        HPUX        UNIX
      SUBROUTINE ASJ45(P,T,D,HA)
c J CERNICHARO
C
C	ATMOSFERA STANDARD U.S. 1962 MES DE JULIO 45 GRADOS LATITUD NORTE
C	HA === ALTURA EN KM DONDE SE DESEA CALCULAR P,T,D
C	P  === PRESION EN MILIBARES
C	D  === DENSIDAD EN GR/M**3
C	T  === TEMPERATURA EN K
C
      REAL*4 P,T,D,HA
      REAL*4 PR(91),TE(91),DEN(91)
      INTEGER I0,I1,I2,J
      REAL*4 X1,X2,X3,Y1,Y2,Y3,A,B,C
*
      DATA TE/294.15,289.65,285.154,279.156,273.168,267.174,261.186,
     $254.702,248.215,241.734,235.254,228.773,222.299,215.825,
     $215.650,215.650,215.650,215.650,216.789,217.982,219.174,
     $220.367,221.559,222.750,223.942,225.132,226.323,227.513,
     $229.492,231.573,233.654,235.735,237.814,240.225,242.697,
     $245.170,247.642,250.115,252.585,255.055,257.525,259.992,
     $262.460,264.927,267.395,269.860,272.325,274.787,275.650,
     $275.650,275.650,275.650,275.650,274.245,271.787,269.330,
     $266.872,264.417,261.960,259.507,257.052,254.600,252.147,
     $248.931,244.521,240.111,235.701,231.295,226.890,222.484,
     $218.079,213.678,209.277,204.880,200.484,196.087,191.691,
     $187.299,182.907,178.515,174.150,174.150,174.150,174.150,
     $174.150,174.150,174.150,174.150,174.150,174.150,174.150/
      DATA PR
     $/.101350E+04,.902198E+03,.801594E+03,.710433E+03,.628063E+03,
     $.553615E+03,.486632E+03,.426404E+03,.372351E+03,.324024E+03,
     $.280904E+03,.242550E+03,.208581E+03,.178571E+03,.152506E+03,
     $.130245E+03,.111251E+03,.950271E+02,.812140E+02,.694686E+02,
     $.594726E+02,.509578E+02,.437053E+02,.375159E+02,.322292E+02,
     $.277141E+02,.238504E+02,.205416E+02,.177123E+02,.152931E+02,
     $.132216E+02,.114455E+02,.992199E+01,.861267E+01,.748802E+01,
     $.651947E+01,.568409E+01,.496251E+01,.433891E+01,.379863E+01,
     $.332989E+01,.292305E+01,.256908E+01,.226070E+01,.199169E+01,
     $.175697E+01,.155167E+01,.137207E+01,.121431E+01,.107476E+01,
     $.951246,.842033,.745358,.659666,.583315,.515225,.454565,
     $.400635,.352641,.310106,.272332,.238893,.209295,.183112,
     $.159876,.139244,.120964,.104821,.905822E-01,.780538E-01,
     $.670582E-01,.574425E-01,.490475E-01,.417458E-01,.354072E-01,
     $.299215E-01,.251894E-01,.211250E-01,.176426E-01,.146698E-01,
     $.121421E-01,.100283E-01,.828248E-02,.684195E-02,.565195E-02,
     $.466893E-02,.358689E-02,.318670E-02,.263297E-02,.217545E-02,
     $.179779E-02/
      DATA DEN/
     $.119194E04,.107953E04,.975726E03,.884585E03,.799782E03,
     $.721113E03,.648694E03,.582972E03,.522462E03,.466877E03,
     $.415928E03,.369346E03,.326869E03,.288235E03,.246363E03,
     $.210402E03,.179719E03,.153510E03,.130507E03,.111022E03,
     $.945290E02,.805568E02,.687200E02,.586725E02,.501363E02,
     $.428845E02,.367119E02,.314533E02,.268872E02,.230062E02,
     $.197129E02,.169141E02,.145345E02,.124899E02,.107483E02,
     $.926367E01,.799602E01,.691194E01,.598427E01,.518837E01,
     $.450452E01,.391664E01,.340998E01,.297272E01,.259482E01,
     $.226811E01,.198495E01,.173947E01,.153465E01,.135829E01,
     $.120219E01,.106417E01,.941988,.837960,.747673,.666424,.593377,
     $.527834,.468961,.416292,.369075,.326876,.289163,.256257,
     $.227775,.202024,.178786,.157877,.139080,.122217,.107121,
     $.936510E-01,.816457E-01,.709824E-01,.615248E-01,.531584E-01,
     $.457778E-01,.392916E-01,.336024E-01,.286279E-01,.242889E-01,
     $.200605E-01,.165682E-01,.136866E-01,.113061E-01,.933969E-02,
     $.717528E-02,.637464E-02,.526696E-02,.435175E-02,.359628E-02/
C
C	INTERPOLACION PARABOLICA ENTRE TRES PUNTOS
C
C	CALL POLI2
C
      I1=HA+1
      I0=I1-1
      I2=I1+1
      IF(I0.GT.0)GO TO 10
      I1=I1+1
      I2=I2+1
      I0=I0+1
10    X1=I0-1
      X2=I1-1
      X3=I2-1
      DO 3 J=1,3
         GO TO (5,6,7)J
5        Y1=PR(I0)
         Y2=PR(I1)
         Y3=PR(I2)
         CALL POLI2(X1,X2,X3,Y1,Y2,Y3,A,B,C)
         P=A+B*HA+C*HA**2
         GO TO 3
6        Y1=TE(I0)
         Y2=TE(I1)
         Y3=TE(I2)
         CALL POLI2(X1,X2,X3,Y1,Y2,Y3,A,B,C)
         T=A+B*HA+C*HA*HA
         GO TO 3
7        Y1=DEN(I0)
         Y2=DEN(I1)
         Y3=DEN(I2)
         CALL POLI2(X1,X2,X3,Y1,Y2,Y3,A,B,C)
         D=A+B*HA+C*HA**2
3     CONTINUE
      RETURN
      END
* Last processed by NICE on
* Customized for :  EEEI        HPUX        UNIX
      FUNCTION KO2(T,P,V,IL)
C----------------------------------------------------------------------
C	OPACIDAD DE LA ATMOSFERA DEBIDA AL OXIGENO (O2)
C	T .... ES LA TEMPERATURA (K)
C	P .... ES LA PRESION EN MILIBARES
C	V .... ES LA FRECUENCIA A LA CUAL SE DESEA CALCULAR KO2
C	IL ... =0  PERFIL CINETICO
C	IL ... =1  PERFIL DE VAN VLECK & WEISSKOPF
C	KO2 .. OPACIDAD EN CM-1
C
C	ANCHURA DE LAS RAYAS DE REBER... BUENA APROXIMACION EN LAS ALAS
C
C J.Cernicharo
C----------------------------------------------------------------------
      REAL*4 KO2,T,P,V
      INTEGER IL,L,J
      REAL*4 FMEN(20),FMAS(20),RN(20),FDEB(6),B1(6),B2(6),B3(6)
*
      REAL*4 TA,V2,SUM,E0,DV1,DV,DV2,A1,A2,A3,E,PI,RD,RD0,GG,B,RR
      REAL*4 FLIN,FVVW
C
C	FMEN ... FRECUENCIAS EN GHZ  N-
C
      DATA FMEN/118.750343,62.486255,60.306044,59.164215,
     $58.323885,57.612488,56.968180,56.363393,55.783819,
     $55.221372,54.671145,54.1302,53.5959,53.0669,52.5424,
     $52.0214,51.50302,50.9873,50.4736,49.9618/
C
C	FMAS ... FRECUENCIAS EN GHZ N+
C
      DATA FMAS/56.264766,58.446580,59.590978,60.434776,
     $61.15057,61.800169,62.411223,62.997991,63.56852,
     $64.127777,64.678914,65.22412,65.764744,66.30206,
     $66.83677,67.36951,67.90073,68.4308,68.9601,69.4887/
C
C	N    ... NUMERO CUANTICO DE ROTACION
C
      DATA RN/1.,3.,5.,7.,9.,11.,13.,15.,17.,19.,21.,23.,
     $25.,27.,29.,31.,33.,35.,37.,39./
C
C	RAYAS CON DN=2
C
      DATA FDEB/368.499,424.7638,487.25,715.3944,773.841,834.147/
      DATA B1/6.79E-6,6.43E-05,2.39E-5,9.79E-6,5.71E-5,1.83E-5/
      DATA B2/.202,.0112,.0112,.0891,.0798,.0798/
      DATA B3/15.6E-4,14.7E-4,14.7E-4,14.4E-4,14E-4,14E-4/
      DATA PI/3.141592654/
      KO2=1.44E-05*P*V/T/T/T
      TA=300./T
      V2=V**2
      SUM=0.
      E0=2.07/T
      DV1=1.41E-03*P*300./T
      DV=DV1
      IF(DV1.GT.0.0527)DV=DV/3.+0.03513
      DV2=DV*DV
      DO 1 L=1,20
         A1=(RN(L)**2+RN(L)+1.)*(2.*RN(L)+1.)/RN(L)/(RN(L)+1.)
         E=E0*RN(L)*(RN(L)+1.)
         A1=A1*2.*V*DV/PI/(V2+DV2)
         A2=RN(L)*(2.*RN(L)+3.)/(RN(L)+1.)
         IF(IL.EQ.0)A2=A2*FLIN(V,FMAS(L),DV)*FMAS(L)
         IF(IL.EQ.1)A2=A2*FVVW(V,FMAS(L),DV)*FMAS(L)
         A3=(RN(L)+1.)*(2.*RN(L)-1.)/RN(L)
         B=DV
         IF(L.EQ.1)B=DV1
         IF(IL.EQ.0.)A3=A3*FMEN(L)*FLIN(V,FMEN(L),B)
         IF(IL.EQ.1.)A3=A3*FMEN(L)*FVVW(V,FMEN(L),B)
1     SUM=SUM+(A1+A2+A3)*EXP(-E)
      KO2=SUM*KO2
C
C	RAYAS CON DN=2
C
      RD=P*TA**3*4.193E-07*V
      RD0=0.
      DO 10 J=1,6
C	IF(V.LE.FDEB(J)+200..AND.V.GE.FDEB(J)-200.)GO TO 15
C	GO TO 10
15       RR=B1(J)*EXP(B2(J)*(1.-TA))
         GG=B3(J)*P*TA**.9
         IF(IL.EQ.0)RR=RR*FLIN(V,FDEB(J),GG)
         IF(IL.EQ.1)RR=RR*FVVW(V,FDEB(J),GG)
         RD0=RD0+RR
10    CONTINUE
      RD=RD*RD0
      KO2=KO2+RD
      RETURN
      END
* Last processed by NICE on
* Customized for :  EEEI        HPUX        UNIX
      FUNCTION KH2O(RHO,T,P,V,IL)
C----------------------------------------------------------------------
C	COEFICIENTE DE ABSORCION DEL VAPOR DE AGUA ATMOSFERICO
C	T ... ES LA TEMPERATURA (K)
C	P ... ES LA PRESION EN MILIBARES
C	RHO . ES LA CONCENTRACION DE VAPOR DE AGUA EN GR/M**3
C	V ... ES LA FRECUENCIA EN GHZ
C	KH2O . COEFICIENTE DE ABSORCION EN CM-1
C	IL=0   PERFIL CINETICO
C	IL=1   PERFIL DE VAN VLECK & WEISSKPOF
C
C J.Cernicharo
C----------------------------------------------------------------------
      REAL*4 RHO,T,P,V
      INTEGER IL
      REAL*4 FRE(19),GL(19),FLM(19),EL(19),DV0(19),DVLM(19),
     $X(19),B1(9),B2(9),B3(9),FDEB(9)
      REAL*4 KH2O,SUM,TA,TKK,TK,SOM,DV,RD,RD0,TT,GG,PI,FV
      INTEGER L,J
*
      REAL*4 FVVW,FLIN
C
C	FRE ... FRECUENCIAS DE LAS TRANSICIONES DEL VAPOR DE AGUA (GHZ)
C
      DATA FRE/22.23507985,183.3100906,321.225644,325.152919,380.197372,
     $390.14,437.346667,439.150812,443.018295,448.001075,
     $470.888947,474.689127,488.491133,556.936002,620.700807,
     $752.033227,916.62,970.31,987.94/
C
C	GL ... DEGENERACION DE LOS NIVELES
C
      DATA GL/3.,1.,3.,1.,3.,1.,1.,3.,3.,3.,1.,1.,1.,3.,3.,
     $1.,1.,1.,1./
C
C	FLM ... CUADRADO DEL ELEMENTO DE MATRIZ LM
C
      DATA FLM/.057,.102,.089,.091,.123,.068,.088,.0101,.088,
     $.132,.102,.118,.036,1.5,.122,2.073,.161,.262,.7557/
C
C	EM ... ENERGIAS EN CM-1 DEL NIVEL SUPERIOR EN LA TRANS L=>M
C	EL ...     "     "  "    "    "   INFERIOR  "  "   "     "
C
C	DATA EM/447.3,142.27,1293.8,326.62,224.84,1538.31,1059.63,
C	1 756.76,1059.90,300.37/
      DATA EL/446.56,136.16,1283.02,315.78,212.16,1525.31,
     $1045.03,742.11,1045.11,285.42,742.074,488.135,586.48,
     $23.794,488.108,70.091,285.217,383.837,37.137/
      DATA DV0/2.85,2.68,2.3,3.03,3.19,2.11,1.5,1.94,1.51,
     $2.47,1.89,2.07,2.58,3.33,2.28,3.13,2.59,2.48,3.09/
      DATA DVLM/13.68,14.49,12.04,15.21,15.84,11.42,7.94,10.44,
     $8.13,14.24,10.56,11.95,14.77,14.66,12.78,13.93,14.06,14.16,
     $15.20/
C
C	X  ... EXPONENTE DE LA TEMPERATURA
C
      DATA X/.626,.649,.42,.619,.63,.33,.29,.36,.332,.51,
     $.380,.38,.57,.645,.6,.69,.676,.56,.66/
      DATA FDEB/68.052,503.56,504.46,658.34,841.01,859.81,899.38,
     $903.28,906.21/
      DATA B1/1.8E-3,3.5E-3,1.2E-3,4.6E-2,1.2E-3,1.5E-3,9.1E-3,
     $6.4E-3,1.79E-2/
      DATA B2/8.75,6.69,6.69,7.76,8.11,7.99,7.84,8.35,5.04/
      DATA B3/2.8E-3,1.27E-3,1.3E-3,3.28E-3,1.7E-3,2.7E-3,3E-3,
     $2.8E-3,2.04E-3/
      DATA PI/3.141592654/,TK/.69503096/
      KH2O=1.44*RHO*V/SQRT(T**3)
      SUM=0.
      TA=300./T
      TKK=TK*T
      DO 1 L=1,19
C	IF(V.LE.FRE(L)+200..AND.V.GE.FRE(L)-200.)GO TO 5
C	GO TO 1
5        SOM=GL(L)*FLM(L)*EXP(-EL(L)/TKK)*(1.-EXP(-FRE(L)/TKK/29.97925))
         DV=DV0(L)*P/1013./((T/300.)**X(L))
         DV=DV*(1.+(4.6E-03*RHO*T/P)*(DVLM(L)/DV0(L)-1.))
         IF(IL.EQ.0)FV=FLIN(V,FRE(L),DV)
         IF(IL.EQ.1)FV=FVVW(V,FRE(L),DV)
         SUM=SUM+SOM*FV
1     CONTINUE
      KH2O=KH2O*SUM
C
C	TERMINO CORRECTOR EMPIRICO (POSIBLE CONTRIBUCION DE LOS DIMEROS
C				    DE H2O)
C
      KH2O=KH2O+1.08E-11*RHO*V*V*P/1000.*(TA)**2.1
C
C	RAYAS DEBILES
C
      RD=1.937E-9*V*RHO*T
      RD0=0.
      DO 10 J=1,9
C	IF(V.LE.FDEB(J)+100..AND.V.GE.FDEB(J)-100.)GO TO 15
C	GO TO 10
15       TT=B2(J)*(1.-TA)
         TT=B1(J)*EXP(TT)*(TA)**3.5
         GG=TA**0.6*P*B3(J)
         IF(IL.EQ.0)FV=FLIN(V,FDEB(J),GG)
         IF(IL.EQ.1)FV=FVVW(V,FDEB(J),GG)
         RD0=RD0+TT*FV
10    CONTINUE
      RD=RD*RD0
      KH2O=KH2O+RD
      RETURN
      END
* Last processed by NICE on
* Customized for :  EEEI        HPUX        UNIX
      FUNCTION FVVW(V,VL,DV)
C----------------------------------------------------------------------
C	PERFIL DE VAN VLECK & WEISSKOPF
C
C J CERNICHARO
C----------------------------------------------------------------------
      REAL*4 FVVW,V,DV,VL,PI,DV2,A1,A2
      DATA PI/3.141592654/
*
      FVVW=DV*V/VL/PI
      DV2=DV*DV
      A1=DV2+(V-VL)**2
      A2=DV2+(V+VL)**2
      FVVW=FVVW*(1./A1+1./A2)
      END
* Last processed by NICE on
* Customized for :  EEEI        HPUX        UNIX
      FUNCTION FLIN(V,VL,DV)
C----------------------------------------------------------------------
C J Cernicharo, Model atmosphere.
C
C	FORMA CINETICA DEL PERFIL
C	V... FRECUENCIA
C	VL.. FRECUENCIA DE LA LINEA
C	DV.. ANCHURA DE LA LINEA
C----------------------------------------------------------------------
      REAL*4 FLIN,V,VL,DV,PI,V2
      DATA PI/3.141592654/
*
      FLIN=4.*V*VL*DV/PI
      V2=V*V
      FLIN=FLIN/(4.*V2*DV*DV+(VL*VL-V2)**2)
      END
