copied stored variables and subroutines from GR1D so GR1D can be link…

src/nulibtable.F90

@@ -17,6 +17,8 @@ module nulibtable
   real*8, allocatable,save :: nulibtable_ewidths(:)
   real*8, allocatable,save :: nulibtable_ebottom(:)
   real*8, allocatable,save :: nulibtable_etop(:)
+  real*8, allocatable,save :: nulibtable_logenergies(:)
+  real*8, allocatable,save :: nulibtable_logetop(:)
 
   real*8, allocatable,save :: nulibtable_emissivities(:,:,:,:)
   real*8, allocatable,save :: nulibtable_absopacity(:,:,:,:)
@@ -640,3 +642,331 @@ subroutine nulibtable_epannihil_single_species_range_energy(xtemp,xeta, &
 
 end subroutine nulibtable_epannihil_single_species_range_energy
 
+
+!!!! ROUTINES FOR COMPATIBILITY WITH GR1D !!!!
+
+!this takes temp,eta, and return phi0/1 over energy (both in and out) and species range
+subroutine nulibtable_inelastic_range_species_range_energy2(xtemp,xeta, &
+     eas,eas_n1,eas_n2,eas_n3,eas_n4)
+
+  use nulibtable
+  implicit none
+
+  real*8, intent(in) :: xtemp, xeta !inputs
+  real*8 :: xltemp, xleta !log versions
+  integer, intent(in) :: eas_n1,eas_n2,eas_n3,eas_n4
+  real*8, intent(out) :: eas(eas_n1,eas_n2,eas_n3,eas_n4)
+  integer :: ins,ing_in,ing_out
+  real*8 :: xeas(eas_n1*eas_n2*(eas_n2+1)/2)
+  integer :: index
+  real*8 :: energy_conversion = 1.60217733d-6*5.59424238d-55
+
+
+  if(size(eas,1).ne.nulibtable_number_species) then
+     stop "nulibtable_inelastic_range_species_range_energy: supplied array dimensions (1) is not commensurate with table"
+  endif
+  if(size(eas,2).ne.nulibtable_number_groups) then
+     stop "nulibtable_inelastic_range_species_range_energy: supplied array dimensions (2) is not commensurate with table"
+  endif  
+  if(size(eas,3).ne.nulibtable_number_groups) then
+     stop "nulibtable_inelastic_range_species_range_energy: supplied array dimensions (3) is not commensurate with table"
+  endif
+  if(size(eas,4).ne.2) then
+     stop "nulibtable_inelastic_range_species_range_energy: supplied array dimensions (4) is not commensurate with table"
+  endif
+
+  xltemp = log10(xtemp)
+  xleta = log10(xeta)
+
+  if (xltemp.lt.nulibtable_logItemp_min) stop "temp below nulib inelastic table minimum temp"
+  if (xltemp.gt.nulibtable_logItemp_max) stop "temp above nulib inelastic table maximum temp"
+  if (xleta.lt.nulibtable_logIeta_min) stop "eta below nulib inelastic table minimum eta"
+  if (xleta.gt.nulibtable_logIeta_max) stop "eta above nulib inelastic table maximum eta"
+
+  xeas = 0.0d0
+  call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_Itable_Phi0,nulibtable_nItemp, &
+       nulibtable_nIeta,eas_n1*eas_n2*(eas_n2+1)/2,nulibtable_logItemp, &
+       nulibtable_logIeta)
+
+  index = 0
+  do ins=1,nulibtable_number_species
+     do ing_in=1,nulibtable_number_groups
+        do ing_out=1,ing_in
+           index = index + 1
+           eas(ins,ing_in,ing_out,1) = 10.0d0**xeas(index)
+        enddo
+     enddo
+     do ing_in=1,nulibtable_number_groups
+        do ing_out=ing_in+1,nulibtable_number_groups
+           eas(ins,ing_in,ing_out,1) =  &
+                exp(-(nulibtable_energies(ing_out)-nulibtable_energies(ing_in))/ &
+                (xtemp*energy_conversion))*eas(ins,ing_out,ing_in,1) !cernohorsky 94
+        enddo
+     enddo
+  enddo
+
+  xeas = 0.0d0
+  call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_Itable_Phi1,nulibtable_nItemp, &
+       nulibtable_nIeta,eas_n1*eas_n2*(eas_n2+1)/2,nulibtable_logItemp, &
+       nulibtable_logIeta)
+
+  index = 0
+  do ins=1,nulibtable_number_species
+     do ing_in=1,nulibtable_number_groups
+        do ing_out=1,ing_in
+           index = index + 1
+           !this way because we interpolate \phi1/\phi0
+           eas(ins,ing_in,ing_out,2) = xeas(index)*eas(ins,ing_in,ing_out,1)
+        enddo
+     enddo
+     do ing_in=1,nulibtable_number_groups
+        do ing_out=ing_in+1,nulibtable_number_groups
+           eas(ins,ing_in,ing_out,2) =  &
+                exp(-(nulibtable_energies(ing_out)-nulibtable_energies(ing_in))/ &
+                (xtemp*energy_conversion))*eas(ins,ing_out,ing_in,2) !cernohorsky 94
+        enddo
+     enddo
+  enddo  
+
+
+
+end subroutine nulibtable_inelastic_range_species_range_energy2
+
+!this takes temp,eta, and return phi0/1 over energy (both in and out) range for a single species
+subroutine nulibtable_inelastic_single_species_range_energy2(xtemp,xeta, &
+     lns,eas,eas_n1,eas_n2,eas_n3)
+
+  use nulibtable
+  implicit none
+
+  real*8, intent(in) :: xtemp, xeta !inputs
+  real*8 :: xltemp, xleta !log versions
+  integer, intent(in) :: lns,eas_n1,eas_n2,eas_n3
+  real*8, intent(out) :: eas(eas_n1,eas_n2,eas_n3)
+  integer :: ing_in,ing_out
+  real*8 :: xeas(eas_n1*(eas_n1+1)/2)
+  integer :: index
+  real*8 :: energy_conversion = 1.60217733d-6*5.59424238d-55
+  integer :: startindex,endindex
+
+  if(size(eas,1).ne.nulibtable_number_groups) then
+     stop "nulibtable_inelastic_single_species_range_energy: supplied array dimensions (1) is not commensurate with table"
+  endif  
+  if(size(eas,2).ne.nulibtable_number_groups) then
+     stop "nulibtable_inelastic_single_species_range_energy: supplied array dimensions (2) is not commensurate with table"
+  endif
+  if(size(eas,3).ne.2) then
+     stop "nulibtable_inelastic_single_species_range_energy: supplied array dimensions (3) is not commensurate with table"
+  endif
+
+  xltemp = log10(xtemp)
+  xleta = log10(xeta)
+
+  if (xltemp.lt.nulibtable_logItemp_min) stop "temp below nulib inelastic table minimum temp"
+  if (xltemp.gt.nulibtable_logItemp_max) stop "temp above nulib inelastic table maximum temp"
+  if (xleta.lt.nulibtable_logIeta_min) stop "eta below nulib inelastic table minimum eta"
+  if (xleta.gt.nulibtable_logIeta_max) stop "eta above nulib inelastic table maximum eta"
+
+  startindex = (lns-1)*nulibtable_number_groups*(nulibtable_number_groups+1)/2+1
+  endindex = startindex + nulibtable_number_groups*(nulibtable_number_groups+1)/2 - 1
+
+  xeas = 0.0d0
+  call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_Itable_Phi0(:,:,startindex:endindex), &
+       nulibtable_nItemp,nulibtable_nIeta,eas_n1*(eas_n1+1)/2, &
+       nulibtable_logItemp,nulibtable_logIeta)
+
+  index = 0
+  do ing_in=1,nulibtable_number_groups
+     do ing_out=1,ing_in
+        index = index + 1
+        eas(ing_in,ing_out,1) = 10.0d0**xeas(index)
+     enddo
+  enddo
+  do ing_in=1,nulibtable_number_groups
+     do ing_out=ing_in+1,nulibtable_number_groups
+        eas(ing_in,ing_out,1) =  &
+             exp(-(nulibtable_energies(ing_out)-nulibtable_energies(ing_in))/ &
+             (xtemp*energy_conversion))*eas(ing_out,ing_in,1) !cernohorsky 94
+     enddo
+  enddo
+
+  xeas = 0.0d0
+  call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_Itable_Phi1(:,:,startindex:endindex), &
+       nulibtable_nItemp,nulibtable_nIeta,eas_n1*(eas_n1+1)/2, &
+       nulibtable_logItemp,nulibtable_logIeta)
+
+  index = 0
+  do ing_in=1,nulibtable_number_groups
+     do ing_out=1,ing_in
+        index = index + 1
+        !this way because we interpolate \phi1/\phi0
+        eas(ing_in,ing_out,2) = xeas(index)*eas(ing_in,ing_out,1)
+     enddo
+  enddo
+  do ing_in=1,nulibtable_number_groups
+     do ing_out=ing_in+1,nulibtable_number_groups
+        eas(ing_in,ing_out,2) =  &
+             exp(-(nulibtable_energies(ing_out)-nulibtable_energies(ing_in))/ &
+             (xtemp*energy_conversion))*eas(ing_out,ing_in,2) !cernohorsky 94
+     enddo
+  enddo
+
+
+
+end subroutine nulibtable_inelastic_single_species_range_energy2
+
+!this takes temp,eta, and return phi0/1 for ep-annihilation over energy (both neutrinos) and species range
+subroutine nulibtable_epannihil_range_species_range_energy2(xtemp,xeta, &
+     eas,eas_n1,eas_n2,eas_n3,eas_n4)
+
+  use nulibtable
+  implicit none
+
+  real*8, intent(in) :: xtemp, xeta !inputs
+  real*8 :: xltemp, xleta !log versions
+  integer, intent(in) :: eas_n1,eas_n2,eas_n3,eas_n4
+  real*8, intent(out) :: eas(eas_n1,eas_n2,eas_n3,eas_n4)
+  integer :: ins,ing_this,ing_that
+  real*8 :: xeas(eas_n1*eas_n2*eas_n3*2)
+  integer :: index
+  real*8 :: energy_conversion = 1.60217733d-6*5.59424238d-55
+
+
+  if(size(eas,1).ne.nulibtable_number_species) then
+     stop "nulibtable_epannihil_range_species_range_energy: supplied array dimensions (1) is not commensurate with table"
+  endif
+  if(size(eas,2).ne.nulibtable_number_groups) then
+     stop "nulibtable_epannihil_range_species_range_energy: supplied array dimensions (2) is not commensurate with table"
+  endif  
+  if(size(eas,3).ne.nulibtable_number_groups) then
+     stop "nulibtable_enannihil_range_species_range_energy: supplied array dimensions (3) is not commensurate with table"
+  endif
+  if(size(eas,4).ne.4) then
+     stop "nulibtable_epannihil_range_species_range_energy: supplied array dimensions (4) is not commensurate with table"
+  endif
+
+  xltemp = log10(xtemp)
+  xleta = log10(xeta)
+
+  if (xltemp.lt.nulibtable_logItemp_min) stop "temp below nulib epannihil table minimum temp"
+  if (xltemp.gt.nulibtable_logItemp_max) stop "temp above nulib epannihil table maximum temp"
+  if (xleta.lt.nulibtable_logIeta_min) stop "eta below nulib epannihil table minimum eta"
+  if (xleta.gt.nulibtable_logIeta_max) stop "eta above nulib epannihil table maximum eta"
+
+  xeas = 0.0d0
+  call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_epannihiltable_Phi0,nulibtable_nItemp, &
+       nulibtable_nIeta,eas_n1*eas_n2*eas_n2*2,nulibtable_logItemp, &
+       nulibtable_logIeta)
+
+  index = 0
+  do ins=1,nulibtable_number_species
+     do ing_this=1,nulibtable_number_groups
+        do ing_that=1,nulibtable_number_groups
+           index = index + 1
+           eas(ins,ing_this,ing_that,1) = 10.0d0**xeas(index)
+           index = index + 1
+           eas(ins,ing_this,ing_that,2) = 10.0d0**xeas(index)
+!           eas(ins,ing_this,ing_that,1) = eas(ins,ing_this,ing_that,2)* &
+!                exp(-(nulibtable_energies(ing_this)+nulibtable_energies(ing_that))/(xtemp*energy_conversion))
+
+        enddo
+     enddo
+  enddo
+
+  xeas = 0.0d0
+  call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_epannihiltable_Phi1,nulibtable_nItemp, &
+       nulibtable_nIeta,eas_n1*eas_n2*eas_n2*2,nulibtable_logItemp, &
+       nulibtable_logIeta)
+
+  index = 0
+  do ins=1,nulibtable_number_species
+     do ing_this=1,nulibtable_number_groups
+        do ing_that=1,nulibtable_number_groups
+           !this way because we interpolate \phi1/\phi0
+           index = index + 1
+           eas(ins,ing_this,ing_that,3) = xeas(index)*eas(ins,ing_this,ing_that,1)
+           index = index + 1
+           eas(ins,ing_this,ing_that,4) = xeas(index)*eas(ins,ing_this,ing_that,2)
+!           eas(ins,ing_this,ing_that,3) = eas(ins,ing_this,ing_that,4)* &
+!                exp(-(nulibtable_energies(ing_this)+nulibtable_energies(ing_that))/(xtemp*energy_conversion))
+        enddo
+     enddo
+  enddo  
+
+
+
+end subroutine nulibtable_epannihil_range_species_range_energy2
+
+!this takes temp,eta, and return epannihilation phi0/1 over energy (both nus) range for a single species
+subroutine nulibtable_epannihil_single_species_range_energy2(xtemp,xeta, &
+     lns,eas,eas_n1,eas_n2,eas_n3)
+
+  use nulibtable
+  implicit none
+
+  real*8, intent(in) :: xtemp, xeta !inputs
+  real*8 :: xltemp, xleta !log versions
+  integer, intent(in) :: lns,eas_n1,eas_n2,eas_n3
+  real*8, intent(out) :: eas(eas_n1,eas_n2,eas_n3)
+  integer :: ing_this,ing_that
+  real*8 :: xeas(eas_n1*eas_n2*2)
+  integer :: index
+  real*8 :: energy_conversion = 1.60217733d-6*5.59424238d-55
+  integer :: startindex,endindex
+
+  if(size(eas,1).ne.nulibtable_number_groups) then
+     stop "nulibtable_epannihil_single_species_range_energy: supplied array dimensions (1) is not commensurate with table"
+  endif  
+  if(size(eas,2).ne.nulibtable_number_groups) then
+     stop "nulibtable_epannihil_single_species_range_energy: supplied array dimensions (2) is not commensurate with table"
+  endif
+  if(size(eas,3).ne.4) then
+     stop "nulibtable_epannihil_single_species_range_energy: supplied array dimensions (3) is not commensurate with table"
+  endif
+
+  xltemp = log10(xtemp)
+  xleta = log10(xeta)
+
+  if (xltemp.lt.nulibtable_logItemp_min) stop "temp below nulib epannihil table minimum temp"
+  if (xltemp.gt.nulibtable_logItemp_max) stop "temp above nulib epannihil table maximum temp"
+  if (xleta.lt.nulibtable_logIeta_min) stop "eta below nulib epannihil table minimum eta"
+  if (xleta.gt.nulibtable_logIeta_max) stop "eta above nulib epannihil table maximum eta"
+
+  startindex = (lns-1)*nulibtable_number_groups*nulibtable_number_groups*2+1
+  endindex = startindex + nulibtable_number_groups*nulibtable_number_groups*2 - 1
+
+  xeas = 0.0d0
+  call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_epannihiltable_Phi0(:,:,startindex:endindex), &
+       nulibtable_nItemp,nulibtable_nIeta,eas_n1*eas_n2*2,nulibtable_logItemp,nulibtable_logIeta)
+
+  index = 0
+  do ing_this=1,nulibtable_number_groups
+     do ing_that=1,nulibtable_number_groups
+        index = index + 1
+        eas(ing_this,ing_that,1) = 10.0d0**xeas(index)
+        index = index + 1
+        eas(ing_this,ing_that,2) = 10.0d0**xeas(index)
+!        eas(ing_this,ing_that,1) = eas(ing_this,ing_that,2)* &
+!             (exp(-(nulibtable_energies(ing_this)+nulibtable_energies(ing_that))/(xtemp*energy_conversion)))
+     enddo
+  enddo
+
+  xeas = 0.0d0
+  call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_epannihiltable_Phi1(:,:,startindex:endindex), &
+       nulibtable_nItemp,nulibtable_nIeta,eas_n1*eas_n2*2,nulibtable_logItemp,nulibtable_logIeta)
+
+  index = 0
+  do ing_this=1,nulibtable_number_groups
+     do ing_that=1,nulibtable_number_groups
+        index = index + 1
+        !this way because we interpolate \phi1/\phi0
+        eas(ing_this,ing_that,3) = xeas(index)*eas(ing_this,ing_that,1)
+        index = index + 1
+        eas(ing_this,ing_that,4) = xeas(index)*eas(ing_this,ing_that,2)
+!        eas(ing_this,ing_that,3) = eas(ing_this,ing_that,4)*&
+!             (exp(-(nulibtable_energies(ing_this)+nulibtable_energies(ing_that))/(xtemp*energy_conversion)))
+     enddo
+  enddo
+
+end subroutine nulibtable_epannihil_single_species_range_energy2
+

src/nulibtable_reader.F90

@@ -85,7 +85,9 @@ subroutine nulibtable_reader(filename,include_Ielectron,include_epannihil_kernel
   allocate(nulibtable_ewidths(nulibtable_number_groups))
   allocate(nulibtable_ebottom(nulibtable_number_groups))
   allocate(nulibtable_etop(nulibtable_number_groups))
-
+  allocate(nulibtable_logenergies(nulibtable_number_groups))
+  allocate(nulibtable_logetop(nulibtable_number_groups))
+
   call h5dopen_f(file_id, "nrho",dset_id, error)
   call h5dread_f(dset_id, H5T_NATIVE_INTEGER, nulibtable_nrho, dims1, error)
   call h5dclose_f(dset_id, error)