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_relax_env
21 : use star_def, only: star_info
22 : use utils_lib, only: is_bad, mesa_error
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
28 :
29 : implicit none
30 :
31 : private
32 : public :: EOP, RELAX_ENV
33 :
34 : contains
35 :
36 0 : subroutine EOP(s,k,T,P,V,E,CP,QQ,SVEL,OP,ierr)
37 : use rsp_eval_eos_and_kap, only: &
38 : eval1_mesa_Rho_given_PT, eval1_gamma_PT_getRho
39 : type (star_info), pointer :: s
40 : integer, intent(out) :: ierr
41 0 : real(dp) :: T,V,P1,DPV,E,OP,F2,P,SVEL,CP,QQ
42 0 : real(dp) :: DU1,DU2,DU3,DU4,DU5,DU6,DU7,DU8,DU9
43 : integer :: k,i
44 :
45 0 : real(dp) :: R,A,PRECEQ,rho_guess,rho,Prad,P_test
46 : data R,A/8.317d7,7.5641d-15/
47 : data PRECEQ/5.d-13/
48 :
49 0 : if(P<=0.d0) GOTO 100 !negative pressure, stop
50 :
51 0 : ierr = 0
52 0 : Prad = Radiation_Pressure(T)
53 0 : if (P <= Prad) then
54 0 : write(*,*) 'P <= Prad', k, P, Prad, T
55 0 : ierr = -1
56 0 : return
57 : call mesa_error(__FILE__,__LINE__,'EOP')
58 : end if
59 0 : rho_guess = eval1_gamma_PT_getRho(s, 0, P, T, ierr)
60 0 : if (ierr /= 0) then
61 0 : write(*,*) 'eval1_gamma_PT_getRho failed', k
62 0 : write(*,*) 'P', P
63 0 : write(*,*) 'T', T
64 0 : ierr = -1
65 0 : return
66 : call mesa_error(__FILE__,__LINE__,'EOP')
67 : end if
68 0 : call eval1_mesa_Rho_given_PT(s, 0, P, T, rho_guess, rho, ierr)
69 0 : if (ierr /= 0) then
70 0 : write(*,*) 'eval1_mesa_Rho_given_PT failed', k
71 0 : write(*,*) 'P', P
72 0 : write(*,*) 'T', T
73 0 : write(*,*) 'rho_guess', rho_guess
74 0 : ierr = -1
75 0 : return
76 : call mesa_error(__FILE__,__LINE__,'EOP')
77 : end if
78 0 : V = 1d0/rho ! initial guess to be improved below
79 :
80 : !write(*,*) 'eval1_mesa_Rho_given_PT', k
81 : !write(*,*) 'P', P
82 : !write(*,*) 'T', T
83 : !write(*,*) 'rho_guess', rho_guess
84 : !write(*,*) 'V', 0, V
85 :
86 0 : I=0
87 : ! NEWTON-RAPHSON ITERATION (TO MAKE P1->P)
88 0 : 1 I=I+1
89 : call mesa_eos_kap(s,0, &
90 : T,V,P1,DPV,DU1,E,DU2,DU3,CP,DU4,DU5,QQ,DU6,DU7,OP &
91 0 : ,DU8,DU9,ierr)
92 0 : if (ierr /= 0) return
93 0 : if(I>25) GOTO 2 !no convergence, stop
94 0 : F2=P1/P-1.d0
95 0 : V=V-F2*P/DPV
96 0 : if (is_bad(V)) then
97 0 : write(*,*) 'V', V
98 0 : write(*,*) 'F2', F2
99 0 : write(*,*) 'DPV', DPV
100 0 : write(*,*) 'P', P
101 0 : write(*,*) 'EOP I', I
102 0 : ierr = -1
103 0 : return
104 : stop
105 : end if
106 :
107 : !write(*,*) 'V', I, V
108 :
109 0 : if(abs(F2)<PRECEQ) GOTO 3
110 0 : GOTO 1
111 : 3 continue
112 : if (abs(P - P_test) > 1d-10*P) then
113 : end if
114 0 : return
115 :
116 0 : 2 write(*,*) 'NO CONVERGENCE IN EOP; T,P= ',T,P
117 0 : ierr = -1
118 0 : return
119 0 : stop
120 0 : 100 write(*,*) 'NEGATIVE OR ZERO PRESSURE IN EOP'
121 0 : ierr = -1
122 0 : return
123 : stop
124 0 : end subroutine EOP
125 :
126 :
127 0 : subroutine RELAX_ENV(s, L0, TH0, TE, NZT, NZN, &
128 0 : M, DM, DM_BAR, R, Vol, T, Et, ierr)
129 0 : use star_utils, only: rand
130 : type (star_info), pointer :: s
131 : integer, intent(in) :: NZN, NZT
132 : real(dp), intent(in) :: L0, TH0, TE
133 : real(dp), intent(inout), dimension(:) :: &
134 : M, DM, DM_BAR, R, Vol, T, Et
135 : integer, intent(out) :: ierr
136 :
137 : logical, parameter :: RSP_eddi = .true. ! use Eddington approx at surface
138 :
139 : real(dp), allocatable, dimension(:) :: &
140 0 : E, P, Lr, Lc, Hp_face, Y_face, K, CPS, QQS, &
141 0 : dC_dr_in,dC_dr_00,dC_dr_out,dC_dT_00,dC_dT_out, &
142 0 : dC_dw_00,dC_dr_in2,dC_dT_in, &
143 0 : DPV,dP_dT_00,DEV,dE_dT_00,dK_dV_00, &
144 0 : dK_dT_00,DVR,DVRM, &
145 0 : CPV, dCp_dT_00, QQV,dQQ_dT_00, &
146 0 : dP_dr_00, dCp_dr_00, dQQ_dr_00, &
147 0 : dP_dr_in,dCp_dr_in,dQQ_dr_in, &
148 0 : dHp_dr_in,dHp_dr_out,dHp_dr_00,dHp_dT_00,dHp_dT_out, &
149 0 : dY_dr_in,dY_dr_00,dY_dr_out,dY_dT_00, &
150 0 : dY_dT_out, &
151 0 : PII,dPII_dr_in,dPII_dr_00,dPII_dr_out, &
152 0 : dPII_dT_00,dPII_dT_out,DPIIZ0, &
153 0 : dsrc_dr_in,dsrc_dr_00,dsrc_dr_out,dsrc_dT_00, &
154 0 : dsrc_dT_out,dsrc_dw_00,SOURCE, &
155 0 : d_damp_dr_in,d_damp_dr_00,d_damp_dr_out,d_damp_dT_00, &
156 0 : d_damp_dT_out,d_damp_dw_00,DAMP, &
157 0 : dsrc_dr_in2,dsrc_dT_in,d_damp_dr_in2,d_damp_dT_in, &
158 0 : d_dampR_dr_in,d_dampR_dr_00,d_dampR_dr_out,d_dampR_dT_00, &
159 0 : d_dampR_dT_out,d_dampR_dw_00,DAMPR, &
160 0 : d_dampR_dr_in2,d_dampR_dT_in, &
161 0 : DLCXM,DLCX0,DLCXP,DLCY0, &
162 0 : DLCYP,DLCZ0,DLCZP, &
163 0 : DLTXM,DLTX0,DLTXP,DLTY0, &
164 0 : DLTYP,DLTZ0,DLTZP,Lt, &
165 0 : PTURB,dPtrb_dw_00,dPtrb_dr_00,dPtrb_dr_in
166 :
167 0 : real(dp) :: T1,DLR,DLRP,DLRM,DLT,DLTP,DLMR,DLMRP,DLMRM,DLMT,DLMTP
168 0 : real(dp) :: W_00,W_out,BW,BK,T2,T3,DLK,DLKP,SVEL,T_0
169 0 : real(dp) :: POM3,POM2,POM
170 : integer :: I,J,IW,IR,IC, IMAXR, IMAXW, IMAXC
171 0 : real(dp) :: MAXR,MAXW,MAXC,PB,PPB,ELB,ELMB
172 :
173 : integer :: INFO,II,ITROUBC,ITROUBT,IZIP
174 0 : real(dp) :: XXR,XXC,XXT,EZH,DXH,DXXC,DXXT,DXKC,DXKT, &
175 0 : IGR1,IGR1XM,IGR1X0,IGR1XP,IGR1Y0,IGR1YP, &
176 0 : IGR2,IGR2XM,IGR2X0,IGR2XP,IGR2Y0,IGR2YP
177 :
178 0 : real(dp) :: FFXM,FFX0,FFXP,FFY0,FFYP,FF
179 0 : real(dp) :: GGXM,GGX0,GGXP,GGY0,GGYP,GG,GPF
180 0 : real(dp) :: EFL02
181 :
182 0 : real(dp) :: POM4,DXXR,TEM1,TEMI,TEMM
183 :
184 0 : real(dp) :: TT,dmN,dmNL,TNL,DDT,MSTAR,PRECR
185 0 : real(dp) :: RINNER,ROUTER,TINNER,MENVEL
186 : integer :: IOP,NEGFLU
187 : logical :: ending
188 :
189 0 : real(dp) :: SUMM,AONE,HAHA
190 0 : real(dp) :: AAA(100),AAP(100),AAT(100)
191 0 : real(dp) :: AALFA, AALFAP, AALFAT
192 : integer :: ICAA, ICAP, ICAT
193 : integer :: NDIVAA, NDIVAP, NDIVAT, dmN_cnt, max_dmN_cnt
194 : logical :: ok_to_adjust_mass, ok_to_adjust_Tsurf
195 0 : real(dp) :: EDFAC, Psurf, CFIDDLE, FSUB, EMR, ELR
196 0 : real(dp) :: PREC1
197 :
198 0 : real(dp) :: XX_max, XX_max_val, XX_max_dx
199 : integer :: i_XX_max, var_XX_max, n, op_err
200 :
201 : !write(*,*) 'RELAX_ENV'
202 0 : EFL02 = EFL0*EFL0
203 0 : FSUB = s% RSP_dq_1_factor
204 0 : EMR = s% RSP_mass
205 0 : ELR = s% RSP_L
206 :
207 0 : n = NZN+1
208 : allocate( &
209 : E(n), P(n), Lr(n), Lc(n), Hp_face(n), Y_face(n), K(n), CPS(n), QQS(n), &
210 : dC_dr_in(n),dC_dr_00(n),dC_dr_out(n),dC_dT_00(n),dC_dT_out(n), &
211 : dC_dw_00(n),dC_dr_in2(n),dC_dT_in(n), &
212 : DPV(n),dP_dT_00(n),DEV(n),dE_dT_00(n),dK_dV_00(n), &
213 : dK_dT_00(n),DVR(n),DVRM(n), &
214 : CPV(n), dCp_dT_00(n), QQV(n),dQQ_dT_00(n), &
215 : dP_dr_00(n), dCp_dr_00(n), dQQ_dr_00(n), &
216 : dP_dr_in(n),dCp_dr_in(n),dQQ_dr_in(n), &
217 : dHp_dr_in(n),dHp_dr_out(n),dHp_dr_00(n),dHp_dT_00(n),dHp_dT_out(n), &
218 : dY_dr_in(n),dY_dr_00(n),dY_dr_out(n),dY_dT_00(n), &
219 : dY_dT_out(n), &
220 : PII(n),dPII_dr_in(n),dPII_dr_00(n),dPII_dr_out(n), &
221 : dPII_dT_00(n),dPII_dT_out(n),DPIIZ0(n), &
222 : dsrc_dr_in(n),dsrc_dr_00(n),dsrc_dr_out(n),dsrc_dT_00(n), &
223 : dsrc_dT_out(n),dsrc_dw_00(n),SOURCE(n), &
224 : d_damp_dr_in(n),d_damp_dr_00(n),d_damp_dr_out(n),d_damp_dT_00(n), &
225 : d_damp_dT_out(n),d_damp_dw_00(n),DAMP(n), &
226 : dsrc_dr_in2(n),dsrc_dT_in(n),d_damp_dr_in2(n),d_damp_dT_in(n), &
227 : d_dampR_dr_in(n),d_dampR_dr_00(n),d_dampR_dr_out(n),d_dampR_dT_00(n), &
228 : d_dampR_dT_out(n),d_dampR_dw_00(n),DAMPR(n), &
229 : d_dampR_dr_in2(n),d_dampR_dT_in(n), &
230 : DLCXM(n),DLCX0(n),DLCXP(n),DLCY0(n), &
231 : DLCYP(n),DLCZ0(n),DLCZP(n), &
232 : DLTXM(n),DLTX0(n),DLTXP(n),DLTY0(n), &
233 : DLTYP(n),DLTZ0(n),DLTZP(n),Lt(n), &
234 0 : PTURB(n),dPtrb_dw_00(n),dPtrb_dr_00(n),dPtrb_dr_in(n))
235 :
236 : ! STORE SOME VALUES FOR FURTHER COMPARISON
237 0 : MSTAR = M(NZN)
238 0 : TINNER = T(1)
239 0 : RINNER = s% R_center ! R0
240 0 : ROUTER = R(NZN)
241 0 : MENVEL = MSTAR-M(1)+dm(1)
242 :
243 0 : TNL = 0d0
244 :
245 0 : if (s% RSP_use_Prad_for_Psurf) then
246 : if(.not.RSP_eddi) then ! EXACT GREY RELATION
247 : T_0= pow(sqrt(3.d0)/4.d0,0.25d0)*TE !0.811194802d0*TE
248 : else ! EDDINGTON APPROXIMATION
249 0 : T_0= pow(0.5d0, 0.25d0)*TE ! T_0= pow(1.0d0/2.d0,0.25d0)*TE
250 : end if
251 0 : Psurf = crad*T_0*T_0*T_0*T_0/3d0
252 : else
253 : Psurf = 0d0
254 : end if
255 :
256 0 : ending=.false.
257 0 : if (s% RSP_trace_RSP_build_model) &
258 0 : write(*,*) '*** relax envelope ***'
259 :
260 0 : NEGFLU=0
261 0 : NDIVAA = 20 ! s% RSP_NDIVAA
262 0 : NDIVAP = 20 ! s% RSP_NDIVAP
263 0 : NDIVAT = 20 ! s% RSP_NDIVAT
264 0 : ok_to_adjust_mass = .true. ! s% RSP_ok_to_adjust_mass
265 0 : ok_to_adjust_Tsurf = .true. ! s% RSP_ok_to_adjust_Tsurf
266 0 : dmN_cnt = 0
267 0 : max_dmN_cnt = s% RSP_relax_max_tries
268 0 : PREC1 = s% RSP_relax_dm_tolerance
269 :
270 : ! SUMM IS ZONE OF THE ENVELOPE UP TO ANCHOR
271 : ! (NOT CHANGED IN THE ITERATIONS)
272 0 : SUMM=0.d0
273 0 : do I=1,NZN-NZT+1
274 0 : SUMM=SUMM+dm(I)
275 : end do
276 :
277 : ! IF BOTH ALFAP NAD ALFAT /= 0, IT IS NECESSARY TO ITERATE WITHOUT
278 : ! TURBULENT FLUX, WITH TURBULENT PRESSURE ONLY, AND AFTER
279 : ! CONVERGENCE TURBULENT FLUX IS ITERATED
280 0 : AALFA = -1 ! turn off ALFA relax
281 0 : if (AALFA <= 0d0) AALFA = ALFA
282 0 : AALFAT = ALFAT
283 0 : AALFAP = ALFAP
284 0 : ALFAT = 0.d0
285 0 : ALFAP = 0.d0
286 :
287 : ! SET ALFA TO ITERATE
288 0 : do I=1,NDIVAA
289 0 : AAA(I)=ALFA+(AALFA-ALFA)*I/dble(NDIVAA)
290 : end do
291 0 : ICAA=1
292 :
293 : ! SET ITERATIONS FOR ALFAP
294 0 : if(AALFAP/=0.d0)then
295 0 : do I=1,NDIVAP
296 0 : AAP(I)=AALFAP*I/dble(NDIVAP)
297 : end do
298 : ICAP=1
299 : end if
300 :
301 : ! SET ITERATIONS FOR ALFAT
302 0 : if(AALFAT/=0.d0)then
303 0 : do I=1,NDIVAT
304 0 : AAT(I)=AALFAT*I/dble(NDIVAT)
305 : end do
306 : ICAT=1
307 : end if
308 : !-
309 0 : PRECR = 1.d-10 !PRECISION FOR NEWTON-RHAPSON ITERATIONS
310 0 : DXH = 0.01d0 !UNDERCORRECTION FOR NEWTON-RHAPSON CORRECTIONS
311 0 : DDT = -dm(NZN)/1000.d0 !INITIAL CHANGE IN OUTER ZONE MASS
312 :
313 0 : IOP = 0
314 :
315 : 999 continue
316 :
317 0 : IZIP = 0
318 0 : if(IOP==0) GOTO 100
319 :
320 : !------ MASS TRICKS ---
321 : ! dmN = dm(NZN) ! THIS IS NOT WORKING IF FSUB/=1
322 0 : dmN = dm(NZN-1) !THIS IS WORKING FOR ALL FSUB VALUES
323 :
324 0 : TT = T(NZN-NZT+1)-TH0
325 :
326 0 : if(IOP==1) GOTO 24
327 :
328 0 : DDT = TT*(dmN-dmNL)/(TT-TNL)
329 : !if (is_bad(DDT)) then
330 : ! write(*,*) 'DDT', DDT
331 : ! write(*,*) 'TT-TNL', TT-TNL
332 : ! write(*,*) 'dmN-dmNL', dmN-dmNL
333 : ! write(*,*) 'TT', TT
334 : ! write(*,*) 'TNL', TNL
335 : ! write(*,*) 'dmNL', dmNL
336 : ! write(*,*) 'dmN', dmN
337 : ! call mesa_error(__FILE__,__LINE__,'failed in RELAX_ENV')
338 : !end if
339 0 : CFIDDLE=0.1d0
340 0 : if(abs(DDT/dmN)>CFIDDLE) DDT=CFIDDLE*dmN*(DDT/abs(DDT))
341 :
342 : ! CHECK IF ALFA ITERATION IS FINISHED
343 0 : if(abs(DDT/dmN)<PREC1.and.ALFA/=AALFA) then
344 0 : ALFA=AAA(ICAA)
345 0 : ICAA=ICAA+1
346 : ! APPLY ARTIFICIAL ZONE MASS CHANGE TO ALLOW ITERATIONS
347 0 : DDT = -dm(NZN)/100000.d0
348 0 : write(*,*) 'ALFA =',ALFA
349 0 : GOTO 24
350 : end if
351 0 : if(abs(DDT/dmN)<PREC1.and.ICAA==NDIVAA+1)then
352 0 : write(*,*) '***** ALFA ITERATION FINISHED *****'
353 0 : write(*,*) 'final ALFA =', ALFA, s% RSP_alfa
354 0 : dmN_cnt = 0
355 0 : NDIVAA=-99
356 : end if
357 :
358 0 : if (s% RSP_relax_alfap_before_alfat) then
359 :
360 : ! CHECK IF ALFAP ITERATION IS FINISHED
361 0 : if(abs(DDT/dmN)<PREC1.and.ALFAP/=AALFAP) then
362 0 : ALFAP=AAP(ICAP)
363 0 : ICAP=ICAP+1
364 : ! APPLY ARTIFICIAL ZONE MASS CHANGE TO ALLOW ITERATIONS
365 0 : DDT = -dm(NZN)/100000.d0
366 0 : write(*,*) 'ALFAP=',ALFAP
367 0 : GOTO 24
368 : end if
369 0 : if(abs(DDT/dmN)<PREC1.and.ICAP==NDIVAP+1)then
370 0 : write(*,*) '***** ALFAP ITERATION FINISHED *****'
371 0 : dmN_cnt = 0
372 0 : NDIVAP=-99
373 : end if
374 :
375 : ! CHECK IF ALFAT ITERATION IS FINISHED
376 0 : if(abs(DDT/dmN)<PREC1.and.ALFAT/=AALFAT) then
377 0 : ALFAT=AAT(ICAT)
378 0 : ICAT=ICAT+1
379 : ! APPLY ARTIFICIAL ZONE MASS CHANGE TO ALLOW ITERATIONS
380 0 : DDT = -dm(NZN)/100000.d0
381 0 : write(*,*) 'ALFAT =',ALFAT
382 0 : GOTO 24
383 : end if
384 0 : if(abs(DDT/dmN)<PREC1.and.ICAT==NDIVAT+1)then
385 0 : write(*,*) '***** ALFAT ITERATION FINISHED *****'
386 0 : dmN_cnt = 0
387 0 : NDIVAT=-99
388 : end if
389 :
390 : else
391 :
392 : ! CHECK IF ALFAT ITERATION IS FINISHED
393 0 : if(abs(DDT/dmN)<PREC1.and.ALFAT/=AALFAT) then
394 0 : ALFAT=AAT(ICAT)
395 0 : ICAT=ICAT+1
396 : ! APPLY ARTIFICIAL ZONE MASS CHANGE TO ALLOW ITERATIONS
397 0 : DDT = -dm(NZN)/100000.d0
398 0 : write(*,*) 'ALFAT =',ALFAT
399 0 : GOTO 24
400 : end if
401 0 : if(abs(DDT/dmN)<PREC1.and.ICAT==NDIVAT+1)then
402 0 : write(*,*) '***** ALFAT ITERATION FINISHED *****'
403 0 : dmN_cnt = 0
404 0 : NDIVAT=-99
405 : end if
406 :
407 : ! CHECK IF ALFAP ITERATION IS FINISHED
408 0 : if(abs(DDT/dmN)<PREC1.and.ALFAP/=AALFAP) then
409 0 : ALFAP=AAP(ICAP)
410 0 : ICAP=ICAP+1
411 : ! APPLY ARTIFICIAL ZONE MASS CHANGE TO ALLOW ITERATIONS
412 0 : DDT = -dm(NZN)/100000.d0
413 0 : write(*,*) 'ALFAP=',ALFAP
414 0 : GOTO 24
415 : end if
416 0 : if(abs(DDT/dmN)<PREC1.and.ICAP==NDIVAP+1)then
417 0 : write(*,*) '***** ALFAP ITERATION FINISHED *****'
418 0 : dmN_cnt = 0
419 0 : NDIVAP=-99
420 : end if
421 :
422 : end if
423 :
424 0 : if(abs(DDT/dmN)<PREC1) then
425 0 : if(NEGFLU>=1)then
426 0 : if (NEGFLU==2)then
427 : !write(*,*) '*** ENVELOPE IS RELAXED ***'
428 : ending=.true.
429 : GOTO 100
430 : end if
431 : !write(*,*) 'NEGATIVE FLUX IS ON'
432 0 : DDT = -dm(NZN)/100000.d0
433 0 : NEGFLU=2
434 0 : GOTO 24
435 : end if
436 : !write(*,*) 'NEGATIVE FLUX WITHOUT Z-DERIV. IS ON'
437 0 : DDT = -dm(NZN)/100000.d0
438 0 : NEGFLU=1
439 0 : GOTO 24
440 : end if
441 :
442 : !write(*,*) 'abs(DDT/dmN), PREC1', DDT, dmN, abs(DDT/dmN), PREC1
443 :
444 : if(abs(DDT/dmN)<PREC1) then
445 : !write(*,*) '*** ENVELOPE IS RELAXED ***'
446 : ending=.true.
447 : GOTO 100
448 : end if
449 :
450 0 : if (dmN_cnt >= max_dmN_cnt) then
451 0 : write(*,*) 'RELAX_ENV has reached max num allowed tries for outer dm', max_dmN_cnt
452 0 : ierr = -1
453 0 : return
454 : call mesa_error(__FILE__,__LINE__)
455 : end if
456 :
457 0 : 24 TNL = TT
458 0 : dmNL = dmN
459 0 : dmN_cnt = dmN_cnt + 1
460 :
461 0 : if (s% RSP_relax_adjust_inner_mass_distribution) then
462 :
463 0 : dmN = dmN-DDT
464 :
465 : if(IOP>=1) then
466 : ! CALCULATE HAHA - ZONE MASS RATIO BELOW ANCHOR
467 : ! TOTAL MASS BELOW ANCHOR: SUMM, FIRST ZONE MASS: AONE
468 0 : AONE=dmN
469 0 : HAHA=1.01d0
470 0 : do I=1,100
471 0 : POM=1.d0/(NZN-NZT+1)*dlog10(1.d0-SUMM/AONE*(1.d0-HAHA))
472 0 : if(dabs(HAHA-10.d0**POM)<1d-10) GOTO 22
473 0 : HAHA=10.d0**POM
474 : end do
475 0 : write(*,*) 'NO CONVERGENCE IN RELAX_ENV ITERATION FOR H'
476 0 : stop
477 : 22 continue
478 0 : HAHA=10.d0**POM
479 : ! SET NEW MASS DISTRIBUTION
480 0 : do I=NZN,1,-1
481 0 : if(I==NZN) then
482 : ! SPECIAL DEFINITIONS FOR THE OUTER ZONE
483 0 : M(I)=MSTAR
484 0 : dm(I)=dmN*FSUB
485 0 : dm_bar(I)=dm(I)*0.5d0
486 : else
487 0 : if(I>=NZN-NZT+1) dm(I)=dmN !DEFINITION DOWN TO ANCHOR
488 0 : if(I<NZN-NZT+1.and.ok_to_adjust_mass) &
489 0 : dm(I)=AONE*pow(HAHA,(NZN-NZT+1)-I)
490 0 : dm_bar(I)=(dm(I)+dm(I+1))*0.5d0
491 0 : M(I)=M(I+1)-dm(I+1)
492 : end if
493 : end do
494 : end if
495 :
496 : end if
497 :
498 : 100 continue
499 :
500 0 : do II=1,8000
501 :
502 : ! LOOP 1 .. EOS
503 0 : !$OMP PARALLEL DO PRIVATE(I,T1,POM,EDFAC,op_err) SCHEDULE(dynamic,2)
504 : do I=2,NZN
505 : T1=P43/dm(I)
506 : Vol(I)=T1*(R(I)**3-R(I-1)**3)
507 : DVRM(I) = -3.d0*T1*R(I-1)**2
508 : DVR(I) = 3.d0*T1*R(I)**2
509 :
510 : if(I==NZN.and.ok_to_adjust_Tsurf)then
511 : if(RSP_eddi)then
512 : EDFAC=1.d0/2.d0
513 : else
514 : EDFAC=sqrt(3.d0)/4.d0
515 : end if
516 : POM=(EDFAC*L0/(4.d0*PI*SIG*R(NZN)**2))**0.25d0
517 : T(NZN)=POM
518 : end if
519 :
520 : call mesa_eos_kap(s,0, &
521 : T(I),Vol(I),P(I),DPV(I),dP_dT_00(I),E(I), &
522 : DEV(I),dE_dT_00(I),CPS(I),CPV(I),dCp_dT_00(I), &
523 : QQS(I),QQV(I),dQQ_dT_00(I),K(I),dK_dV_00(I),dK_dT_00(I),op_err)
524 : if (op_err /= 0) ierr = op_err
525 : if (ierr /= 0) cycle
526 : if (is_bad(dQQ_dT_00(I))) then
527 : write(*,*) 'dQQ_dT_00(I)', i, dQQ_dT_00(I)
528 : write(*,*) 'QQS(I)', i, QQS(I)
529 : write(*,*) 'T(I)', i, T(I)
530 : write(*,*) 'Vol(I)', i, Vol(I)
531 : stop
532 : end if
533 :
534 : dP_dr_00(I) = DPV(I)*DVR(I)
535 : dP_dr_in(I) = DPV(I)*DVRM(I)
536 : dCp_dr_00(I) = CPV(I)*DVR(I)
537 : dCp_dr_in(I) = CPV(I)*DVRM(I)
538 : dQQ_dr_00(I) = QQV(I)*DVR(I)
539 : dQQ_dr_in(I) = QQV(I)*DVRM(I)
540 : end do
541 : !$OMP END PARALLEL DO
542 0 : if (ierr /= 0) return
543 :
544 : ! FIRST ZONE, CALCULATION OF R0 AND EOS
545 0 : P(1)=P(2)+G*M(1)*dm_bar(1)/(P4*R(1)**4)
546 : ! WITH GIVEN T AND P ITERATE FOR V AND CALCULATE ITS DERIVATIVES
547 0 : call EOP(s,0,T(1),P(1),Vol(1),E(1),CPS(1),QQS(1),SVEL,K(1),ierr)
548 0 : if (ierr /= 0) return
549 0 : s% R_center=pow(R(1)**3-3.d0*Vol(1)*dm(1)/P4,1.d0/3.d0)
550 0 : T1=P43/dm(1)
551 0 : DVRM(1) = -3.d0*T1*s% R_center**2
552 0 : DVR(1) = 3.d0*T1*R(1)**2
553 : ! CALCULATE EOS FOR THE FIRST ZONE TO OBTAIN DERIVATIVES
554 : call mesa_eos_kap(s,0, &
555 : T(1),Vol(1),P(1),DPV(1),dP_dT_00(1),E(1), &
556 : DEV(1),dE_dT_00(1),CPS(1),CPV(1),dCp_dT_00(1), &
557 0 : QQS(1),QQV(1),dQQ_dT_00(1),K(1),dK_dV_00(1),dK_dT_00(1),ierr)
558 0 : if (ierr /= 0) return
559 : !write(*,*) '1st zone T V P',T(1),Vol(1),P(1)
560 :
561 0 : dP_dr_00(1) = DPV(1)*DVR(1)
562 0 : dP_dr_in(1) = DPV(1)*DVRM(1)
563 0 : dCp_dr_00(1) = CPV(1)*DVR(1)
564 0 : dCp_dr_in(1) = CPV(1)*DVRM(1)
565 0 : dQQ_dr_00(1) = QQV(1)*DVR(1)
566 0 : dQQ_dr_in(1) = QQV(1)*DVRM(1)
567 :
568 : ! SET E_T (NOW =w) BELOW AND ABOVE BOUNDARIES
569 : ! w = EFL02 BELOW IBOTOM, INSREAD OF =0 BETTER, BECAUSE OF
570 : ! TURBULENT FLUX, NEAR THE IBOTOM
571 0 : Et(NZN) = 0.d0
572 0 : do I=1,IBOTOM
573 0 : Et(I) = 0.d0
574 : end do
575 :
576 0 : do I=1,NZN-1
577 0 : POM= (R(I)**2)/(2.d0*G*M(I))
578 0 : Hp_face(I)=POM*(P(I)*Vol(I)+P(I+1)*Vol(I+1))
579 0 : dHp_dr_in(I)=POM*(P(I)*DVRM(I)+Vol(I)*dP_dr_in(I))
580 : dHp_dr_00(I)=2.d0*Hp_face(I)/R(I)+POM*(P(I)*DVR(I)+Vol(I)*dP_dr_00(I) &
581 0 : +P(I+1)*DVRM(I+1)+Vol(I+1)*dP_dr_in(I+1))
582 0 : dHp_dr_out(I)=POM*(P(I+1)*DVR(I+1)+Vol(I+1)*dP_dr_00(I+1))
583 0 : dHp_dT_00(I)=POM*Vol(I)*dP_dT_00(I)
584 0 : dHp_dT_out(I)=POM*Vol(I+1)*dP_dT_00(I+1)
585 : end do
586 0 : POM=(R(NZN)**2)/(2.d0*G*M(NZN))
587 0 : Hp_face(NZN)=POM*P(NZN)*Vol(NZN)
588 0 : dHp_dr_in(NZN)=POM*(P(NZN)*DVRM(NZN)+Vol(NZN)*dP_dr_in(NZN))
589 : dHp_dr_00(NZN)=2.d0*Hp_face(NZN)/R(NZN)+POM* &
590 0 : (P(NZN)*DVR(NZN)+Vol(NZN)*dP_dr_00(NZN))
591 0 : dHp_dr_out(NZN)=0.d0
592 0 : dHp_dT_00(NZN)=POM*Vol(NZN)*dP_dT_00(NZN)
593 0 : dHp_dT_out(NZN)=0.d0
594 :
595 0 : do I=1,NZN-1
596 0 : POM=0.5d0*(QQS(I)/CPS(I)+QQS(I+1)/CPS(I+1))
597 0 : POM2=0.5d0*(P(I+1)-P(I))
598 : IGR1=0.5d0*(QQS(I)/CPS(I)+QQS(I+1)/CPS(I+1))*(P(I+1)-P(I)) &
599 0 : -(dlog(T(I+1))-dlog(T(I)))
600 : IGR1X0=POM2*((dQQ_dr_00(I)-QQS(I)/CPS(I)*dCp_dr_00(I))/CPS(I)+ &
601 : (dQQ_dr_in(I+1)-QQS(I+1)/CPS(I+1)*dCp_dr_in(I+1))/CPS(I+1)) &
602 0 : +POM*(dP_dr_in(I+1)-dP_dr_00(I))
603 : IGR1XM=POM2*(dQQ_dr_in(I)-QQS(I)/CPS(I)*dCp_dr_in(I))/CPS(I) &
604 0 : +POM*(-dP_dr_in(I))
605 : IGR1XP=POM2*(dQQ_dr_00(I+1)-QQS(I+1)/CPS(I+1)*dCp_dr_00(I+1))/CPS(I+1) &
606 0 : +POM*(dP_dr_00(I+1))
607 : IGR1Y0=POM2*(dQQ_dT_00(I)-QQS(I)/CPS(I)*dCp_dT_00(I))/CPS(I) &
608 0 : +POM*(-dP_dT_00(I))+1.d0/T(I)
609 : IGR1YP=POM2*(dQQ_dT_00(I+1)-QQS(I+1)/CPS(I+1)*dCp_dT_00(I+1))/CPS(I+1) &
610 0 : +POM*(dP_dT_00(I+1))-1.d0/T(I+1)
611 :
612 0 : POM=2.d0/(Vol(I)+Vol(I+1))
613 0 : POM2=8.d0*PI*(R(I)**2)/dm_bar(I)*Hp_face(I)
614 0 : IGR2=4.d0*PI*(R(I)**2)*Hp_face(I)*POM/dm_bar(I)
615 : IGR2X0=2.d0*IGR2/R(I)+IGR2/Hp_face(I)*dHp_dr_00(I) &
616 0 : -POM2/(Vol(I)+Vol(I+1))**2*(DVR(I)+DVRM(I+1))
617 : IGR2XM=-POM2/(Vol(I)+Vol(I+1))**2*DVRM(I) &
618 0 : +IGR2/Hp_face(I)*dHp_dr_in(I)
619 : IGR2XP=-POM2/(Vol(I)+Vol(I+1))**2*DVR(I+1) &
620 0 : +IGR2/Hp_face(I)*dHp_dr_out(I)
621 0 : IGR2Y0=IGR2/Hp_face(I)*dHp_dT_00(I)
622 0 : IGR2YP=IGR2/Hp_face(I)*dHp_dT_out(I)
623 :
624 0 : Y_face(I)=IGR1*IGR2
625 0 : dY_dr_00(I)=IGR1*IGR2X0+IGR2*IGR1X0
626 0 : dY_dr_in(I)=IGR1*IGR2XM+IGR2*IGR1XM
627 0 : dY_dr_out(I)=IGR1*IGR2XP+IGR2*IGR1XP
628 0 : dY_dT_00(I)=IGR1*IGR2Y0+IGR2*IGR1Y0
629 0 : dY_dT_out(I)=IGR1*IGR2YP+IGR2*IGR1YP
630 : end do
631 :
632 : ! PROTECT FROM 'ALWAYS ZERO' SOLUTION WHEN (Y>0)
633 : ! (SEEMS UNNECESSARY)
634 0 : do I=IBOTOM+1,NZN-1
635 : ! if((Y_face(I)+Y_face(I-1))>0.d0.and.Et(I)==0.d0)
636 : ! x Et(I)=1.d+6!1.d-6
637 0 : if(ALFA==0.d0) Et(I)=0.d0
638 : end do
639 :
640 0 : do I=1,NZN-1
641 : POM=sqrt(2.d0/3.d0)*0.5d0
642 : FF=POM*((E(I) +P(I) *Vol(I) )/T(I) &
643 0 : +(E(I+1)+P(I+1)*Vol(I+1))/T(I+1))
644 0 : FFY0=POM*(dE_dT_00(I)+dP_dT_00(I)*Vol(I)-(E(I)+P(I)*Vol(I))/T(I))/T(I)
645 : FFYP=POM*(dE_dT_00(I+1)+dP_dT_00(I+1)*Vol(I+1)-(E(I+1)+P(I+1)*Vol(I+1)) &
646 0 : /T(I+1))/T(I+1)
647 0 : FFXP=POM* ((DEV(I+1)+P(I+1))*DVR(I+1)+Vol(I+1)*dP_dr_00(I+1))/T(I+1)
648 : FFX0=POM*(((DEV(I) +P(I) )*DVR(I) +Vol(I) *dP_dr_00(I) )/T(I)+ &
649 : ((DEV(I+1)+P(I+1))*DVRM(I+1)+Vol(I+1)*dP_dr_in(I+1)) &
650 0 : /T(I+1))
651 0 : FFXM=POM* ((DEV(I) +P(I) )*DVRM(I) +Vol(I)*dP_dr_in(I))/T(I)
652 :
653 0 : POM=ALFAS*ALFA
654 0 : POM2=0.5d0*(CPS(I)+CPS(I+1))
655 0 : GG=POM*POM2*Y_face(I)
656 0 : GGXM=POM*(POM2*dY_dr_in(I)+Y_face(I)*0.5d0*dCp_dr_in(I))
657 0 : GGX0=POM*(POM2*dY_dr_00(I)+Y_face(I)*0.5d0*(dCp_dr_00(I)+dCp_dr_in(I+1)))
658 0 : GGXP=POM*(POM2*dY_dr_out(I)+Y_face(I)*0.5d0*dCp_dr_00(I+1))
659 0 : GGY0=POM*(POM2*dY_dT_00(I)+Y_face(I)*0.5d0*dCp_dT_00(I))
660 0 : GGYP=POM*(POM2*dY_dT_out(I)+Y_face(I)*0.5d0*dCp_dT_00(I+1))
661 0 : GPF=GG/FF
662 :
663 : ! corelation PI defined without e_t
664 0 : POM=1.d0
665 :
666 0 : PII(I)=POM*GG
667 0 : DPIIZ0(I)=0.d0
668 0 : dPII_dr_00(I)=POM*GGX0
669 0 : dPII_dr_in(I)=POM*GGXM
670 0 : dPII_dr_out(I)=POM*GGXP
671 0 : dPII_dT_00(I)=POM*GGY0
672 0 : dPII_dT_out(I)=POM*GGYP
673 : end do
674 :
675 0 : do I=IBOTOM+1,NZN-1
676 : ! SOURCE TERM
677 0 : POM=0.5d0*(PII(I)/Hp_face(I)+PII(I-1)/Hp_face(I-1))
678 0 : POM2=T(I)*P(I)*QQS(I)/CPS(I)
679 0 : POM3=sqrt(Et(I))
680 0 : SOURCE(I)=POM*POM2*POM3
681 0 : TEM1=POM2*POM3*0.5d0
682 0 : TEMI=-PII(I)/Hp_face(I)**2
683 0 : TEMM=-PII(I-1)/Hp_face(I-1)**2
684 :
685 : ! X -> w
686 : ! Y -> T
687 : ! Z -> R
688 :
689 : dsrc_dr_out(I)= TEM1*(dPII_dr_out(I)/Hp_face(I) &
690 0 : +TEMI*dHp_dr_out(I))
691 : dsrc_dr_00(I)= TEM1*(dPII_dr_00(I)/Hp_face(I)+dPII_dr_out(I-1)/Hp_face(I-1) &
692 0 : +TEMI*dHp_dr_00(I)+TEMM*dHp_dr_out(I-1))
693 : dsrc_dr_in(I)= TEM1*(dPII_dr_in(I)/Hp_face(I)+dPII_dr_00(I-1)/Hp_face(I-1) &
694 0 : +TEMI*dHp_dr_in(I)+TEMM*dHp_dr_00(I-1))
695 : dsrc_dr_in2(I)=TEM1*(dPII_dr_in(I-1)/Hp_face(I-1) &
696 0 : +TEMM*dHp_dr_in(I-1))
697 :
698 : dsrc_dT_out(I)= TEM1*(dPII_dT_out(I)/Hp_face(I) &
699 0 : +TEMI*dHp_dT_out(I))
700 : dsrc_dT_00(I)= TEM1*(dPII_dT_00(I)/Hp_face(I)+dPII_dT_out(I-1)/Hp_face(I-1) &
701 0 : +TEMI*dHp_dT_00(I)+TEMM*dHp_dT_out(I-1))
702 : dsrc_dT_in(I)= TEM1*(dPII_dT_00(I-1)/Hp_face(I-1) &
703 0 : +TEMM*dHp_dT_00(I-1))
704 :
705 0 : dsrc_dw_00(I)= POM*POM2*0.5d0/sqrt(Et(I))
706 :
707 0 : POM=POM*POM3
708 :
709 : dsrc_dT_00(I)=dsrc_dT_00(I)+ &
710 : POM/CPS(I)*(P(I)*QQS(I)+T(I)*QQS(I)*dP_dT_00(I)+T(I)*P(I)*dQQ_dT_00(I) &
711 0 : -T(I)*P(I)*QQS(I)/CPS(I)*dCp_dT_00(I))
712 : dsrc_dr_00(I)=dsrc_dr_00(I)+ &
713 : POM*T(I)/CPS(I)*(QQS(I)*dP_dr_00(I)+P(I)*dQQ_dr_00(I) &
714 0 : -P(I)*QQS(I)/CPS(I)*dCp_dr_00(I))
715 : dsrc_dr_in(I)=dsrc_dr_in(I)+ &
716 : POM*T(I)/CPS(I)*(QQS(I)*dP_dr_in(I)+P(I)*dQQ_dr_in(I) &
717 0 : -P(I)*QQS(I)/CPS(I)*dCp_dr_in(I))
718 :
719 : ! SOURCE TERM ALWAYS POSITIVE
720 : ! THIS IS FORMULATION USED IN BUDAPEST-FLORIDA CODE
721 : ! IT REQUIRES E_T AS MAIN VARIABLE - OTHERWISE CONVERGENCE
722 : ! IS EXTREMELY SLOW, AND POSSIBLE ONLY IF LOW ACCURACY
723 : ! FOR CONVERGENCE CONDITION, DXXC< 1d-2 - 1.d-4, IS SET
724 : ! if(SOURCE(I)<=0.d0)then
725 : ! SOURCE(I) = 0.d0
726 : ! dsrc_dr_out(I) = 0.d0
727 : ! dsrc_dr_00(I) = 0.d0
728 : ! dsrc_dr_in(I) = 0.d0
729 : ! dsrc_dr_in2(I)= 0.d0
730 : ! dsrc_dT_out(I) = 0.d0
731 : ! dsrc_dT_in(I) = 0.d0
732 : ! dsrc_dT_00(I) = 0.d0
733 : ! dsrc_dw_00(I) = 0.d0
734 : ! end if
735 :
736 : ! DAMP TERM
737 0 : POM=(CEDE/ALFA)*(Et(I)**1.5d0-EFL02**1.5d0)
738 0 : POM2=0.5d0*(Hp_face(I)+Hp_face(I-1))
739 0 : DAMP(I)=POM/POM2
740 0 : TEM1=-0.5d0*POM/POM2**2
741 0 : d_damp_dr_out(I) =TEM1*dHp_dr_out(I)
742 0 : d_damp_dr_00(I) =TEM1*(dHp_dr_00(I)+dHp_dr_out(I-1))
743 0 : d_damp_dr_in(I) =TEM1*(dHp_dr_in(I)+dHp_dr_00(I-1))
744 0 : d_damp_dr_in2(I)=TEM1*(dHp_dr_in(I-1))
745 0 : d_damp_dT_00(I) =TEM1*(dHp_dT_00(I)+dHp_dT_out(I-1))
746 0 : d_damp_dT_out(I) =TEM1*dHp_dT_out(I)
747 0 : d_damp_dT_in(I) =TEM1*dHp_dT_00(I-1)
748 0 : d_damp_dw_00(I) =1.5d0*(CEDE/ALFA)/POM2*sqrt(Et(I))
749 :
750 : ! RADIATIVE DAMP TERM
751 0 : if(GAMMAR==0.d0)then
752 0 : DAMPR(I) = 0.d0
753 0 : d_dampR_dr_out(I) = 0.d0
754 0 : d_dampR_dr_00(I) = 0.d0
755 0 : d_dampR_dr_in(I) = 0.d0
756 0 : d_dampR_dr_in2(I)= 0.d0
757 0 : d_dampR_dT_out(I) = 0.d0
758 0 : d_dampR_dT_in(I) = 0.d0
759 0 : d_dampR_dT_00(I) = 0.d0
760 0 : d_dampR_dw_00(I) = 0.d0
761 : else
762 0 : POM=(GAMMAR**2)/(ALFA**2)*4.d0*SIG
763 0 : POM2=T(I)**3*Vol(I)**2/(CPS(I)*K(I))
764 0 : POM3=Et(I)
765 0 : POM4=0.5d0*(Hp_face(I)**2+Hp_face(I-1)**2)
766 0 : DAMPR(I)=POM*POM2*POM3/POM4
767 0 : TEM1=-DAMPR(I)/POM4
768 0 : d_dampR_dr_out(I) =TEM1*(Hp_face(I)*dHp_dr_out(I))
769 : d_dampR_dr_00(I) =TEM1*(Hp_face(I)*dHp_dr_00(I) &
770 0 : +Hp_face(I-1)*dHp_dr_out(I-1))
771 : d_dampR_dr_in(I) =TEM1*(Hp_face(I)*dHp_dr_in(I) &
772 0 : +Hp_face(I-1)*dHp_dr_00(I-1))
773 0 : d_dampR_dr_in2(I)=TEM1*(Hp_face(I-1)*dHp_dr_in(I-1))
774 0 : d_dampR_dT_out(I) =TEM1*(Hp_face(I)*dHp_dT_out(I))
775 : d_dampR_dT_00(I) =TEM1*(Hp_face(I)*dHp_dT_00(I) &
776 0 : +Hp_face(I-1)*dHp_dT_out(I-1))
777 0 : d_dampR_dT_in(I) =TEM1*(Hp_face(I-1)*dHp_dT_00(I-1))
778 :
779 0 : d_dampR_dw_00(I)=POM*POM2/POM4
780 :
781 0 : TEM1=POM*POM3/POM4
782 : d_dampR_dr_00(I)=d_dampR_dr_00(I) &
783 : +TEM1*T(I)**3*(2.d0*Vol(I)*DVR(I) &
784 : -Vol(I)**2*( 1.d0/CPS(I)*dCp_dr_00(I) &
785 : +1.d0/K(I)*dK_dV_00(I)*DVR(I))) &
786 0 : /(CPS(I)*K(I))
787 :
788 : d_dampR_dr_in(I)=d_dampR_dr_in(I) &
789 : +TEM1*T(I)**3*(2.d0*Vol(I)*DVRM(I) &
790 : -Vol(I)**2*( 1.d0/CPS(I)*dCp_dr_in(I) &
791 : +1.d0/K(I)*dK_dV_00(I)*DVRM(I))) &
792 0 : /(CPS(I)*K(I))
793 :
794 : d_dampR_dT_00(I)=d_dampR_dT_00(I) &
795 : +TEM1*Vol(I)**2*(3.d0*T(I)**2 &
796 : -T(I)**3*( 1.d0/CPS(I)*dCp_dT_00(I) &
797 : +1.d0/K(I)*dK_dT_00(I))) &
798 0 : /(CPS(I)*K(I))
799 :
800 : end if
801 :
802 0 : dC_dr_00(I) =dsrc_dr_00(I) -d_damp_dr_00(I) -d_dampR_dr_00(I)
803 0 : dC_dr_out(I) =dsrc_dr_out(I) -d_damp_dr_out(I) -d_dampR_dr_out(I)
804 0 : dC_dr_in(I) =dsrc_dr_in(I) -d_damp_dr_in(I) -d_dampR_dr_in(I)
805 0 : dC_dr_in2(I)=dsrc_dr_in2(I)-d_damp_dr_in2(I)-d_dampR_dr_in2(I)
806 0 : dC_dT_in(I) =dsrc_dT_in(I) -d_damp_dT_in(I) -d_dampR_dT_in(I)
807 0 : dC_dT_00(I) =dsrc_dT_00(I) -d_damp_dT_00(I) -d_dampR_dT_00(I)
808 0 : dC_dT_out(I) =dsrc_dT_out(I) -d_damp_dT_out(I) -d_dampR_dT_out(I)
809 0 : dC_dw_00(I) =dsrc_dw_00(I) -d_damp_dw_00(I) -d_dampR_dw_00(I)
810 : end do
811 :
812 0 : do I=1,NZN
813 : ! CONVECTIVE LUMINOSITY
814 0 : if(I<IBOTOM.or.I>=NZN.or.alfa==0d0) then
815 0 : Lc(I)=0.d0
816 0 : DLCX0(I)=0.d0
817 0 : DLCXM(I)=0.d0
818 0 : DLCXP(I)=0.d0
819 0 : DLCY0(I)=0.d0
820 0 : DLCYP(I)=0.d0
821 0 : DLCZ0(I)=0.d0
822 0 : DLCZP(I)=0.d0
823 : else
824 : POM=4.d0*PI*(R(I)**2)*(T(I)/Vol(I)+T(I+1)/Vol(I+1))*0.5d0* &
825 0 : (ALFAC/ALFAS)
826 0 : POM3=(sqrt(Et(I))+sqrt(Et(I+1)))*0.5d0
827 0 : Lc(I)=POM*PII(I)*POM3
828 0 : DLCZ0(I)=0d0
829 0 : DLCZP(I)=0d0
830 0 : if (Et(I) /= 0d0) DLCZ0(I)=POM*PII(I)*0.25d0/sqrt(Et(I))
831 0 : if (Et(I+1) /= 0d0) DLCZP(I)=POM*PII(I)*0.25d0/sqrt(Et(I+1))
832 0 : POM2=4.d0*PI*(R(I)**2)*PII(I)*0.5d0*(ALFAC/ALFAS)*POM3
833 0 : POM=POM*POM3
834 : DLCX0(I)=POM*dPII_dr_00(I)+2.d0*Lc(I)/R(I) &
835 : -POM2*(T(I) /(Vol(I)**2 )*DVR(I)+ &
836 0 : T(I+1)/(Vol(I+1)**2)*DVRM(I+1))
837 0 : DLCXM(I)=POM*dPII_dr_in(I)-POM2*T(I) /(Vol(I)**2 )*DVRM(I)
838 0 : DLCXP(I)=POM*dPII_dr_out(I)-POM2*T(I+1)/(Vol(I+1)**2)*DVR(I+1)
839 0 : DLCY0(I)=POM*dPII_dT_00(I)+POM2/Vol(I)
840 0 : DLCYP(I)=POM*dPII_dT_out(I)+POM2/Vol(I+1)
841 :
842 0 : if(I==IBOTOM) DLCZ0(I)=0.d0
843 0 : if(I==NZN-1) DLCZP(I)=0.d0
844 :
845 0 : if(NEGFLU==0.and.PII(I)<0.d0)then
846 0 : Lc(I)=0.d0
847 0 : DLCX0(I)=0.d0
848 0 : DLCXM(I)=0.d0
849 0 : DLCXP(I)=0.d0
850 0 : DLCY0(I)=0.d0
851 0 : DLCYP(I)=0.d0
852 0 : DLCZ0(I)=0.d0
853 0 : DLCZP(I)=0.d0
854 : end if
855 0 : if(NEGFLU==1.and.PII(I)<0.d0)then
856 0 : if(II>300)then
857 0 : DLCZ0(I)=0.d0
858 0 : DLCZP(I)=0.d0
859 : end if
860 : end if
861 : end if
862 :
863 : ! TURBULENT LUMINOSITY
864 : if(I<IBOTOM.or.I>=NZN.or. &
865 0 : ALFAT==0.d0.or.ALFA==0.d0)then
866 0 : Lt(I)=0.d0
867 0 : DLTX0(I)=0.d0
868 0 : DLTXM(I)=0.d0
869 0 : DLTXP(I)=0.d0
870 0 : DLTY0(I)=0.d0
871 0 : DLTYP(I)=0.d0
872 0 : DLTZ0(I)=0.d0
873 0 : DLTZP(I)=0.d0
874 : else
875 0 : POM=-2.d0/3.d0*ALFA*ALFAT*(4.d0*PI*(R(I)**2))**2
876 0 : POM2=Hp_face(I)*(1.d0/Vol(I)**2+1.d0/Vol(I+1)**2)*0.5d0
877 0 : POM3=(Et(I+1)**1.5d0-Et(I)**1.5d0)/dm_bar(I)
878 0 : Lt(I)=POM*POM2*POM3
879 0 : if (is_bad(Lt(I))) then
880 0 : write(*,*) 'Lt(I)', I, Lt(I)
881 0 : write(*,*) 'POM', I, POM
882 0 : write(*,*) 'POM2', I, POM2
883 0 : write(*,*) 'POM3', I, POM3
884 0 : stop
885 : end if
886 : DLTX0(I)=4.d0*Lt(I)/R(I) &
887 : +Lt(I)/POM2*Hp_face(I)*(-1.d0/Vol(I)**3*DVR(I) &
888 : -1.d0/Vol(I+1)**3*DVRM(I+1)) &
889 0 : +Lt(I)/Hp_face(I)*dHp_dr_00(I)
890 :
891 : DLTXM(I)=Lt(I)/POM2*Hp_face(I)*(-1.d0/Vol(I)**3*DVRM(I)) &
892 0 : +Lt(I)/Hp_face(I)*dHp_dr_in(I)
893 : DLTXP(I)=Lt(I)/POM2*Hp_face(I)*(-1.d0/Vol(I+1)**3*DVR(I+1)) &
894 0 : +Lt(I)/Hp_face(I)*dHp_dr_out(I)
895 0 : DLTY0(I)=Lt(I)/Hp_face(I)*dHp_dT_00(I)
896 0 : DLTYP(I)=Lt(I)/Hp_face(I)*dHp_dT_out(I)
897 0 : DLTZ0(I)=-POM*POM2*1.5d0*sqrt(Et(I) )/dm_bar(I)
898 0 : DLTZP(I)= POM*POM2*1.5d0*sqrt(Et(I+1))/dm_bar(I)
899 : end if
900 : end do
901 :
902 : ! TURBULENT PRESSURE (ZONE)
903 0 : do I=IBOTOM+1,NZN-1
904 0 : if(ALFAP==0.d0)then
905 0 : PTURB(I) = 0.d0
906 0 : dPtrb_dw_00(I) = 0.d0
907 0 : dPtrb_dr_00(I) = 0.d0
908 0 : dPtrb_dr_in(I) = 0.d0
909 : else
910 0 : PTURB(I) = ALFAP*Et(I)/Vol(I)
911 0 : dPtrb_dw_00(I) = ALFAP/Vol(I)
912 0 : TEM1=-ALFAP*Et(I)/Vol(I)**2
913 0 : dPtrb_dr_00(I) = TEM1*DVR(I)
914 0 : dPtrb_dr_in(I) = TEM1*DVRM(I)
915 : end if
916 : end do
917 :
918 0 : do I=1,IBOTOM
919 0 : PTURB(I)= 0.d0
920 0 : dPtrb_dw_00(I)= 0.d0
921 0 : dPtrb_dr_00(I)= 0.d0
922 0 : dPtrb_dr_in(I)= 0.d0
923 0 : dC_dr_00(I) = 0.d0
924 0 : dC_dr_out(I) = 0.d0
925 0 : dC_dr_in(I) = 0.d0
926 0 : dC_dr_in2(I)= 0.d0
927 0 : dC_dT_in(I) = 0.d0
928 0 : dC_dT_00(I) = 0.d0
929 0 : dC_dT_out(I) = 0.d0
930 0 : dC_dw_00(I) = 0.d0
931 : end do
932 0 : do I=1,IBOTOM-1
933 0 : DLCX0(I)= 0.d0
934 0 : DLCXM(I)= 0.d0
935 0 : DLCXP(I)= 0.d0
936 0 : DLCY0(I)= 0.d0
937 0 : DLCYP(I)= 0.d0
938 0 : DLCZ0(I)= 0.d0
939 0 : DLCZP(I)= 0.d0
940 0 : DLTX0(I)= 0.d0
941 0 : DLTXM(I)= 0.d0
942 0 : DLTXP(I)= 0.d0
943 0 : DLTY0(I)= 0.d0
944 0 : DLTYP(I)= 0.d0
945 0 : DLTZ0(I)= 0.d0
946 0 : DLTZP(I)= 0.d0
947 : end do
948 0 : PTURB(NZN)= 0.d0
949 0 : dPtrb_dw_00(NZN)= 0.d0
950 0 : dPtrb_dr_00(NZN)= 0.d0
951 0 : dPtrb_dr_in(NZN)= 0.d0
952 0 : DLCX0(NZN)= 0.d0
953 0 : DLCXM(NZN)= 0.d0
954 0 : DLCXP(NZN)= 0.d0
955 0 : DLCY0(NZN)= 0.d0
956 0 : DLCYP(NZN)= 0.d0
957 0 : DLCZ0(NZN)= 0.d0
958 0 : DLCZP(NZN)= 0.d0
959 0 : DLTX0(NZN)= 0.d0
960 0 : DLTXM(NZN)= 0.d0
961 0 : DLTXP(NZN)= 0.d0
962 0 : DLTY0(NZN)= 0.d0
963 0 : DLTYP(NZN)= 0.d0
964 0 : DLTZ0(NZN)= 0.d0
965 0 : DLTZP(NZN)= 0.d0
966 0 : dC_dr_00(NZN) = 0.d0
967 0 : dC_dr_out(NZN) = 0.d0
968 0 : dC_dr_in(NZN) = 0.d0
969 0 : dC_dr_in2(NZN)= 0.d0
970 0 : dC_dT_in(NZN) = 0.d0
971 0 : dC_dT_00(NZN) = 0.d0
972 0 : dC_dT_out(NZN) = 0.d0
973 0 : dC_dw_00(NZN) = 0.d0
974 :
975 0 : if(ALFA==0.d0) then !RADIATIVE CASE
976 0 : SOURCE(I) = 0.d0
977 0 : DAMP(I) = 0.d0
978 0 : DAMPR(I) = 0.d0
979 0 : dC_dr_00(I) = 0.d0
980 0 : dC_dr_out(I) = 0.d0
981 0 : dC_dr_in(I) = 0.d0
982 0 : dC_dr_in2(I) = 0.d0
983 0 : dC_dT_in(I) = 0.d0
984 0 : dC_dT_00(I) = 0.d0
985 0 : dC_dT_out(I) = 0.d0
986 0 : dC_dw_00(I) = 0.d0
987 : end if
988 :
989 :
990 : ! INITIALIZE HD(11,3*NZN)
991 0 : do I=1,3*NZN
992 0 : do J=1,11
993 0 : HD(J,I)=0.d0
994 : end do
995 : end do
996 :
997 : ! LOOP 2 .. LUM PLUGS
998 : DLR = 0.d0
999 0 : DLRP = 0.d0
1000 0 : DLRM = 0.d0
1001 0 : DLT = 0.d0
1002 0 : DLTP = 0.d0
1003 0 : DLR = 0.d0 !!-1.d0
1004 0 : do 5 I=1,NZN
1005 0 : IR = 3+3*(I-1)
1006 0 : IC = 2+3*(I-1)
1007 0 : IW = 1+3*(I-1)
1008 : ! SET LUM(I-1)
1009 0 : DLMR = DLR
1010 0 : DLMRP = DLRP
1011 0 : DLMRM = DLRM
1012 0 : DLMT = DLT
1013 0 : DLMTP = DLTP
1014 0 : if(I==NZN) GOTO 6
1015 : ! Lr(I)=Eq. A.4, Stellingwerf 1975, Appendix A
1016 : ! CALC LUM(I)
1017 0 : W_00=T(I)**4
1018 0 : W_out=T(I+1)**4
1019 0 : BW=dlog(W_out/W_00)
1020 0 : BK=dlog(K(I+1)/K(I))
1021 0 : T1=-CL*R(I)**4/dm_bar(I)
1022 0 : T2=(W_out/K(I+1)-W_00/K(I))/(1.d0-BK/BW)
1023 0 : Lr(I)=T1*T2
1024 0 : T3=T1/(BW-BK)
1025 0 : DLK= (T3/K(I)) *(W_00*BW/K(I) -T2) !dL(i)/dK(i)
1026 0 : DLKP=-(T3/K(I+1))*(W_out*BW/K(I+1)-T2) !dL(i)/dK(i+1)
1027 0 : DLRP= DLKP*dK_dV_00(I+1)*DVR(I+1)
1028 0 : DLRM= DLK *dK_dV_00(I) *DVRM(I)
1029 0 : DLR= 4.d0*T1*T2/R(I)+DLK*dK_dV_00(I)*DVR(I)+DLKP*dK_dV_00(I+1)*DVRM(I+1)
1030 0 : DLTP=4.d0*(T3/T(I+1))*(W_out*BW/K(I+1)-T2*BK/BW)+DLKP*dK_dT_00(I+1)
1031 0 : DLT=-4.d0*(T3/T(I))*(W_00*BW/K(I)-T2*BK/BW)+DLK*dK_dT_00(I)
1032 0 : GOTO 7
1033 : ! OUTER LUM BOUNDARY CONDITION
1034 : 6 continue
1035 0 : Lr(I)=L0
1036 0 : DLT = 4.d0*L0/T(I) !L=4piR^2sigT^4
1037 0 : DLR = 2.d0*L0/R(I) !L=4piR^2sigT^4
1038 0 : DLRM = 0.d0
1039 0 : DLRP = 0.d0
1040 0 : DLTP = 0.d0
1041 : 7 continue
1042 : ! CALC ENERGY EQUATION(I)
1043 0 : HR(IW) = -(Lr(I)/L0+Lc(I)/L0+Lt(I)/L0-1.d0)
1044 : !write(*,*) 'energy', I, HR(IW)
1045 0 : if (is_bad(HR(IW))) then
1046 0 : write(*,*) 'L0', I, L0
1047 0 : write(*,*) 'Lt(I)', I, Lt(I)
1048 0 : write(*,*) 'Lc(I)', I, Lc(I)
1049 0 : write(*,*) 'Lr(I)', I, Lr(I)
1050 0 : stop
1051 : end if
1052 0 : HD(6,IW) = (DLT +DLCY0(I)+DLTY0(I))/L0 !Y(i)
1053 0 : HD(9,IW) = (DLTP+DLCYP(I)+DLTYP(I))/L0 !Y(i+1)
1054 0 : HD(5,IW) = (DLRM+DLCXM(I)+DLTXM(I))/L0 !X(i-1)
1055 0 : HD(8,IW) = (DLR +DLCX0(I)+DLTX0(I))/L0 !X(i)
1056 0 : HD(11,IW) = (DLRP+DLCXP(I)+DLTXP(I))/L0 !X(i+1)
1057 0 : HD(7,IW) = ( DLCZ0(I)+DLTZ0(I))/L0 !Z(i)
1058 0 : HD(10,IW) = ( DLCZP(I)+DLTZP(I))/L0 !Z(i+1)
1059 : ! CALC MOMEMTUM EQUATION(I)
1060 0 : T1=-P4*(R(I)**2)/dm_bar(I)
1061 0 : if(I==NZN) then
1062 0 : HR(IR) = -T1*(Psurf-P(I)-PTURB(I))+G*M(I)/R(I)**2
1063 : else
1064 0 : HR(IR)=-T1*(P(I+1)-P(I)+PTURB(I+1)-PTURB(I))+G*M(I)/R(I)**2
1065 : end if
1066 : !write(*,*) 'momentum', I, HR(IR)
1067 0 : HD(3,IR) = +T1*(-dP_dr_in(I)-dPtrb_dr_in(I)) !X(i-1)
1068 0 : HD(4,IR) = +T1*(-dP_dT_00(I)) !Y(i)
1069 :
1070 0 : HD(2,IR) = 0.d0 !Z(i-1)
1071 0 : HD(5,IR) = +T1*(-dPtrb_dw_00(I)) !Z(i)
1072 0 : HD(8,IR) = +T1*(dPtrb_dw_00(I+1)) !Z(i+1)
1073 :
1074 0 : if(I==NZN) GOTO 111
1075 : HD(6,IR) = +T1*(dP_dr_in(I+1)-dP_dr_00(I)+dPtrb_dr_in(I+1)-dPtrb_dr_00(I))+ &
1076 0 : 4.d0*G*M(I)/(R(I)**3) !X(i)
1077 0 : HD(7,IR) = +T1*(dP_dT_00(I+1)) !Y(i+1)
1078 0 : HD(9,IR) = +T1*(dP_dr_00(I+1)+dPtrb_dr_00(I+1)) !X(i+1)
1079 0 : GOTO 112
1080 0 : 111 HD(7,IR)= 0.d0
1081 0 : HD(9,IR)= 0.d0
1082 0 : HD(6,IR)= +T1*(-dP_dr_00(I))+4.d0*G*M(I)/(R(I)**3)
1083 : 112 continue
1084 :
1085 : ! CALC CONVECTIVE ENERGY EQUATION
1086 0 : if(I<=IBOTOM.or.I==NZN.or.ALFA==0.d0)then
1087 0 : do J=1,11
1088 0 : HD(J,IC)=0.d0
1089 : end do
1090 0 : HD(6,IC)=1.d0
1091 0 : HR(IC)=0.d0
1092 : else
1093 0 : T1=-1.d0/dm(I)
1094 0 : HR(IC)= -(T1*(Lt(I)-Lt(I-1))+SOURCE(I)-DAMP(I)-DAMPR(I))
1095 : !write(*,*) 'conv energy', I, HR(IC)
1096 0 : if (is_bad(HR(IC))) then
1097 0 : write(*,*) 'Lt(I', I, Lt(I)
1098 0 : write(*,*) 'Lt(I-1)', I-1, Lt(I-1)
1099 0 : write(*,*) 'SOURCE(I)', I, SOURCE(I)
1100 0 : write(*,*) 'DAMP(I)', I, DAMP(I)
1101 0 : write(*,*) 'DAMPR(I)', I, DAMPR(I)
1102 0 : stop
1103 : end if
1104 0 : HD(3,IC) = +T1*( -DLTZ0(I-1)) !Z(i-1)
1105 0 : HD(6,IC) = dC_dw_00(I) +T1*(DLTZ0(I)-DLTZP(I-1)) !Z(i)
1106 0 : HD(9,IC) = +T1*(DLTZP(I)) !Z(i+1)
1107 0 : HD(1,IC) = dC_dr_in2(I)+T1*( -DLTXM(I-1)) !X(i-2)
1108 0 : HD(4,IC) = dC_dr_in(I) +T1*(DLTXM(I)-DLTX0(I-1)) !X(i-1)
1109 0 : HD(7,IC) = dC_dr_00(I) +T1*(DLTX0(I)-DLTXP(I-1)) !X(i)
1110 0 : HD(10,IC) = dC_dr_out(I) +T1*(DLTXP(I)) !X(i+1)
1111 0 : HD(2,IC) = dC_dT_in(I) +T1*( -DLTY0(I-1)) !Y(i-1)
1112 0 : HD(5,IC) = dC_dT_00(I) +T1*(DLTY0(I)-DLTYP(I-1)) !Y(i)
1113 0 : HD(8,IC) = dC_dT_out(I) +T1*(DLTYP(I)) !Y(i+1)
1114 : end if
1115 :
1116 : ! SET ANCHOR T(NZN)=CONST
1117 0 : if(I==NZN) then
1118 0 : do j=1,11
1119 0 : HD(J,IW) = 0.d0
1120 : end do
1121 : ! SET dT(NZN)=0
1122 0 : HD(6,IW) = 1.d0
1123 0 : HR(IW) = 0.d0
1124 : ! SET DERIVATIVES d(*)/dT(NZN)=0 IN ZONE/ITERFACE NZN
1125 0 : HD(4,IR) = 0.d0
1126 0 : HD(5,IC) = 0.d0
1127 : end if
1128 0 : if(I==NZN-1)then
1129 : ! SET DERIVATIVES d(*)/dT(NZN)=0 IN ZONE/ITERFACE NZN-1
1130 0 : HD(9,IW) = 0.d0
1131 0 : HD(7,IR) = 0.d0
1132 0 : HD(8,IC) = 0.d0
1133 : end if
1134 :
1135 0 : 5 continue
1136 :
1137 0 : if(ending) GOTO 889
1138 :
1139 0 : do J=1,5 !translate hd into band storage of LAPACK/LINPACK
1140 0 : do I=1,3*NZN+1
1141 0 : ABB(J,I)=0.0d0
1142 : end do
1143 : end do
1144 0 : do J=1,5
1145 0 : do I=1,3*NZN+1-J
1146 0 : ABB(11-J,I+J)=HD(6+J,I) !upper diagonals
1147 0 : ABB(11+J,I)=HD(6-J,I+J) !lower diagonals
1148 : end do
1149 : end do
1150 0 : do I=1,3*NZN+1
1151 0 : ABB(11,I)=HD(6,I)
1152 : end do
1153 :
1154 : !- LAPACK
1155 0 : call DGBTRF(3*NZN,3*NZN,5,5,ABB,LD_ABB,IPVT,INFO)
1156 0 : if(INFO/=0) then
1157 0 : write(*,*) 'hyd: LAPACK/dgbtrf problem no., iter.',INFO,II
1158 0 : open(61,file='x.dat',status='unknown')
1159 0 : write(61,'(F6.3,tr2,f8.2,tr2,f7.2,tr2,d9.3)')EMR,ELR,TE-0.5d0, &
1160 0 : gam
1161 0 : close(61)
1162 0 : open(61,file='logic.dat',status='unknown')
1163 0 : write(61,'(I1)') 3
1164 0 : close(61)
1165 0 : stop
1166 : end if
1167 0 : call DGBTRS('n',3*NZN,5,5,1,ABB,LD_ABB,IPVT,HR,3*NZN,INFO)
1168 0 : if(INFO/=0) then
1169 0 : write(*,*) 'hyd: LAPACK/dgbtrs problem no., iter.',INFO,II
1170 0 : stop
1171 : end if
1172 : !-
1173 0 : do I=1,3*NZN,1
1174 0 : DX(I) = HR(I)
1175 : end do
1176 :
1177 0 : EZH = 1.d0
1178 0 : XXR = 0.d0
1179 0 : XXT = 0.d0
1180 0 : XXC = 0.d0
1181 0 : XX_max = 0
1182 0 : XX_max_val = 0
1183 0 : XX_max_dx = 0
1184 0 : i_XX_max = 0
1185 0 : var_XX_max = 0
1186 0 : do I=2,NZN
1187 0 : IR = 3+3*(I-1)
1188 0 : IC = IR-1
1189 0 : IW = IR-2
1190 0 : if(Et(I)>(1.d-20)*EFL02) XXC=dabs(DX(IC)/Et(I))
1191 0 : XXR=((DX(IR-3)-DX(IR))/(R(I)-R(I-1)))/DXH
1192 0 : XXT=dabs(DX(IW)/T(I))/DXH
1193 :
1194 : if (XXC > XX_max) then
1195 : XX_max = XXC
1196 : XX_max_val = Et(I)
1197 : XX_max_dx = DX(IC)
1198 : i_XX_max = i
1199 : var_XX_max = 1
1200 : end if
1201 : if (XXR > XX_max) then
1202 : XX_max = XXR
1203 : XX_max_val = R(I)
1204 : XX_max_dx = DX(IR)
1205 : i_XX_max = i
1206 : var_XX_max = 2
1207 : end if
1208 : if (XXT > XX_max) then
1209 : XX_max = XXT
1210 : XX_max_val = T(I)
1211 : XX_max_dx = DX(IW)
1212 : i_XX_max = i
1213 : var_XX_max = 3
1214 : end if
1215 :
1216 0 : EZH=1.d0/dmax1(1.d0/EZH,XXR,XXT,XXC)
1217 : end do
1218 :
1219 0 : DXXR = 0.d0
1220 0 : DXXT = 0.d0
1221 0 : DXXC = 0.d0
1222 0 : ITROUBT = 0
1223 0 : ITROUBC = 0
1224 : ! It seems that for BL Her models fixed undercorrection factor works best
1225 0 : do I=1,NZN
1226 0 : IW=1+3*(I-1)
1227 0 : IR=IW+2
1228 0 : IC=IW+1
1229 0 : T(I) = T(I) +EZH*DX(IW)/2d0
1230 0 : R(I) = R(I) +EZH*DX(IR)/2d0
1231 0 : Et(I) = Et(I)+EZH*DX(IC)/2d0
1232 0 : if(Et(I)<=0.d0) then
1233 0 : Et(I)=EFL02*1d-4*rand(s)
1234 : end if
1235 0 : DXKT=DXXT
1236 0 : DXKC=DXXC
1237 0 : DXXR=dmax1(DXXR,dabs(DX(IR)/R(I)))
1238 0 : DXXT=dmax1(DXXT,dabs(DX(IW)/T(I)))
1239 0 : if(Et(I)>(1.d-20)*EFL02) &
1240 0 : DXXC=dmax1(DXXC,dabs(DX(IC)/Et(I)))
1241 0 : if(DXXC>DXKC) ITROUBC=I
1242 0 : if(DXXT>DXKT) ITROUBT=I
1243 : end do
1244 : !write(*,*) 'apply corr', II,dabs(DXXC),dabs(DXXR),dabs(DXXT),ITROUBC, &
1245 : ! Et(ITROUBC),EZH
1246 : !call mesa_error(__FILE__,__LINE__,'debug')
1247 : !write(*,*) 'DXXC/PRECR', ITROUBC, DXXC/PRECR, DXXC
1248 0 : if(dabs(DXXC)<PRECR.and.dabs(DXXR)<PRECR.and. &
1249 0 : dabs(DXXT)<PRECR) then
1250 0 : IZIP=1
1251 0 : IOP=IOP+1
1252 0 : GOTO 999
1253 : end if
1254 :
1255 : end do
1256 :
1257 0 : write(*,'(A)')
1258 0 : write(*,*) ' NO CONVERGENCE IN RELAX_ENV, ITERATION: ',II
1259 0 : write(*,*) ' try increasing RSP_relax_dm_tolerance', s% RSP_relax_dm_tolerance
1260 0 : write(*,'(A)')
1261 0 : ierr = -1
1262 0 : return
1263 0 : stop
1264 :
1265 : 889 continue
1266 :
1267 :
1268 : !MAXFLUXD=-2d10
1269 : !do I=1,NZN
1270 : ! if(abs(HR(1+3*(I-1)))>MAXFLUXD) &
1271 : ! MAXFLUXD=abs(HR(1+3*(I-1)))
1272 : !end do
1273 : !write(*,*) 'MAXFLUXD= ',MAXFLUXD
1274 :
1275 0 : POM=MSTAR-M(1)+dm(1)
1276 0 : if (s% RSP_trace_RSP_build_model) then
1277 0 : write(*,*) ' inner temper. rel. diff',(T(1)-TINNER)/TINNER
1278 0 : write(*,*) ' inner radius rel. diff',(R(1)-RINNER)/RINNER
1279 0 : write(*,*) ' outer radius rel. diff',(R(NZN)-ROUTER)/ROUTER
1280 0 : write(*,*) ' envelope mass rel. diff',(POM-MENVEL)/MENVEL
1281 : end if
1282 : 222 format(A,d14.8)
1283 : if (.false.) then
1284 : MAXR=-1d0
1285 : MAXW=-1d0
1286 : MAXC=-1d0
1287 : write(*,*) 'static structure check'
1288 : !open(15,file='ss_lin.dat',status='unknown')
1289 : do J=1,NZN
1290 : IR=4+3*(J-1)
1291 : IW=2+3*(J-1)
1292 : PB=P(J)+PTURB(J)
1293 : if (J == NZN) then
1294 : PPB = 0
1295 : else
1296 : PPB=P(J+1)+PTURB(J+1)
1297 : end if
1298 : T1=P4*(R(J)**2)
1299 : T2=(PPB-PB)/dm_bar(J)
1300 : T3=G*M(J)/(R(J)**2)
1301 : if (dabs((T3+T1*T2)/T3)>MAXR)then
1302 : MAXR=dabs((T3+T1*T2)/T3)
1303 : IMAXR=J
1304 : end if
1305 : ELB=Lc(J)+Lr(J)+Lt(J)
1306 : if (J == 1) then
1307 : ELMB=0
1308 : else
1309 : ELMB=Lc(J-1)+Lr(J-1)+Lt(J-1)
1310 : end if
1311 : if (dabs((ELB-L0)/L0)>MAXW)then
1312 : MAXW=dabs((ELB-L0)/L0)
1313 : IMAXW=J
1314 : end if
1315 : if (J > IBOTOM .and. J < NZN) then
1316 : ELB=Lt(J)
1317 : ELMB=Lt(J-1)
1318 : T1=(ELB-ELMB)/dm(J)
1319 : T2=SOURCE(J)-DAMP(J)-DAMPR(J)
1320 : if ((T2/=0.d0).and.(T1/=0.d0))then
1321 : if (dabs(T1-T2)/max(T1,T2)>MAXC)then
1322 : MAXC=dabs(T1-T2)/max(T1,T2)
1323 : IMAXC=J
1324 : end if
1325 : end if
1326 : end if
1327 : end do
1328 : if (MAXR/=-1d0)write(*,*) 'MAX DIFFERENCE R: ',MAXR,' ZONE: ',IMAXR
1329 : if (MAXW/=-1d0)write(*,*) 'MAX DIFFERENCE W: ',MAXW,' ZONE: ',IMAXW
1330 : if (MAXC/=-1d0)write(*,*) 'MAX DIFFERENCE C: ',MAXC,' ZONE: ',IMAXC
1331 : end if
1332 : return
1333 0 : end subroutine RELAX_ENV
1334 :
1335 : end module rsp_relax_env
|
