C+
C NAME:
C	FillMap
C PURPOSE:
C	Fills maps at a given longitude with an average of the data from
C 	the other longitudes. Found a big mistake 10/30/00 BVJ. IP and IM
C	were taking data from the wrong place.
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
	subroutine FillMap(Mode,XCbeg,XCend,nLng,nLat,Const,DMap)
C INPUTS:
C	Mode		   integer	0 - only bad values are processed
C					1 - all values are processed
C					2 - only good values are processed, all weighted equally
C					3 - make all 90 deg. longitudes the same value
C					4 - make all valid 90 deg. longitudes the same value
C					5 - make all 90 deg. longitudes the same value if equatorial longitudes are valid
C	XCbeg		   real		Beginning rotation values
C	XCend		   real		Ending rotation value
C	nLng		   integer	# Longitude points
C	nLat		   integer	# Latitude points
C	Const		   real		Longitude filter constant (in terms of the numbers of intervals in one rotation)
C	DMap(nLng,nLat) real		Input map
C OUTPUTS:
C	DMap(nLng,nLat) real		Output map
C CALLS:
C	BadR4
C PROCEDURE:
C	Bad values (indicated by BadR4()) are processed
C MODIFICATION HISTORY:
C	MAY, 1999 B. Jackson (STEL,UCSD)
C	JUN-2001, Paul Hick (UCSD/CASS)
C	    Renamed the subroutine from FillMapL to FillMap.
C	    Both tomography programs (ipsd and ipsdt call this subroutine).
C-
	real		DMap(nLng,nLat)		! Map

	! Average the other non-Bad longitude values.

	Bad = BadR4()
	ExpPWC = 50.

	! This part Filters or smooths the data

	NXC  = (XCend-XCbeg) + 0.1
	Lmod = nLng/NXC
	JJ   = nLat/2 + 1			! middle latitude

	do I=1,nLng
	    do J=1,nLat
		DMNT   = 0.
		ANDMNT = 0.
		Ithr = 0
		if(DMap(I,J).eq.Bad.and.Mode.eq.0) Ithr = 1
		if(Mode.eq.1) Ithr = 1
		if(Ithr.eq.1) then
		    do II=1,NXC
			IP = I + II*Lmod
			IM = I - II*Lmod
			if(IP.le.nLng) then
			    if(DMap(IP,J).ne.Bad) then
				ExpPW = ((II*Lmod)/Const)**2
				Expot = 999999.
				if(ExpPW.gt.ExpPWC) Expot = 0.
				if(Expot.ne.0.) Expot = Exp(-ExpPW)
				DMNT   = DMNT   + DMap(IP,J)*Expot
				ANDMNT = ANDMNT + Expot
			    end if
			end if

			if(IM.ge.1) then
			    if(DMap(IM,J).ne.Bad) then
				ExpPW = ((II*Lmod)/Const)**2
				Expot = 999999.
				if(ExpPW.gt.ExpPWC) Expot = 0.
				if(Expot.ne.0.) Expot = Exp(-ExpPW)
				DMNT   = DMNT   + DMap(IM,J)*Expot
				ANDMNT = ANDMNT + Expot
			    end if
			end if
		    end do
		end if

		! This part when Mode = 2, blends the data so that there is no seam

		if(DMap(I,J).ne.Bad.and.Mode.eq.2) then
		    do II=1,NXC
			IP = I + II*Lmod
			IM = I - II*Lmod
			if(IP.le.nLng) then
			    if(DMap(IP,J).ne.Bad) then
				DMNT   = DMNT   + DMap(IP,J)
				ANDMNT = ANDMNT + 1.
			    end if
			end if
			if(IM.ge.1) then
			    if(DMap(IM,J).ne.Bad) then
				DMNT   = DMNT   + DMap(IM,J)
				ANDMNT = ANDMNT + 1.
			    end if
			end if
		    end do
		end if

		! This part when Mode = 3, sums all valid data at + or - 90 deg.

		if(Mode.ge.3.and.I.eq.1) then
		    if(J.eq.1.or.J.eq.nLat) then
			do II=1,nLng
			    if(DMap(II,J).ne.Bad) then
				DMNT   = DMNT   + DMap(II,J)
				ANDMNT = ANDMNT + 1.
			    end if
			end do
		    end if
		end if

		if(ANDMNT.ne.0.) then

		    if(Mode.eq.0)then
			DMap(I,J) = DMNT/ANDMNT
		    end if

		    if(Mode.eq.1)then
			if(DMap(I,J).eq.Bad)then
			    DMap(I,J) = DMNT/ANDMNT
			else
			    DMap(I,J) = (DMap(I,J) + DMNT)/(ANDMNT + 1.)
			end if
		    end if

		    if(Mode.eq.2)then
			DMap(I,J) = (DMap(I,J) + DMNT)/(ANDMNT + 1.)
		    end if

		    if(Mode.eq.3)then
			do II=1,nLng
			    DMap(II,J) = DMNT/ANDMNT
			end do
			!print *, 'Mode = 3', I, J, DMNT, ANDMNT
		    end if

		    if(Mode.eq.4)then
			do II=1,nLng
			    if(DMap(II,J).ne.Bad)then
				DMap(II,J) = DMNT/ANDMNT
			    end if
			end do
			!print *, 'Mode = 4', I, J, DMNT, ANDMNT
		    end if

		    if(Mode.eq.5)then
			do II=1,nLng
			    if(DMap(II,JJ).ne.Bad)then
				DMap(II,J) = DMNT/ANDMNT
			    end if
			end do
			!print *, 'Mode = 5', I, J, DMNT, ANDMNT
		    end if

		end if

	    end do
	end do

	return
	end