diff options
Diffstat (limited to 'src/nuc_eos/nuc_eos.F90')
| -rw-r--r-- | src/nuc_eos/nuc_eos.F90 | 730 |
1 files changed, 0 insertions, 730 deletions
diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90 deleted file mode 100644 index 65b0562..0000000 --- a/src/nuc_eos/nuc_eos.F90 +++ /dev/null @@ -1,730 +0,0 @@ -! ######################################################### -! -! Copyright C. D. Ott and Evan O'Connor, July 2009 -! -! UNITS: density g/cm^3 -! temperature MeV -! ye number fraction per baryon -! energy erg/g -! pressure dyn/cm^2 -! chemical potentials MeV -! entropy k_B / baryon -! cs2 cm^2/s^2 (not relativistic) -! -! keyerr --> error output; should be 0 -! rfeps --> root finding relative accuracy, set around 1.0d-10 -! keytemp: 0 -> coming in with eps -! 1 -> coming in with temperature -! 2 -> coming in with entropy -! - -#include "cctk.h" - -subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& - xdpderho,xdpdrhoe,xxa,xxh,xxn,xxp,xabar,xzbar,xmu_e,xmu_n,xmu_p, & - xmuhat,keytemp,keyerr,rfeps) - - use eosmodule - implicit none - - real*8, intent(inout) :: xrho,xye - real*8, intent(inout) :: xtemp,xenr,xent - real*8, intent(out) :: xprs,xcs2,xdedt - real*8, intent(out) :: xdpderho,xdpdrhoe,xxa,xxh,xxn,xxp - real*8, intent(out) :: xabar,xzbar,xmu_e,xmu_n,xmu_p,xmuhat - real*8, intent(in) :: rfeps - integer, intent(in) :: keytemp - integer, intent(out) :: keyerr - - - ! local variables - real*8 :: lr,lt,y,xx,xeps,leps,xs,xpressure - real*8 :: d1,d2,d3 - real*8 :: ff(nvars) - integer :: keyerrt - integer :: keyerrr - - keyerrt = 0 - keyerrr = 0 - - if(xrho.gt.eos_rhomax) then - keyerr = 103 - return - endif - - if(xrho.lt.eos_rhomin) then - keyerr = 104 - return - endif - - if(xye.gt.eos_yemax) then - keyerr = 101 - return - endif - - if(xye.lt.eos_yemin) then - keyerr = 102 - return - endif - - if(keytemp.eq.1) then - if(xtemp.gt.eos_tempmax) then - keyerr = 105 - return - endif - - if(xtemp.lt.eos_tempmin) then - keyerr = 106 - return - endif - endif - - lr = log10(xrho) - lt = log10(max(xtemp,eos_tempmin)) - y = xye - xeps = xenr + energy_shift - - keyerr = 0 - - if(keytemp.eq.0) then - !need to find temperature based on xeps - if(xeps .gt. 0.0d0) then - leps = log10(xeps) - call findtemp(lr,lt,y,leps,keyerrt,rfeps) - else - keyerrt = 667 - endif - if(keyerrt.ne.0) then - call CCTK_WARN (0, "Did not find temperature") - endif - xtemp = 10.0d0**lt - - elseif(keytemp.eq.2) then - !need to find temperature based on xent - xs = xent - call findtemp_entropy(lr,lt,y,xs,keyerrt,rfeps) - xtemp = 10.0d0**lt - - elseif(keytemp.eq.3) then - !need to find rho based on xprs - xpressure = log10(xprs) - call findrho_press(lr,lt,y,xpressure,keyerrr,rfeps) - if(keyerrr.ne.0) then - write(*,*) "Problem in findrho_press:", keyerr - keyerr = keyerrr - return - endif - xrho = 10.0d0**lr - - endif - - ! have temperature, proceed: - call findall(lr,lt,y,ff) - - !unless we want xprs to be constant (keytemp==3), reset xprs - if(.not.keytemp.eq.3) then - xprs = 10.0d0**ff(1) - endif - - if(.not.keytemp.eq.0) then - xenr = 10.0d0**ff(2) - energy_shift - endif - - if(.not.keytemp.eq.2) then - xent = ff(3) - endif - - xcs2 = ff(5) - -! derivatives - xdedt = ff(6) - - xdpdrhoe = ff(7) - - xdpderho = ff(8) - -! chemical potentials - xmuhat = ff(9) - - xmu_e = ff(10) - - xmu_p = ff(11) - - xmu_n = ff(12) - -! compositions - xxa = ff(13) - - xxh = ff(14) - - xxn = ff(15) - - xxp = ff(16) - - xabar = ff(17) - - xzbar = ff(18) - - -end subroutine nuc_eos_full - -subroutine nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) - - implicit none - real*8 xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe - real*8,parameter:: polyK1 = 5.0d14 - real*8,parameter :: polyGamma = 1.33333333333333d0 - integer keytemp - - xprs=polyK1*xrho**(polyGamma) - xenr=xprs/xrho/(polyGamma-1.d0) - xcs2=polyGamma*xprs/xrho - xdpderho=0.0d0 - xdpdrhoe=0.0d0 - - xenr=max(8.0d18,xenr) - xenr=min(3.0d20,xenr) - -end subroutine nuc_poly_eos - -subroutine nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) - - implicit none - real*8 xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe - real*8,parameter :: idealK1 = 4.93483302030614d14 -! real*8,parameter :: idealK1 = 1.2435d15 * (0.5d0**(4.d0/3.d0)) - real*8,parameter :: idealgamma = 1.41d0 - integer keytemp - - if(keytemp.eq.1) then -! energy wanted - xprs=idealK1*xrho**(idealgamma) - xenr=xprs/xrho/(idealgamma-1.d0) - xcs2=idealgamma*xprs/xrho - endif - - xprs = (idealgamma - 1.0d0) * xrho * xenr - - xdpderho = (idealgamma - 1.0d0 ) * xrho - xdpdrhoe = (idealgamma - 1.0d0 ) * xenr - - xcs2= xdpdrhoe+xdpderho*xprs/xrho**2 - -end subroutine nuc_low_eos - -subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& - xdpderho,xdpdrhoe,xmunu,keytemp,keyerr,rfeps) - - use eosmodule - implicit none - - real*8, intent(inout) :: xrho,xye - real*8, intent(inout) :: xtemp,xenr,xent - real*8, intent(in) :: rfeps - real*8, intent(out) :: xprs,xmunu,xcs2,xdedt - real*8, intent(out) :: xdpderho,xdpdrhoe - integer, intent(in) :: keytemp - integer, intent(out) :: keyerr - - ! local variables - real*8 :: lr,lt,y,xt,xx,xeps,leps,xs,xpressure - real*8 :: d1,d2,d3,ff(8) - integer :: keyerrt - integer :: keyerrr - - keyerrt = 0 - keyerrr = 0 - - if(xrho.gt.eos_rhomax) then - keyerr = 103 - return - endif - - if(xrho.lt.eos_rhomin*1.126d0) then - call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) - xent = 0.0d0 - return - endif - - if(xye.gt.eos_yemax) then - keyerr = 101 - return - endif - - if(xye.lt.eos_yemin) then - keyerr = 102 - return - endif - - if(keytemp.eq.1) then - if(xtemp.gt.eos_tempmax) then - call CCTK_WARN (0, "nuc_eos: temp > tempmax") - endif - - if(xtemp.lt.eos_tempmin) then - keyerr = 106 - return - endif - endif - - lr = log10(xrho) - lt = log10(max(xtemp,eos_tempmin)) - y = xye - - keyerr = 0 - - if(keytemp.eq.0) then - !need to find temperature based on xeps - xeps = xenr + energy_shift - if(xeps .gt. 0.0d0) then - leps = log10(xeps) - call findtemp(lr,lt,y,leps,keyerrt,rfeps) - else - keyerrt = 667 - endif - if(keyerrt.eq.667) then - ! turns out, we can still converge, but - ! too a bogus, possibly negative temperature - xt = xtemp - call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps) - if(keyerrt.eq.0) then - keyerrt = 668 - xtemp = xt - else - keyerrt = 667 - keyerr = keyerrt - return - endif - keyerr = keyerrt - else - xtemp = 10.0d0**lt - endif - elseif(keytemp.eq.2) then - !need to find temperature based on xent - xs = xent - call findtemp_entropy(lr,lt,y,xs,keyerrt,rfeps) - if(keyerrt.ne.0) then - keyerr = keyerrt - return - endif - xtemp = 10.0d0**lt - - elseif(keytemp.eq.3) then - !need to find rho based on xprs - xpressure = log10(xprs) - call findrho_press(lr,lt,y,xpressure,keyerrr,rfeps) - if (keyerrr.ne.0) then - keyerr = keyerrr - write(*,*) "Problem in findrho_press:", keyerr - return - endif - xrho = 10.0d0**lr - endif - - ! have rho,temp,ye; proceed: - if(keyerr.ne.668) then - call findall_short(lr,lt,y,ff) - else - call findall_short_lowT(lr,xt,y,ff) - endif - - !unless we want xprs to be constant (keytemp==3), reset xprs - if(.not.keytemp.eq.3) then - xprs = 10.0d0**ff(1) - endif - - !unless we want xenr to be constant (keytemp==0), reset xenr - if(.not.keytemp.eq.0) then - xenr = 10.0d0**ff(2) - energy_shift - endif - - !unless we want xent to be constant (keytemp==2), reset xent - if(.not.keytemp.eq.2) then - xent = ff(3) - endif - - xmunu = ff(4) - - xcs2 = ff(5) - - xdedt = ff(6) - - xdpdrhoe = ff(7) - - xdpderho = ff(8) - -end subroutine nuc_eos_short - - -subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& - keytemp,keyerr,rfeps) - - use eosmodule - implicit none - - real*8, intent(in) :: xrho,xye - real*8, intent(inout) :: xtemp,xenr - real*8, intent(in) :: rfeps - real*8, intent(out) :: xprs - integer, intent(in) :: keytemp - integer, intent(out) :: keyerr - - ! local variables - real*8 :: xcs2,xdpderho,xdpdrhoe - real*8 :: lr,lt,y,xt,xx,xeps,leps,xs - real*8 :: d1,d2,d3,ff(8) - integer :: keyerrt - character(len=256) :: warnline - - keyerrt = 0 - - if(xrho.gt.eos_rhomax) then - keyerr = 103 - return - endif - - if(xrho.lt.eos_rhomin*1.126d0) then - call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) - return - endif - - if(xye.gt.eos_yemax) then - keyerr = 101 - return - endif - - if(xye.lt.eos_yemin) then - keyerr = 102 - return - endif - - if(keytemp.eq.1) then - if(xtemp.gt.eos_tempmax) then - call CCTK_WARN (0, "nuc_eos: temp > tempmax") - endif - - if(xtemp.lt.eos_tempmin) then - keyerr = 106 - return - endif - endif - - keyerr = 0 - - if(keytemp.gt.1) then - call CCTK_WARN (0, "eos_nuc_press does not support keytemp other than 0 and 1") - endif - - lr = log10(xrho) - lt = log10(max(xtemp,eos_tempmin)) - y = xye - if(keytemp.eq.0) then - !need to find temperature based on xeps - xeps = xenr + energy_shift - if(xeps .gt. 0.0d0) then - leps = log10(xeps) - call findtemp(lr,lt,y,leps,keyerrt,rfeps) - else - keyerrt = 667 - endif -! if(xenr .lt. -9.00725347549139d20) then -! if(xenr.lt.-2.28659899d20) then -! write(warnline,"(A4,i5,1P10E15.6)") "UGGx",keyerrt,lr,lt,xt,xtemp,y,xeps,xenr -! call CCTK_WARN(1,warnline) -! endif - if(keyerrt.eq.667) then - ! turns out, we can still converge, but - ! too a bogus, possibly negative temperature - xt = xtemp - call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps) - if(keyerrt.eq.0) then - keyerrt = 668 - if(xeps.lt.0.0d0) then - keyerrt = 668 ! impossible to do anything - endif - xtemp = xt - else - keyerrt = 667 - keyerr = keyerrt - return - endif -! if(xenr .lt. -9.00725347549139d20) then -! if(xenr.lt.-2.28659899d20) then -! write(warnline,"(A4,i5,1P10E15.6)") "UGG",keyerrt,lr,lt,xt,xtemp,y,xeps,xenr -! call CCTK_WARN(0,warnline) -! endif - keyerr = keyerrt - else - xtemp = 10.0d0**lt - endif - endif - ! have temperature, proceed: - if(keyerr.ne.668) then - call findall_press_eps(lr,lt,y,ff) - else - call findall_press_eps_lowT(lr,xt,y,ff) - endif - xprs = 10.0d0**ff(1) - xenr = 10.0d0**ff(2) - energy_shift - -end subroutine nuc_eos_press_eps - -subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpderho,& - xdpdrhoe,& - keytemp,keyerr,rfeps) - - use eosmodule - implicit none - - real*8, intent(in) :: xrho,xye - real*8, intent(inout) :: xtemp,xenr - real*8, intent(in) :: rfeps - real*8, intent(out) :: xdpdrhoe,xdpderho - integer, intent(in) :: keytemp - integer, intent(out) :: keyerr - - ! local variables - real*8 :: xcs2,xprs - real*8 :: lr,lt,y,xt,xx,xeps,leps,xs - real*8 :: d1,d2,d3,ff(2) - integer :: keyerrt - - keyerrt = 0 - - if(xrho.gt.eos_rhomax) then - keyerr = 103 - return - endif - - if(xrho.lt.eos_rhomin*1.126d0) then - call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) - return - endif - - if(xye.gt.eos_yemax) then - keyerr = 101 - return - endif - - if(xye.lt.eos_yemin) then - keyerr = 102 - return - endif - - if(keytemp.eq.1) then - if(xtemp.gt.eos_tempmax) then - call CCTK_WARN (0, "nuc_eos: temp > tempmax") - endif - - if(xtemp.lt.eos_tempmin) then - keyerr = 106 - return - endif - endif - - lr = log10(xrho) - lt = log10(max(xtemp,eos_tempmin)) - y = xye - - keyerr = 0 - if(keytemp.eq.0) then - !need to find temperature based on xeps - xeps = xenr + energy_shift - if(xeps .gt. 0.0d0) then - leps = log10(xeps) - call findtemp(lr,lt,y,leps,keyerrt,rfeps) - else - keyerrt = 667 - endif - if(keyerrt.eq.667) then - ! turns out, we can still converge, but - ! too a bogus, possibly negative temperature - xt = xtemp - call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps) - if(keyerrt.eq.0) then - keyerrt = 668 - xtemp = xt - else - keyerrt = 667 - keyerr = keyerrt - return - endif - keyerr = keyerrt - else - xtemp = 10.0d0**lt - endif - elseif(keytemp.gt.1) then - call CCTK_WARN (0, "eos_nuc_press does not support keytemp > 1") - endif - - ! have temperature, proceed: - if(keyerr.ne.668) then - call findall_dpdr_dpde(lr,lt,y,ff) - else - call findall_dpdr_dpde_lowT(lr,xtemp,y,ff) - endif - - xdpdrhoe = ff(1) - xdpderho = ff(2) - - if(keytemp.eq.1.and.keyerr.ne.668) then - call findthis(lr,lt,y,xeps,alltables(:,:,:,2),d1,d2,d3) - xeps = 10.0d0**xeps - energy_shift - xenr = xeps - endif - -end subroutine nuc_eos_dpdr_dpde - - -subroutine findthis(lr,lt,y,value,array,d1,d2,d3) - - use eosmodule - - implicit none - - integer rip,rim - integer tip,tim - integer yip,yim - - real*8 lr,lt,y,value,d1,d2,d3 - real*8 array(*) - -! Ewald's interpolator - call intp3d(lr,lt,y,value,1,array,nrho,ntemp,nye,logrho,logtemp,ye,d1,d2,d3) - -end subroutine findthis - - -subroutine findall(lr,lt,y,ff) - - use eosmodule - implicit none - - real*8 ff(nvars) - real*8 ffx(nvars,1) - real*8 lr,lt,y - integer i - -! Ewald's interpolator - call intp3d_many(lr,lt,y,ffx,1,alltables,& - nrho,ntemp,nye,nvars,logrho,logtemp,ye) - ff(:) = ffx(:,1) - - -end subroutine findall - - -subroutine findall_short(lr,lt,y,ff) - - use eosmodule - implicit none - - real*8 ffx(8,1) - real*8 ff(8) - real*8 lr,lt,y - integer i - integer :: nvarsx = 8 - - -! Ewald's interpolator - call intp3d_many(lr,lt,y,ffx,1,alltables(:,:,:,1:8), & - nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) - ff(:) = ffx(:,1) - -end subroutine findall_short - -subroutine findall_short_lowT(lr,t,y,ff) - - use eosmodule - implicit none - - real*8 ffx(8,1) - real*8 ff(8) - real*8 lr,t,y - integer i - integer :: nvarsx = 8 - - -! Ewald's interpolator - call intp3d_many_linearTlow(lr,t,y,ffx,1,alltables(:,:,:,1:8), & - nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) - ff(:) = ffx(:,1) - -end subroutine findall_short_lowT - - -subroutine findall_press_eps(lr,lt,y,ff) - - use eosmodule - implicit none - - real*8 ffx(2,1) - real*8 ff(2) - real*8 lr,lt,y - integer i - integer :: nvarsx = 2 - - -! Ewald's interpolator - call intp3d_many(lr,lt,y,ffx,1,alltables(:,:,:,1:2), & - nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) - ff(:) = ffx(:,1) - -end subroutine findall_press_eps - -subroutine findall_press_eps_lowT(lr,t,y,ff) - - use eosmodule - implicit none - - real*8 ffx(2,1) - real*8 ff(2) - real*8 lr,t,y - integer i - integer :: nvarsx = 2 - -! Ewald's interpolator - call intp3d_many_linearTlow(lr,t,y,ffx,1,alltables(:,:,:,1:2), & - nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) - ff(:) = ffx(:,1) - -end subroutine findall_press_eps_lowT - - -subroutine findall_dpdr_dpde(lr,lt,y,ff) - - use eosmodule - implicit none - - real*8 ffx(2,1) - real*8 ff(2) - real*8 lr,lt,y - integer i - integer :: nvarsx = 2 - - -! Ewald's interpolator - call intp3d_many(lr,lt,y,ffx,1,alltables(:,:,:,7:8), & - nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) - ff(:) = ffx(:,1) - -end subroutine findall_dpdr_dpde - -subroutine findall_dpdr_dpde_lowT(lr,t,y,ff) - - use eosmodule - implicit none - - real*8 ffx(2,1) - real*8 ff(2) - real*8 lr,t,y - integer i - integer :: nvarsx = 2 - - -! Ewald's interpolator - call intp3d_many_linearTlow(lr,t,y,ffx,1,alltables(:,:,:,7:8), & - nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) - ff(:) = ffx(:,1) - -end subroutine findall_dpdr_dpde_lowT |