!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