C+
C NAME:
C	MkDMaptdn
C PURPOSE:
C	Change the density map from the line of sight (L.O.S.) values of G-level 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 MkDMaptdn(LON,LAT,ITIM,GWTij,PW,FIX,NL,NLOS,NLOSP1,
C		nLng,nLat,nT,DMAP,ANLIMIT,GSSIG,SIGB,NS,NBADSG,SIGMAP,
C		ANMAP,DMA,FIXM,WTSM,DMEANF,DFIXM2,WTP,Alng,ConstL,CONR,CONT,aNday,Ztmp)
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 times at reference surface
C	GWTij(NLOS,NL)		real		m Weights of each point along the L.O.S.
C	PW			real		Power of G to N
C	FIX(NL)			real		Ratio observed/model G for total L.O.S.
C	NL			integer		Number of G-level data points
C	NLOS			integer		Number of L.O.S. distance segments
C	NLOSP1			integer		Number of L.O.S. distance segments + 1
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	DMAP(nLng,nLat,nT)	real		Density map
C	SIGMAP			real		Density map convergence criteria
C	ANLIMIT			real		Number of points constituting deconvolution  
C	GSSIG(NL)		real		Deviation of source from mean in sigma
C	SIGB			real		G-level sigma cut-off
C	NS			integer		1 = Activate sigma cut-off, 0 = do not
C	NBADSG(NL)		integer		Bad source G-value 0 - bad, 1 - good
C OUTPUTS:
C	DMAP(nLng,nLat,nT)	real		Density map
C	BADD			real		Value of map where there is no projection
C	NBADSG(NL)		integer		Bad source G-value 0 - bad, 1 - good
C	SIGMAP			real		Density map convergence criteria
C
C SCRATCH ARRAYS:
C	ANMAP(nLng,nLat,nT)	real
C	DMA(nLng,nLat,nT)	real
C	FIXM(nLng,nLat,nT)	real
C	WTSM(nLng,nLat,nT)	real
C	DMEANF(nLng,nLat,nT)	real
C	DFIXM2(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 MkDMaptdn(LON,LAT,ITIM,GWTij,PW,FIX,NL,NLOS,NLOSP1,
     &		nLng,nLat,nT,DMAP,ANLIMIT,GSSIG,SIGB,NS,NBADSG,SIGMAP,
     &		ANMAP,DMA,FIXM,WTSM,DMEANF,DFIXM2,WTP,Alng,ConstL,CONR,CONT,aNday,Ztmp)

	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
     &			NBADSG(NL)		! Bad source velocity 0 - bad, 1 - good

	real		GWTij(NLOS,NL),		! m weights of each point along the L.O.S.
     &			FIX(NL),		! Ratio observed/model G for total L.O.S.
     &			GSSIG(NL),		! Deviation of source from mean in sigma
     &			DMAP(nLng,nLat,nT)	! Density map

C	Scratch arrays

	real		ANMAP(nLng,nLat,nT),	! Number of points per bin
     &			DMA(nLng,nLat,nT),
     &			FIXM(nLng,nLat,nT),		
     &			WTSM(nLng,nLat,nT),	
     &			DMEANF(nLng,nLat,nT),	
     &			DFIXM2(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,ANMAP)
	call ArrR4Bad(nLngLatnT,DMA)
	call ArrR4Zero(nLngLatnT,FIXM)
	call ArrR4Bad(nLngLatnT,WTSM)
	call ArrR4Zero(nLngLatnT,DMEANF)
	call ArrR4Zero(nLngLatnT,DFIXM2)
	call ArrR4Zero(nLOSnL,WTP)

	NSS = 0
	do I=1,NL
	    if (abs(GSSIG(I)) .gt. SIGB) then		! Source above limit
		NSS = NSS+1
		if (NS .eq. 1) NBADSG(I) = 0		! Mark source as bad
	    end if

	    if (NBADSG(I) .ne. 0) then			! Only use valid sources
		LNlast = 0
		LTlast = 0
		ITlast = 0
		do J=1,NLOS
		    LN = LON(J,I)			! map longitude
		    LT = LAT(J,I)			! map latitude
                    IT = ITIM(J,I)                      ! integer time
		    WTSj = GWTij(J,I)
C		    FWS = WTSj*(FIX(I)**(1.0/(2.0*PW)))
C		    if(FIXM(LN,LT,IT).eq.Bad) FIXM(LN,LT,IT) = 0.
		    FIXM(LN,LT,IT) = FIXM(LN,LT,IT)+WTSj*FIX(I)
		    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)
C	print *, LN, LT, IT, FIXM(LN,LT,IT), WTSj, FIX(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,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
			DMEANF(LN,LT,IT) = FIXM(LN,LT,IT)/WTSM(LN,LT,IT)
			DMA(LN,LT,IT) = DMEANF(LN,LT,IT)*DMAP(LN,LT,IT)
C	print *, LN,LT,IT,FIXM(LN,LT,IT),WTSM(LN,LT,IT),DMEANF(LN,LT,IT),DMAP(LN,LT,IT),DMA(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 (NBADSG(I) .ne. 0) then			! Only use 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 (DMA(LN,LT,IT) .ne. Bad) DFIXM2(LN,LT,IT)=DFIXM2(LN,LT,IT)
     &		    +WTP(J,I)*((FIX(I)-DMEANF(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 schedule(static) reduction(+:CONVER,NCONVER)
	do K=1,nT
	    do J=1,nLat
		do I=1,nLng
		    if (ANMAP(I,J,K) .lt. ANLIMIT) DMA(I,J,K) = Bad ! data not deconvolved
		    if (DMA(I,J,K) .ne. Bad) then
			CONVER = CONVER+DFIXM2(I,J,K)
			NCONVER = NCONVER+1
		    end if
		    DMAP(I,J,K) = DMA(I,J,K)
C	if(DMAP(I,J,K).gt.1.0e7.and.DMAP(I,J,K).lt.1.0e22) print *, I, J, DMAP(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)') 'density 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,F8.4,A,F7.3,A)')
     &		'#D map convergence =',SIGMAP,', from last =',abs(1-SIGMAP/SIGOLD)*100,'%'

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

	return
	end