!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)
      use precision_set
      use contants
      include 'module_implicit.f'
!     Paramenters have to precede their first use.
      integer, parameter :: nel=38   ! Only to A=38 Sr is allowed at present. 
      integer, parameter :: nion=5   ! Ionization stages, -1,0,1,2,5 are calculated---except for
!                                    !   for those that do not have them, of course.
!                                    !   Stage 5 is treated as having infinite ionization energy
!                                    !   represented by 9999.999 which is just a placeholder.
!                                    !   The stage -1 is calculated for only H:  i.e., H-.
!                                    !   The ih_minus_not = 1 turns off the H-
!                                    !   which is useful for comparison to the
!                                    !   the analytic pure H xne_h which does not use stage H-.
      include 'function_type_declaration.f' ! A file that declares the function types.
      character, intent(inout) :: aname(nel)*2 ! Output:  Atomic name symbols.
      integer :: i,j,k,kk,L,m,n ! Usually indexes for doloops.  However, following 
!                               !   Mihalas-1978-112, j is also ionization stage and k is element sometimes.
      integer, intent(in) :: iacceleration   ! The parameter controling the acceleration of the iteration
!                                            !  for electron density.
!                                            !  0 for none, 1 for midpoint, 2 for linear, 3 for quadratic. 
!                                            !  See discussion of acceleration in 
!                                            !  file xne_acceleration.f. 
!                                            !  iacceleration=3 recommeded circa 2020jan06.
!          https://stackoverflow.com/questions/37723973/fortran-2003-2008-elegant-default-arguments
!          According to this site, f95 does NOT allow default values for dummy arguments!!!
!          MRC-84 says optional attribute allowed.  But how do you tell what to do when
!          you do not supply a value?
!          The website above has no elegant, universal workaround to suggest.
      integer, intent(in) :: idiagnostic   ! Input:  0 for no diagnostics of the saha solver.  
!                                          !   1 for diagnostics for saha_solver.f 
!                                          !   2 for diagnostics for saha_solver.f 
!                                          !     and 'partition_function_read.f' 
      integer, intent(in) :: ih_minus_not  ! Input:  Set to 0 for H- used.
!                                           !   Set to 1 to turn off H- which is useful for comparison to the
!                                           !   the analytic solution for pure H which does not have 
!                                           !   the H- stage. 
      integer :: ih_minus_not_1             ! Equals ih_minus_not, unless ixne_initial=3 in which case it is 1.
      integer, intent(in) :: imassnumfrac  ! Input:  1 for mass fraction, 2 for number fraction.
      integer, save :: ipartfunc = 0 ! 0 initially causes a read-in of the partition function table.  
!                                    !   1 after the partition function has been read.
      integer, intent(out) :: iter ! Output:  Count of iterations to convergence.
      integer, intent(in) :: ixne_initial   ! Input:  For the initial electron density.
!                                           !   1 the anlytical hydrogen formula in which case the
!                                           !     input xne is ignored.
!                                           !   2 uses the input xne as the initial value which is a good
!                                           !     idea if you are doing a for series of states where each
!                                           !     succeeding state is close to the last one, and so
!                                           !     the converged last state xne is a good initial xne for
!                                           !     the next state xne.
!                                           !   3 is for testing with pure H (with no H-) where one supplies
!                                           !     reduced ion density xx_reduced
!                                           !     and an initial reduced electron density
!                                           !     see if on converges to the exact result.
!                                           !     The initial xne is ignored in this case.
      integer :: jmin,jmax   ! Minimum and maximum ionization stages.  Which vary with the element.
      integer, intent(in) :: kkel ! Input:  Number of elements to be calculated for.
!                                 !   Except for ixne_initial=3 in which case only H is solved for without H-.
      integer :: kkel_1           ! Equals kkel, unless ixne_initial=3 in which case it is 1.
      integer, intent(in) :: kel(kkel) ! Input:  Elements by atomic number to be calculated for.
      real (kind=np), intent(inout) :: aweight(nel)  ! Output:  Atomic weights from the read-in 
                                                         !   partition function table.
!                                                        !   They are rounded off to whole numbers from
!                                                        !   terrestrial averages (I think).  The exact
!                                                        !   isotope weighted atomic weights hard 
!                                                        !   to know a priori for general astrophysical systems.
      real (kind=np) :: aweight_1(kkel) ! Atomic weights for the input elements. 
      real (kind=np) :: denom  ! A denominator value used briefly.
      real (kind=np), intent(out) :: fraction(-1:nion,kkel)  ! Output:  Ionization fractions.
      real (kind=np), save :: partfunc(-1:nion,nel) ! Initially the partition functions which are currently 
!                                                   !   just the ground state statistical weights from 
!                                                   !   the read-in table.  But these are then converted
!                                                   !   to the ratios of the partition functions since
!                                                   !   those are actually used.  See Mihalas-1978-113. 
      real (kind=np) :: pair(2,3)               ! Used for the linear acceleration of the iteration:
!                                               !   slope=(pair(2,2)-pair(2,1))/(pair(1,2)-pair(1,1)) .
      real (kind=np) :: phi(-1:nion,kkel)        ! Mihalas-1978-113's phi-tilde function. 
      real (kind=np) :: phi_prod(-1:nion,kkel)   ! Mihalas-1978-114's phi-tilde function products. 
      real (kind=np), intent(in) :: r_dilution  ! Input:  If r_dilution < 0, does nothing.
!                                               !   Otherwise tem(used)=tem(input)*W**(1/4),
!                                               !   where w=(1/2)(1-sqrt(1-r_dilution**(-2)))
!                                               !   is the dilution factor (Mihalas-1978-120)
!                                               !   for photospheric approximation for temperature.
!                                               !   r_dilution = radius/radius_photosphere.
      real (kind=np) :: relerr  ! Relative change-in or relative error-in. 
      real (kind=np), intent(in) :: rho         ! Input:  Matter density in g/cm**3. 
      real (kind=np) :: sumel                   ! Sum of the element abundances.
      real (kind=np), intent(in) :: tem         ! Input:  Temperature in kelvin.
      real (kind=np) :: tem_use                 ! Temperature in kelvin used in calculation.
!                                               !   It is tem if r_dilution < 0 and the dilution
!                                               !   temperature if r_dilution > 0.
      real (kind=np) :: w_dilution  ! If r_dilution > 0, it is the diltuion factor.
      real (kind=np), intent(out) :: xne_atom_ratio  ! Output:  Ratio (electron density/atom density).
      real (kind=np), intent(in) :: xne_converged ! Input:  The convergence criterion.
!                                                 ! The iteration stops when
!                                                 !   relative changes in electron density from
!                                                 !   xne to xne_out are less/equal to
!                                                 !   xne_converged.  Probaly xne_converged=1/100 
!                                                 !   is good enough for
!                                                 !   actual use of the saha solver and not just testing.
!                                                 !   We know errors in the converged value to 1 percent
!                                                 !   are not the limiting errors.
      real (kind=np), optional, intent(in) :: x_reduced,xx_reduced ! Input.  Values used for
!                                                        !   testing in the ixne_initial=3 case.
      real (kind=np), intent(in) :: xel(kkel)   ! Input:  Mass/number fraction elements.
      real (kind=np), save :: xion_en(-1:nion,nel) ! Ionization energies in eV from the read-in 
!                                                  !   partition function table.
      real (kind=np) :: xnum(kkel)  ! Number density of elements.
      real (kind=np) :: xnum_total  ! Total number density of elements collectively
      real (kind=np), intent(inout) :: xne  ! Input/Output:  Electron density, input is inital value 
!                                           ! for iteration of ixne_initial=2 and output is the 
!                                           ! converged value
      real (kind=np) :: xne_accel ! The accelerated electron density for the next iteration
!                                 !  of the Saha electron density equation.
      real (kind=np), intent(out) :: xne_h  ! Output:  Electron density for the pure hydrogen (without H-)
!                                           !   is used as the initial value for the iteration if xne < 0.
!                                           !   The analytic solution is exact (within numerical error)
!                                           !   for partition functions used.
      real (kind=np) :: xne_L,xne_q ! The linear and quadratic accelerated xne values for next
!                                   !   Saha equation evaluation.  
!                                   !   iacceleration=(0,1,2,3) chooses which of
!                                   !   xne_out,xne_mid,xne_L,xne_q to use. 
      real (kind=np) :: xne_min,xne_max ! The current allowed range for the converged electron density.
      real (kind=np) :: xne_mid ! The midpoint value for current allowed range for the 
!                               !   converged electron density.
      real (kind=np) :: xne_out ! The electron density value produced by the Saha electron density equation.
      real (kind=np) :: xne_j ! The unnormalized electron density contribution of ionization stage j 
!                             !  of element k.
!
      if(ipartfunc .eq. 0)then
         ipartfunc=1
         call partition_function_read(nel,nion,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. 2) then
             do405 : do k=1,nel      ! k is atomic number as for Mihalas-1978-112.
               jmax=min(k,nion)
               write(*,fmt=910) aname(k),k,aweight(k),                               &
     &                          xion_en(-1:jmax,k)
               write(*,fmt=920) partfunc(-1:jmax,k)
  910          format(1x,a2,i4,f6.0,18x,7f9.3)
  920          format(31x,7f9.3)
             end do do405
         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
!
      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
        phi(-1:nion-1,kk)=partfunc(-1:nion-1,kel(kk))                        &  ! Mihalas-1978-113
     &    *saha_con*tem_use**(-1.5_np)                                       &
     &    *exp(xion_en(-1:nion-1,kel(kk))/(xk_boltzmann_ev*tem_use))
!        print*,'phi assigned, etc.'
!        print*,partfunc(0,kel(kk)),saha_con,tem_use
!        print*,xion_en(0,kel(kk)),
!     &         xion_en(0,kel(kk))/(xk_boltzmann_ev*tem_use),phi(0,1)
      end do do408
!
      do410 : do kk=1,kkel
        if(kel(kk) .gt. 1) then
            jmin=0
          else
            jmin=-1
        end if
        jmax=min(kel(kk),nion)
        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) stop 
!          print*,'For ixne_initial=2 you need xne > 0'
!          print*,'Your input value is xne= ',xne
        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. 1) then
              jmin=0
            else
              jmin=-1   ! The case where H- is included.
          end if
          jmax=min(kel(kk),nion)
          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))
          fraction(j,kk)=fraction(j,kk)/denom
          xne_out=xne_out+xnum(kk)*xne_j/denom
        end do do520
        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 .ge. 1) 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
