!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.
!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. 2)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=1,nion  ! We allow for C6+ = CVII mainly for testing purposes. 
           jel_ion(2,k)=k
         end do do402 
         jel_ion(2,nion+1:nel)=5 ! Standard upper limit on ionization stages.
!
         call partition_function_read(nel,nion,jel_ion,idiagnostic,                          &
     &                  aname,aweight,xion_en,partfunc)
         partfunc(-1:nion-1,:)=partfunc(-1:nion-1,:)/partfunc(0:nion,:)  
!                              ! The ratios are actually what is needed (Mihalas-1978-113).
         if(idiagnostic .eq. 3) then
             do405 : do k=1,nel      ! k is atomic number as for Mihalas-1978-112.
               jmax=jel_ion(2,k)
               write(*,fmt=910) aname(k),k,aweight(k),                               &
     &                          xion_en(-1:jmax,k)
               write(*,fmt=920) partfunc(-1:jmax-1,k)
  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 .eq. 0) then
          print*,'You have kkel=',kkel
          print*,'kkel must have greater than zero.'
          stop 'saha_solver.f'
      end if
!
      do406 : do kk=1,kkel    ! kk is the calculation index number for the element.
        aweight_1(kk)=aweight(kel(kk))
      end do do406
!
      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
      do408 : do kk=1,kkel
        jmax=jel_ion(2,kel(kk))-1    ! This is no phi function for the highest ionization stage.
!                                    !   So the jmax's for phi are 1 less than for ions.
!        print*,kk,kel(kk),jmax
        phi(-1:jmax,kk)=partfunc(-1:jmax,kel(kk))                        &  ! Mihalas-1978-113
     &    *saha_con*tem_use**(-1.5_np)                                       &
     &    *exp(xion_en(-1:jmax,kel(kk))/(xk_boltzmann_ev*tem_use))
        print*,'phi assigned, etc.'
        print*,kk,kel(kk),saha_con,tem_use
        write(*,fmt=922) partfunc(0:jmax,kel(kk))
        write(*,fmt=922) xion_en(0:jmax,kel(kk))
        write(*,fmt=922) xion_en(0:jmax,kel(kk))                             &
     &                   /(xk_boltzmann_ev*tem_use)
        write(*,fmt=922) phi(0:jmax,kk)
      end do do408
  922 format(1p,8e15.7)
      if(kel(1) .eq. 6) stop
!
      do410 : do kk=1,kkel
        jmin=jel_ion(1,kel(kk))
        jmax=jel_ion(2,kel(kk))
        phi_prod(jmax,kk)=1.0_np
        do430 : do j=jmax-1,jmin,-1
          phi_prod(j,kk)=phi(j,kk)*phi_prod(j+1,kk)   ! Mihalas-1978-114
        end do do430
      end do do410
!
      kkel_1=kkel 
      ih_minus_not_1=ih_minus_not
      if(ixne_initial .eq. 1) then
          xne_h=func_xne_h_analytic_general(xnum_total*phi(0,1))           &
     &          /phi(0,1)
           if(idiagnostic .ge. 1) then
              print*,'xne_h,xnum_total,phi(0,1)'
              write(*,'(1p,3e15.7)') xne_h,xnum_total,phi(0,1)
           end if
          xne=xne_h
        else if(ixne_initial .eq. 2) then
          if(xne .le. 0.0_np) 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/phi(0,1)
         xnum_total=xx_reduced/phi(0,1)
         kkel_1=1    ! This is a pure hydrogen case no matter what the other parameters are.
         ih_minus_not_1=1  ! No H- in this case.
      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 kk=1,kkel_1
          if(kel(kk) .gt. 1  .or.  ih_minus_not_1 .eq. 0) then
              jmin=jel_ion(1,kel(kk))
            else
              jmin=0  ! The case where H- is excluded. 
          end if
          jmax=jel_ion(2,kel(kk))
          xne_j=0.0_np
          do530 : do j=jmin,jmax
            fraction(j,kk)=xne**(jmax-j)*phi_prod(j,kk)   ! Mihalas-1978-114--115
            xne_j=xne_j+real(j,kind=np)*fraction(j,kk)   ! Mihalas-1978-114--115
!                       See MRC-162 for real function.
!            print*,fraction(j,kk),xne_j
          end do do530
          denom=sum(fraction(jmin:jmax,kk))  ! What order sum sum in?
!                                            !   One should add from smallest
!                                            !   to largest for each sign first.
!                                            !   It may be too finicky to worry
!                                            !   About question in this procedure. 
          fraction(j,kk)=fraction(j,kk)/denom
          xne_out=xne_out+xnum(kk)*xne_j/denom
        end do do520
!        print*,xne,xne_out
        relerr=(xne_out-xne)/xne
        if(abs(relerr) .le. xne_converged) 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 
            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. 100) stop 'The iteration failing iter=100'
        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
