Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018-2019 Radek Smolec & The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module rsp_build
21 : use star_def, only: star_info
22 : use utils_lib, only: is_bad
23 : use const_def, only: dp, crad
24 : use eos_lib, only: Radiation_Pressure
25 : use rsp_def
26 : use rsp_eval_eos_and_kap, only: X, Y, Z
27 : use rsp_lina, only: mesa_eos_kap, do_LINA
28 : use rsp_relax_env, only: EOP, RELAX_ENV
29 :
30 : implicit none
31 :
32 : private
33 : public :: do_rsp_build
34 :
35 : real(dp) :: PREC,FSUB,TIN,CFIDDLE,ALF, &
36 : HHFAC,DdmFAC,EFL02,EMR,ELR, &
37 : E_0,E_1,T_0,T_1,V_0,V_1,P_0,P_1,QQ_0,QQ_1, &
38 : CP_0,CP_1,OP_0,OP_1,R_1,M_0,dm_bar_0
39 : real(dp), dimension(15) :: PERS,ETO
40 : real(dp), pointer, dimension(:) :: &
41 : M, DM, DM_BAR, R, Vol, T, w, Et, E, P, Lr, Lc, Hp_face, Y_face, K, CPS, QQS
42 : logical, parameter :: RSP_eddi = .true. ! use Eddington approx at surface
43 :
44 : contains
45 :
46 0 : subroutine do_rsp_build(s,ierr)
47 : !use rsp_create_env, only: do_rsp_create_env
48 : type (star_info), pointer :: s
49 : integer, intent(out) :: ierr
50 : integer :: NZT ! warst. z joniz./il. warstw
51 0 : real(dp) :: TH0 ! temp. warstwy jonizacji
52 0 : real(dp) :: Mass,L
53 0 : real(dp) :: H,dmN
54 : integer :: NMODES ! ilosc modow rozwazanych (N=NMODES)
55 : integer :: NDIM1,NDIM2 ! maks. ilosc modow/warstw
56 : real(dp) :: VEL0(15)
57 : integer :: I,J,kk,NSEQ
58 0 : real(dp) :: GEFF,MBOL
59 : character (len=250) :: FILENAME
60 :
61 : integer :: IO,II
62 0 : real(dp) :: SS, AA, BB
63 0 : real(dp), allocatable :: TA(:), VEL(:,:), TEMP(:)
64 0 : real(dp) :: TAUTEFF, TAUATTEFF
65 : logical :: RELAX
66 0 : real(dp) :: amix1, amix2
67 :
68 0 : ierr = 0
69 :
70 0 : if (.not. associated(M)) then
71 : allocate( &
72 : M(NZN+1), DM(NZN+1), DM_BAR(NZN+1), T(NZN+1), E(NZN+1), w(NZN+1), &
73 : Et(NZN+1), P(NZN+1), R(NZN+1), Vol(NZN+1), &
74 : Lr(NZN+1), Lc(NZN+1), Hp_face(NZN+1), Y_face(NZN+1), K(NZN+1), &
75 0 : CPS(NZN+1), QQS(NZN+1))
76 : end if
77 :
78 0 : allocate(TA(NZN+1),VEL(NZN+1,15),TEMP(NZN))
79 :
80 0 : TH0 = s% RSP_T_anchor
81 0 : TIN = s% RSP_T_inner
82 0 : NZT = s% RSP_nz_outer
83 0 : FSUB = s% RSP_dq_1_factor
84 0 : ddmfac = 1d0 ! s% RSP_ddmfac
85 0 : hhfac = 1.02d0 ! s% RSP_hhfac
86 0 : nmodes = s% RSP_nmodes
87 0 : RELAX = s% RSP_relax_initial_model
88 0 : EFL02 = EFL0*EFL0
89 :
90 : ! START MAIN CYCLE
91 :
92 : ! INITIALIZE ARRAYS
93 0 : NDIM1=15 !max. number of modes
94 0 : NDIM2=NZN+1
95 0 : do J=1,NDIM1
96 0 : do I=1,NDIM2
97 0 : VEL(I,J)=0.d0
98 : end do
99 : end do
100 0 : do I=1,NDIM1
101 0 : PERS(I)= 0.d0
102 : VEL0(I)= 0.d0
103 0 : ETO(I) = 0.d0
104 : !SGR(I) = 0.d0
105 : end do
106 0 : do I=1,NDIM2
107 0 : R(I) = 0.d0
108 0 : P(I) = 0.d0
109 0 : Vol(I) = 0.d0
110 0 : E(I) = 0.d0
111 0 : T(I) = 0.d0
112 0 : K(I)= 0.d0
113 0 : M(I)= 0.d0
114 0 : dm(I)= 0.d0
115 0 : dm_bar(I)= 0.d0
116 0 : w(I)=0.d0
117 0 : Hp_face(I) = 0.d0
118 0 : Y_face(I) = 0.d0
119 0 : CPS(I) = 0.d0
120 0 : QQS(I) = 0.d0
121 0 : Lc(I) = 0.d0
122 0 : Lr(I) = 0.d0
123 0 : w(I) = 0.d0
124 : end do
125 :
126 : ! SET STANDARD PARAMETERS
127 : ! NSEQ = SEQUENCE NUMBER (DUMMY)
128 0 : NSEQ = 1
129 :
130 0 : CFIDDLE = 0.02d0 !0.02
131 0 : ALF = 1.0d-6
132 : ! PRECISIONS
133 0 : PREC = 1d-10
134 :
135 0 : EMR = s% RSP_mass
136 0 : ELR = s% RSP_L
137 0 : TE = s% RSP_Teff
138 0 : Mass=EMR*SUNM
139 0 : L=ELR*SUNL
140 :
141 0 : if (s% RSP_trace_RSP_build_model) then
142 0 : write(*,*) '*** build initial model ***'
143 0 : write(*,'(a9,f15.5)') 'M/Msun', EMR
144 0 : write(*,'(a9,f15.5)') 'L/Lsun', ELR
145 0 : write(*,'(a9,f15.5)') ' Teff', TE
146 : end if
147 :
148 0 : call STAH(s,Mass,L,TE,H,dmN,TH0,NZT,NZN,ierr)
149 0 : if (ierr /= 0) return
150 :
151 : ! INITIAL GUESS FOR w=E_T (NOW DEFINED AT THE ZONE)
152 0 : TEMP(1)=0.d0
153 0 : do I=2,NZN
154 : !TEMP(I)=0.5d0*(w(I)**2+w(I-1)**2)
155 : !TEMP(I)=sqrt(w(I)**2*w(I-1)**2)
156 : TEMP(I)=(w(I)**2*dm(I)+w(I-1)**2*dm(I-1))/ &
157 0 : (dm(I)+dm(I-1))
158 : end do
159 0 : do I=1,NZN
160 0 : w(I)=TEMP(I)
161 : ! LINE BELOW MUST BE PRESENT IF MAIN VARIABLE IS E_T
162 : ! (DERIVATIVES ~1/sqrt(E_T) AND WITH w=0, LAPACK
163 : ! PROBLEM ARISES; NO SUCH PROBLEM WHEN E_T^2 IS USED
164 : ! AND LINE BELOW IS NOT NECESSARY THEN)
165 0 : if(w(I)<=EFL02) w(I)=EFL02
166 : !write(*,*) NZN-I+1, sqrt(w(i))
167 : end do
168 :
169 : ! NOTE: w(I) now holds Et = w**2
170 : ! watch out
171 :
172 0 : if(RELAX) then
173 : call RELAX_ENV(s, L, TH0, TE, NZT, NZN, &
174 0 : M, DM, DM_BAR, R, Vol, T, w, ierr)
175 0 : if (ierr /= 0) return
176 : end if
177 :
178 : ! optical depth. just for output.
179 0 : TA(NZN) = s% tau_factor*s% tau_base
180 0 : SS=TA(NZN)
181 0 : do I=NZN,2,-1
182 0 : SS=SS+K(I)*(R(I)-R(I-1))/Vol(I)
183 0 : TA(I-1)=SS
184 : end do
185 0 : II=0; IO=0
186 0 : do I=NZN,1,-1
187 0 : if(TA(I)>2.d0/3.d0)then
188 0 : II=I
189 0 : IO=I+1
190 0 : GOTO 77
191 : end if
192 : end do
193 : 77 continue
194 0 : if (IO == 0) then
195 0 : write(*,*) 'failed to find photosphere'
196 0 : ierr = -1
197 0 : return
198 : end if
199 0 : AA=(T(IO)-T(II))/(TA(IO)-TA(II))
200 0 : BB=T(IO)-AA*TA(IO)
201 0 : TAUTEFF=AA*2.d0/3.d0+BB
202 :
203 0 : do I=NZN,1,-1
204 : if(T(I)>TE)then
205 : II=I
206 : IO=I+1
207 : GOTO 78
208 : end if
209 : end do
210 : 78 continue
211 0 : AA=(T(IO)-T(II))/(TA(IO)-TA(II))
212 0 : BB=T(IO)-AA*TA(IO)
213 0 : TAUATTEFF=(TE-BB)/AA
214 :
215 0 : call cleanup_for_LINA(s, M, DM, DM_BAR, R, Vol, T, w, P, ierr)
216 :
217 : 1 format(1X,1P,5D26.16)
218 :
219 0 : GEFF=G*Mass/R(NZN)**2
220 0 : MBOL=-2.5d0*dlog10(ELR)+4.79d0
221 :
222 0 : if(NMODES==0) GOTO 11 ! jesli masz liczyc tylko static envelope
223 :
224 0 : if (.not. (s% use_RSP_new_start_scheme .or. s% use_other_RSP_linear_analysis)) then
225 0 : if (s% RSP_trace_RSP_build_model) write(*,*) '*** linear analysis ***'
226 0 : do I=1,NZN ! LINA changes Et, so make a work copy for it
227 0 : Et(I) = w(I)
228 : end do
229 : call do_LINA(s, L, NZN, NMODES, VEL, PERS, ETO, &
230 0 : M, DM, DM_BAR, R, Vol, T, Et, Lr, ierr)
231 0 : if (ierr /= 0) then
232 : return
233 : end if
234 0 : FILENAME=trim(s% log_directory) // '/' // 'LINA_period_growth.data'
235 0 : open(15,file=trim(FILENAME),status='unknown')
236 0 : write(*,'(a)') ' P(days) growth'
237 0 : do I=1,NMODES
238 0 : s% rsp_LINA_periods(i) = PERS(I)
239 0 : s% rsp_LINA_growth_rates(i) = ETO(I)
240 0 : write(*,'(I3,2X,99e16.5)') I-1, &
241 0 : PERS(I)/86400.d0,ETO(I)
242 0 : write(15,'(I3,2X,99e16.5)') I-1, &
243 0 : PERS(I)/86400.d0,ETO(I)
244 : end do
245 0 : close(15)
246 0 : s% RSP_have_set_velocities = .true.
247 : else
248 0 : PERS(1:NMODES) = 0d0
249 0 : VEL0(1:NMODES) = 0d0
250 0 : do I=1,NZN
251 0 : do j=1,nmodes
252 0 : VEL(I,J) = 0d0
253 : end do
254 : end do
255 : end if
256 :
257 : 11 continue
258 :
259 : 5568 format(1X,1P,5E15.6)
260 :
261 : 444 format(F6.3,tr2,f8.2,tr2,f7.2,tr2,d9.3)
262 0 : if (s% RSP_trace_RSP_build_model) then
263 0 : write(*,*) '*** done creating initial model ***'
264 0 : write(*,'(A)')
265 : end if
266 : ! recall that w is actually Et = w**2 at this point
267 0 : call set_build_vars(s,M,DM,DM_BAR,R,Vol,T,w,Lr,Lc)
268 :
269 0 : IWORK=0
270 0 : TEFF = TE
271 0 : RSTA = R(NZN)
272 0 : do I=1,NZN
273 0 : kk = NZN+1-i
274 0 : s% L(kk)=L
275 0 : s% L_start(kk)=0.0d0 !SMOLEC!
276 0 : s% v(kk)=0d0
277 : end do
278 0 : s% L_center=L
279 0 : if(ALFA==0.d0) EFL0=0.d0
280 0 : s% rsp_period=s% RSP_default_PERIODLIN
281 0 : if (is_bad(s% rsp_period)) then
282 0 : write(*,1) 'rsp_period', s% rsp_period
283 0 : call mesa_error(__FILE__,__LINE__,'rsp_build read_model')
284 : end if
285 0 : amix1 = s% RSP_fraction_1st_overtone
286 0 : amix2 = s% RSP_fraction_2nd_overtone
287 0 : if((AMIX1+AMIX2)>1.d0) write(*,*) 'AMIX DO NOT ADD UP RIGHT'
288 0 : if (.not. s% use_RSP_new_start_scheme) then
289 0 : PERIODLIN=PERS(s% RSP_mode_for_setting_PERIODLIN+1)
290 0 : s% rsp_period=PERIODLIN
291 0 : s% v_center = 0d0
292 0 : do I=1,NZN
293 : s% v(NZN+1-i)=1.0d5*s% RSP_kick_vsurf_km_per_sec* &
294 0 : ((1.0d0-AMIX1-AMIX2)*VEL(I,1)+AMIX1*VEL(I,2)+AMIX2*VEL(I,3))
295 : end do
296 : end if
297 :
298 0 : end subroutine do_rsp_build
299 :
300 :
301 0 : subroutine STAH(s,MX,L,TE,H0,dmN0,TH0,NZT,NZN,ierr)
302 : ! INTEGRATE STATIC ENVELOPE
303 : ! CONVECTIVE FLUX INCLUDED (WUCHTERL & FEUCHTINGER 1998)
304 : ! TURBULENT PRESSERE AND OVERSHOOTING NEGLECTED
305 : ! (ALPHA_P = ALPHA_T = 0)
306 : ! DIFFUSION APPROXIMATION FOR RADIATIVE TRANSFER
307 : ! HYDROGEN ZONE DEPTH (NZT, IN ZONES), AND TEMPERATURE(TH0) FIXED
308 : ! ARGUMENTS.. M(GM), L(CGS), TE,TH0(DEG K)
309 :
310 : type (star_info), pointer :: s
311 : real(dp), intent(in) :: MX,L,TE,TH0
312 : real(dp), intent(out) :: H0,dmN0
313 : integer, intent(in) :: NZT
314 : integer, intent(inout) :: NZN
315 : integer, intent(out) :: ierr
316 :
317 0 : real(dp) :: dmN,dm_0,H,Psurf,DDT
318 0 : real(dp) :: GPF
319 0 : real(dp) :: F2,F1,D,HH,TT,dmL
320 0 : real(dp) :: T4_1,RM,T4_0,WE,dmNL
321 0 : real(dp) :: FACQ,HH1,HH2
322 : integer :: N,N1,I,ITIN,dmN_cnt,NCHANG,IG,H_cnt
323 0 : real(dp) :: HP_0,HP_1,IGR_0,IGR_1,w_0
324 0 : real(dp) :: Lr_0,Lc_0,SVEL_0,HSTART,tau_sum,TH0_tol,TIN_tol, &
325 0 : dmN_too_large, dmN_too_small, H_too_large, H_too_small
326 : logical :: adjusting_dmN, in_photosphere, in_outer_env, &
327 : have_dmN_too_large, have_dmN_too_small, &
328 : have_H_too_large, have_H_too_small, have_T
329 :
330 : include 'formats'
331 : ierr = 0
332 0 : IG = 0
333 0 : HH1 = 1.005d0 !1.005
334 0 : HH2 = 1.005d0 !1.005
335 0 : NCHANG = 0
336 0 : HSTART = 1.01d0 ! s% RSP_hstart
337 0 : FACQ = 1.005d0 ! 1.005
338 0 : ITIN = 1
339 0 : tau_sum = 0
340 0 : adjusting_dmN = .true.
341 0 : in_photosphere = .true.
342 0 : in_outer_env = .true.
343 0 : have_dmN_too_large = .false.
344 0 : have_dmN_too_small = .false.
345 0 : have_H_too_large = .false.
346 0 : have_H_too_small = .false.
347 :
348 0 : TH0_tol = s% RSP_T_anchor_tolerance
349 0 : TIN_tol = s% RSP_T_inner_tolerance
350 :
351 0 : call ZNVAR(s,H,dmN,L,TE,MX,ierr)
352 0 : if (ierr /= 0) return
353 0 : dmN_cnt = 1
354 0 : H_cnt = 1
355 :
356 0 : start_from_top_loop: do
357 0 : if (s% RSP_trace_RSP_build_model) write(*,*) 'call setup_outer_zone'
358 0 : call setup_outer_zone(ierr)
359 0 : if (ierr /= 0) return
360 0 : if (s% RSP_trace_RSP_build_model) write(*,*) 'call store_N'
361 0 : call store_N
362 0 : zone_loop: do
363 0 : R_1=pow(R_1**3-3.d0*V_0*dm_0/P4,1.d0/3.d0)
364 0 : N=N-1
365 0 : if (s% RSP_trace_RSP_build_model) write(*,*) 'zone_loop', N, T_0, TIN
366 0 : if (N==0 .or. T_0 >= TIN) then
367 0 : if (s% RSP_trace_RSP_build_model) write(*,*) 'call next_H'
368 0 : call next_H ! sets HH
369 0 : if (N==0 .and. abs(T(1)-TIN)<TIN*TIN_tol) then
370 0 : s% M_center = M_0 - dm_0
371 0 : s% star_mass = s% RSP_mass
372 0 : s% mstar = s% star_mass*SUNM
373 0 : s% xmstar = s% mstar - s% M_center
374 0 : s% M_center = s% mstar - s% xmstar ! this is how it is set when read file
375 0 : s% L_center = s% RSP_L*SUNL
376 0 : s% R_center = pow(r(1)**3 - Vol(1)*dm(1)/P43, 1d0/3d0)
377 0 : s% v_center = 0
378 0 : if (s% RSP_trace_RSP_build_model) &
379 0 : write(*,*) ' inner dm growth scale', HH
380 : exit start_from_top_loop ! done
381 : end if
382 0 : if (H_cnt >= s% RSP_max_inner_scale_tries) then
383 0 : write(*,*) 'failed to find inner dm scaling to satisfy tolerance for T_inner'
384 0 : write(*,*) 'you might try increasing RSP_T_inner_tolerance'
385 0 : ierr = -1
386 0 : return
387 : !call mesa_error(__FILE__,__LINE__)
388 : end if
389 0 : if (s% RSP_trace_RSP_build_model) write(*,*) 'call prepare_for_new_H'
390 0 : call prepare_for_new_H
391 0 : cycle zone_loop
392 : end if
393 0 : if (s% RSP_trace_RSP_build_model) write(*,*) 'call setup_next_zone'
394 0 : call setup_next_zone
395 : ! pick T to make Lr + Lc = L
396 0 : if (s% RSP_trace_RSP_build_model) write(*,*) 'call get_T'
397 0 : have_T = get_T(ierr)
398 0 : if (ierr /= 0) return
399 0 : if (.not. have_T) then
400 0 : call failed
401 0 : DDT = dmN/1.d3
402 0 : dmN = dmN-DDT
403 0 : cycle start_from_top_loop
404 : end if
405 0 : if (s% RSP_trace_RSP_build_model) write(*,*) 'call get_V'
406 0 : call get_V(ierr)
407 0 : if (ierr /= 0) return
408 : if((NZN-N+1==NZT .or. T_0 >= TH0) &
409 0 : .and. adjusting_dmN) then
410 0 : if (s% RSP_trace_RSP_build_model) &
411 0 : write(*,*) 'call next_dmN', dmN_cnt, NZN-N+1, T_0, TH0, abs(T_0-TH0), TH0_tol*TH0
412 0 : if (dmN_cnt >= s% RSP_max_outer_dm_tries) then
413 0 : write(*,*) 'failed to find outer dm to satisfy tolerance for T_anchor'
414 0 : write(*,*) 'you might try increasing RSP_T_anchor_tolerance'
415 0 : ierr = -1
416 0 : return
417 : !call mesa_error(__FILE__,__LINE__)
418 : end if
419 0 : call next_dmN
420 0 : if(NZN-N+1==NZT.and.abs(T_0-TH0)<TH0_tol*TH0) then
421 0 : adjusting_dmN = .false.
422 0 : if (s% RSP_trace_RSP_build_model) &
423 0 : write(*,*) ' outer dm/Msun', dmN/SUNM
424 : end if
425 : ! one last repeat with final dmN
426 : !if (s% RSP_trace_RSP_build_model) write(*,*) 'cycle start_from_top_loop', dmN_cnt, dmN/SUNM
427 : cycle start_from_top_loop
428 : end if
429 0 : if (s% RSP_trace_RSP_build_model) write(*,*) 'call store_N'
430 0 : call store_N
431 : end do zone_loop
432 : end do start_from_top_loop
433 :
434 0 : s% R_center=R_1; H0=H; dmN0=dmN
435 0 : if(N/=0) call change_NZN
436 :
437 : contains
438 :
439 : subroutine report_location_of_photosphere
440 : real(dp) :: tau, dtau
441 : integer :: i
442 : include 'formats'
443 : tau = 0d0 ! surface currently at tau=0
444 : do i=NZN,1,-1
445 : dtau = dm(i)*K(i)/(4d0*pi*r(i)**2)
446 : tau = tau + dtau
447 : if (tau >= 2d0/3d0) then
448 : write(*,2) 'cells from surface tau=0 to tau=2/3', &
449 : NZN+1-i, tau-dtau, 2d0/3d0, tau, T(i), TE
450 : return
451 : end if
452 : end do
453 : end subroutine report_location_of_photosphere
454 :
455 0 : subroutine setup_outer_zone(ierr)
456 : use rsp_eval_eos_and_kap, only: get_surf_P_T_kap
457 : integer, intent(out) :: ierr
458 0 : real(dp) :: tau_surf, kap_guess, T_surf, kap_surf, Teff_atm
459 : include 'formats'
460 0 : dm_0=dmN*FSUB
461 0 : M_0=MX
462 0 : dm_bar_0=(dm_0/2.d0)
463 : if(.not.RSP_eddi) then ! EXACT GREY RELATION
464 : WE=TE**4
465 : T4_0=WE*sqrt(3.d0)/4.d0 !0.4330127018d0
466 : T_0= pow(sqrt(3.d0)/4.d0,0.25d0)*TE !0.811194802d0*TE
467 : else ! EDDINGTON APPROXIMATION
468 0 : WE=TE**4
469 0 : T4_0=WE*0.5d0 ! T4_0=WE*1.0d0/2.d0
470 0 : T_0=pow(0.5d0, 0.25d0)*TE ! T_0= pow(1.0d0/2.d0,0.25d0)*TE
471 : end if
472 0 : RM=sqrt(L/(P4*SIG*WE))
473 0 : R_1=RM
474 0 : if (s% RSP_use_atm_grey_with_kap_for_Psurf) then
475 0 : tau_surf = s% RSP_tau_surf_for_atm_grey_with_kap
476 0 : kap_guess = 1d-2
477 : call get_surf_P_T_kap(s, &
478 : MX, RM, L, tau_surf, kap_guess, &
479 0 : T_surf, Psurf, kap_surf, Teff_atm, ierr)
480 : !write(*,*) 'Psurf from atm, Prad_surf', Psurf, crad*T_0*T_0*T_0*T_0/3d0
481 0 : if (ierr /= 0) then
482 0 : write(*,*) 'failed in get_surf_P_T_kap'
483 0 : return
484 : end if
485 0 : else if (s% RSP_use_Prad_for_Psurf) then
486 0 : Psurf = crad*T_0*T_0*T_0*T_0/3d0
487 : else
488 0 : Psurf = 0d0
489 : end if
490 0 : Psurf_from_atm = Psurf
491 0 : P_0=Psurf+G*M_0*dm_bar_0/(P4*R_1**4)
492 : call EOP(s,-1,T_0,P_0,V_0, &
493 0 : E_0,CP_0,QQ_0,SVEL_0,OP_0,ierr)
494 0 : if (ierr /= 0) return
495 0 : w_0=0.d0
496 0 : Lc_0=0.d0
497 0 : Lr_0=L
498 0 : N=NZN
499 0 : end subroutine setup_outer_zone
500 :
501 0 : subroutine get_V(ierr)
502 : integer, intent(out) :: ierr
503 : ierr = 0
504 0 : T_0=sqrt(sqrt(T4_0))
505 : call EOP(s,N,T_0,P_0,V_0, &
506 0 : E_0,CP_0,QQ_0,SVEL_0,OP_0,ierr)
507 0 : if (ierr /= 0) return
508 0 : if(N/=NZN)then
509 0 : call CFLUX(HP_0,IGR_0,Lc_0,w_0,GPF,N)
510 0 : if(Lc_0>=L) then
511 0 : write(*,*) 'trouble!',I
512 0 : stop
513 : end if
514 : end if
515 0 : end subroutine get_V
516 :
517 0 : subroutine next_dmN
518 : real(dp), parameter :: search_factor = 2d0
519 0 : dmN_cnt = dmN_cnt+1
520 0 : dmNL=dmN
521 0 : if (T_0 < TH0) then ! dmN is too small
522 0 : dmN_too_small = dmN
523 0 : have_dmN_too_small = .true.
524 0 : if (.not. have_dmN_too_large) then
525 0 : dmN = dmN * search_factor
526 0 : DDT = dmN - dmNL
527 0 : return
528 : end if
529 : else ! T_0 > TH0, dmN is too large
530 0 : dmN_too_large = dmN
531 0 : have_dmN_too_large = .true.
532 0 : if (.not. have_dmN_too_small) then
533 0 : dmN = dmN / search_factor
534 0 : DDT = dmN - dmNL
535 0 : return
536 : end if
537 : end if
538 : ! search using bounds
539 : ! just bisect since for dmN too large, stop short of target cell
540 0 : dmN = 0.5d0*(dmN_too_large + dmN_too_small)
541 0 : DDT = dmN - dmNL
542 : end subroutine next_dmN
543 :
544 0 : subroutine next_H ! same scheme as next_dmN. bound and bisect.
545 : real(dp), parameter :: search_factor = 1.05d0
546 : real(dp) :: HH_prev
547 0 : HH_prev = HH
548 0 : H_cnt = H_cnt+1
549 : !write(*,*) 'next_H have_H_too_large', have_H_too_large, H_too_large
550 : !write(*,*) 'next_H have_H_too_small', have_H_too_small, H_too_small
551 0 : if (T_0 < TIN) then ! H is too small
552 0 : H_too_small = H
553 0 : have_H_too_small = .true.
554 0 : if (.not. have_H_too_large) then
555 0 : HH = H * search_factor
556 : !write(*,*) 'too small next_H HH, HH_prev', HH, HH_prev
557 0 : return
558 : end if
559 : else ! T_0 > TIN, H is too large
560 0 : H_too_large = H
561 0 : have_H_too_large = .true.
562 0 : if (.not. have_H_too_small) then
563 0 : HH = H / search_factor
564 : !write(*,*) 'too large next_H HH, HH_prev', HH, HH_prev
565 0 : return
566 : end if
567 : end if
568 : ! search using bounds. keep it simple.
569 : ! just bisect since for H too large, stop short of target cell
570 0 : HH = 0.5d0*(H_too_large + H_too_small)
571 : !write(*,*) 'next_H HH, HH_prev', HH, HH_prev
572 : !if (abs(HH - HH_prev) < 1d-6*HH) call mesa_error(__FILE__,__LINE__,'next_H')
573 : end subroutine next_H
574 :
575 0 : subroutine store_N
576 0 : real(dp) :: dtau
577 0 : R(N) = R_1
578 0 : P(N) = P_0
579 0 : Vol(N) = V_0
580 0 : E(N) = E_0
581 0 : T(N) = T_0
582 0 : K(N) = OP_0
583 0 : M(N) = M_0
584 0 : dm(N) = dm_0
585 0 : dm_bar(N) = dm_bar_0
586 0 : Hp_face(N) = HP_0
587 0 : Y_face(N)= IGR_0
588 0 : CPS(N) = CP_0
589 0 : QQS(N) = QQ_0
590 0 : Lc(N) = Lc_0
591 0 : Lr(N) = Lr_0
592 0 : w(N) = w_0
593 0 : dtau = dm(N)*K(N)/(P4*R(N)**2)
594 0 : if (in_photosphere .and. T(N) >= TE) then
595 0 : if (s% RSP_testing) &
596 0 : write(*,*) 'nz phot, tau_sum, T', &
597 0 : NZN-N, tau_sum, tau_sum+dtau, T(N+1), TE, T(N)
598 0 : in_photosphere = .false.
599 : end if
600 0 : tau_sum = tau_sum + dtau
601 0 : end subroutine store_N
602 :
603 0 : subroutine setup_next_zone
604 0 : P_1 = P_0
605 0 : V_1 = V_0
606 0 : OP_1 = OP_0
607 0 : T4_1 = T4_0
608 0 : M_0 = M_0-dm_0
609 0 : dmL = dm_0 !dmL IS dm IN A PREVIOUS ZONE
610 0 : HP_1 = HP_0
611 0 : IGR_1 = IGR_0
612 0 : CP_1 = CP_0
613 0 : QQ_1 = QQ_0
614 0 : T_1 = sqrt(sqrt(T4_1))
615 0 : E_1 = E_0
616 : ! RESET dm FOR THE INNER ZONES
617 0 : if(N==NZN-1) then
618 0 : dm_0=dmN
619 0 : else if (T_1 > TH0) then
620 0 : if (in_outer_env) then
621 0 : in_outer_env = .false.
622 0 : if (s% RSP_testing) write(*,*) 'nz outer', NZN-N, T_1, TH0
623 : end if
624 0 : dm_0=dm_0*H
625 : end if
626 0 : dm_bar_0=(dm_0+dmL)/2.d0
627 0 : P_0=P_1+G*M_0*dm_bar_0/(P4*R_1**4)
628 0 : end subroutine setup_next_zone
629 :
630 0 : real(dp) function eval_T_residual(ierr)
631 : integer, intent(out) :: ierr
632 : call EOP(s,0, &
633 0 : T_0,P_0,V_0,E_0,CP_0,QQ_0,SVEL_0,OP_0,ierr)
634 0 : if (ierr /= 0) return
635 0 : call CFLUX(HP_0,IGR_0,Lc_0,w_0,GPF,N)
636 0 : TT=4.d0*SIG*P4**2*R_1**4/(3.d0*dm_bar_0*L)
637 0 : T4_0 = T_0**4
638 : Lr_0=TT*(T4_0/OP_0-T4_1/OP_1)/ &
639 0 : (1.d0-dlog(OP_0/OP_1)/dlog(T4_0/T4_1))*L
640 0 : eval_T_residual = (Lr_0 + Lc_0)/L - 1d0
641 0 : end function eval_T_residual
642 :
643 0 : real(dp) function get_T_residual(lnT, dfdx, lrpar, rpar, lipar, ipar, ierr)
644 : ! returns with ierr = 0 if was able to evaluate f and df/dx at x
645 : ! if df/dx not available, it is okay to set it to 0
646 : use const_def, only: dp
647 : integer, intent(in) :: lrpar, lipar
648 : real(dp), intent(in) :: lnT
649 : real(dp), intent(out) :: dfdx
650 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
651 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
652 : integer, intent(out) :: ierr
653 0 : ierr = 0
654 0 : dfdx = 0
655 0 : T_0 = exp(lnT)
656 0 : get_T_residual = eval_T_residual(ierr)
657 0 : end function get_T_residual
658 :
659 0 : logical function get_T_estimate()
660 : use num_lib, only: safe_root_with_brackets
661 0 : real(dp) :: Tmax, epsx, epsy, residual, lnT, &
662 0 : lnT_min, lnT_max, resid_T_min, resid_T_max, dfdx
663 : integer :: n, ierr
664 : integer, parameter :: lrpar=1, lipar=1, imax=50
665 0 : real(dp), pointer :: rpar(:)
666 0 : integer, pointer :: ipar(:)
667 : include 'formats'
668 0 : Tmax = 0.99d0*(3d0*P_0/crad)**0.25d0 ! Prad must be < P_0
669 0 : lnT_min = log(T_1)
670 0 : lnT_max = log(Tmax)
671 : resid_T_min = get_T_residual(lnT_min, &
672 0 : dfdx, lrpar, rpar, lipar, ipar, ierr)
673 : resid_T_max = get_T_residual(lnT_max, &
674 0 : dfdx, lrpar, rpar, lipar, ipar, ierr)
675 0 : epsx = 1d-5
676 0 : epsy = 1d-5
677 : ierr = 0
678 : lnT = safe_root_with_brackets( &
679 : get_T_residual, lnT_min, lnT_max, &
680 : resid_T_min, resid_T_max, &
681 : imax, epsx, epsy, &
682 0 : lrpar, rpar, lipar, ipar, ierr)
683 0 : if (ierr /= 0) then
684 0 : write(*,2) 'T_1', N, T_1
685 0 : write(*,2) 'Tmax', N, Tmax
686 0 : write(*,2) 'resid_T_min', N, resid_T_min
687 0 : write(*,2) 'resid_T_max', N, resid_T_max
688 0 : return
689 : call mesa_error(__FILE__,__LINE__,'get_T failed in root find')
690 : end if
691 0 : T_0 = exp(lnT)
692 : residual = get_T_residual(lnT, &
693 0 : dfdx, lrpar, rpar, lipar, ipar, ierr)
694 0 : get_T_estimate = .true.
695 0 : end function get_T_estimate
696 :
697 0 : logical function get_T(ierr)
698 : integer, intent(out) :: ierr
699 0 : real(dp) :: Prad
700 0 : ierr = 0
701 0 : if (get_T_estimate()) then
702 : ! use T_0 from get_T_estimate as initial guess
703 0 : T4_0 = T_0*T_0*T_0*T_0
704 : else
705 0 : T4_0=T4_1*1.2d0
706 0 : T_0=sqrt(sqrt(T4_0))
707 0 : if (s% RSP_testing) then
708 0 : do
709 0 : Prad = Radiation_Pressure(T_0)
710 0 : if (Prad < P_0) exit
711 0 : write(*,*) '1 reduce T_0 in get_T', N
712 0 : T_0 = 0.99d0*T_0
713 : end do
714 : end if
715 : end if
716 :
717 0 : Lc_loop1: do ! reduce T if Lc >= L
718 : call EOP(s,N, &
719 0 : T_0,P_0,V_0,E_0,CP_0,QQ_0,SVEL_0,OP_0,ierr)
720 0 : if (ierr /= 0) return
721 0 : call CFLUX(HP_0,IGR_0,Lc_0,w_0,GPF,N)
722 0 : if(Lc_0<L) exit Lc_loop1
723 0 : T_0=T_0-10.d0*(Lc_0/L)
724 0 : if (T_0 <= 0) then
725 0 : ierr = -1
726 0 : return
727 : call mesa_error(__FILE__,__LINE__,'T_0 <= 0 in Lc_loop1')
728 : end if
729 0 : T4_0=T_0**4
730 : end do Lc_loop1
731 :
732 0 : TT=4.d0*SIG*P4**2*R_1**4/(3.d0*dm_bar_0*L)
733 : Lr_0=TT*(T4_0/OP_0-T4_1/OP_1)/ &
734 0 : (1.d0-dlog(OP_0/OP_1)/dlog(T4_0/T4_1))*L
735 0 : F1=(Lr_0+Lc_0)/L-1.d0
736 :
737 0 : D=T4_0/1.d3
738 0 : I=0
739 0 : T1_loop: do ! adjust T to make Lr + Lc = L
740 0 : I=I+1
741 0 : if(I>10000) then
742 0 : get_T = .false.
743 0 : if (s% RSP_testing) then
744 0 : write(*,*) 'failed get_T', N, T_0
745 0 : call mesa_error(__FILE__,__LINE__,'get_T')
746 : end if
747 0 : return
748 : end if
749 0 : do while (abs(D/T4_0)>0.5d0)
750 0 : D=(T4_0/2.d0)*(D/abs(D))
751 : end do
752 0 : T4_0 = T4_0-D
753 0 : T_0=sqrt(sqrt(T4_0))
754 0 : Lc_loop: do ! reduce T if Lc >= L
755 : call EOP(s,N, &
756 0 : T_0,P_0,V_0,E_0,CP_0,QQ_0,SVEL_0,OP_0,ierr)
757 0 : if (ierr /= 0) return
758 0 : call CFLUX(HP_0,IGR_0,Lc_0,w_0,GPF,N)
759 0 : if (Lc_0<L) exit Lc_loop
760 0 : T_0=T_0-10.d0*(Lc_0/L)
761 0 : if (T_0 <= 0) then
762 0 : ierr = -1
763 0 : return
764 : call mesa_error(__FILE__,__LINE__,'T_0 <= 0 in Lc_loop')
765 : end if
766 0 : T4_0=T_0**4
767 : end do Lc_loop
768 : Lr_0=TT*(T4_0/OP_0-T4_1/OP_1)/ &
769 0 : (1.d0-dlog(OP_0/OP_1)/dlog(T4_0/T4_1))*L
770 0 : F2=(Lr_0+Lc_0)/L-1.d0
771 0 : if(ABS(F2)<PREC .or. F2==F1) exit T1_loop
772 0 : D=F2*(T4_0-T4_0)/(F2-F1)
773 0 : F1=F2
774 0 : T4_0=T4_0
775 : end do T1_loop
776 :
777 0 : get_T = .true.
778 0 : if (s% RSP_testing) write(*,*) 'done get_T', N, T_0
779 0 : end function get_T
780 :
781 0 : subroutine failed
782 0 : write(*,*) 'NO CONVERGENCE IN STA INNER LOOP',I
783 0 : if((.not.adjusting_dmN) .OR. dmN_cnt>1 .or. IG > 54) then
784 0 : write(*,*) 'zone ',N,'IGR= ',IGR_0
785 0 : stop
786 : end if
787 0 : dmN=dmN/4.d0
788 0 : write(*,*) 'PHOENIX CONDITION'
789 0 : IG=IG+1
790 0 : end subroutine failed
791 :
792 0 : subroutine prepare_for_new_H
793 0 : ITIN = ITIN+1
794 : ! STRATOWE WARTOSCI DLA ITERACJI PONIZEJ TH0
795 0 : N1 = NZN-NZT+1
796 0 : R_1 = R(N1)
797 0 : V_0 = Vol(N1)
798 0 : P_0 = P(N1)
799 0 : T_0 = T(N1)
800 0 : T4_0 = T_0**4
801 0 : M_0 = M(N1)
802 0 : dm_0 = dm(N1)
803 0 : OP_0 = K(N1)
804 0 : N = N1
805 0 : H = HH
806 0 : HP_0 = Hp_face(N1)
807 0 : IGR_0 = Y_face(N1)
808 0 : QQ_0 = QQS(N1)
809 0 : CP_0 = CPS(N1)
810 0 : Lc_0 = Lc(N1)
811 0 : E_0 = E(N1)
812 0 : end subroutine prepare_for_new_H
813 :
814 0 : subroutine change_NZN
815 0 : NZN=NZN-N
816 0 : do I=1,NZN
817 0 : R(I)=R(I+N)
818 0 : P(I)=P(I+N)
819 0 : Vol(I)=Vol(I+N)
820 0 : E(I)=E(I+N)
821 0 : T(I)=T(I+N)
822 0 : K(I)=K(I+N)
823 0 : M(I)=M(I+N)
824 0 : dm(I)=dm(I+N)
825 0 : dm_bar(I)=dm_bar(I+N)
826 0 : Hp_face(I) = Hp_face(I+N)
827 0 : Y_face(I) = Y_face(I+N)
828 0 : CPS(I) = CPS(I+N)
829 0 : QQS(I) = QQS(I+N)
830 0 : Lc(I) = Lc(I+N)
831 0 : Lr(I) = Lr(I+N)
832 0 : w(I) = w(I+N)
833 : end do
834 0 : end subroutine change_NZN
835 :
836 : end subroutine STAH
837 :
838 :
839 0 : subroutine ZNVAR(s,H,dmN,L,TE,M,ierr)
840 : type (star_info), pointer :: s
841 : integer, intent(out) :: ierr
842 0 : real(dp) :: H,dmN,L,TE,M,Psurf
843 0 : real(dp) :: T0,R,TAU0,P,V, &
844 0 : dtau, kap, alfa, G_M_dtau_div_R2, Prad, Pgas_0, &
845 0 : dP_dV, dkap_dV, xx, residual, d_residual_dlnV, &
846 0 : dkap_dlnV, dlnV, dP_dlnV, lnV
847 : integer :: I
848 0 : H=HHFAC
849 0 : ierr = 0
850 : if(.not.RSP_eddi) then ! EXACT GREY RELATION
851 : T0= pow(sqrt(3.d0)/4.d0,0.25d0)*TE !0.811194802d0*TE
852 : else ! EDDINGTON APPROXIMATION
853 0 : T0= pow(0.5d0, 0.25d0)*TE ! T0= pow(1.0d0/2.d0,0.25d0)*TE
854 : end if
855 : if (s% RSP_use_Prad_for_Psurf) then
856 : Psurf = crad*T0*T0*T0*T0/3d0
857 : else
858 0 : Psurf = 0d0
859 : end if
860 0 : R=sqrt(L/(4.d0*PI*SIG))/TE**2
861 0 : TAU0=0.003d0
862 :
863 : ! dtau is optical depth we want to middle of 1st cell; ~ 0.003
864 : ! dtau = dm*kap/(4*pi*R^2)
865 : ! mass of 1st cell will be 2*dm, so dm is mass above center of cell.
866 : ! R is derived from L and Teff.
867 : ! T0 is derived from Teff. it is temperature at center of 1st cell.
868 : ! Psurf is pressure just outside 1st cell; either 0 or, better, Prad(T0).
869 : ! need to find P, pressure at center of 1st cell
870 : ! P = Psurf + G*M*dm/(4*pi*R^4)
871 : ! substitute for dm to get
872 : ! P = Psurf + G*M*dtau/(R^2*kap)
873 : ! kap depends on P, so need to solve this implicit equation iteratively.
874 : ! once find P and kap, can evaluate dm = 4*pi*R^2*dtau/kap
875 :
876 : ! note that since T0 is fixed, we can only change P by changing Pgas.
877 : ! for most cases this is not a problem since Prad << Pgas.
878 : ! but for massive blue stars, we can have Prad >> Pgas.
879 : ! to handle both cases well, we need to iterate on Pgas instead of P
880 : ! that will let us make sure we don't create guesses that imply Pgas < 0.
881 :
882 0 : dtau = TAU0 ! s% RSP_outer_dtau_target
883 : ! make rough initial guess for opacity based on T0
884 0 : if (T0 < 4700d0) then
885 0 : kap = 1d-3
886 0 : else if (T0 > 5100d0) then
887 0 : kap = 1d-2*(1d0 + X)
888 : else
889 0 : alfa = (T0 - 4700)/(5100 - 4700)
890 0 : kap = alfa*1d-2*(1d0 + X) + (1d0-alfa)*1d-3
891 : end if
892 0 : G_M_dtau_div_R2 = G*M*dtau/R**2
893 0 : Prad = crad*T0*T0*T0*T0/3d0
894 0 : Pgas_0 = G_M_dtau_div_R2/kap
895 0 : P = Pgas_0 + Prad ! initial guess for P
896 0 : call EOP(s,0,T0,P,V,xx,xx,xx,xx,xx,ierr) ! initial V
897 0 : if (ierr /= 0) return
898 : !write(*,*) 'init T0,P,V', T0,P,V
899 0 : do I=1,25
900 : call mesa_eos_kap(s,-4, &
901 0 : T0,V,P,dP_dV,xx,xx,xx,xx,xx,xx,xx,xx,xx,xx,kap,dkap_dV,xx,ierr)
902 0 : if (ierr /= 0) return
903 0 : residual = P - (Prad + G_M_dtau_div_R2/kap)
904 0 : if (abs(residual) < 1d-6*P) exit ! done
905 0 : dP_dlnV = dP_dV*V
906 0 : lnV = log(V)
907 0 : dkap_dlnV = dkap_dV*V
908 0 : d_residual_dlnV = dP_dlnV + G_M_dtau_div_R2*dkap_dlnV/kap**2
909 0 : dlnV = - residual/d_residual_dlnV
910 0 : V = exp(lnV + dlnV)
911 : !write(*,*) 'T0, new V, P, kap, residual, dP_dV', i, T0, V, P, kap, residual, dP_dV
912 : end do
913 :
914 : !write(*,*) 'V, P, kap, residual', i, V, P, kap, residual
915 :
916 0 : dmN = 4*pi*R**2*dtau/kap
917 0 : if (s% RSP_testing) write(*,*) 'initial dmN', dmN/SUNM
918 : !stop
919 : end subroutine ZNVAR
920 :
921 :
922 0 : subroutine CFLUX(HP_0,IGR_0,Lc_0,OMEGA_0,GPF,N)
923 :
924 0 : real(dp) :: POM,POM2,HP_0,IGR_0,OMEGA_0,PII,Lc_0
925 0 : real(dp) :: FF,GG,GPF,ENT
926 0 : real(dp) :: AA,BB,CC,DELTA
927 : integer :: N
928 : !-
929 0 : if(ALFA==0.d0) then
930 0 : OMEGA_0=0.d0
931 0 : Lc_0=0.d0
932 0 : return
933 : end if
934 :
935 : ! PRESSURE SCALE HEIGHT
936 0 : HP_0=R_1**2/(G*M_0)*(P_0*V_0+P_1*V_1)/2.d0
937 0 : POM=P4*R_1**2*HP_0/dm_bar_0/(V_0+V_1)*2.d0
938 :
939 : ! SUPERADIABATIC GRADIENT, Y
940 : IGR_0=POM*((QQ_0/CP_0+QQ_1/CP_1)/2.d0*(P_1-P_0) &
941 0 : -(dlog(T_1)-dlog(T_0))) !hyt!
942 :
943 0 : POM=sqrt(2.d0/3.d0)*0.5d0
944 0 : ENT=(E_0+P_0*V_0)/T_0+(E_1+P_1*V_1)/T_1
945 0 : FF=POM*ENT
946 0 : POM=ALFAS*ALFA
947 0 : POM2=0.5d0*(CP_0+CP_1)
948 0 : GG=POM*POM2*IGR_0
949 0 : GPF=GG/FF
950 :
951 :
952 : ! BOTTOM BOUNDARY CONDITION FOR CONVECTION
953 0 : if(N<=IBOTOM)then
954 0 : Lc_0=0.d0
955 0 : OMEGA_0=0.d0
956 0 : return
957 : end if
958 0 : if(IGR_0>0.d0)then
959 : if(.true. .or. GAMMAR==0.d0)then ! gammar breaks Cep 11.5M model BP
960 : OMEGA_0=sqrt(ALFA/CEDE*FF*GPF* &
961 0 : (T_0*P_0*QQ_0/CP_0+T_1*P_1*QQ_1/CP_1)*0.5d0)
962 : else
963 : AA=CEDE/ALFA
964 : POM=(GAMMAR**2/ALFA**2)*4.d0*SIG
965 : BB=(POM/HP_0)* &
966 : (T_0**3*V_0**2/CP_0/OP_0 &
967 : +T_1**3*V_1**2/CP_1/OP_1)*0.5d0
968 : CC=-(T_0*P_0*QQ_0/CP_0 &
969 : +T_1*P_1*QQ_1/CP_1)*0.5d0*FF*GPF
970 : DELTA=BB**2-4.d0*AA*CC
971 : if(DELTA<=0.d0) then
972 : write(*,*) 'CFLUX: Error! : Y>0, but no solution found'
973 : stop
974 : end if
975 : OMEGA_0=(-BB+sqrt(DELTA))/(2.d0*AA)
976 : end if
977 0 : PII=FF*OMEGA_0*GPF
978 0 : Lc_0=P4*R_1**2*(T_0/V_0+T_1/V_1)*0.5d0*PII*(ALFAC/ALFAS)
979 : ! (REPLACES ALFAS BY ALFAC)
980 : else
981 0 : OMEGA_0=0.d0
982 0 : Lc_0=0.d0
983 : end if
984 : end subroutine CFLUX
985 :
986 : end module rsp_build
|
