Skip to content

Commit 6312a35

Browse files
DebraheemDebraheem
committed
[ci skip] add fixed surface atm options
1 parent c25c7f2 commit 6312a35

File tree

1 file changed

+158
-26
lines changed

1 file changed

+158
-26
lines changed

‎star/private/hydro_vars.f90‎

Lines changed: 158 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -781,13 +781,19 @@ subroutine get_surf_PT( &
781781
real(dp) :: T_surf
782782
real(dp) :: P_surf_atm
783783
real(dp) :: P_surf
784+
real(dp) :: Pextra
785+
real(dp) :: kap_surf
786+
real(dp) :: M_surf
787+
784788

785789
include 'formats'
786790

787791
! Set up stellar surface parameters
788792

789793
L_surf = s% L(1)
790794
R_surf = s% r(1)
795+
kap_surf = s% opacity(1)
796+
M_surf = s% m(1)
791797
Teff = s% Teff
792798

793799
! Initialize partials
@@ -823,38 +829,164 @@ subroutine get_surf_PT( &
823829
endif
824830

825831
else
826-
827832
! Evaluate temperature and pressure based on atm_option
833+
! (yes, we might want to translate atm_option into an integer,
834+
! but these string comparisons are cheap)
835+
836+
! The first few are special, 'trivial-atmosphere' options
837+
838+
select case (s% atm_option)
839+
840+
case ('fixed_Teff')
841+
842+
! set Tsurf from Eddington T-tau relation
843+
! for current surface tau and Teff = `atm_fixed_Teff`.
844+
! set Psurf = Radiation_Pressure(Tsurf)
845+
846+
Teff = s% atm_fixed_Teff
847+
Teff4 = Teff*Teff*Teff*Teff
848+
T_surf = pow(3._dp/4._dp*Teff4*(tau_surf + 2._dp/3._dp), 0.25_dp)
849+
lnT_surf = log(T_surf)
850+
lnP_surf = Radiation_Pressure(T_surf)
851+
852+
if (.not. skip_partials) then
853+
dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
854+
dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
855+
endif
828856

829-
if (.false. .and. 1 == s% solver_test_partials_k .and. &
830-
s% solver_iter == s% solver_test_partials_iter_number) then
831-
star_debugging_atm_flag = .true.
832-
end if
857+
case ('fixed_Tsurf')
833858

834-
call get_atm_PT( &
835-
s, tau_surf, L_surf, R_surf, s% m(1), s% cgrav(1), skip_partials, &
836-
Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
837-
lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
838-
ierr)
839-
if (ierr /= 0) then
840-
if (s% report_ierr) then
841-
write(*,1) 'tau_surf', tau_surf
842-
write(*,1) 'L_surf', L_surf
843-
write(*,1) 'R_surf', R_surf
844-
write(*,1) 'Teff', Teff
845-
write(*,1) 's% m(1)', s% m(1)
846-
write(*,1) 's% cgrav(1)', s% cgrav(1)
847-
write(*,*) 'failed in get_atm_PT'
859+
! set Teff from Eddington T-tau relation for given
860+
! Tsurf and tau=2/3 set Psurf =
861+
! Radiation_Pressure(Tsurf)
862+
863+
T_surf = s% atm_fixed_Tsurf
864+
lnT_surf = log(T_surf)
865+
T_surf4 = T_surf*T_surf*T_surf*T_surf
866+
Teff = pow(4._dp/3._dp*T_surf4/(tau_surf + 2._dp/3._dp), 0.25_dp)
867+
lnP_surf = Radiation_Pressure(T_surf)
868+
869+
if (.not. skip_partials) then
870+
dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
871+
dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
872+
endif
873+
874+
case ('fixed_Psurf')
875+
876+
! set Teff from L and R using L = 4*pi*R^2*boltz_sigma*T^4.
877+
! set Tsurf using Eddington T-tau relation
878+
879+
if (L_surf < 0._dp) then
880+
if (s% report_ierr) then
881+
write(*,2) 'get_surf_PT: L_surf <= 0', s% model_number, L_surf
882+
call mesa_error(__FILE__,__LINE__)
883+
end if
884+
s% retry_message = 'L_surf < 0'
885+
ierr = -1
886+
return
848887
end if
849-
return
850-
end if
851888

852-
if (.false. .and. 1 == s% solver_test_partials_k .and. &
853-
s% solver_iter == s% solver_test_partials_iter_number) then
854-
s% solver_test_partials_val = atm_test_partials_val
855-
s% solver_test_partials_dval_dx = atm_test_partials_dval_dx
856-
end if
889+
lnP_surf = safe_log(s% atm_fixed_Psurf)
890+
if (R_surf <= 0._dp) then
891+
T_surf4 = 1._dp
892+
T_surf = 1._dp
893+
lnT_surf = 0._dp
894+
if (.not. skip_partials) then
895+
dlnT_dlnR = 0._dp
896+
dlnT_dL = 0._dp
897+
endif
898+
Teff = s% T(1)
899+
else
900+
Teff = pow(L_surf/(4._dp*pi*R_surf*R_surf*boltz_sigma), 0.25_dp)
901+
T_surf4 = 3d0/4d0*pow(Teff, 4.d0)*(tau_surf + 2._dp/3._dp)
902+
T_surf = pow(T_surf4, 0.25_dp)
903+
lnT_surf = log(T_surf)
904+
if (.not. skip_partials) then
905+
dlnT_dlnR = -0.5_dp
906+
dlnT_dL = 1._dp/(4._dp*L_surf)
907+
endif
908+
end if
909+
910+
if (.not. skip_partials) then
911+
dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
912+
dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
913+
endif
914+
915+
case ('fixed_Psurf_and_Tsurf')
916+
917+
lnP_surf = safe_log(s% atm_fixed_Psurf)
918+
T_surf = s% atm_fixed_Tsurf
919+
lnT_surf = log(T_surf)
920+
T_surf4 = T_surf*T_surf*T_surf*T_surf
921+
Teff = pow(4d0/3d0*T_surf4/(tau_surf + 2d0/3d0), 0.25d0)
922+
923+
if (.not. skip_partials) then
924+
dlnT_dL = 0; dlnT_dlnR = 0; dlnT_dlnM = 0; dlnT_dlnkap = 0
925+
dlnP_dL = 0; dlnP_dlnR = 0; dlnP_dlnM = 0; dlnP_dlnkap = 0
926+
endif
927+
928+
case default
929+
930+
! Everything else -- the 'non-trivial atmospheres' ---
931+
! gets passed to atm_support
932+
933+
if (.false. .and. 1 == s% solver_test_partials_k .and. &
934+
s% solver_iter == s% solver_test_partials_iter_number) then
935+
star_debugging_atm_flag = .true.
936+
end if
857937

938+
call get_atm_PT( &
939+
s, tau_surf, L_surf, R_surf, s% m(1), s% cgrav(1), skip_partials, &
940+
Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
941+
lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
942+
ierr)
943+
if (ierr /= 0) then
944+
if (s% report_ierr) then
945+
write(*,1) 'tau_surf', tau_surf
946+
write(*,1) 'L_surf', L_surf
947+
write(*,1) 'R_surf', R_surf
948+
write(*,1) 'Teff', Teff
949+
write(*,1) 's% m(1)', s% m(1)
950+
write(*,1) 's% cgrav(1)', s% cgrav(1)
951+
write(*,*) 'failed in get_atm_PT'
952+
end if
953+
return
954+
end if
955+
956+
if (.false. .and. 1 == s% solver_test_partials_k .and. &
957+
s% solver_iter == s% solver_test_partials_iter_number) then
958+
s% solver_test_partials_val = atm_test_partials_val
959+
s% solver_test_partials_dval_dx = atm_test_partials_dval_dx
960+
end if
961+
end select
962+
end if
963+
964+
! if using fixed surface, calculate Pextra.
965+
if (s% atm_option == 'fixed_Tsurf' .or. s% atm_option == 'fixed_Psurf_and_Tsurf' &
966+
.or. s% atm_option == 'fixed_Psurf' .or. s% atm_option == 'fixed_Teff') then
967+
! add extra pressure for fixed atmosphere options
968+
if (s% Pextra_factor /= 0._dp) then
969+
P_surf_atm = exp(lnP_surf)
970+
Pextra = s% Pextra_factor*(kap_surf/tau_surf)*(L_surf/M_surf)/(6._dp*pi*clight*s% cgrav(1))
971+
P_surf = P_surf_atm + Pextra
972+
if (P_surf < 1E-50_dp) then
973+
lnP_surf = -50*ln10
974+
if (.not. skip_partials) then
975+
dlnP_dL = 0._dp
976+
dlnP_dlnR = 0._dp
977+
dlnP_dlnM = 0._dp
978+
dlnP_dlnkap = 0._dp
979+
endif
980+
else
981+
lnP_surf = log(P_surf)
982+
if (.not. skip_partials) then
983+
dlnP_dL = dlnP_dL*P_surf_atm/P_surf
984+
dlnP_dlnR = dlnP_dlnR*P_surf_atm/P_surf
985+
dlnP_dlnM = dlnP_dlnM*P_surf_atm/P_surf
986+
dlnP_dlnkap = dlnP_dlnkap*P_surf_atm/P_surf
987+
endif
988+
end if
989+
end if
858990
end if
859991

860992
! Check outputs

0 commit comments

Comments
 (0)