64 integer :: ai, ni, zi, snfile, fstat
65 character :: nl, buf(30)
67 character*10,
dimension(track_nuclei_nr) :: tmp_nucleinames
68 character*1500 :: head_toplist
69 character*5 :: help_char
70 info_entry(
"analysis_init")
87 '# 1:iteration 2:time[s], 3:T[GK], 4:rho[g/cm3], 5:Ye '//nl &
88 ,
'# 6:R[km], 7:Y_n, 8:Y_p, 9:Y_alpha, 10:Y_lights '//nl &
89 ,
'# 11:Y_heavies 12:<Z> 13:<A> 14:entropy [kB/baryon] (15:Sn [MeV])'
96 allocate(
sn(net_size))
100 read(snfile, *, iostat=fstat) zi, ni, ai, val
113 write(
tsfile,
'((A1,(A22,3x),*(A23,3x)))') &
114 "#",
"time[s]",
"Temperature [GK]",
"tau_ga [s]",
"tau_ag [s]",
"tau_ng [s]",
"tau_gn [s]", &
115 "tau_pg [s]",
"tau_gp [s]",
"tau_np [s]",
"tau_pn [s]", &
116 "tau_an [s]",
"tau_na [s]",
"tau_ap [s]",
"tau_pa [s]",
"tau_beta [s]",
"tau_alpha [s]", &
117 "tau_nfiss [s]",
"tau_sfiss [s]",
"tau_bfiss [s]"
128 write(
track_unit,
'((A,20X,*(A,18X)))')
"# time[s]",tmp_nucleinames(:)
136 write (
toplist_unit,
"('#',31X,A,434X,A,413X,'|',5X,A,185X,A)")
"Top nuclear energy producers for each reaction given in erg/g/s", &
137 "Top nuclear energy consumers for each reaction given in erg/g/s", &
138 "Top nuclear energy producers for each isotope given in erg/g/s", &
139 "Top nuclear energy consumers for each isotope given in erg/g/s"
142 head_toplist =
"# Time[s] Etot[erg/g/s] Reaction_1"
144 head_toplist = trim(adjustl(head_toplist))//
" Reaction_"//trim(adjustl(
int_to_str(i)))
146 head_toplist = trim(adjustl(head_toplist))//
" Ereac_1"
148 head_toplist = trim(adjustl(head_toplist))//
" Ereac_"//trim(adjustl(
int_to_str(i)))
150 head_toplist = trim(adjustl(head_toplist))//
" Reaction_1"
152 head_toplist = trim(adjustl(head_toplist))//
" Reaction_"//trim(adjustl(
int_to_str(i)))
154 head_toplist = trim(adjustl(head_toplist))//
" Ereac_1"
156 head_toplist = trim(adjustl(head_toplist))//
" Ereac_"//trim(adjustl(
int_to_str(i)))
158 head_toplist = trim(adjustl(head_toplist))//
" | Iso1"
160 help_char =
"Iso"//trim(adjustl(
int_to_str(i)))
161 head_toplist = trim(adjustl(head_toplist))//
" "//adjustr(help_char)
163 head_toplist = trim(adjustl(head_toplist))//
" Eiso_1"
165 head_toplist = trim(adjustl(head_toplist))//
" Eiso_"//trim(adjustl(
int_to_str(i)))
167 head_toplist = trim(adjustl(head_toplist))//
" Iso1"
169 help_char =
"Iso"//trim(adjustl(
int_to_str(i)))
170 head_toplist = trim(adjustl(head_toplist))//
" "//adjustr(help_char)
172 head_toplist = trim(adjustl(head_toplist))//
" Eiso_1"
174 head_toplist = trim(adjustl(head_toplist))//
" Eiso_"//trim(adjustl(
int_to_str(i)))
184 '# Newton-Raphson diagnostic output'//nl &
185 ,
'# 1 : global iteration count (cnt)'//nl &
186 ,
'# 2 : global time [s]'//nl &
187 ,
'# 3 : temperature T9 [GK]'//nl &
188 ,
'# 4 : adapt stepsize loop counter (k)'//nl &
189 ,
'# 5 : NR loop counter (nr_count)'//nl &
190 ,
'# 6 : global timestep [s]'//nl &
191 ,
'# 7 : mass conservation error'//nl &
192 ,
'# 8 : maximal abundance change (eps)'//nl &
193 ,
'# 9 : most rapidly evolving isotope (epsl)'//nl &
194 ,
'# 10: abundance of the isotope epsl, y(epsl)'
201 write(
nuc_heat_file,
'((A,5X))')
"# Generated Energy for each nuclear reaction type given in erg/g/s"
203 "# time[s]",
"Engen(Total)",
"S_src",
"Engen(weak)",
"Engen(n,g)",&
204 "Engen(p,g)",
"Engen(a,g)",
"Engen(n,p)",
"Engen(n,a)", &
205 "Engen(p,a)",
"Engen(fiss)"
212 write(
nu_loss_gain_unit,
'(A)')
"# File containing the neutrino loss/gain for the nuclear heating."
213 write(
nu_loss_gain_unit,
'(A)')
"# Negative values mean that energy is added to the system."
216 "#",
"time[s]",
"T[GK]",
"rho[g/cm3]",
"R[km]",
"nu_total",
"nu_beta",&
217 "nu_nuheat",
"nu_thermal"
221 info_exit(
"analysis_init")
234 real(
r_kind),
intent(in) :: ctime
235 real(
r_kind),
intent(in) :: t9
236 real(
r_kind),
intent(in) :: rho_b
237 real(
r_kind),
intent(in) :: entropy
238 real(
r_kind),
intent(in) :: rkm
239 real(
r_kind),
dimension(net_size),
intent(in) :: y
240 real(
r_kind),
dimension(0:net_size),
intent(inout) :: pf
242 info_entry(
'output_initial_step')
244 info_exit(
'output_initial_step')
261 integer,
intent(in) :: cnt
262 real(
r_kind),
intent(in) :: ctime
263 real(
r_kind),
intent(in) :: t9
264 real(
r_kind),
intent(in) :: rho_b
265 real(
r_kind),
intent(in) :: entropy
266 real(
r_kind),
intent(in) :: rkm
267 real(
r_kind),
dimension(net_size),
intent(in) :: y
268 real(
r_kind),
dimension(0:net_size),
intent(inout) :: pf
270 info_entry(
'output_final_step')
277 if (timescales_every.gt.0) timescales_every= 1
279 print
'("---------------")'
281 print
'("============================")'
284 print
'(A)',
"End of trajectory file reached."
286 print
'(A)',
"Final time reached."
288 print
'(A)',
"Final temperature reached."
290 print
'(A)',
"Final density reached."
292 print
'(A)',
"Freeze-out reached."
295 info_exit(
'output_final_step')
328 integer,
intent(in) :: cnt
329 real(
r_kind),
intent(in) :: ctime
330 real(
r_kind),
intent(in) :: t9
331 real(
r_kind),
intent(in) :: rho_b
332 real(
r_kind),
intent(in) :: entropy
333 real(
r_kind),
intent(in) :: rkm
334 real(
r_kind),
dimension(net_size),
intent(in) :: y
335 real(
r_kind),
dimension(0:net_size),
intent(inout) :: pf
337 real(
r_kind) :: ysum,yasum,abar,ye
339 info_entry(
"output_iteration")
345 print
'(A)',
"#-- 1:it 2:time[s] 3:T9[GK] 4:rho[g/cm3] &
346 5:entropy[kB/baryon]"
349 print
'(I8,X,ES14.7,X,F10.7,X,ES14.7,X,F12.7)', &
350 cnt,ctime,t9,rho_b,entropy
354 if (verbose_level.ge.2)
then
426 ysum= sum(y(1:net_size))
427 yasum= sum(y(1:net_size)*
isotope(1:net_size)%mass)
428 ye = sum(
isotope(1:net_size)%p_nr*y(1:net_size))
468 ysum= sum(y(1:net_size))
469 yasum= sum(y(1:net_size)*
isotope(1:net_size)%mass)
470 ye = sum(
isotope(1:net_size)%p_nr*y(1:net_size))
495 if (timescales_every .gt. 0)
then
496 if ((mod(cnt,timescales_every).eq.0) .and. (
evolution_mode .ne. em_nse))
then
526 info_exit(
"output_iteration")
550 real(
r_kind),
intent(in) :: ctime
551 real(
r_kind),
dimension(net_size),
intent(in) :: y
552 real(
r_kind),
dimension(track_nuclei_nr) :: tmp_y
555 info_entry(
"output_track_nuclei")
559 if (tmp_y(i) .lt. 1d-99) tmp_y(i) = 0.
562 write(
track_unit,
'(*(es23.16,5x))') ctime,tmp_y(:)
563 info_exit(
"output_track_nuclei")
592 real(
r_kind),
intent(in) :: ctime
593 real(
r_kind),
intent(in) :: temp
594 real(
r_kind),
intent(in) :: rho
595 real(
r_kind),
intent(in) :: rkm
630 integer,
intent(in) :: cnt
631 real(
r_kind),
intent(in) :: ctime
632 real(
r_kind),
intent(in) :: t9
633 real(
r_kind),
intent(in) :: rho_b
634 real(
r_kind),
intent(in) :: entropy
635 real(
r_kind),
intent(in) :: rkm
636 real(
r_kind),
dimension(net_size),
intent(in) :: y
638 real(
r_kind) :: ye,yneu,ypro,yhe4,ylight,yheavies,ysum,yasum,yzsum,abar,zbar
646 ye = sum(
isotope(1:net_size)%p_nr*y(1:net_size))
647 ylight= max(1.0d-99,sum(y(ipro+1:ihe4-1)))
649 if (ihe4.ge.1) yheavies= max(1.0d-99,sum(y(ihe4+1:net_size)))
651 ysum = sum(y(1:net_size))
652 yasum = sum(y(1:net_size)*
isotope(1:net_size)%mass)
653 yzsum = sum(y(1:net_size)*
isotope(1:net_size)%p_nr)
659 if (ineu.ge.1) yneu= max(1.0d-99,y(ineu))
662 if (ipro.ge.1) ypro= max(1.0d-99,y(ipro))
665 if (ihe4.ge.1) yhe4= max(1.0d-99,y(ihe4))
668 write (
mainout_unit,
'(i8,X,100(es23.16,2X))') cnt,ctime,t9,rho_b,ye,rkm, &
669 yneu,ypro,yhe4,ylight,yheavies,zbar,abar,entropy,
sn_ave
671 write (
mainout_unit,
'(i8,X,100(es23.16,2X))') cnt,ctime,t9,rho_b,ye,rkm, &
672 yneu,ypro,yhe4,ylight,yheavies,zbar,abar,entropy
698 real(
r_kind),
intent(in) :: ctime
699 real(
r_kind),
intent(in) :: t9
700 real(
r_kind),
intent(in) :: rho_b
701 real(
r_kind),
dimension(net_size),
intent(in):: y
719 real(r_kind),
intent(in) :: ctime
720 real(r_kind),
intent(in) :: T9
721 real(r_kind),
intent(in) :: rho_b
722 real(r_kind),
dimension(net_size),
intent(in):: Y
757 integer,
intent(in) :: cnt
758 integer,
intent(in) :: k
759 integer,
intent(in) :: nr_c
760 real(
r_kind),
intent(in) :: ctime
761 real(
r_kind),
intent(in) :: t9
762 real(
r_kind),
intent(in) :: stepsize
763 real(
r_kind),
intent(in) :: mtot
764 real(
r_kind),
dimension(net_size) :: y_p,y
770 info_entry(
"output_nr_diagnostic")
783 if(y_p(i).lt.1.d-10) cycle
784 if(dabs(y_p(i)-y(i)).eq.0.d0) cycle
785 epst = dabs(y(i)/y_p(i)-1.d0)
786 if (epst.gt.eps)
then
794 '(i7,2(1x,es19.11),2(1x,i4),3(1x,es19.11),a5,x,es19.11)') &
795 cnt,ctime,t9,k,nr_c,stepsize,dabs(mtot-1d0),eps,
isotope(epsl)%name, &
799 info_exit(
"output_nr_diagnostic")
821 integer,
intent(in) :: cnt
822 real(
r_kind),
intent(in) :: ctime
826 info_entry(
'output_final_stats')
836 print
'(A,ES19.12)',
'Final time: ', ctime
837 print
'(A,I10)',
'Total number of iterations: ', cnt
838 print
'(A,ES12.5)',
'Elapsed time [s]: ', simtime
840 info_exit(
'output_final_stats')
862 real(
r_kind),
intent(in) :: t,t9,rho_b
863 real(
r_kind),
dimension(net_size),
intent(in) :: y
864 character (len=21) :: snapfile
865 integer :: snapunit, z
868 info_entry(
"printsnap")
870 write(snapfile,
'("snaps/snapsh_",i4.4,".dat")')
snapcnt
873 write(snapunit,
'(3a8)')
'time',
'temp',
'dens'
874 write(snapunit,
'(3es23.14)')t,t9,rho_b
875 write(snapunit,
'(4(a8))')
'nin',
'zin',
'y',
'x'
878 write(snapunit,
'(2(i3,2x),2(es23.14))') &
887 info_exit(
"printsnap")
902 real(
r_kind),
dimension(net_size),
intent(in) :: y
904 real(
r_kind) :: ysum,yasum,yzsum,abar_ions,zbar_ions
907 info_entry(
"calc_nseparation_energy")
911 do i=max(1,ihe4),net_size
912 abar_ions = abar_ions +
isotope(i)%mass * y(i)
913 zbar_ions = zbar_ions +
isotope(i)%p_nr * y(i)
916 ysum = sum(y(max(1,ihe4):net_size))
917 abar_ions = abar_ions / ysum
918 zbar_ions = zbar_ions / ysum
921 do i=max(1,ihe4),net_size
934 info_exit(
"calc_nseparation_energy")
969 logical,
intent(in) :: hdf5_mode
971 logical,
intent(in) :: write_toplist
977 integer :: help_parts
980 real(
r_kind) :: engen_q_value
981 real(
r_kind) :: engen_q_weak
982 real(
r_kind) :: engen_q_ng,engen_q_gn
983 real(
r_kind) :: engen_q_pg,engen_q_gp
984 real(
r_kind) :: engen_q_ga,engen_q_ag
985 real(
r_kind) :: engen_q_pa,engen_q_ap
986 real(
r_kind) :: engen_q_pn,engen_q_np
987 real(
r_kind) :: engen_q_na,engen_q_an
988 real(
r_kind) :: engen_fiss
989 real(
r_kind) :: engen_tot_ncap, engen_tot_pcap, engen_tot_acap, &
990 engen_tot_npcap,engen_tot_nacap, engen_tot_pacap
991 real(
r_kind) :: engen_tot_exc
992 real(
r_kind) :: entropy_src_term
993 real(
r_kind) :: mic2,gi,twopi_h2c2,mui
995 real(
r_kind),
dimension(0:net_size)::pf
996 real(
r_kind) :: etaele,etapos
999 real(
r_kind) :: tmp_energy
1000 real(
r_kind),
dimension(net_size) ::det_bet_engen
1001 real(
r_kind),
dimension(net_size) ::det_ag_engen
1002 real(
r_kind),
dimension(net_size) ::det_pg_engen
1003 real(
r_kind),
dimension(net_size) ::det_ng_engen
1004 real(
r_kind),
dimension(net_size) ::det_pn_engen
1005 real(
r_kind),
dimension(net_size) ::det_ap_engen
1006 real(
r_kind),
dimension(net_size) ::det_an_engen
1007 real(
r_kind),
dimension(net_size) ::det_other_engen
1008 real(
r_kind),
dimension(net_size) ::det_fiss_engen
1009 real(
r_kind),
dimension(net_size) ::det_tot_engen
1011 integer,
dimension(nreac) ::toplist_idx
1012 real(
r_kind),
dimension(nreac) ::toplist_energy
1013 logical,
dimension(nreac) ::toplist_mask
1014 character*30 ::reac_dummy
1015 character*30,
dimension(toplist_amount)::tplist_str
1016 real(
r_kind),
dimension(toplist_amount)::tplist
1017 character*30,
dimension(toplist_amount)::btlist_str
1018 real(
r_kind),
dimension(toplist_amount)::btlist
1019 real(
r_kind),
dimension(net_size) ::toplist_energy_iso
1020 integer,
dimension(net_size) ::toplist_idx_iso
1021 character*5,
dimension(toplist_amount) ::tplist_str_iso
1022 real(
r_kind),
dimension(toplist_amount)::tplist_iso
1023 character*5,
dimension(toplist_amount) ::btlist_str_iso
1024 real(
r_kind),
dimension(toplist_amount)::btlist_iso
1026 info_entry(
"calc_nuclear_heating")
1047 det_bet_engen(:) = 0.0
1048 det_ag_engen(:) = 0.0
1049 det_pg_engen(:) = 0.0
1050 det_ng_engen(:) = 0.0
1051 det_pn_engen(:) = 0.0
1052 det_ap_engen(:) = 0.0
1053 det_an_engen(:) = 0.0
1054 det_other_engen(:) = 0.0
1055 det_fiss_engen(:) = 0.0
1056 det_tot_engen(:) = 0.0
1061 if (write_toplist)
then
1062 toplist_energy(:) = 0d0
1063 toplist_energy_iso(:) = 0d0
1064 toplist_mask(:) = .false.
1075 if ((
evolution_mode .eq. em_nse) .and. (.not. rr_tmp%is_weak)) cycle
1078 if (rat .eq. 0) cycle
1081 select case (rr_tmp%group)
1084 tmp_energy = rat*
y(rr_tmp%parts(1)) * rr_tmp%q_value
1086 engen_q_value = engen_q_value + tmp_energy
1090 if (rr_tmp%is_weak)
then
1091 engen_q_weak = engen_q_weak + tmp_energy
1093 det_bet_engen(rr_tmp%parts(1)) = det_bet_engen(rr_tmp%parts(1))+tmp_energy
1097 if(rr_tmp%reac_type.eq.
rrt_gn)
then
1098 engen_q_gn = engen_q_gn + tmp_energy
1101 if (rr_tmp%parts(2) .eq.
ineu)
then
1102 help_parts = rr_tmp%parts(3)
1104 help_parts = rr_tmp%parts(2)
1106 det_ng_engen(help_parts) = det_ng_engen(help_parts)+tmp_energy
1111 if(rr_tmp%reac_type.eq.
rrt_gp)
then
1112 engen_q_gp = engen_q_gp + tmp_energy
1115 if (rr_tmp%parts(2) .eq.
ipro)
then
1116 help_parts = rr_tmp%parts(3)
1118 help_parts = rr_tmp%parts(2)
1120 det_pg_engen(help_parts) = det_pg_engen(help_parts)+tmp_energy
1124 if(rr_tmp%reac_type.eq.
rrt_ga)
then
1125 engen_q_ga = engen_q_ga + tmp_energy
1128 if (rr_tmp%parts(2) .eq.
ihe4)
then
1129 help_parts = rr_tmp%parts(3)
1131 help_parts = rr_tmp%parts(2)
1133 det_ag_engen(help_parts) = det_ag_engen(help_parts)+tmp_energy
1138 if(rr_tmp%reac_type.eq.
rrt_o)
then
1140 det_other_engen(rr_tmp%parts(1)) = det_other_engen(rr_tmp%parts(1))+tmp_energy
1146 tmp_energy = rat*
y(rr_tmp%parts(1)) *
y(rr_tmp%parts(2)) &
1149 engen_q_value = engen_q_value + tmp_energy
1153 if(rr_tmp%reac_type.eq.
rrt_ng)
then
1154 engen_q_ng = engen_q_ng + tmp_energy
1157 if (rr_tmp%parts(2) .eq.
ineu)
then
1158 help_parts = rr_tmp%parts(1)
1160 help_parts = rr_tmp%parts(2)
1162 det_ng_engen(help_parts) = det_ng_engen(help_parts)+tmp_energy
1166 if(rr_tmp%reac_type.eq.
rrt_pg)
then
1167 engen_q_pg = engen_q_pg + tmp_energy
1170 if (rr_tmp%parts(2) .eq.
ipro)
then
1171 help_parts = rr_tmp%parts(1)
1173 help_parts = rr_tmp%parts(2)
1175 det_pg_engen(help_parts) = det_pg_engen(help_parts)+tmp_energy
1179 if(rr_tmp%reac_type.eq.
rrt_ag)
then
1180 engen_q_ag = engen_q_ag + tmp_energy
1183 if (rr_tmp%parts(2) .eq.
ihe4)
then
1184 help_parts = rr_tmp%parts(1)
1186 help_parts = rr_tmp%parts(2)
1188 det_ag_engen(help_parts) = det_ag_engen(help_parts)+tmp_energy
1192 if(rr_tmp%reac_type.eq.
rrt_pa)
then
1193 engen_q_pa = engen_q_pa + tmp_energy
1196 if (rr_tmp%parts(3) .eq.
ihe4)
then
1197 help_parts = rr_tmp%parts(4)
1199 help_parts = rr_tmp%parts(3)
1201 det_ap_engen(help_parts) = det_ap_engen(help_parts)+tmp_energy
1205 if(rr_tmp%reac_type.eq.
rrt_ap)
then
1206 engen_q_ap = engen_q_ap + tmp_energy
1209 if (rr_tmp%parts(2) .eq.
ihe4)
then
1210 help_parts = rr_tmp%parts(1)
1212 help_parts = rr_tmp%parts(2)
1214 det_ap_engen(help_parts) = det_ap_engen(help_parts)+tmp_energy
1218 if(rr_tmp%reac_type.eq.
rrt_np)
then
1219 engen_q_np = engen_q_np + tmp_energy
1222 if (rr_tmp%parts(3) .eq.
ipro)
then
1223 help_parts = rr_tmp%parts(4)
1225 help_parts = rr_tmp%parts(3)
1227 det_pn_engen(help_parts) = det_pn_engen(help_parts)+tmp_energy
1231 if(rr_tmp%reac_type.eq.
rrt_pn)
then
1232 engen_q_pn = engen_q_pn + tmp_energy
1235 if (rr_tmp%parts(1) .eq.
ipro)
then
1236 help_parts = rr_tmp%parts(2)
1238 help_parts = rr_tmp%parts(1)
1240 det_pn_engen(help_parts) = det_pn_engen(help_parts)+tmp_energy
1244 if(rr_tmp%reac_type.eq.
rrt_an)
then
1245 engen_q_an = engen_q_an + tmp_energy
1248 if (rr_tmp%parts(1) .eq.
ihe4)
then
1249 help_parts = rr_tmp%parts(2)
1251 help_parts = rr_tmp%parts(1)
1253 det_an_engen(help_parts) = det_an_engen(help_parts)+tmp_energy
1257 if(rr_tmp%reac_type.eq.
rrt_na)
then
1258 engen_q_na = engen_q_na + tmp_energy
1261 if (rr_tmp%parts(3) .eq.
ihe4)
then
1262 help_parts = rr_tmp%parts(4)
1264 help_parts = rr_tmp%parts(3)
1266 det_an_engen(help_parts) = det_an_engen(help_parts)+tmp_energy
1270 if(rr_tmp%reac_type.eq.
rrt_o)
then
1272 det_other_engen(rr_tmp%parts(1)) = det_other_engen(rr_tmp%parts(1))+tmp_energy
1278 tmp_energy = rat*
y(rr_tmp%parts(1)) *
y(rr_tmp%parts(2)) &
1279 *
y(rr_tmp%parts(3)) * rr_tmp%q_value
1282 engen_q_value = engen_q_value + tmp_energy
1284 if(rr_tmp%reac_type.eq.
rrt_o)
then
1286 det_other_engen(rr_tmp%parts(1)) = det_other_engen(rr_tmp%parts(1))+tmp_energy
1292 tmp_energy = rat*
y(rr_tmp%parts(1)) *
y(rr_tmp%parts(2)) &
1293 *
y(rr_tmp%parts(3)) *
y(rr_tmp%parts(4))* rr_tmp%q_value
1295 engen_q_value = engen_q_value + tmp_energy
1297 if(rr_tmp%reac_type.eq.
rrt_o)
then
1299 det_other_engen(rr_tmp%parts(1)) = det_other_engen(rr_tmp%parts(1))+tmp_energy
1308 if (write_toplist)
then
1311 if (tmp_energy .ne. 0d0)
then
1312 toplist_energy(i) = tmp_energy
1313 toplist_mask(i) = .true.
1321 fission:
do i=1,
nfiss
1324 if (rat .eq. 0) cycle
1325 if ((fr_tmp%mode.eq.2).or.(fr_tmp%mode.eq.3))
then
1326 engen_q_value = engen_q_value + rat*
y(
fissrate(i)%fissparts(1)) *
fissrate(i)%averageQ
1327 engen_fiss = engen_fiss + rat*
y(
fissrate(i)%fissparts(1)) *
fissrate(i)%averageQ
1329 else if (fr_tmp%mode.eq.1)
then
1336 det_fiss_engen(
fissrate(i)%fissparts(1)) = det_fiss_engen(
fissrate(i)%fissparts(1))+tmp_energy
1346 if (
y(i).gt.1d-25)
then
1347 tmp_energy = -(isotope(i)%mass_exc)*(
y(i)-
y_p(i))
1348 engen_tot_exc = engen_tot_exc + tmp_energy
1350 if (write_toplist)
then
1351 toplist_energy_iso(i) = tmp_energy
1355 det_tot_engen(i) = -(isotope(i)%mass_exc)*(
y(i)-
y_p(i))
1361 entropy_src_term = 0.d0
1365 if (
rhob .gt. 1d-10)
then
1368 mue = etaele *
unit%k_mev *
t9*1d9
1372 mue = mue +
unit%mass_e
1374 kbt =
t9*1d9 *
unit%k_mev
1376 if (
y(i).gt.1d-25)
then
1377 mic2 = isotope(i)%mass*
unit%mass_u - isotope(i)%p_nr *
unit%mass_e + isotope(i)%mass_exc
1378 gi = 2d0*isotope(i)%spin + 1d0
1380 mui = -kbt*log( gi*pf(i) / (
rhob*
y(i)*
unit%n_a) * &
1381 (twopi_h2c2 * mic2 * kbt)**1.5d0 )
1383 entropy_src_term = entropy_src_term -(mic2 + mui + isotope(i)%p_nr*mue)*(
y(i)-
y_p(i))
1389 engen_q_value = engen_q_value /
unit%hix
1390 engen_q_weak = engen_q_weak /
unit%hix
1391 engen_q_ng = engen_q_ng /
unit%hix
1392 engen_q_gn = engen_q_gn /
unit%hix
1393 engen_q_pg = engen_q_pg /
unit%hix
1394 engen_q_gp = engen_q_gp /
unit%hix
1395 engen_q_ga = engen_q_ga /
unit%hix
1396 engen_q_ag = engen_q_ag /
unit%hix
1397 engen_q_np = engen_q_np /
unit%hix
1398 engen_q_pn = engen_q_pn /
unit%hix
1399 engen_q_na = engen_q_na /
unit%hix
1400 engen_q_an = engen_q_an /
unit%hix
1401 engen_q_pa = engen_q_pa /
unit%hix
1402 engen_q_ap = engen_q_ap /
unit%hix
1403 engen_fiss = engen_fiss /
unit%hix
1407 det_bet_engen = det_bet_engen /
unit%hix
1408 det_ag_engen = det_ag_engen /
unit%hix
1409 det_pg_engen = det_pg_engen /
unit%hix
1410 det_ng_engen = det_ng_engen /
unit%hix
1411 det_pn_engen = det_pn_engen /
unit%hix
1412 det_ap_engen = det_ap_engen /
unit%hix
1413 det_an_engen = det_an_engen /
unit%hix
1414 det_other_engen = det_other_engen /
unit%hix
1415 det_fiss_engen = det_fiss_engen /
unit%hix
1419 engen_tot_ncap = engen_q_ng + engen_q_gn
1420 engen_tot_pcap = engen_q_pg + engen_q_gp
1421 engen_tot_acap = engen_q_ag + engen_q_ga
1422 engen_tot_npcap= engen_q_np + engen_q_pn
1423 engen_tot_nacap= engen_q_na + engen_q_an
1424 engen_tot_pacap= engen_q_pa + engen_q_ap
1426 if (write_toplist)
then
1429 toplist_energy(:) = toplist_energy(:) /
unit%hix
1430 toplist_energy_iso(:) = toplist_energy_iso(:) /
stepsize/
unit%hix
1435 if (toplist_energy(1) .lt. toplist_energy(2))
then
1436 call raise_exception(
"There was a bug in the toplist. Energy is not sorted.", &
1437 "calc_nuclear_heating")
1442 tplist_str(:) =
"---------"
1446 if (j .gt.
nreac)
exit
1448 if (toplist_mask(toplist_idx(j)))
then
1451 tplist(i+1) = toplist_energy(j)
1460 if (tplist(1) .lt. tplist(2))
then
1461 call raise_exception(
"There was a bug in the toplist. Energy is not sorted.", &
1462 "calc_nuclear_heating")
1467 btlist_str(:) =
"---------"
1471 if (j .gt.
nreac)
exit
1473 if (toplist_mask(toplist_idx(
nreac+1-j)))
then
1476 btlist(i+1) = toplist_energy(
nreac+1-j)
1487 call quicksort(0,net_size,toplist_energy_iso,toplist_idx_iso)
1491 tplist_str_iso(:) =
"-----"
1495 if (j .gt. net_size)
exit
1497 if (toplist_energy_iso(j) .ne. 0d0)
then
1498 tplist_str_iso(i+1) =
net_names(toplist_idx_iso(j))
1499 tplist_iso(i+1) = toplist_energy_iso(j)
1510 btlist_str_iso(:) =
"-----"
1514 if (j .gt. net_size)
exit
1516 if (toplist_energy_iso(net_size+1-j) .ne. 0d0)
then
1517 btlist_str_iso(i+1) =
net_names(toplist_idx_iso(net_size+1-j))
1518 btlist_iso(i+1) = toplist_energy_iso(net_size+1-j)
1531 "(A,3X), 2X,"//int_to_str(
toplist_amount)//
"(es12.4,3X),' | ',"//&
1535 time,engen_tot_exc, tplist_str, tplist, btlist_str,btlist, tplist_str_iso, tplist_iso, &
1536 btlist_str_iso, btlist_iso
1567 engen_tot_pcap, engen_tot_acap,engen_tot_npcap, &
1568 engen_tot_nacap, engen_tot_pacap, engen_q_weak, &
1570 det_bet_engen, det_ag_engen, det_pg_engen, &
1571 det_ng_engen, det_pn_engen, det_ap_engen, &
1572 det_an_engen, det_other_engen, det_fiss_engen, &
1579 info_exit(
"calc_nuclear_heating")
1613 real(
r_kind),
intent(in) :: ctime
1614 real(
r_kind),
intent(in) :: t9
1615 logical,
optional,
intent(in) :: hdf5
1616 real(
r_kind),
dimension(net_size),
intent(in) :: y
1634 real(
r_kind) :: tau_alpha
1635 real(
r_kind) :: tau_nfiss
1636 real(
r_kind) :: tau_sfiss
1637 real(
r_kind) :: tau_bfiss
1640 real(
r_kind) :: ng_sum,gn_sum,pg_sum,gp_sum,beta_sum
1641 real(
r_kind) :: np_sum,pn_sum,an_sum,na_sum,ag_sum,ga_sum
1642 real(
r_kind) :: ap_sum,pa_sum
1643 real(
r_kind) :: alpha_sum, nfiss_sum, sfiss_sum, bfiss_sum
1649 real(
r_kind) :: q_sm = 1.e-20
1651 info_entry(
"calc_av_timescales")
1674 ysum = sum(y(max(1,
ihe4):net_size))
1682 if ((rr_tmp%group.eq.1).or.(rr_tmp%group.eq.2).or.(rr_tmp%group.eq.3) .or.(rr_tmp%group.eq.11))
then
1683 if ((rr_tmp%reac_type .eq.
rrt_betm).or.(rr_tmp%reac_type .eq.
rrt_betp).and.(rr_tmp%source.ne.
'ms99'))
then
1685 nuc = rr_tmp%parts(1)
1686 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1687 beta_sum = beta_sum + rat*y(nuc)
1693 if ((rr_tmp%reac_type.eq.
rrt_alpd))
then
1694 nuc = rr_tmp%parts(1)
1695 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1696 alpha_sum = alpha_sum + rat*y(nuc)
1701 if(rr_tmp%reac_type.eq.
rrt_ap)
then
1702 nuc = rr_tmp%parts(2)
1703 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1704 ap_sum = ap_sum + rat*y(nuc)*y(
ihe4)
1709 if(rr_tmp%reac_type.eq.
rrt_pa)
then
1710 nuc = rr_tmp%parts(2)
1711 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1712 pa_sum = pa_sum + rat*y(nuc)*y(
ipro)
1717 if(rr_tmp%reac_type.eq.
rrt_ng)
then
1718 nuc = rr_tmp%parts(2)
1719 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1720 ng_sum = ng_sum + rat*y(nuc)*y(
ineu)
1725 if(rr_tmp%reac_type.eq.
rrt_gn)
then
1726 nuc = rr_tmp%parts(1)
1727 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1728 gn_sum = gn_sum + rat*y(nuc)
1733 if(rr_tmp%reac_type.eq.
rrt_pg)
then
1734 nuc = rr_tmp%parts(2)
1735 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1736 pg_sum = pg_sum + rat*y(nuc)*y(
ipro)
1741 if(rr_tmp%reac_type.eq.
rrt_gp)
then
1742 nuc = rr_tmp%parts(1)
1743 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1744 gp_sum = gp_sum + rat*y(nuc)
1749 if(rr_tmp%reac_type.eq.
rrt_ag)
then
1750 nuc = rr_tmp%parts(2)
1751 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1752 ag_sum = ag_sum + rat*y(nuc)*y(
ihe4)
1757 if(rr_tmp%reac_type.eq.
rrt_ga)
then
1758 nuc = rr_tmp%parts(1)
1759 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1760 ga_sum = ga_sum + rat*y(nuc)
1765 if(rr_tmp%reac_type.eq.
rrt_np)
then
1766 nuc = rr_tmp%parts(2)
1767 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1768 np_sum = np_sum + rat*y(nuc)*y(
ineu)
1773 if(rr_tmp%reac_type.eq.
rrt_pn)
then
1774 nuc = rr_tmp%parts(2)
1775 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1776 pn_sum = pn_sum + rat*y(nuc)*y(
ipro)
1781 if(rr_tmp%reac_type.eq.
rrt_an)
then
1782 nuc = rr_tmp%parts(2)
1783 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1784 an_sum = an_sum + rat*y(nuc)*y(
ihe4)
1789 if(rr_tmp%reac_type.eq.
rrt_na)
then
1790 nuc = rr_tmp%parts(2)
1791 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1792 na_sum = na_sum + rat*y(nuc)*y(
ineu)
1798 fission:
do i=1,
nfiss
1802 if(fr_tmp%reac_type.eq.
rrt_nf)
then
1803 nuc = fr_tmp%fissparts(2)
1804 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1805 nfiss_sum = nfiss_sum + rat*y(nuc)*y(
ineu)
1810 if(fr_tmp%reac_type.eq.
rrt_sf)
then
1811 nuc = fr_tmp%fissparts(1)
1812 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1813 sfiss_sum = sfiss_sum + rat*y(nuc)
1818 if(fr_tmp%reac_type.eq.
rrt_bf)
then
1819 nuc = fr_tmp%fissparts(1)
1820 if(rat.gt.1.d-20 .and. y(nuc).gt.1.d-20)
then
1821 bfiss_sum = bfiss_sum + rat*y(nuc)
1827 tau_ng = ysum / ng_sum
1828 tau_gn = ysum / gn_sum
1829 tau_pg = ysum / pg_sum
1830 tau_gp = ysum / gp_sum
1831 tau_ag = ysum / ag_sum
1832 tau_ga = ysum / ga_sum
1833 tau_np = ysum / np_sum
1834 tau_pn = ysum / pn_sum
1835 tau_an = ysum / an_sum
1836 tau_na = ysum / na_sum
1837 tau_pa = ysum / pa_sum
1838 tau_ap = ysum / ap_sum
1839 tau_beta = ysum / beta_sum
1840 tau_alpha = ysum / alpha_sum
1841 tau_nfiss = ysum / nfiss_sum
1842 tau_sfiss = ysum / sfiss_sum
1843 tau_bfiss = ysum / bfiss_sum
1845 if (.not.
present(hdf5))
then
1847 write(
tsfile,
'(*(es23.16,3x))') &
1848 ctime, t9, tau_ga,tau_ag, tau_ng, tau_gn, tau_pg, tau_gp, tau_np, &
1849 tau_pn, tau_an, tau_na,tau_pa,tau_ap, tau_beta, tau_alpha, &
1850 tau_nfiss, tau_sfiss, tau_bfiss
1855 tau_gp, tau_np, tau_pn, tau_an, tau_na, &
1856 tau_ap, tau_pa, tau_beta, tau_alpha, &
1857 tau_nfiss, tau_sfiss, tau_bfiss)
1864 info_exit(
'calc_av_timescales')
1882 real(
r_kind),
dimension(net_size),
intent(in) :: y
1883 integer,
dimension(net_size) :: list
1884 real(
r_kind),
dimension(:),
allocatable :: xa,yz
1885 integer :: finfile,finsfile,finzfile
1886 integer :: mval,zmval
1888 integer :: mass,protonnum
1893 write(finfile,
'(A)')
"# 1:A 2:Z 3:N 4:Yi 5:Xi"
1899 if (y(i).le.1.d-20) cycle
1900 write(finfile,
'(3(i5,2x),2es23.14)')
isotope(i)%mass, &
1917 if(y(i).gt.1.d-20)
then
1918 xa(mass) = xa(mass) + y(i)*mass
1919 if (protonnum.gt.0) yz(protonnum) = yz(protonnum) + y(i)
1924 write(finsfile,
'(A)')
"# 1:A 2:Y(A) 3:X(A)"
1927 write(finsfile,
'(i5,2x,es23.14,2x,es23.14)')i, xa(i)/i ,xa(i)
1930 write(finzfile,
'(A)')
"# 1:Z 2:Y(Z)"
1932 write(finzfile,
'(i5,2x,es23.14)')i, yz(i)
1938 end subroutine finab
1956 integer,
intent(in) :: nsize
1957 integer,
dimension(nsize),
intent(out) :: list
1958 type(sorttype),
dimension(nsize) :: tmp
1960 integer :: amin,zmin
1963 info_entry(
'finab_sort')
1971 tmp(i)%chk = .false.
1976 if(tmp(j)%chk) cycle
1977 select case (tmp(j)%a-amin)
1983 if (tmp(j)%z.lt.zmin)
then
1993 tmp(nxt)%chk = .true.
1998 info_exit(
'finab_sort')
2019 info_entry(
'analysis_finalize')
2041 info_exit(
'analysis_finalize')