C+
C NAME:
C	MkGModeltd_in
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:
C	call MkGModeltd_in(XCbe,XCtbeg,XCtend,XCtpr,XLON,XLAT,RP,
C		WTS2,DFAC,PW,DEN1AU,XCtpr_in,XLON_in,XLAT_in,RP_in,WTS1_in,DFAC_in,
C		NL,N_IN,NLOS,NLOSP1,DMAP,
C		nLng,nLat,nT,RRS,GM2,GWTij,GWTij_in,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*8		Projected time
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 FUNCTIONS/SUBROUTINES:
C	Get3DTval
C PROCEDURE:
C MODIFICATION HISTORY:
C	NOV, 1995 B. Jackson (STEL,UCSD)
C-
	subroutine MkGModeltd_in(XCbe,XCtbeg,XCtend,XCtpr,XLON,XLAT,RP,
     &		WTS2,DFAC,PW,DEN1AU,XCtpr_in,XLON_in,XLAT_in,RP_in,WTS1_in,DFAC_in,
     &		NL,N_IN,NLOS,NLOSP1,DMAP,
     &		nLng,nLat,nT,RRS,GM2,GWTij,GWTij_in,GMWTij,GMWTi)
     
	real		XLON(NLOSP1,NL),	! Projected Carr. rot. value of point on L.O.S.  
     &			XCbe(2,nT),		! Beginning and ending Carr. variables
     &			XLAT(NLOSP1,NL),	! Heliographic lat. (in deg.) point on L.O.S.
     &			RP  (NLOSP1,NL),	! Distance above Sun of point on L.O.S.
     &			WTS2(NLOS  ,NL),	! Weights of each point along the L.O.S.
     &			DFAC(NLOSP1,NL),	! Density factors for each L.O.S. point
     &			XLON_in(N_IN),		! Projected Carr. rot. value of in-situ point
     &			XLAT_in(N_IN),		! Projected heliographic lat. (in deg.) of in-situ point
     &			RP_in  (N_IN),		! Distance above Sun of in-situ point
     &			WTS1_in(N_IN),		! Weights of each in-situ point
     &			DFAC_in(1,N_IN),	! Density factors for each in-situ point
     &			DMAP(nLng,nLat,nT),	! Density map
     &			GM2(NL+N_IN),		! Model G-levels^2 for a given source
     &			GWTij(NLOS ,NL),	! L.O.S. weights
     &			GWTij_in(N_IN)		! in-situ point weights

	real*8		XCtpr(NLOSP1,NL)	! Projected time value on surface
	real*8		XCtpr_in(N_IN)		! Projected in-situ point time value on surface

	real*8		XCtbeg,XCtend

	real		GMWTij(NLOS,NL),	! Scratch arrays
     &			GMWTi(NL+N_IN)

	BADD = BadR4()
	R1AU = 1.
	NDENj = 0
	do I=1,NL+N_IN
	    GMWTi(I) = 0
            if(I.le.NL) then
              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) stop 'Bad density in MkGModeltd_in'
		DENj = DENj*DFAC(J,I)*((RRS/RP(J,I))**2)
C	if(DENj.le.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.) then
		  NDENj = NDENj + 1
		  DENj = 0.001
	        end if
C		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)
C		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)
C		GMWTij(J,I) = (DENMj**PW)*WTS2(J,I)*cosd(XLAT(J,I))
		GMWTi(I)    = GMWTi(I) + GMWTij(J,I)
	      end do
	    else
              call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,DMAP,1,
     &                  XLON_in(I-NL),XLAT_in(I-NL),XCtpr_in(I-NL),DENj)
	      if (DENj .eq. BADD) stop 'Bad density in MkGModeltd_in'
	      if(DENj.le.0.) then
		NDENj = NDENj + 1
		DENj = 0.001
	      end if
	      DENj = DENj*DFAC_in(1,I-NL)*((RRS/RP_in(I-NL))**2)
	      GWTij_in(I-NL) = (DENj**(2.0*PW))*WTS1_in(I-NL)
	        
	      DENMj = DEN1AU*((R1AU/RP_in(I-NL))**2)
	      GMWTij1 = (DENMj**(2.0*PW))*WTS1_in(I-NL)
C	      GMWTij1 = (DENMj**PW)*WTS1_in(I-NL)*cosd(XLAT_in(I-NL))
	      GMWTi(I)  = GMWTi(I) + GMWTij1

	    end if
C	if(GMWTi(I).eq.0) print *, I, DENMj, PW, WTS2(1,I), WTS2(NLOS,I)
	end do

	if(NDENj.ne.0) print *, NDENj, ' L.O.S. density values in MkGmodeltdn are below zero'
	do I=1,NL+N_IN
	    GM2(I) = 0.
            if(I.le.NL) then
            
C	Scintillation index line of sight integration for both the modeled
C	density weighting (GWTij) and the mean density weighting (GMWTij).

	      do J=1,NLOS
		GWTij(J,I)   = GWTij(J,I)/GMWTi(I)
		GM2(I) 	     = GM2(I) + GWTij(J,I)
	      end do
	    else
	      GWTij_in(I-NL) = GWTij_in(I-NL)/GMWTi(I)
	      GM2(I) 	     = GM2(I) + GWTij_in(I-NL)
	    end if
	end do

	do I=1,NL+N_IN
          if(I.le.NL) then	  
	    do J=1,NLOS
	      GWTij(J,I)   = GWTij(J,I)/GM2(I)
	    end do
	  else
	    GWTij_in(I-NL) = GWTij_in(I-NL)/GM2(I)
	  end if
	end do

	return
	end