C+
C NAME:
C	MkGModeltdn
C PURPOSE:
C	Make model line of sight G-levels from a density map at a
C	given height. Following, a latex file containing the scintillation
C	index value m and the weighting function.
C	Line-of-sight integration
C
C	Scintillation index \\
C	   $$ m^2 = \int dz W(z) $$ 
C	Weighting function \\
C	   $$ W(z) = (\Delta N_e(z))^2 \int dq 
C	          \sin^2(\frac{q^2\lambda z}{4\pi})
C	          \exp(-\frac{\theta_o^2 q^2 z^2}{2})
C	          q^{-3} $$
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
	subroutine MkGModeltdn(XCbe,XCtbeg,XCtend,XCtpr,XLON,XLAT,RP,
     &		WTS2,DFAC,PW,DEN1AU,NL,NLOS,NLOSP1,DMAP,
     &		nLng,nLat,nT,RRS,GM2,GWTij,GMWTij,GMWTi)
C INPUTS:
C       XCbe(2,nT)      real    	Beginning and ending Carrington variables
C	XCtbeg		real		Beginning time interval
C	XCtend		real		Ending time variable
C	XCtpr(NLOSP1,NL)real		Projected time (in XC)
C	XLON(NLOSP1,NL)	integer		Projected Carr. rot. value of point on L.O.S.
C	XLAT(NLOSP1,NL)	integer		Heliographic lat. (in deg.) point on L.O.S.
C	RP  (NLOSP1,NL)	real		Distance above Sun of point on L.O.S.
C	WTS2(NLOS  ,NL)	real		Weights of each point along the L.O.S.
C	DFAC(NLOSP1,NL)	real		Density factors for each L.O.S. point
C	PW		real		Density power
C	DEN1AU		real		Average density at 1 AU
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	DMAP(nLng,nLat,nT)	real	Density map
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	dRR		real		Interval between height maps
C	RRS		real		Height of deconvolution surfaces
C OUTPUTS:
C	GM2(NL)		real		Model G-levels for a given source
C	GWTij(NLOS ,NL)	real		Value weights
C SCRATCH ARRAYS:
C	GMWTij(NLOS,NL) real
C	GMWTi(NL)	real
C CALLS:
C	Get3DTval
C PROCEDURE:
C MODIFICATION HISTORY:
C	NOV-1995, B. Jackson (STEL,UCSD)
C-
	real		XLON(NLOSP1,NL)		! Projected Carr. rot. value of point on L.O.S.  
	real		XCbe(2,nT)		! Beginning and ending Carr. variables
	real		XLAT(NLOSP1,NL)		! Heliographic lat. (in deg.) point on L.O.S.
	real		RP  (NLOSP1,NL)		! Distance above Sun of point on L.O.S.
	real		WTS2(NLOS  ,NL)		! Weights of each point along the L.O.S.
	real		DFAC(NLOSP1,NL)		! Density factors for each L.O.S. point
	real		DMAP(nLng,nLat,nT)	! Density map
	real		GM2(NL)			! Model G-levels^2 for a given source
	real		GWTij(NLOS ,NL)

	real		GMWTij(NLOS,NL)		! Scratch arrays
	real		GMWTi(NL)

	double precision XCtbeg
	double precision XCtend
	double precision XCtpr(NLOSP1,NL)	! Projected time value on surface

	BADD = BadR4()

	R1AU  = 1.0
	NDENj = 0

!$omp parallel private(I,J,GMWTiT,GM2T,DENj,DENMj)
!$omp do reduction(+:NDENj) schedule(static)
	do I=1,NL
	    GMWTiT = 0
            do J=1,NLOS
                call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,DMAP,1,
     &                  XLON(J,I),XLAT(J,I),XCtpr(J,I),DENj)
		if (DENj .eq. BADD) call stop_here('Bad density in MkGModeltdn')
		DENj = DENj*DFAC(J,I)*((RRS/RP(J,I))**2)
		!if (DENj .le. 0.0) print *, 'Bad density in MkGmodeltdn',J,I,XLON(J,I),XLAT(J,I),XCtpr(J,I),DENj,DFAC(J,I),RRS,RP(J,I),WTS2(J,I),PW
		if (DENj .le. 0.0) NDENj = NDENj + 1
		if (DENj .le. 0.0) DENj = 0.001
		!GWTij(J,I) = (DENj**(2.0*PW))*WTS2(J,I)/(RP(J,I)**2) ! Asai's
		GWTij(J,I)  = (DENj**(2.0*PW))*WTS2(J,I)
		!GWTij(J,I) = (DENj**PW)*WTS2(J,I)*cosd(XLAT(J,I))

		DENMj = DEN1AU*((R1AU/RP(J,I))**2)
		GMWTij(J,I) = (DENMj**(2.0*PW))*WTS2(J,I)
		!GMWTij(J,I) = (DENMj**PW)*WTS2(J,I)*cosd(XLAT(J,I))
		GMWTiT = GMWTiT+GMWTij(J,I)
	    end do
	    GMWTi(I) = GMWTiT
	    if (GMWTiT .eq. 0) print *, I, DENMj, PW, WTS2(1,I), WTS2(NLOS,I)
	end do
!$omp end do

!$omp master
	if (NDENj .ne. 0) print *, NDENj, ' L.O.S. density values in MkGmodeltdn are below zero'
!$omp end master

	! Scintillation index line of sight integration for both the modeled
	! density weighting (GWTij) and the mean density weighting (GMWTij).

!$omp do schedule(static)
	do I=1,NL

	    !call ArrR4DivideByConstant(NLOS,GWTij(1,I),GMWTi(I),GWTij(1,I))
	    !GM2(I) = ArrR4Total(NLOS,GWTij(1,I))
	    !call ArrR4DivideByConstant(NLOS,GWTij(1,I),GM2(I),GWTij(1,I))

	    do J=1,NLOS
		GWTij(J,I) = GWTij(J,I)/GMWTi(I)
	    end do

	    GM2T = 0.0

	    do J=1,NLOS
		GM2T = GM2T+GWTij(J,I)
	    end do

	    do J=1,NLOS
		GWTij(J,I) = GWTij(J,I)/GM2T
	    end do

	    GM2(I) = GM2T

	end do
!$omp end do nowait
!$omp end parallel

	return
	end