C+
C NAME:
C	MkVMap
C PURPOSE:
C	Change the velocity map from the line of sight (L.O.S.) values of velocity and
C	weights by the ratio of observed to model values for each line of sight. 
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
	subroutine MkVMap(LON,LAT,RP,XEV,VWTij,FIX,NL,NLOS,NLOSP1,DMAP,DREF,
     &		XLON,XLAT,DFAC,nLng,nLat,PW,VMAP,ANLIMIT,VLIM,VSSIG,VM,SIGB,NS,
     &		NBADSV,SIGMAP,ANMAP,VMA,FIXM,WTSM,VMEANF,VFIXM2,WTP,CONR)
C INPUTS:
C	LON(NLOSP1,NL)		integer		Projected Carr. rot. value of point on L.O.S.
C	LAT(NLOSP1,NL)		integer		Heliographic lat. point on L.O.S.
C	RP (NLOSP1,NL)		real		Distance above Sun of point on L.O.S.
C	XEV(NL)			real		Source Elongation
C	VWTij(NLOS,NL)		real		Weights of each point along the L.O.S.
C	FIX(NL)			real		Ratio observed/model V for total L.O.S.
C	NL			integer		Number of velocity data points
C	NLOS			integer		Number of L.O.S. distance segments
C	NLOSP1			integer		Number of L.O.S. distance segments + 1
C	DMAP(nLng,nLat)		real		Density map
C	DREF			real		Density map reference height
C	XLON(NLOSP1,NL)		real		Projected Carr. rot. value of point on L.O.S.
C	XLAT(NLOSP1,NL)		real		Latitude (in deg.) of point on L.O.S.
C	DFAC(NLOSP1,NL)		real		Interpolated DF value between maps at RP
C	nLng			integer		Number of longitude points in DMAP
C	nLat			integer		Number of latitude points in DMAP
C	PW			real		Power of density relationship to g-level
C	VMAP(nLng,nLat)		real		Velocity map
C	SIGMAP			real		Velocity map convergence criteria
C	ANLIMIT			real		Number of points constituting deconvolution  
C	VLIM			real		Limit on velocity (maximum)
C	VSSIG(NL)		real		Deviation of source from mean in sigma
C	VM(NL)			real		Model value
C	SIGB			real		Velocity sigma cut-off
C	NS			integer		1 = Activate sigma cut-off, 0 = do not
C	NBADSV(NL)		integer		Bad source velocity 0 - bad, 1 - good
C	CONR			real		Gridsphere filter factor
C SCRATCH FILES:
C	ANMAP(nLng,nLat)	real		Scratch array
C	VMA(nLng,nLat)		real		Scratch array
C	FIXM(nLng,nLat)		real		Scratch array
C	WTSM(nLng,nLat)		real		Scratch array
C	VMEANF(nLng,nLat)	real		Scratch array
C	VFIXM2(nLng,nLat)	real		Scratch array
C	WTP(NLOS,NL)		real		Scratch array 
C OUTPUTS:
C	VMAP(nLng,nLat)		real		Velocity map
C	NBADSV(NL)		integer		Bad source velocity 0 - bad, 1 - good
C	SIGMAP			real		Velocity map convergence criteria
C CALLS:
C	BadR4, ArrR4Bad, ArrR4Zero, GridSphere2D, Say
C PROCEDURE:
C MODIFICATION HISTORY:
C	NOV, 1995 B. Jackson (STEL,UCSD)
C-
	integer		LON(NLOSP1,NL),		! Projected Carr. rot. value of point on L.O.S.  
     &			LAT(NLOSP1,NL),		! Heliographic lat. (in deg.) point on L.O.S.
     &			NBADSV(NL)		! Bad source velocity 0 - bad, 1 - good
	real		RP  (NLOSP1,NL),	! Distance above Sun of point on L.O.S.
     &			XLON(NLOSP1,NL),	! Projected Carr. rot. value of point on L.O.S.
     &			XLAT(NLOSP1,NL),	! Latitude (in deg.) of point on L.O.S.
     &			VWTij(NLOS,NL),		! Weights of each point along the L.O.S.
     &			DFAC(NLOSP1,NL),	! Interpolated DF value between maps at RP
     &			DMAP(nLng,nLat),	! Density map
     &			VMAP(nLng,nLat),	! Velocity map
     &			FIX  (NL),		! Ratio observed/model G for total L.O.S.
     &			XEV  (NL),		! Source elongation 
     &			VSSIG(NL),		! Deviation of source from mean in sigma
     &			VM   (NL)		! Model value
	real		ANMAP(nLng,nLat),
     &			VMA(nLng,nLat),
     &			FIXM(nLng,nLat),
     &			WTSM(nLng,nLat),
     &			VMEANF(nLng,nLat),
     &			VFIXM2(nLng,nLat),
     &			WTP(NLOS,NL)

	character	cStr*150

	include		'mapangles.h'

	nLngLat = nLng*nLat
	nLOSnL  = NLOS*NL
	Bad = BadR4()
	nLat1 = nLat - 1
	daLat2 = 90./nLat1
	call ArrR4Bad	(nLngLat,VMA)
	call ArrR4Zero	(nLngLat,FIXM)
	call ArrR4Bad	(nLngLat,WTSM)
	call ArrR4Zero	(nLngLat,VMEANF)
	call ArrR4Zero	(nLngLat,VFIXM2)
	call ArrR4Bad	(nLngLat,ANMAP)
	call ArrR4Zero	(nLOSnL,WTP)
	NSS = 0
	do I=1,NL
	    if (NBADSV(I) .ne. 0) then
		if (abs(VSSIG(I)) .gt. SIGB) then	! source above limit
		    NSS = NSS+1
		    if (NS .eq. 1) NBADSV(I) = 0	! Mark source as bad
	    	end if
	    end if
	    if (NBADSV(I) .ne. 0) then
		LNlast = 0
		LTlast = 0
		sinE = sind(abs(XEV(I)))
		do J=1,NLOS
		    LN = LON(J,I)			! map longitude
		    LT = LAT(J,I)			! map latitude
c		    VSN = sinE/RP(J,I)
		    WTSj = VWTij(J,I)
		    Fws = WTSj*FIX(I)
		    if (VMAP(LN,LT) .ge. VLIM .and. FIX(I) .gt. 1.0) Fws = WTSj
		    FIXM(LN,LT) = FIXM(LN,LT)+Fws
		    WTP(J,I) = WTSj
		    if(WTSM(LN,LT).eq.Bad) WTSM(LN,LT) = 0.
		    WTSM(LN,LT) = WTSM(LN,LT)+WTP(J,I)
		    if (LN .ne. LNlast .or. LT .ne. LTlast) then
			if(ANMAP(LN,LT).eq.Bad) ANMAP(LN,LT) = 0.
			ANMAP(LN,LT) = ANMAP(LN,LT)+1
		    end if
		    LNlast = LN
		    LTlast = LT
 		end do
	    end if
	end do
	call GridSphere2D(ALng,nLng,nLat,1,ANMAP,CONR,3,0.0,90.0)

C  Now modify the latitude map count to approximate the cos latitude filter effect

	do J=1,nLat
	    alat = abs(XLdeg(J))
	    if(alat.eq.90.) alat = alat - dalat2
	    sqonebycos = SQRT(1./cosd(alat))
	    do I=1,nLng
		if(ANMAP(I,J).ne.Bad) then
		    ANMAP(I,J) = ANMAP(I,J)*sqonebycos
		end if
	    end do		    
	end do

	do LN=1,nLng
	    do LT=1,nLat
		if (FIXM(LN,LT) .gt. 0. .and. WTSM(LN,LT) .gt. 0.) then
		    VMEANF(LN,LT) = FIXM(LN,LT)/WTSM(LN,LT) 
		    VMA(LN,LT) = VMEANF(LN,LT)*VMAP(LN,LT)
		else
		    ANMAP(LN,LT) = 0.
		end if
	    end do
	end do
C  calculation of the Chi square distribution for each latitude and longitude
	do I=1,NL
	    if (NBADSV(I) .ne. 0) then
	    do J=1,NLOS
		LN = LON(J,I)			! map longitude
		LT = LAT(J,I)			! map latitude
		if (VMA(LN,LT) .ne. Bad) VFIXM2(LN,LT) = VFIXM2(LN,LT)+WTP(J,I)*((FIX(I)-VMEANF(LN,LT))**2)
	    end do
	    end if
	end do
C  calculation of the "total" Chi square distribution for all latitudes and longitudes
	CONVER  = 0
	NCONVER = 0
	do I=1,nLng
	    do J=1,nLat
		if (ANMAP(I,J) .lt. ANLIMIT) VMA(I,J) = Bad ! data not deconvolved
		if (VMA(I,J) .ne. Bad) then
		    CONVER  = CONVER +VFIXM2(I,J)
		    NCONVER = NCONVER+1
		end if
		VMAP(I,J) = VMA(I,J)
	    end do
	end do
	SIGOLD = SIGMAP
	SIGMAP = CONVER/NCONVER
	write (cStr,'(A,F4.1,A,F5.2,A)') 'velocity sources above',SIGB,' sigma:',100.*NSS/NL,'%'
	if (NS .eq. 1) write (cStr(itrim(cStr)+1:),'(A,I4,A)') '#you just removed',NSS,' sources'
	write (cStr(itrim(cStr)+1:),'(A,F7.3,A,F7.3,A)') '#V map convergence =',SIGMAP,
     &		', from last =',abs(1-SIGMAP/SIGOLD)*100,'%'
	call Say('MkVMap','I','Info',cStr)
	return
	end