C+
C NAME:
C	FixModeltdn_mhdn
C PURPOSE:
C	To determine the change needed for the various lines of sight (LOS)
C	[g = m/<m> ==> dm = <m> dg]
C	The velocity and density (Thomson scattering) analyses are handled the same but differently 
C	than g-levels according to the mode input.
C	For the MHD version of FixModeltdn, more strict model values are retained, documented, and removed by setting NBS to 0.
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call FixModeltdn_mhdn(ID,NL,OBSV,MODV,PW,NBS,NBSF,SIG,FIX,SSIG,RAT)
C INPUTS:
C	ID		integer	1: Fix V model
C				2: Fix G model (Only this mode does anything different. For ipshtd this mode is not used.
C				3: Fix Thomson scattering model
C	NL		integer	# data points
C	OBSV(NL)	real	Observed values
C	MODV(NL)	real	Model values
C	SIG		real	previous convergence criteria
C	PW		real    Power of G
C	NBS(NL)		integer Bad source if <1
C	NBSF		integer The number of sources removed because of a bad model so far
C OUTPUTS:
C	SIG		real	new convergence criteria
C	FIX(NL)		real	Ratio of observed to model values	????
C	SSIG(NL)	real	deviation from mean			????
C	RAT		real	Ratio by which to multiply DEN1AU (G only)
C FUNCTIONS/SUBROUTINES:
C PROCEDURE:
C MODIFICATION HISTORY:
C	NOV, 1995 B. Jackson (STEL,UCSD)
C-
	subroutine FixModeltdn_mhdn(ID,NL,OBSV,MODV,PW,NBS,NBSF,SIG,FIXR,OFIXR,SSIG,RAT)
	real		OBSV(NL),	! Observed G-levels
     &			MODV(NL),	! Modeled G-levels
     &			FIXR(NL),	! Ratio of observed to modeled G-values
     &			OFIXR(NL),	! Old value of observed to modeled G-values. 1.0 if this is the first iteration
     &			SSIG(NL)	! G-source deviation from mean
	integer		NBS(NL)		! Bad source value if < 1
	character	cStr*100,
     &			cID(3)		/'V','G','D'/

	SIGOLD = SIG			! Only used in Say statement
	SIG = 0.0
	RAT = 0.			! Only useful for ID=1 (G level)
	NLL = 0
	PWD = 2.0
	NBSFS = NBSF			! the starting number of model sources removed on this iteration.
	NBSFLO = 0
	do I=1,NL
C	    if(NBS(I).ne.0) then		! Old way
	    if(NBS(I).gt.0) then
	        OBS = OBSV(I)
		AMO = MODV(I)
	        if (ID .eq. 2) then
		  OBS = OBSV(I)*OBSV(I) 	! OK for g-value linearization 7/14/00 B. Jackson
		  AMO = MODV(I)*MODV(I)
		end if
		if(AMO.le.1.0E-4) then
		  NBSFLO = NBSFLO + 1
		  if(ISNaN(AMO)) then
		    if(ID.eq.1) then
		      print *, 'Bad velocity model value. NBSV set to zero.'
		    end if
		    if(ID.eq.2) then
		      print *, 'Bad g-value model value. NBSG set to zero.'
		    end if
		  end if
		end if
		if(AMO.le.1.0E-4.or.ISNaN(AMO)) then
		  NBS(I) = 99
		  NBSF = NBSF + 1
		end if
	        R = OBS/AMO
	        FIXR(I) = 1.0
C	        if(NBS(I).ne.0) then			! Old way
C	        if(NBS(I).gt.0) then
		if(NBS(I).ne.99) then
		  if (R.gt.10.0) R = 10.0		! Fix in 4/27/2011 BVJ
		  if (R.lt.0.1)  R = 0.1
		else
		  R = OFIXR(I)
		end if
		FIXR(I) = R
		OFIXR(I) = R
		NLL = NLL + 1
C		RAT = RAT +  SQRT(1/R)	! Out 7/9/00 B. Jackson
		RAT = RAT + (1./R)	! In 7/9/00 B. Jackson
		R = R-1
		SSIG(I) = R
		SIG = SIG + R*R

C		if(ID.eq.1) then 
C		  write(*,'(A,I2,2I6,I2,2F8.2,F7.3,F8.1,F8.2)') 'ID,I,NLL,OBS,AMO,R,RAT,SIG',ID,I,NLL,NBS(I),OBS,AMO,R,RAT,SIG 
C		end if

C		if(ID.eq.2) then
C		  write(*,'(A,I2,2I6,I2,2F8.4,F7.3,F8.1,F8.2)') 'ID,I,NLL,OBS,AMO,R,RAT,SIG',ID,I,NLL,NBS(I),OBS,AMO,R,RAT,SIG
C		end if

	    end if
	end do
	NBSFN = NBSF - NBSFS
	if(NBSFN.ne.0) then 
	  if(ID.eq.1) then
	    write (*,'(A,I4,A,I4)') '%Fixmodeltdn_mhd-I-Info, # low v models',NBSFLO,' # bad v models',NBSFN-NBSFLO 
	  end if	  
	  if(ID.eq.2) then
	    write (*,'(A,I4,A,I4)') '%Fixmodeltdn_mhd-I-Info, # low g models',NBSFLO,' # bad g models',NBSFN-NBSFLO 
	  end if
	end if

	if(ID.eq.1) RAT =  RAT/NLL
	if(ID.eq.2) RAT = (RAT/NLL)**(1./PW)	! In 7/14/00 B. Jackson
	if(ID.eq.3) RAT =  RAT/NLL
C	if(ID.eq.2) RAT = (RAT/NLL)**2.0	! Out 7/14/00 B. Jackson

C	if(NLL.eq.0) print *, 'RAT, SIG, NLL, NL, PW', RAT, SIG, NLL, NL, PW

	SIG = sqrt(SIG/NLL)
	do I=1,NL
	    if(NBS(I).gt.0) SSIG(I) = SSIG(I)/SIG
	end do
	write (cStr,'(A,F6.3,A,F7.3,A)') cID(ID)//' convergence = ',SIG,
     &		', from last =',abs(1-SIG/SIGOLD)*100,'%'
	if (ID .eq. 1) write (cStr(itrim(cStr)+1:),'(A,F7.3)') '#Summed model to observed velocity ratio ',RAT
	if (ID .ge. 2) write (cStr(itrim(cStr)+1:),'(A,F7.3)') '#Summed model to observed density ratio ',RAT
	call Say('FixModeltdn','I','Info',cStr)
	return
	end