C+
C NAME:
C	MkVMaptdN
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:
C	call MkVMaptdn(LON,LAT,ITIM,RP,XEV,VWTij,FIX,NL,NLOS,NLOSP1,XLON,XLAT,
C		nLng,nLat,nT,VMAP,ANLIMIT,VLIM,VSSIG,VM,SIGB,NS,NBADSV,SIGMAP,
C		ANMAP,VMA,FIXM,WTSM,VMEANF,VFIXM2,WTP,Alng,ConstL,CONR,CONT,aNday,Ztmp)
C     
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	ITIM(NLOSP1,NL)		integer		Projected time values
C	RP (NLOSP1,NL,NT)	real		Distance above Sun of point on L.O.S.
C	XEV(NL)			real		Source Elongation
C	VWTij(NLOS,NL,NT)	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	XLON(NLOSP1,NL,NT)	real		Projected Carr. rot. value of point on L.O.S.
C	XLAT(NLOSP1,NL,NT)	real		Latitude (in deg.) of point on L.O.S.
C	nLng			integer		Number of longitude points in DMAP
C	nLat			integer		Number of latitude points in DMAP
C	nT			integer		Number of time intervals
C	VMAP(nLng,nLat,nT)	real		Velocity map
C	SIGMAP			real		Velocity map convergence criteria
C	ANLIMIT			integer		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 OUTPUTS:
C	VMAP(nLng,nLat,nT)	real		Velocity map
C	NBADSV(NL)		integer		Bad source velocity 0 - bad, 1 - good
C	SIGMAP			real		Velocity map convergence criteria
C SCRATCH ARRAYS:
C	ANMAP(nLng,nLat,nT)	integer
C	VMA(nLng,nLat,nT)	real
C	FIXM(nLng,nLat,nT)	real
C	WTSM(nLng,nLat,nT)	real
C	VMEANF(nLng,nLat,nT)	real
C	VFIXM2(nLng,nLat,nT)	real
C	WTP(NLOS,NL)		real
C	Ztmp(nLng,nLat,nT)	real
C INCLUDE:
C	include		'mapangles.h'
C PROCEDURE:
C MODIFICATION HISTORY:
C	NOV, 1995 B. Jackson (STEL,UCSD)
C-
	subroutine MkVMaptdn(LON,LAT,ITIM,RP,XEV,VWTij,FIX,NL,NLOS,NLOSP1,XLON,XLAT,
     &		nLng,nLat,nT,VMAP,ANLIMIT,VLIM,VSSIG,VM,SIGB,NS,NBADSV,SIGMAP,
     &		ANMAP,VMA,FIXM,WTSM,VMEANF,VFIXM2,WTP,Alng,ConstL,CONR,CONT,aNday,Ztmp)
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.
     &			ITIM(NLOSP1,NL),	! Projected time values
     &			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.
     &			VMAP(nLng,nLat,nT),	! 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

C 	Scratch files

	real		ANMAP(nLng,nLat,nT),	! Number of points per bin
     &			VMA(nLng,nLat,nT),
     &			FIXM(nLng,nLat,nT),		
     &			WTSM(nLng,nLat,nT),	
     &			VMEANF(nLng,nLat,nT),	
     &			VFIXM2(nLng,nLat,nT),	
     &			WTP(NLOS,NL),
     &			Ztmp(nLng,nLat,nT)

	character	cStr*150

	include		'mapangles.h'

	nLngLat	= nLng*nLat
	nLngLatnT = nLng*nLat*nT
	nLOSnL	  = NLOS*NL
	Bad = BadR4()
	nLat1 = nLat - 1
	daLat2 = 90./nLat1

	call ArrR4Bad (nLngLatnT,VMA)
	call ArrR4Zero(nLngLatnT,FIXM)
	call ArrR4Bad(nLngLatnT,WTSM)
	call ArrR4Zero(nLngLatnT,VMEANF)
	call ArrR4Zero(nLngLatnT,VFIXM2)
	call ArrR4Zero(nLOSnL,WTP)
	call ArrR4Bad(nLngLatnT,ANMAP)

	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			! Use only valid sources
		LNlast = 0
		LTlast = 0
		ITlast = 0
		sinE = sind(abs(XEV(I)))

		do J=1,NLOS
		    LN = LON(J,I)			! map longitude
		    LT = LAT(J,I)			! map latitude
		    IT = ITIM(J,I)			! integer time
C		    VSN = sinE/RP(J,I)
		    WTSj = VWTij(J,I)
		    Fws = WTSj*FIX(I)
		    if (VMAP(LN,LT,IT) .ge. VLIM .and. FIX(I) .gt. 1.0) Fws = WTSj
C		    if (VMAP(LN,LT,IT) .ge. VLIM) Fws = WTSj
C		    if(FIXM(LN,LT,IT).eq.Bad) FIXM(LN,LT,IT) = 0.
		    FIXM(LN,LT,IT) = FIXM(LN,LT,IT) + Fws
		    WTP(J,I) = WTSj
		    if(WTSM(LN,LT,IT).eq.Bad) WTSM(LN,LT,IT) = 0.
		    WTSM(LN,LT,IT) = WTSM(LN,LT,IT)+WTP(J,I)

		    if(LN.ne.LNlast.or.LT.ne.LTlast.or.IT.ne.ITlast) then
			if(ANMAP(LN,LT,IT).eq.Bad) ANMAP(LN,LT,IT) = 0.
     			ANMAP(LN,LT,IT) = ANMAP(LN,LT,IT)+1.			
		    end if

		    LNlast = LN
		    LTlast = LT
		    ITlast = IT
 		end do
	    end if
	end do

!$omp parallel private(I,J,N,amin,amax,alat,sqonebycos)
!$omp do schedule(static)
	do N = 1, nT
C	    call ArrR4getminmax(nLngLat,FIXM(1,1,N),amin,amax)
C	    if(amax.ne.Bad) then
C		call GridSphere2D(ALng,nLng,nLat,1,FIXM(1,1,N),CONR,3,0.0,90.0)
C	    end if
C	    call ArrR4getminmax(nLngLat,WTSM(1,1,N),amin,amax)
C	    if(amax.ne.Bad) then
C		call GridSphere2D(ALng,nLng,nLat,1,WTSM(1,1,N),CONR,4,0.0,90.0)
C	    end if
	    call ArrR4getminmax(nLngLat,ANMAP(1,1,N),amin,amax)
	    if(amax.ne.Bad) then
		call GridSphere2D(ALng,nLng,nLat,1,ANMAP(1,1,N),CONR,3,0.0,90.0)
C
C  Now modify the latitude map count to approximate the cos latitude filter effect
C
		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,N).ne.Bad) then
			    ANMAP(I,J,N) = ANMAP(I,J,N)*sqonebycos
			end if
		    end do		    
		end do
C
	    end if
	end do
!$omp end do nowait
!$omp end parallel

C	call Timesmooth(nLng,nLat,nT,FIXM, 0,1.0*CONT/aNday,0.,Ztmp)
C	call Timesmooth(nLng,nLat,nT,WTSM, 0,1.0*CONT/aNday,0.,Ztmp)

	call Timesmooth(nLng,nLat,nT,ANMAP,0,1.0*CONT/aNday,0.,Ztmp)

C	do N=1,nT
C	    call ArrR4getminmax(nLngLat,WTSM(1,1,N),amin,amax)
C	    a70 = WTSM(30,8,N)
C	    print *, 'Final map wt min, max, +70 deg.', N, amin, amax, a70
C	    call ArrR4getminmax(nLngLat,ANMAP(1,1,N),amin,amax)
C	    a70 = ANMAP(30,8,N)
C	    print *, 'Final map # min, max, +70 deg.', N, amin, amax, a70 
C	end do

!$omp parallel private(LN,LT,IT)
!$omp do schedule(static)
	do IT=1,nT
	    do LT=1,nLat
		do LN=1,nLng
		    if (FIXM(LN,LT,IT) .gt. 0. .and. WTSM(LN,LT,IT) .gt. 0.) then
			VMEANF(LN,LT,IT) = FIXM(LN,LT,IT)/WTSM(LN,LT,IT) 
			VMA(LN,LT,IT) = VMEANF(LN,LT,IT)*VMAP(LN,LT,IT)
		     else
			ANMAP(LN,LT,IT) = 0.
		    end if
		end do
	    end do
	end do
!$omp end do nowait
!$omp end parallel

C  calculation of the Chi square distribution for each latitude and longitude

	do I=1,NL
	    if (NBADSV(I) .ne. 0) then			! Use only valid sources
	    do J=1,NLOS
		LN = LON(J,I)			! map longitude
		LT = LAT(J,I)			! map latitude
		IT = ITIM(J,I)			! map time
		if (VMA(LN,LT,IT) .ne. Bad) VFIXM2(LN,LT,IT) = VFIXM2(LN,LT,IT)
     &		    +WTP(J,I)*((FIX(I)-VMEANF(LN,LT,IT))**2)
	    end do
	    end if
	end do

C  calculation of the "total" Chi square distribution for all latitudes and longitudes

	CONVER  = 0
	NCONVER = 0
!$omp parallel private(I,J,K)
!$omp do reduction(+:CONVER,NCONVER) schedule(static)
	do K=1,nT
	    do J=1,nLat
		do I=1,nLng
		    if (ANMAP(I,J,K) .lt. ANLIMIT) VMA(I,J,K) = Bad ! data not deconvolved
		    if (VMA(I,J,K) .ne. Bad) then
			CONVER  = CONVER +VFIXM2(I,J,K)
			NCONVER = NCONVER+1
 		    end if
		    VMAP(I,J,K) = VMA(I,J,K)
		end do
	    end do
	end do
!$omp end do nowait
!$omp end parallel

	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,F9.4,A,F7.3,A)') '#V map convergence =',SIGMAP,
     &		', from last =',abs(1-SIGMAP/SIGOLD)*100,'%'

	call Say('MkVMap','I','Info',cStr)

	return
	end