program fpreact_test ! test code for fpreact -- fortran interface to preact module ! this code tests a sampling of reaction tables; it tests ! all tables use adas_mod use fpreact_calls use periodic_table_mod implicit NONE ! note that the fpreact calls are generally vectorized, so we ! create vectors here -- even if only of length 1. integer, parameter :: jZimp_lbound = 3 , jZimp_ubound = 74 !3 <= Z impurity <=74 and 58,59,60 has some !problems with neutral radiation numerics, !so they are forbidden. integer :: jZimp ! Z for impurity integer i ! loop variable real*8 enr8(1) ! energy or energy/nucleon (KeV or KeV/AMU) real*8 tr8(1) ! temperature or T/nuclean (KeV or KeV/AMU) real*8 ansr(1),ansrr(1) ! result returned (usually , m**3/sec) real*8 var8(1) ! variation parameter (for beam-beam fusion table) real enr4(1),tr4(1),ansr4(1),var4(1),ansrr4(1) ! real versions of the above real*8, dimension(:,:), allocatable :: rioniz8, rrecom8, rcoron8, rrad8 ! impurity rates real*4, dimension(:,:), allocatable :: rioniz4, rrecom4, rcoron4, rrad4 ! impurity rates integer :: izbeam,iztarg integer ierr,ihandle,idum1,idum2 integer, dimension(2) :: izimp =(/6,8/) !SV reaction is available only for C+6, and O+8 integer, dimension(3) :: izneut = (/1,2,3/) !SV reaction. Neutrals are H0, He0 and LiO integer :: iimp,ineut real*8 ranger8(2) real ranger4(2) real*8 aa,aa1,aa2,aa3 real a integer iz,ic real :: zz character*16 tmp_r4,tmp_r8 character(len=4) :: fusion_type REAL*8,parameter:: ps_mp = 1.6726e-27_R8 ! proton mass, kg REAL*8, parameter :: elex_ah = 1.0079_r8 iz=2 a=0. aa=0. do iz=1,109 write(*,*)'tr_species_convert ',iz call tr_species_convert(iz,1,iz*2,aa,aa1,aa2,ic) call tr_species_convert(iz,1,int(aa2),aa,aa1,aa2,ic) zz= aa2/ps_mp if(iz.eq.3) then zz=6.9610_r8 call tr_species_convert(iz,1,int(zz),aa,aa1,aa2,ic) endif call tr_species_convert(iz,1,zz,aa,aa1,aa3,ic) if(aa2.ne.aa3) write(*,*)'tr_species_convert',aa2/ps_mp,aa3/ps_mp enddo tmp_r8=to_periodic_table(3,6.88_r8,1,1) write(*,*)'to_pseudo_periodic_table ',tmp_r8 tmp_r8=to_periodic_table(3,6.01_r8,1,1) write(*,*)'to_pseudo_periodic_table ',tmp_r8 tmp_r8=to_pseudo_periodic_table(1,3,6.88_r8,1,1) write(*,*)'to_pseudo_periodic_table ',tmp_r8 tmp_r8=to_pseudo_periodic_table(1,3,6.01_r8,1,1) write(*,*)'to_pseudo_periodic_table ',tmp_r8 tmp_r8=to_pseudo_periodic_table(1,3,6.96_r8,1,1) write(*,*)'to_pseudo_periodic_table ',tmp_r8 tmp_r8=to_pseudo_periodic_table(1,3,6.88,1,1) write(*,*)'to_pseudo_periodic_table ',tmp_r8 tmp_r8=to_pseudo_periodic_table(1,3,6.01,1,1) write(*,*)'to_pseudo_periodic_table ',tmp_r8 tmp_r8=to_pseudo_periodic_table(1,3,6.96,1,1) write(*,*)'to_pseudo_periodic_table ',tmp_r8 stop enr8 = 10.00_R8 ! KeV enr4 = real(enr8) tr8 = 3.00_R8 ! KeV tr4 = real(tr8) write(6,*) ' ' write(6,*) ' -->fpreact_test: atomic reaction tests...............' write(6,*) 'IMPURITY STOPPING CROSS_SECTION' !------------------- ! test: SV impurity stopping cross-section ! iimp=5 do ineut=1,size(izneut) do iimp=jZimp_lbound, jZimp_ubound call zstop_sigv(izneut(ineut),iimp,enr8,1,ansr) call zstop_sigv(izneut(ineut),iimp,enr4,1,ansr4) call zstop_sigv(izneut(ineut),dble(iimp),enr8,1,ansrr) call zstop_sigv(izneut(ineut),real(iimp),enr4,1,ansrr4) write(6,1001) izneut(ineut),iimp,enr8,ansr,ansr4,ansrr,ansrr4 call zstop_handle(ineut,iimp,ihandle) call sv_coldtarget_erange(ihandle,ranger8(1),ranger8(2),idum1,idum2) enddo enddo 1001 format(/ & ' ** impurity stopping sigma*v ** '/ & ' Z(neutral atom) = ',i1,' Z(impurity) = ',i2,' E/A = ',f10.2, & ' KeV/AMU'/ & ' -->real*8 sigma*v = ',1pe20.13,/' real sigma*v = ',1pe13.6/& ' -->real*8 sigma*v = ',1pe20.13,/' real sigma*v = ',1pe13.6/ ) 1021 format(/ & ' ** impurity stopping table energy limits'/ & ' Z(neutral atom) = ',i1,' Z(impurity) = ',i2,' E/A_min = ',1pe20.13, & ' KeV/AMU'/ & ' E/A_max = ',1pe20.13, & ' KeV/AMU'/ ) !------------------- ! test: electron impact ionization ! write(6,*) ' ' write(6,*) ' ' write(6,*) 'ELECTRON IMPACT IONOZATION' do ineut=1,size(izneut) call sigvte_ioniz(izneut(ineut),tr8,1,ansr) call sigvte_ioniz(izneut(ineut),tr4,1,ansr4) write(6,1002) izneut(ineut),tr8,ansr,ansr4 call sigvte_ioniz_handle(izneut(ineut),ihandle) call ei_sigvte_trange(ihandle,ranger8(1),ranger8(2),idum1,idum2) enddo 1002 format(/ & ' ** electron impact ionization ** '/ & ' Z(neutral atom) = ',i1,' electron temperature = ',f10.2,' KeV'/ & ' -->real*8 = ',1pe20.13,/' real = ',1pe13.6) 2001 format(/' electron ionization Z=',i1,' table limits (Te, KeV):'/ & ' (preact ei_sigvte_trange call, real*8): ',2(1x,1pe20.13)) !------------------- !------------------- enr8 = 100.00_R8 ! KeV (beam energy/nucleon) enr4 = real(enr8) tr8 = 10.00_R8 ! KeV (target temperature/nucleon) tr4 = real(tr8) !------------------- ! test: impact ionization ! write(6,*)' ' write(6,*)' ' write(6,*)'IMPACT IONIZATION NEUTRAL_BEAM ' do ineut=1,size(izneut) do iimp=1,size(izneut) izbeam=ineut iztarg=iimp call hatom_btsigv(fpreact_II,fpreact_NEUTRAL_BEAM, & enr8,tr8,1,izbeam,iztarg, ansr,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') call hatom_btsigv(fpreact_II,fpreact_NEUTRAL_BEAM, & enr4,tr4,1,izbeam,iztarg, ansr4,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') write(6,1003) ' ** ion impact ionization ** ', & izbeam,iztarg,enr8,tr8,ansr,ansr4 call hatom_handle(2,izbeam,iztarg,ihandle) call ii_warmtarget_erange(ihandle,ranger8(1),ranger8(2),idum1,idum2) call ii_coldtarget_erange(ihandle,ranger8(1),ranger8(2),idum1,idum2) enddo enddo 1023 format(/ & ' ** II warm target table energy limits'/ & ' Z(neutral atom) = ',i1,' Z(target) = ',i2,' E/A_min = ',1pe20.13, & ' KeV/AMU'/ & ' E/A_max = ',1pe20.13, & ' KeV/AMU'/ ) 1024 format(/ & ' ** II cold target table energy limits'/ & ' Z(neutral atom) = ',i1,' Z(target) = ',i2,' E/A_min = ',1pe20.13, & ' KeV/AMU'/ & ' E/A_max = ',1pe20.13, & ' KeV/AMU'/ ) write(6,*)' ' write(6,*)' ' write(6,*)'IMPACT IONIZATION ION_BEAM ' do ineut=1,size(izneut) do iimp=1,size(izneut) izbeam=ineut iztarg=iimp call hatom_btsigv(fpreact_II,fpreact_ION_BEAM, & enr8,tr8,1,izbeam,iztarg, ansr,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') call hatom_btsigv(fpreact_II,fpreact_ION_BEAM, & enr4,tr4,1,izbeam,iztarg, ansr4,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') write(6,1003) ' ** ion impact ionization ** ', & izbeam,iztarg,enr8,tr8,ansr,ansr4 enddo enddo !------------------- ! test: charge exchange rate coefficient write(6,*)' ' write(6,*)' ' write(6,*)'CHARGE EXCHANGE RATE COEFFICIENT NEUTRAL_BEAM' ! do ineut=1,size(izneut) do iimp=1,ineut izbeam=ineut iztarg=iimp call hatom_btsigv(fpreact_CX,fpreact_NEUTRAL_BEAM, & enr8,tr8,1,izbeam,iztarg, ansr,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') call hatom_btsigv(fpreact_CX,fpreact_NEUTRAL_BEAM, & enr4,tr4,1,izbeam,iztarg, ansr4,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') write(6,1003) ' ** charge exchange ** ', & izbeam,iztarg,enr8,tr8,ansr,ansr4 call hatom_handle(1,izbeam,iztarg,ihandle) call cx_warmtarget_erange(ihandle,ranger8(1),ranger8(2),idum1,idum2) call cx_coldtarget_erange(ihandle,ranger8(1),ranger8(2),idum1,idum2) enddo enddo 1025 format(/ & ' ** CX warm target table energy limits'/ & ' Z(neutral atom) = ',i1,' Z(target) = ',i2,' E/A_min = ',1pe20.13, & ' KeV/AMU'/ & ' E/A_max = ',1pe20.13, & ' KeV/AMU'/ ) 1026 format(/ & ' ** CX cold target table energy limits'/ & ' Z(neutral atom) = ',i1,' Z(target) = ',i2,' E/A_min = ',1pe20.13, & ' KeV/AMU'/ & ' E/A_max = ',1pe20.13, & ' KeV/AMU'/ ) write(6,*)' ' write(6,*)' ' write(6,*)'CHARGE EXCHANGE RATE COEFFICIENT ION_BEAM' ! do ineut=1,size(izneut) do iimp=1,ineut izbeam=iimp iztarg=ineut call hatom_btsigv(fpreact_CX,fpreact_ION_BEAM, & enr8,tr8,1,izbeam,iztarg, ansr,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') call hatom_btsigv(fpreact_CX,fpreact_ION_BEAM, & enr4,tr4,1,izbeam,iztarg, ansr4,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') write(6,1003) ' ** charge exchange ** ', & izbeam,iztarg,enr8,tr8,ansr,ansr4 enddo enddo !------------------- ! test: charge exchange Eav rate coefficient ! write(6,*)' ' write(6,*)' ' write(6,*)'CHARGE EXCHANGE Eav RATE COEFFICIENT NEUTRAL_BEAM' do ineut=1,size(izneut) do iimp=1,ineut izbeam=ineut iztarg=iimp call hatom_btsigv(fpreact_EAV,fpreact_NEUTRAL_BEAM, & enr8,tr8,1,izbeam,iztarg, ansr,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') call hatom_btsigv(fpreact_EAV,fpreact_NEUTRAL_BEAM, & enr4,tr4,1,izbeam,iztarg, ansr4,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') write(6,1003) ' ** charge exchange of target reagent ** ', & izbeam,iztarg,enr8,tr8,ansr,ansr4 enddo enddo write(6,*)' ' write(6,*)' ' write(6,*)'CHARGE EXCHANGE Eav RATE COEFFICIENT ION_BEAM' do ineut=1,size(izneut) do iimp=1,ineut izbeam=iimp iztarg=ineut call hatom_btsigv(fpreact_EAV,fpreact_ION_BEAM, & enr8,tr8,1,izbeam,iztarg, ansr,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') call hatom_btsigv(fpreact_EAV,fpreact_ION_BEAM, & enr4,tr4,1,izbeam,iztarg, ansr4,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') write(6,1003) ' ** charge exchange of target reagent ** ', & izbeam,iztarg,enr8,tr8,ansr,ansr4 enddo enddo !------------------- ! test: charge exchange rate coefficient ! write(6,*)' ' write(6,*)' ' write(6,*)'CHARGE EXCHANGE Fd RATE COEFFICIENT NEUTRAL_BEAM' do ineut=1,size(izneut) do iimp=1,ineut izbeam=ineut iztarg=iimp call hatom_btsigv(fpreact_FD,fpreact_NEUTRAL_BEAM, & enr8,tr8,1,izbeam,iztarg, ansr,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') call hatom_btsigv(fpreact_FD,fpreact_NEUTRAL_BEAM, & enr4,tr4,1,izbeam,iztarg, ansr4,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') write(6,1003) ' ** charge exchange (normalized ) ** ', & izbeam,iztarg,enr8,tr8,ansr,ansr4 enddo enddo write(6,*)' ' write(6,*)' ' write(6,*)'CHARGE EXCHANGE Fd RATE COEFFICIENT ION_BEAM ' do ineut=1,size(izneut) do iimp=1,ineut izbeam=iimp iztarg=ineut call hatom_btsigv(fpreact_FD,fpreact_ION_BEAM, & enr8,tr8,1,izbeam,iztarg, ansr,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') call hatom_btsigv(fpreact_FD,fpreact_ION_BEAM, & enr4,tr4,1,izbeam,iztarg, ansr4,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_btsigv returned an error code!') write(6,1003) ' ** charge exchange (normalized ) ** ', & izbeam,iztarg,enr8,tr8,ansr,ansr4 enddo enddo 1003 format(/,A,/ & ' Z(neutral beam) = ',i1,' Z(target plasma) = ',i1/ & ' E(beam) = ',f10.2,' KeV/AMU; T(target) = ',f10.2,' KeV/AMU'/ & ' -->real*8 result = ',1pe20.13,/' real result = ',1pe13.6) write(6,*)' ' write(6,*)' ' write(6,*)'CHARGE EXCHANGE sigma*v COEFFICIENT' do ineut=1,size(izneut) !neutrals do iimp=1,ineut !ion species izbeam=ineut iztarg=iimp call hatom_sigv(fpreact_CX,enr8,1,izbeam,iztarg, ansr,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_sigv returned an error code!') call hatom_sigv(fpreact_CX,enr4,1,izbeam,iztarg, ansr4,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_sigv returned an error code!') write(6,1013) ' ** charge exchange sigma*v (m**3/sec) **', & izbeam,iztarg,enr8,ansr,ansr4 enddo enddo write(6,*)' ' write(6,*)' ' write(6,*)' IMPACT IONIZATION sigma*v COEFFICIENT' do ineut=1,size(izneut) !neutrals do iimp=1,size(izneut)!ion species izbeam=ineut iztarg=iimp call hatom_sigv(fpreact_II,enr8,1,izbeam,iztarg, ansr,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_sigv returned an error code!') call hatom_sigv(fpreact_II,enr4,1,izbeam,iztarg, ansr4,ierr) if(ierr.ne.0) call errmsg_exit(' ?hatom_sigv returned an error code!') write(6,1013) ' ** impact ionization sigma*v (m**3/sec) **', & izbeam,iztarg,enr8,ansr,ansr4 enddo enddo 1013 format(/,A,/ & ' Z(neutral beam) = ',i1,' Z(ion species) = ',i1/ & ' E(beam) = ',f10.2,' KeV/AMU'/ & ' -->real*8 result = ',1pe20.13,/' real result = ',1pe13.6) !------------------- ! test: electron, H recombination ! write(6,*)' ' write(6,*)' ' write(6,*)'ELECTRON H RECOMBINATION' call sigvte_reco(tr8,1,ansr) call sigvte_reco(tr4,1,ansr4) write(6,1004) tr8,ansr,ansr4 1004 format(/ & ' ** e-,H+ recombination, for Te = ',f10.2,' KeV:'/ & ' -->real*8 result = ',1pe20.13,/' real result = ',1pe13.6) !------------------- write(6,*) ' ' write(6,*) ' ' write(6,*) ' ' write(6,*) 'COLD TARGET FUSION REACTION' !------------------- ! test: COLD TARGET fusion ! the following have the same reactivity ! DT -- D beam, T cold target 50 KeV D beam ! TD -- T beam, D cold target 75 KeV T beam ! the energy is must be measured in the rest frame of the target! ! (similar issue for D3/3D and T3/3T). ! ! The results will be not the same but close to each other, because each reaction ! has own table. So it is matter of an interpolation. tr8=50.0_R8 tr4=real(tr8) ! DT -- DT fusion: D beam, T target fusion_type=' DT ' call fusn_coldtarget_DT(tr8,ansr,1) call fusn_coldtarget_DT(tr4,ansr4,1) write(6,1015) fusion_type,tr8,ansr,ansr4 ! TD -- TD fusion: T beam, D target tr8=75.0_R8 tr4=real(tr8) fusion_type=' TD ' call fusn_coldtarget_TD(tr8,ansr,1) call fusn_coldtarget_TD(tr4,ansr4,1) write(6,1015) fusion_type,tr8,ansr,ansr4 tr8=50.0_R8 tr4=real(tr8) ! D3 -- D-He3 fusion: D beam, He3 target fusion_type='DHe3' call fusn_coldtarget_D3(tr8,ansr,1) call fusn_coldtarget_D3(tr4,ansr4,1) write(6,1015) fusion_type,tr8,ansr,ansr4 tr8=75.0_R8 tr4=real(tr8) ! 3D -- He3-D fusion: He3 beam, D target fusion_type='He3D' call fusn_coldtarget_3D(tr8,ansr,1) call fusn_coldtarget_3D(tr4,ansr4,1) write(6,1015) fusion_type,tr8,ansr,ansr4 tr8=50.0_R8 tr4=real(tr8) ! T3 -- T-He3 fusion: D beam, T target fusion_type='THe3' call fusn_coldtarget_T3(tr8,ansr,1) call fusn_coldtarget_T3(tr4,ansr4,1) write(6,1015) fusion_type,tr8,ansr,ansr4 ! 3T -- He3-T fusion: T beam, T target fusion_type='He3T' call fusn_coldtarget_3T(tr8,ansr,1) call fusn_coldtarget_3T(tr4,ansr4,1) write(6,1015) fusion_type,tr8,ansr,ansr4 ! DDp -- DD fusion --> p,T fusion_type=' DDp' call fusn_coldtarget_DDp(tr8,ansr,1) call fusn_coldtarget_DDp(tr4,ansr4,1) write(6,1015) fusion_type,tr8,ansr,ansr4 ! DDn -- DD fusion --> n,He3 fusion_type=' DDn' call fusn_coldtarget_DDn(tr8,ansr,1) call fusn_coldtarget_DDn(tr4,ansr4,1) write(6,1015) fusion_type,tr8,ansr,ansr4 ! TT -- TT fusion --> 2n,He4 fusion_type=' TT ' call fusn_coldtarget_TT(tr8,ansr,1) call fusn_coldtarget_TT(tr4,ansr4,1) write(6,1015) fusion_type,tr8,ansr,ansr4 1015 format(/ & ' ** sigma*v[m**3/sec]',1x,a,1x,'cold target fusion, beam energy = ',f10.2,' KeV:'/ & ' -->real*8 result = ',1pe20.13,/' real result = ',1pe13.6) !------------------- write(6,*) ' ' write(6,*) ' ' write(6,*) ' ' write(6,*) 'THERMONUCLEAR FUSION REACTION' !------------------- ! test: thermonuclear fusion ! --thermonuclear rate coefficients (vs. Temperature) ! DT -- DT fusion: D beam, T target fusion_type=' DT ' call fusn_thermo_DT(tr8,ansr,1) call fusn_thermo_DT(tr4,ansr4,1) write(6,1005) fusion_type,tr8,ansr,ansr4 ! TD -- TD fusion: T beam, D target fusion_type=' TD ' call fusn_thermo_TD(tr8,ansr,1) call fusn_thermo_TD(tr4,ansr4,1) write(6,1005) fusion_type,tr8,ansr,ansr4 ! DDp -- DD fusion --> p,T fusion_type=' DDp' call fusn_thermo_DDp(tr8,ansr,1) call fusn_thermo_DDp(tr4,ansr4,1) write(6,1005) fusion_type,tr8,ansr,ansr4 ! DDn -- DD fusion --> n,He3 fusion_type=' DDn' call fusn_thermo_DDn(tr8,ansr,1) call fusn_thermo_DDn(tr4,ansr4,1) write(6,1005) fusion_type,tr8,ansr,ansr4 ! TT -- TT fusion --> 2n,He4 fusion_type=' TT ' call fusn_thermo_TT(tr8,ansr,1) call fusn_thermo_TT(tr4,ansr4,1) write(6,1005) fusion_type,tr8,ansr,ansr4 ! D3 -- D-He3 fusion: D beam, He3 target fusion_type='DHe3' call fusn_thermo_D3(tr8,ansr,1) call fusn_thermo_D3(tr4,ansr4,1) write(6,1005) fusion_type,tr8,ansr,ansr4 ! 3D -- He3-D fusion: He3 beam, D target fusion_type='He3D' call fusn_thermo_3D(tr8,ansr,1) call fusn_thermo_3D(tr4,ansr4,1) write(6,1005) fusion_type,tr8,ansr,ansr4 ! T3 -- T-He3 fusion: D beam, T target fusion_type='THe3' call fusn_thermo_T3(tr8,ansr,1) call fusn_thermo_T3(tr4,ansr4,1) write(6,1005) fusion_type,tr8,ansr,ansr4 ! 3T -- He3-T fusion: T beam, T target fusion_type='He3T' call fusn_thermo_3T(tr8,ansr,1) call fusn_thermo_3T(tr4,ansr4,1) write(6,1005) fusion_type,tr8,ansr,ansr4 1005 format(/ & ' ** <> [m**3/sec]',1x,a,1x,'thermonuclear fusion, Ti = ',f10.2,' KeV:'/ & ' -->real*8 result = ',1pe20.13,/' real result = ',1pe13.6) !------------------- ! test: beam-target fusion ! write(6,*) ' ' write(6,*) ' ' write(6,*) ' ' write(6,*) 'BEAM-TARGET FUSION .' fusion_type=' DT ' call fusn_btsigv_DT(enr8,tr8,ansr,1) call fusn_btsigv_DT(enr4,tr4,ansr4,1) write(6,1006) fusion_type,'D',enr8,'T',tr8,ansr,ansr4 fusion_type=' TD ' call fusn_btsigv_TD(enr8,tr8,ansr,1) call fusn_btsigv_TD(enr4,tr4,ansr4,1) write(6,1006) fusion_type,'T',enr8,'D',tr8,ansr,ansr4 fusion_type=' DDp' call fusn_btsigv_DDp(enr8,tr8,ansr,1) call fusn_btsigv_DDp(enr4,tr4,ansr4,1) write(6,1006) fusion_type,'D',enr8,'D',tr8,ansr,ansr4 fusion_type=' DDn' call fusn_btsigv_DDn(enr8,tr8,ansr,1) call fusn_btsigv_DDn(enr4,tr4,ansr4,1) write(6,1006) fusion_type,'D',enr8,'D',tr8,ansr,ansr4 fusion_type=' TT ' call fusn_btsigv_TT(enr8,tr8,ansr,1) call fusn_btsigv_TT(enr4,tr4,ansr4,1) write(6,1006) fusion_type,'T',enr8,'T',tr8,ansr,ansr4 fusion_type='DHe3' call fusn_btsigv_D3(enr8,tr8,ansr,1) call fusn_btsigv_D3(enr4,tr4,ansr4,1) write(6,1006) fusion_type,'D',enr8,'He3',tr8,ansr,ansr4 fusion_type='He3D' call fusn_btsigv_3D(enr8,tr8,ansr,1) call fusn_btsigv_3D(enr4,tr4,ansr4,1) write(6,1006) fusion_type,'He3',enr8,'D',tr8,ansr,ansr4 fusion_type='THe3' call fusn_btsigv_T3(enr8,tr8,ansr,1) call fusn_btsigv_T3(enr4,tr4,ansr4,1) write(6,1006) fusion_type,'T',enr8,'He3',tr8,ansr,ansr4 fusion_type='He3T' call fusn_btsigv_3T(enr8,tr8,ansr,1) call fusn_btsigv_3T(enr4,tr4,ansr4,1) write(6,1006) fusion_type,'He3',enr8,'T',tr8,ansr,ansr4 1006 format(/ & ' ** ',1x,a,1x,' beam-target fusion **'/ & 1x,a,1x,' beam energy = ',f10.2,'KeV',1x,a,1x,' target Ti = ',f10.2,' KeV:'/ & ' -->real*8 result = ',1pe20.13,/' real result = ',1pe13.6) !------------------- ! test: beam-beam fusion ! write(6,*) ' ' write(6,*) ' ' write(6,*) ' ' write(6,*) ' BEAM-BEAM FUSION ' tr8=0.50_R8 !fractional beam energy variation in the target frame tr4=real(tr8) fusion_type=' DT ' call fusn_bbsigv_DT(enr8,tr8,ansr,1) call fusn_bbsigv_DT(enr4,tr4,ansr4,1) write(6,1007) fusion_type,'D',enr8,'T',tr8,ansr,ansr4 ! fusion_type=' TD ' ! call fusn_bbsigv_TD(enr8,tr8,ansr,1) ! call fusn_bbsigv_TD(enr4,tr4,ansr4,1) ! write(6,1007) fusion_type,'T',enr8,'D',tr8,ansr,ansr4 fusion_type=' DDp' call fusn_bbsigv_DDp(enr8,tr8,ansr,1) call fusn_bbsigv_DDp(enr4,tr4,ansr4,1) write(6,1007) fusion_type,'D',enr8,'D',tr8,ansr,ansr4 fusion_type=' DDn' call fusn_bbsigv_DDn(enr8,tr8,ansr,1) call fusn_bbsigv_DDn(enr4,tr4,ansr4,1) write(6,1007) fusion_type,'D',enr8,'D',tr8,ansr,ansr4 fusion_type=' TT ' call fusn_bbsigv_TT(enr8,tr8,ansr,1) call fusn_bbsigv_TT(enr4,tr4,ansr4,1) write(6,1007) fusion_type,'T',enr8,'T',tr8,ansr,ansr4 fusion_type='DHe3' call fusn_bbsigv_D3(enr8,tr8,ansr,1) call fusn_bbsigv_D3(enr4,tr4,ansr4,1) write(6,1007) fusion_type,'D',enr8,'He3',tr8,ansr,ansr4 ! fusion_type='He3D' ! call fusn_bbsigv_3D(enr8,tr8,ansr,1) ! call fusn_bbsigv_3D(enr4,tr4,ansr4,1) ! write(6,1007) fusion_type,'He3',enr8,'D',tr8,ansr,ansr4 fusion_type='THe3' call fusn_bbsigv_T3(enr8,tr8,ansr,1) call fusn_bbsigv_T3(enr4,tr4,ansr4,1) write(6,1007) fusion_type,'T',enr8,'He3',tr8,ansr,ansr4 ! fusion_type='He3T' ! call fusn_bbsigv_3T(enr8,tr8,ansr,1) ! call fusn_bbsigv_3T(enr4,tr4,ansr4,1) ! write(6,1007) fusion_type,'He3',enr8,'T',tr8,ansr,ansr4 1007 format(/ & ' ** <>',1x,a,1x,'beam-beam fusion **'/ & 1x,a,1x,'beam energy = ',f10.2,'KeV',1x,a,1x,', variation fraction = ',f10.2,':'/ & ' -->real*8 result = ',1pe20.13,/' real result = ',1pe13.6) write(6,*) ' ' write(6,*) ' ' write(6,*) ' ' write(6,*) 'IMPURITY REACTION ' ! ------------------ ! test: impurity radiation rate ! tr8(1) = .1_R8 ! Te in keV tr4(1) = .1 var8(1) = 1.e19_R8 ! Ne in m-3 var4(1) = 1.e19 ! call imptable_def_brem(0) default ! -- ionization -- do jZimp = jZimp_lbound, jZimp_ubound !allowed range forZ impurity impurity: if( jZimp /= 58 .and. jZimp /= 59 .and. jZimp /= 60) then allocate( rioniz8(jZimp+1,1), rioniz4(jZimp+1,1)) allocate( rrecom8(jZimp+1,1), rrecom4(jZimp+1,1)) allocate( rcoron8(jZimp+1,1), rcoron4(jZimp+1,1)) allocate( rrad8(jZimp+1,1), rrad4(jZimp+1,1)) call imptable_ioniz(jZimp, 1, tr8, var8, jZimp+1, rioniz8) call imptable_ioniz(jZimp, 1, tr4, var4, jZimp+1, rioniz4) write(6,1101) ' ** Impurity ionization rate ' // & trim(name_pseudo_periodic_table(0,jZimp,0.)) // ' **', & jZimp, tr8(1), var8(1) do i = 1, jZimp+1 write(6,1102) (i-1), rioniz8(i,1), rioniz4(i,1) end do ! -- recombination -- call imptable_recom(jZimp, 1, tr8, var8, jZimp+1, rrecom8) call imptable_recom(jZimp, 1, tr4, var4, jZimp+1, rrecom4) write(6,1101) ' ** Impurity recombination rate ' // & trim(name_pseudo_periodic_table(0,jZimp,0.)) // ' **', & jZimp, tr8(1), var8(1) do i = 1, jZimp+1 write(6,1102) (i-1), rrecom8(i,1), rrecom4(i,1) end do ! -- bremsstrahlung -- call imptable_brem(jZimp, 1, tr8, var8, jZimp+1, rrad8) call imptable_brem(jZimp, 1, tr4, var4, jZimp+1, rrad4) write(6,1101) ' ** Impurity bremsstrahlung ' // & trim(name_pseudo_periodic_table(0,jZimp,0.)) // ' **', & jZimp, tr8(1), var8(1) do i = 1, jZimp+1 write(6,1102) (i-1), rrad8(i,1), rrad4(i,1) end do rrad8=0.0_R8 rrad4=0.0 ! -- radiation (no bremsstrahlung) -- call imptable_rad(jZimp, 1, tr8, var8, jZimp+1, rrad8) call imptable_rad(jZimp, 1, tr4, var4, jZimp+1, rrad4) write(6,1101) ' ** Impurity radiation rate (no bremsstrahlung) ' // & trim(name_pseudo_periodic_table(0,jZimp,0.)) // ' **', & jZimp, tr8(1), var8(1) do i = 1, jZimp+1 write(6,1102) (i-1), rrad8(i,1), rrad4(i,1) end do ! -- coronal equilibrium -- call coronal(jZimp, 1, rioniz8, rrecom8, jZimp+1, rcoron8) call coronal(jZimp, 1, rioniz4, rrecom4, jZimp+1, rcoron4) write(6,1101) ' ** Impurity coronal equilibrium ' // & trim(name_pseudo_periodic_table(0,jZimp,0.)) // ' **', & jZimp, tr8(1), var8(1) do i = 1, jZimp+1 write(6,1102) (i-1), rcoron8(i,1), rcoron4(i,1) end do ! -- coronal equilibrium radiaiton -- call coronal_rad(jZimp, 1, rcoron8, rrad8, jZimp+1, enr8) call coronal_rad(jZimp, 1, rcoron4, rrad4, jZimp+1, enr4) aa=0.0_r8 write(6,1101) ' ** Impurity coronal equilibrium radiation (no bremsstrahlung) ' // & trim(name_pseudo_periodic_table(0,jZimp,aa)) // ' **', & jZimp, tr8(1), var8(1) write(6,1103) enr8(1), enr4(1) write(6,'(//)') deallocate( rioniz8, rioniz4) deallocate( rrecom8, rrecom4) deallocate( rcoron8, rcoron4) deallocate( rrad8, rrad4) ! call imptable_print(jZimp) else write(6,'(//)') write(6,*)' -->fpreact_test: abnormal exit: Z impurity out of range or .eq. forbidden numbers :58, 59, 60' endif impurity enddo 1101 format(/,A,/ & ' Z(impurity) = ',i2,' Te = ',f10.2,' KeV',& ' Ne = ',es10.3, ' m-3') 1102 format(' charge = ',i2,' -->real*8 result = ',1pe20.13,/12x,' real result = ',1pe13.6) 1103 format(' -->real*8 result = ',1pe20.13,/' real result = ',1pe13.6) a=12.0107 aa=12.0107_r8 tmp_r4=to_periodic_table(6,a,-1,1) tmp_r8=to_periodic_table(6,aa,-1,1) write(*,*) 'to periodic r4 ',tmp_r4 write(*,*) 'to periodic r8 ',tmp_r8 call inv_periodic_table('C', .true., iz, a, ic) write(*,*) 'r4 ',iz, a, ic call inv_periodic_table('C', .true., iz, aa, ic) write(*,*) 'r8 ',iz, aa, ic tmp_r4= name_periodic_table (iz,a) tmp_r8= name_periodic_table (iz,aa) write(*,*) 'name r4 ',tmp_r4 write(*,*) 'name r8 ',tmp_r8 a=1. aa=1._r8 tmp_r4=to_pseudo_periodic_table(1,1,a,-1,1) tmp_r8=to_pseudo_periodic_table(1,1,aa,-1,1) write(*,*) 'to pseudo periodic r4 ',tmp_r4 write(*,*) 'to pseudo periodic r8 ',tmp_r8 call inv_pseudo_periodic_table(1,'H', .true., iz, a, ic) write(*,*) 'pseudo r4 ',iz, a, ic call inv_pseudo_periodic_table(1,'H', .true., iz, aa, ic) write(*,*) 'pseudo r8 ',iz, aa, ic call standard_amu(6,a) write(*,*) 'r4 A=',a call standard_amu(6,aa) write(*,*) 'r8 A=',aa iz=11.8 write(*,*)'iz= ',iz iz=11.8+0.01_r8 write(*,*)'iz= ',iz iz=11.8+0.5_r8 write(*,*)'iz= ',iz iz=11.1 write(*,*)'iz= ',iz iz=11.1+0.01_r8 write(*,*)'iz= ',iz iz=11.1+0.5_r8 write(*,*)'iz= ',iz call tr_species_convert(6,5,12,aa1,aa,aa2,ic) aa2=aa2/1.6726e-27_R8 aa1=aa1/1.6022e-19_R8 aa=aa/1.6022e-19_R8 write(*,*)aa1,aa,aa2,ic call tr_species_convert(3,3,6,aa1,aa,aa2,ic) aa2=aa2/1.6726e-27_R8 aa1=aa1/1.6022e-19_R8 aa=aa/1.6022e-19_R8 write(*,*)aa1,aa,aa2,ic end program fpreact_test