1 | !Notes:
|
---|
2 | !Check if have to add this module to some control list so that it becomes numerically part of the rest of the astrobear program (for instance, is there somehting thst needs to happen so that public variables in other parts of this code are understood by this new module? -- I remember having to add a Bonnor ebert module list/functions/ to some master list -- check this out, should be documented somewhere in the code)
|
---|
3 |
|
---|
4 | !Tests:
|
---|
5 | !check EOS stop check works when EOS dne Ideal gas law. - passed
|
---|
6 | !could plot T(MeV) in visit and Lnuc to check change in energy is consistent with this...?
|
---|
7 |
|
---|
8 | MODULE NeutrinoCooling
|
---|
9 | USE PhysicsDeclarations
|
---|
10 | USE GlobalDeclarations
|
---|
11 | USE EOS
|
---|
12 | ! REAL(KIND=qPrec), PUBLIC :: NeutrinoCoolPar(1:3)=(/1.d-8, 2d18, 1.9d25/)
|
---|
13 |
|
---|
14 | !!! This module requires the use of an ideal gas EOS
|
---|
15 |
|
---|
16 | CONTAINS
|
---|
17 |
|
---|
18 | !Finalize initialization for Neutrino Cooling
|
---|
19 | SUBROUTINE Init_NeutrinoCool
|
---|
20 |
|
---|
21 | IF(iEOS/=0) THEN
|
---|
22 | IF (MPI_ID == 0) THEN
|
---|
23 | PRINT *, '!!!Neutrino cooling not supported for EOS other than ideal gas law!!!'
|
---|
24 | PRINT *, '!!!Consider switching iEOS to 0 in physics.data!!!'
|
---|
25 | PRINT *, '!!!Simulation now stopping!!!'
|
---|
26 | STOP
|
---|
27 | END IF
|
---|
28 | END IF
|
---|
29 |
|
---|
30 | END SUBROUTINE Init_NeutrinoCool
|
---|
31 |
|
---|
32 | SUBROUTINE Neutrino_Cooling(q, dqdt, Temp)
|
---|
33 | REAL(KIND=qPrec) :: q(:),dqdt(:), Temp
|
---|
34 | dqdt(iE) = dqdt(iE) - NeutrinoCoolingStrength(q,Temp) !net previously calculated dqdt(ie)-neutrinocoolingrate
|
---|
35 | END SUBROUTINE Neutrino_Cooling
|
---|
36 |
|
---|
37 | FUNCTION NeutrinoCoolingStrength(q, Temp)
|
---|
38 | REAL(KIND=qPrec) :: Temp, q(:) !note Temp is the physical temperature of the gas -- in Kelvin****
|
---|
39 | REAL(KIND=qPrec) :: NeutrinoCoolingStrength
|
---|
40 | REAL(KIND=qPrec) :: Lnuc, Lpp !physical definitions?
|
---|
41 | REAL(KIND=qPREC) :: Kb_MeV_K !boltzmann constant in MeV/K
|
---|
42 | REAL(KIND=qPREC) :: T !temperature scaling dependence of cooling rate, in MeV
|
---|
43 | REAL(KIND=qPREC) :: rho !density scaling dependence of cooling rate, in g/cc
|
---|
44 |
|
---|
45 | Kb_MeV_K=8.61733035d-11
|
---|
46 | !T=Temp*tempscale !now in K -- this is wrong, Temp is already in Kelvin -- fixed 10/2/18
|
---|
47 | T=Temp*Kb_MeV_K !now in MeV
|
---|
48 | rho=q(1)*rscale !now in cgs
|
---|
49 | Lnuc=2d18*T**6d0 !erg/g/s
|
---|
50 | Lnuc=Lnuc*rho !erg/cm**3/s --- energy density (pressure)/time --- correct units for input into source solver
|
---|
51 | Lnuc=Lnuc*edotscale !edotscale = timescale/pscale; now in computational units
|
---|
52 | Lpp=1.9d25*T**9d0/rho !erg/g/s
|
---|
53 | Lpp=Lpp*rho*edotscale !compuational units
|
---|
54 |
|
---|
55 | NeutrinoCoolingStrength = (Lnuc + Lpp)
|
---|
56 |
|
---|
57 | ! IF (MPI_ID == 0) print *, 'NeutrinoCoolingStrength =', NeutrinoCoolingStrength
|
---|
58 | ! IF (MPI_ID == 0) print *, 'rho=q(1)*rscale =', rho
|
---|
59 | ! IF (MPI_ID == 0) print *, 'T=TempScale*Temp*Kb_MeV_K', T
|
---|
60 | ! IF (MPI_ID == 0) print *, 'Temp=', Temp, 'Temp/Tempscale=', temp/tempscale
|
---|
61 | ! IF (MPI_ID == 0) print *, 'Tempscale=', Tempscale, 'edotscale', edotscale
|
---|
62 |
|
---|
63 |
|
---|
64 | END FUNCTION NeutrinoCoolingStrength
|
---|
65 |
|
---|
66 | END MODULE NeutrinoCooling
|
---|