!<pre>
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
! D.J. Jeffery, 2020jan06
!      include 'saha_solver.f'  ! The Saha solver.
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
      include 'xne_acceleration.f'   ! Electron density solution iteration acceleration.
      include 'partition_function_read.f'  ! Reads in the partition function table.
      include 'saha_solver_functions.f'  ! Includes some saha_solver.f functions.
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
      subroutine saha_solver(kkel,kel,imassnumfrac,xel,                     &
     &                       rho,tem,r_dilution,                            &
     &                       iacceleration,idiagnostic,ih_minus_not,        &
     &                       xne_converged,                                 &
     &                       ixne_initial,x_reduced,xx_reduced,             &
     &                       xne,xne_h,aname,aweight,xne_atom_ratio,        &
     &                       fraction,iter)
      include 'saha_solver_declarations.f'
!
      if(ipartfunc .eq. 0  .or.  idiagnostic .eq. 3)then
         ipartfunc=1
         jel_ion(1,1)=-1      ! We allow for the H- ion.
         jel_ion(1,2:nel)=0   ! All other elements start with the neutral atom.
         do402 : do k_atom=1,nion  ! We allow for C6+ = CVII mainly for testing purposes. 
           jel_ion(2,k_atom)=k_atom
         end do do402 
         jel_ion(2,nion+1:nel)=nion-1 ! Standard upper limit on ionization stages.
!
         call partition_function_read(nel,nion,jel_ion,idiagnostic,                          &
     &                  aname,aweight,xion_en,partfunc)
!          This code uses the partition functions themselves, not their ratios
!          as in Mihalas-1978-113--114.
         if(idiagnostic .eq. 3) then
             do405 : do k_atom=1,nel      ! k_atom is atomic number as for Mihalas-1978-112.
               jmax=jel_ion(2,k_atom)
               write(*,fmt=910) aname(k_atom),k_atom,aweight(k_atom),                         &
     &                          xion_en(-1:jmax,k_atom)
               write(*,fmt=920) partfunc(-1:jmax-1,k_atom)
  910          format(1x,a2,i4,f6.0,18x,8f9.3)
  920          format(31x,8f9.3)
             end do do405
             return    ! See MRC-78 for return statement.
         end if
      end if
!
      if(iacceleration .lt. 0  .or.  iacceleration .gt. 3) then
          print*,'You have given iacceleration=',iacceleration
          print*,'Allowed iacceleration values 0,1,2,3'
          stop 'saha_solver.f'
      end if
      if(kkel .le. 0) then
          print*,'You have kkel=',kkel
          print*,'kkel must be greater than zero.'
          stop 'saha_solver.f'
      end if
!
      do410 : do k_index=1,kkel    
        aweight_1(k_index)=aweight(kel(k_index))
      end do do410
!
      sumel=sum(xel(:))   ! For normalization if that was not set initially.
      if(imassnumfrac .eq. 1) then    ! For mass fractions as input parameters.
          xnum(:)=rho*(xel(:)/sumel)/(aweight_1(:)*amu)
          xnum_total=sum(xnum(:))
        else                          ! For number fractions as input parameters.
          xnum_total=rho/sum( (xel(:)/sumel)*aweight_1(:)*amu )
          xnum(:)=xnum_total*(xel(:)/sumel)
      end if
!
      if(r_dilution .le. 0.0_np) then
         tem_use=tem
        else
         w_dilution=0.5_np*(1.0_np-sqrt(1.0_np-r_dilution**(-2)))  ! Mihalas-1978-120
         tem_use=tem*w_dilution**0.25_np
      end if
!
      tmp=func_phi(partfunc(0,1),partfunc(1,1),                          &
     &    1.0_np,xion_en,tem_use) 
! With ene=1.0_np, this is exactly ithe phi-tilde function (Mihalas-1978-113)
! tmp = 1/n_fiducial_electron_density_sort_of.
      if(ixne_initial .eq. 1) then
          xne_h=func_xne_h_analytic_exact(xnum_total*tmp)/tmp
           if(idiagnostic .ge. 1) then
              print*,'xne_h,xnum_total,tmp'
              write(*,'(1p,3e15.7)') xne_h,xnum_total,tmp
           end if
          xne=xne_h                      ! The hydrogen approximation.
!                                        !   I suspect it is usually an 
!                                        !   underestimate, but not by more
!                                        !   than a factor of 100 for large ionization
!                                        !   since atoms can at most yield about 100 electons each.
!                                        !   For low ionization, since H has a relatively
!                                        !   1st ionization, it might also be low often.
!                                        !   But it is exactly right in the limit of
!                                        !   pure hydrogen (with no H-).
        else if(ixne_initial .eq. 2) then
          if(xne .lt. -xprecision) then
              print*,'For ixne_initial=2 you need xne > = 0'
              print*,'Your input value is xne= ',xne
              stop 'saha_solver.f'
          end if
        else 
         xne=x_reduced/tmp
         xnum_total=xx_reduced/tmp
         xnum(1)=xnum_total
         if(kkel .ne. 1  .or. kel(1) .ne. 1                                &
     &    .or. ih_minus_not .ne. 1                                         &
     &    .or. abs(kel(1)-1.0_np) .gt. xprecision) then
            print*,'This must be a pure hydrogen case'
            print*,'(no H-) since it is for testing.'
            print*,'So one of the following has a wrong value.'
           print*,'kkel,kel(1),ih_minus_not,kel(1)'
           print*,kkel,kel(1),ih_minus_not,kel(1)
           stop 'saha_solver.f'
         end if
      end if
!
      if(idiagnostic .ge. 1) then
        write(*,'(a)') 'The Saha solution for:'
        write(*,fmt=930) iacceleration,rho,tem_use,xne
  930   format('iacceleration=',i5,1p,' rho=',e15.7,                        &
     &    ' T=',e15.7,' xne_initial=',e15.7)
        write(*,*) 'kel(:)',kel(:)
        write(*,*) 'xel(:)',xel(:)
        write(*,fmt=940)
      end if
!
      iter=0
      pair(:,:)=0.0_np
      do510 : do
        iter=iter+1
        xne_out=0.0_np
        do520 : do k_index=1,kkel
          k_atom=kel(k_index)
          if(k_atom .gt. 1  .or.  ih_minus_not .eq. 0) then
              jmin=jel_ion(1,k_atom)
            else
              jmin=0  ! The case where H- is excluded. 
          end if
          jmax=jel_ion(2,k_atom)
          do530 : do j=jmin,jmax-1         
             phi_j=func_phi(partfunc(j,k_atom),partfunc(j+1,k_atom),         &
     &                      xne,xion_en(j,k_atom),tem_use)
!               The phi_j functions always grow with j (I think).
!               The higher the ionization stage, the higher the ionization energy.
!               So if very low ionization, phi_j grows from the start and jfid=jmin.
!               If the very high ionization, then jmax becomes jfid, since N_j/N_j+1
!               is infinite.
             if(phi_j .ge. 1.0_np) exit do530
          end do do530 
          jfid=j     ! MRC-61 confirms this should be j of exit or j=(jmax-1)+1=jmax for normal leave.
!
          fraction(jfid,k_index)=1.0_np   ! The fraction is normalized below.
          xne_j=jfid                      ! The unnormalzed contribution of j=jfid is jfid. 
          product=1.0_np/partfunc(jfid,k_atom)
          do540 : do j=jfid-1,jmin,-1     ! No execution if jfid=jmin
            product=product                                                  &
     &        *func_phi_r(xne,xion_en(j,k_atom),tem_use)
            if(product .lt. xprecision) then
                fraction(jmin:j,k_index)=0.0_np    ! These cannot contribute relative to 1
!                                                  !   and underflows are uncertain (MRC-225).
                exit do540
            end if 
            fraction(j,k_index)=partfunc(j,k_atom)*product
            xne_j=xne_j+real(j,kind=np)*fraction(j,k_index)   ! Mihalas-1978-114--115
          end do do540
!
          product=1.0_np/partfunc(jfid,k_atom)
          do550 : do j=jfid+1,jmax,1       ! No execution is jfid=jmax
            product=product                                                   &
     &        *func_phi_r_inverse(xne,xion_en(j,k_atom),tem_use)
            if(product .lt. xprecision) then
                fraction(j:jmax,k_index)=0.0_np    ! These cannot contribute relative to 1
!                                                  !   and underflows are uncertain (MRC-225).
                exit do550
            end if
            fraction(j,k_index)=partfunc(j,k_atom)*product
            xne_j=xne_j+real(j,kind=np)*fraction(j,k_index)   ! Mihalas-1978-114--115.
          end do do550 
          denom=sum(fraction(jmin:jmax,k_index))  ! What order of summation in sum?
!                                                 !   One should add from smallest
!                                                 !   to largest for each sign first.
!                                                 !   It may be too finicky to worry
!                                                 !   About question in this procedure. 
          if(k_atom .eq. 6) then
              write(*,"('frac',7e15.7)") fraction(jmin:jmax,k_index)
          end if
          fraction(:,k_index)=fraction(:,k_index)/denom
          if(k_atom .eq. 6) then
              write(*,"('frac',7e15.7)") fraction(jmin:jmax,k_index)
          end if
          xne_out=xne_out+xnum(k_index)*xne_j/denom
        end do do520
!
!        print*,xne,xne_out
        relerr=(xne_out-xne)/xne
        if(abs(relerr) .le. xne_converged  .or.                           &
     &     abs(sum(fraction(0,1:kkel))/real(kkel,kind=np)-1.0_np)         &
     &     .le. xprecision) then
            xne_atom_ratio=xne_out/xnum_total
            if(idiagnostic .ge. 1) then
                write(*,fmt=942) iter,xne,xne_out,0.0_np,relerr,          &
     &                                xne_atom_ratio 
                if(kel(1) .eq. 6) then
                   print*,'fraction(3,6),denom'
                   print*,fraction(3,6),denom
                   write(*,"('frac',7e15.7)")                             &
     &               fraction(0:6,1)
                end if 
            end if
            xne=xne_out   ! Convergence achieved.
            exit do510
        end if
        pair(:,1)=pair(:,2)
        pair(:,2)=pair(:,3)
        pair(1,3)=xne
        pair(2,3)=xne_out
        call xne_acceleration(iacceleration,iter,                         &
     &                 pair,xne,xne_out,                                  &
     &                 xne_min,xne_max,xne_mid,xne_L,xne_q,xne_accel)
        if(idiagnostic .eq. 2) then
           write(*,fmt=942) iter,xne,xne_out,xne_accel,relerr,            &
     &                      xne_accel/xnum_total
           if(iter .eq. 50) then
               print*,'The iteration failing iter=500'
               return
           end if 
        end if
        xne=xne_accel   ! This must come after the idiagnostic printout.
      end do do510
!
  940 format(' iter',12x,'xne',8x,'xne_out',6x,'xne_accel',9x,'relerr',   &
     &       5x,'Ionization')
  942 format(i5,1p,5e15.7)
!
      end subroutine saha_solver
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12