
      SUBROUTINE gtbmap(bkgfil, bkgmap, wmsiz1, wmsiz2, wmstrt, wmend, 
     &                  status)

      INTEGER wmsiz1, wmsiz2, wmstrt(2), wmend(2), status

      REAL bkgmap(wmsiz1, wmsiz2)

      CHARACTER*(*) bkgfil

c Subroutine to read a background map file. This routine is clever 
c in as much as it can operate with a map that is a factor of the 
c WMAP size but it is not clever in that it does no checking whether 
c the map was generated by the same instrument. Note that the WMAP
c size is that of the actual input, not the WMAP with extra bits added
c on the ends as used internally in this program.

c Arguments :
c      bkgfil    c      i: Name of background map
c      bkgmap    r      r: Background map
c      wmsize    i      i: WMAP (and background map) size
c      status    i      r: 0 == OK

c  kaa  2/4/97

      INTEGER inmap

      REAL value

      INTEGER ilun, block, baxsze(2), factor(2), wmsize(2)
      INTEGER i, j, k, l, k1, l1, ipt

      CHARACTER contxt*72, comment*72

      LOGICAL qxpand(2), anyf

      INTEGER fcstln
      EXTERNAL fcstln

* Dynamic memory allocation stuff
C  the following MEM common block definition is in the system iraf77.inc file
C
C Get IRAF MEM common into main program.
C
      LOGICAL          MEMB(1)
      INTEGER*2        MEMS(1)
      INTEGER*4        MEMI(1)
      INTEGER*4        MEML(1)
      REAL             MEMR(1)
      DOUBLE PRECISION MEMD(1)
      COMPLEX          MEMX(1)
      EQUIVALENCE (MEMB, MEMS, MEMI, MEML, MEMR, MEMD, MEMX)
      COMMON /MEM/ MEMD

C     datatype gives a symbolic code for the data type, e.g.,  4 = Integer*4
C     1 is boolean
C     3 is short integer
C     4 is integer
C     5 is long integer
C     6 is single precision real
C     7 is double precision real
C     8 complex



      status = 0
      DO i = 1, wmsiz2
         DO j = 1, wmsiz1
            bkgmap(j,i) = 0.
         ENDDO
      ENDDO

c If no background map file is set then return

      IF ( bkgfil .EQ. 'none' ) RETURN

c Open the background map file

      CALL getlun(ilun)
      CALL ftopen(ilun, bkgfil, 0, block, status)
      WRITE(contxt, '(a,a)') 'Error opening background file : ',
     &                       bkgfil(:MIN(fcstln(bkgfil),len(contxt)-32))
      IF ( status .NE. 0 ) GOTO 999

c Get the size of the map

      CALL ftgkyj(ilun, 'NAXIS1', baxsze(1), comment, status)
      contxt = 'Failed to read NAXIS1 keyword'
      IF ( status .NE. 0 ) GOTO 999
      CALL ftgkyj(ilun, 'NAXIS2', baxsze(2), comment, status)
      contxt = 'Failed to read NAXIS2 keyword'
      IF ( status .NE. 0 ) GOTO 999

c Find out whether this map is compatible with the WMAP.

      wmsize(1) = wmsiz1
      wmsize(2) = wmsiz2
      DO i = 1, 2
         IF ( MOD((wmend(i)-wmstrt(i)+1),baxsze(i)) .EQ. 0 ) THEN
            qxpand(i) = .TRUE.
            factor(i) = (wmend(i)-wmstrt(i)+1)/baxsze(i)
         ELSEIF ( MOD(baxsze(i),(wmend(i)-wmstrt(i)+1)) .EQ. 0 ) THEN
            qxpand(i) = .FALSE.
            factor(i) = baxsze(i)/(wmend(i)-wmstrt(i)+1)
         ELSE
            CALL fcecho(
     &       'The input background map is not compatible with the WMAP')
            WRITE(contxt,'(a,i1,a,i5,a,i5,a)') 'Axis ', i, ' has size ',
     &       baxsze(i), ' in the background map and ', 
     &       (wmend(i)-wmstrt(i)+1), ' in the WMAP'
            CALL fcecho(contxt)
            status = -1
            RETURN
         ENDIF
      ENDDO

      IF ( ( qxpand(1) .AND. .NOT.qxpand(2) ) .OR.
     &     ( qxpand(2) .AND. .NOT.qxpand(1) ) ) THEN
         CALL fcecho(
     &   'One axis of the background map is larger than the WMAP while')
         CALL fcecho(
     &   'the other axis is smaller. I cannot handle this case')
         status = -1
         RETURN
      ENDIF

c We seem to be OK so grab the memory for the input map and read it in

      inmap = 0
      CALL udmget(baxsze(1)*baxsze(2), 6, inmap, status)
      contxt = 'Insufficient memory for inmap array'
      IF ( status .NE. 0 ) goto 999

      CALL ftgpve(ilun, 0, 1, baxsze(1)*baxsze(2), -99., MEMR(inmap), 
     &            anyf, status)
      contxt = 'Failed to read the background map array'
      IF ( status .NE. 0 ) goto 999

      CALL ftclos(ilun, status)
      contxt = 'Failed to close the background map file'
      IF ( status .NE. 0 ) goto 999
      CALL frelun(ilun)

c Now transform this input into a map the same size as the WMAP. We either
c compress the map down by adding together bins or we expand it out by
c dividing the values evenly into the new bins.

      IF ( qxpand(1) ) THEN

         ipt = inmap
         DO i = 1, baxsze(2)
            DO j = 1, baxsze(1)

               value = MEMR(ipt) / factor(1) / factor(2)
               DO k = 1, factor(2)
                  k1 = (i-1)*factor(2)+k+wmstrt(2)-1
                  DO l = 1, factor(1)
                     l1 = (j-1)*factor(1)+l+wmstrt(1)-1
                     bkgmap(l1,k1) = value
                  ENDDO
               ENDDO
               ipt = ipt + 1

            ENDDO
         ENDDO

      ELSEIF ( .NOT.qxpand(1) ) THEN

         ipt = inmap
         DO i = 1, baxsze(2)
            k = (i-1)/factor(2) + 1 + wmstrt(2) - 1
            DO j = 1, baxsze(1)
               l = (j-1)/factor(1) + 1 + wmstrt(1) - 1

               bkgmap(l,k) = bkgmap(l,k) + MEMR(ipt)

               ipt = ipt + 1

            ENDDO
         ENDDO

      ENDIF

 999  CONTINUE
      IF ( status .NE. 0 ) THEN
         CALL fcerr(contxt)
         CALL fcerrm(status)
      ENDIF

      RETURN
      END
