C+
C NAME:
C	FillMapOLS
C PURPOSE:
C	To fill bad spots in maps at opposite the region of interest with a smooth average of the 
C	data from the other nearby longitudes. This includes points on both east and west sides 
C	of the volume.
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call FillMapOLS(Mo,XCbe,XCint,nLng,nLat,nT,N,nTmax,RRS,Speed,Const,AEmax,AEmin,AEang,DMap)
C INPUTS:
C	Mo		   integer	1 - Smooth/interpolate region most distant from Earth
C					2 - Limit values farthest from Earth
C	XCbe(2,nT)	   real		Boundary values of time maps
C	XCint(NTmax)	   real		Start and end times of intervals
C	nLng		   integer	# Longitude points
C	nLat		   integer	# Latitude points
C	nT		   integer	# Time points
C	N		   integer	Map time value
C	nTmax		   integer	Total # of map time values
C	RRS		   real		Deconvolution source surface distance
C	Speed		   real		Average speed used for data analysis
C	Const		   real		Longitude filter constant (in terms of the numbers of intervals in one rotation)
C	AEmax		   real		Anti-Earth maximum value
C	AEmin		   real		Anti-Earth minimum value
C	AEang		   real		Angle from Earth to begin longitude limit
C	DMap(nLng,nLat,nT) real		Input map
C OUTPUTS:
C	DMap(nLng,nLat,nT) real		Output map
C FUNCTIONS/SUBROUTINES:
C PROCEDURE:
C	Bad values (indicated by BadR4()) are processed
C MODIFICATION HISTORY:
C	MAY, 1999 B. Jackson (STEL,UCSD), Modified 2008/03/03 BVJ
C-
	subroutine FillMapOLS(Mo,XCbe,XCint,nLng,nLat,nT,N,NTmax,RRS,Speed,Const,AEmax,AEmin,AEang,DMap)

	real		XCbe(2,nT),			! Boundary values of time maps
     &			XCint(nTmax),			! Start and end times of intervals		
     &			DMap(nLng,nLat,nT)		! Map
C
C  Average the other non-Bad longitude values.
C
	
	if(Speed.ne.0.0) then
	    XCmid = ((XCint(N)+XCint(N+1))/2.0) - ((0.8-RRS)*((400./Speed)*4.)/25.38)
	else
	    XCmid = ((XCint(N)+XCint(N+1))/2.0)
	end if
	XCbeg = XCbe(1,N)
	XCend = XCbe(2,N)

	Bad = BadR4()
	ExpPWC = 50.

	NXC  = (XCend-XCbeg) + 0.1			! Usually 3
	Lmod = nLng/NXC					! Usually 18 for the IPS
	IXCopp = (XCend - XCmid)*Lmod + 0.5		! Middle I of region of interest
C
C	This part Filters or smooths the data
C
	if(Mo.eq.1) then
		
	  nlngB1 = IXCopp - Lmod/2
	  nLngE1 = nLngB1 + 1
	  nLngE2 = IXCopp + Lmod/2
	  nLngB2 = nLngE2 - 1
	  NX = 2

C	  print *, 'NXC, Lmod, XCbeg, XCend, XCmid, IXCopp, nLngB, nLngE', NXC, Lmod, XCbeg, XCend, XCmid, IXCopp, nLngB1, nLngE1, nLngB2, nLngE2

	  do J=1,nLat
	     do I=nLngB1,nLngE1
	        DMNT   = 0.
	        ANDMNT = 0.
		do II=1,NX
		    IP = I + II
		    IM = I - II
		    if(IP.le.nLng) then
			if(DMap(IP,J,N).ne.Bad) then
			    ExpPW = (II/Const)**2
			    Expot = 999999.
			    if(ExpPW.gt.ExpPWC) Expot = 0.
			    if(Expot.ne.0.) Expot = Exp(-ExpPW)
			    DMNT   = DMNT   + DMap(IP,J,N)*Expot 			   
			    ANDMNT = ANDMNT + Expot
			end if		
		    end if
		    if(IM.ge.1) then
			if(DMap(IM,J,N).ne.Bad) then
			    ExpPW = (II/Const)**2
			    Expot = 999999.
			    if(ExpPW.gt.ExpPWC) Expot = 0.
			    if(Expot.ne.0.) Expot = Exp(-ExpPW)
			    DMNT   = DMNT   + DMap(IM,J,N)*Expot 			   
			    ANDMNT = ANDMNT + Expot
			end if		
		    end if
		end do
C	        if(ANDMNT.ne.0.) then  			! Changed BVJ 2008/03/03

C		    if(DMap(I,J,N).eq.Bad) then
C			DMap(I,J,N) = DMNT/ANDMNT
C		    end if

C	        end if  		 		  		 				
	     end do	
	     do I=nLngB2,nLngE2
	        DMNT   = 0.
	        ANDMNT = 0.
		do II=1,NX
		    IP = I + II
		    IM = I - II
		    if(IP.le.nLng) then
			if(DMap(IP,J,N).ne.Bad) then
			    ExpPW = (II/Const)**2
			    Expot = 999999.
			    if(ExpPW.gt.ExpPWC) Expot = 0.
			    if(Expot.ne.0.) Expot = Exp(-ExpPW)
			    DMNT   = DMNT   + DMap(IP,J,N)*Expot 			   
			    ANDMNT = ANDMNT + Expot
			end if		
		    end if
		    if(IM.ge.1) then
			if(DMap(IM,J,N).ne.Bad) then
			    ExpPW = (II/Const)**2
			    Expot = 999999.
			    if(ExpPW.gt.ExpPWC) Expot = 0.
			    if(Expot.ne.0.) Expot = Exp(-ExpPW)
			    DMNT   = DMNT   + DMap(IM,J,N)*Expot 			   
			    ANDMNT = ANDMNT + Expot
			end if		
		    end if
		end do		 		  		 				
	     end do
	     
	     do I=nLngB1,nLngE1
	       if(ANDMNT.ne.0.) then

		 if(DMap(I,J,N).eq.Bad) then
		   DMap(I,J,N) = DMNT/ANDMNT
		 end if

	       end if  
	     end do
	     do I=nLngB2,nLngE2
	       if(ANDMNT.ne.0.) then

		 if(DMap(I,J,N).eq.Bad) then
		   DMap(I,J,N) = DMNT/ANDMNT
		 end if

	       end if  
	     end do

	  end do

	end if
C
C	This part limits the Anti-Earth maps beyond an angle AEang
C
	if(Mo.eq.2) then

	  Idist = Lmod*AEang/360. + .5			! Integer distance of anti-Earth angle
	  Ibeg = IXCopp - Idist
	  Iend = IXCopp + Idist

C	  print *, 'XCmid, XCbeg, IXCopp, Idist, Ibeg, Iend', XCmid, XCbeg, IXCopp, Idist, Ibeg, Iend

	  do I=1,nLng

	    if(I.lt.Ibeg .or. I.gt.Iend) then
	    
	      do J=1,nLat

		if(DMap(I,J,N).ne.Bad) then
		    if(Dmap(I,J,N).lt.AEmin) Dmap(I,J,N) = AEmin
		    if(Dmap(I,J,N).gt.AEmax) Dmap(I,J,N) = AEmax
		end if

	      end do
	    end if
	  end do

	end if

	return
	end