Blog: meeting notes: NeutrinoCooling.f90

File NeutrinoCooling.f90, 2.9 KB (added by Erica Kaminski, 6 years ago)
Line 
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
8MODULE 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
16CONTAINS
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
66END MODULE NeutrinoCooling