Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 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_lina
21 : use star_def, only: star_info
22 : use utils_lib, only: is_bad, mesa_error
23 : use const_def, only: dp, i8, crad
24 : use rsp_def
25 :
26 : implicit none
27 :
28 : private
29 : public :: mesa_eos_kap, SORT, do_LINA
30 :
31 : contains
32 :
33 0 : subroutine do_LINA(s, L0, NZN, NMODES, VEL, PERS, ETO, &
34 0 : M, DM, DM_BAR, R, Vol, T, Et, Lr, ierr)
35 : ! LINear Analysis (LINA) USING LAPACK/EISPACK
36 : ! SETS THE MATRIX FOR THE LINEAR EIGENVALUE PROBLEM
37 :
38 : type (star_info), pointer :: s
39 : real(dp), intent(in) :: L0
40 : integer, intent(in) :: NMODES, NZN
41 : real(dp), intent(out) :: VEL(:,:) ! (NZN+1,15)
42 : real(dp), intent(out), dimension(15) :: PERS, ETO
43 : real(dp), intent(inout), dimension(:) :: &
44 : M, DM, DM_BAR, R, Vol, T, Et, Lr
45 : integer, intent(out) :: ierr
46 :
47 0 : real(dp), dimension(15) :: OMEG, EK
48 : real(dp), allocatable, dimension(:) :: &
49 0 : E, P, Lc, Hp_face, Y_face, K, CPS, QQS, &
50 0 : MX10,MX00,MX01,MY00,MY01, &
51 0 : MU10,MU00,MU01,MZ00,MZ01, &
52 0 : EVUU0,EVUUM, &
53 0 : EX20,EX10,EX00,EX01, &
54 0 : EY10,EY00,EY01, &
55 0 : EU10,EU00,EZ10,EZ00,EZ01, &
56 0 : CX10,CX00,CX01,CY00,CY01, &
57 0 : CZ00,CX20,CY10,CU00,CU10, &
58 0 : CZ01,CZ10, &
59 0 : dC_dr_in, dC_dr_00, dC_dr_out, dC_dT_out, &
60 0 : dC_dT_00, dC_dT_in,dC_dw_00,dC_dr_in2, &
61 0 : DPV,dP_dT_00,DEV,dE_dT_00,dK_dV_00, &
62 0 : dK_dT_00,DVR,DVRM, &
63 0 : ELRP,ELRM,ELR,ELTP,ELT, &
64 0 : ELZ0,ELZP, &
65 0 : CPV,dCp_dT_00,QQV,dQQ_dT_00 , &
66 0 : dP_dr_00, dCp_dr_00, dQQ_dr_00, &
67 0 : dP_dr_in,dCp_dr_in,dQQ_dr_in, &
68 0 : dHp_dr_in,dHp_dr_out,dHp_dr_00,dHp_dT_00,dHp_dT_out, &
69 0 : dY_dr_in,dY_dr_00,dY_dr_out,dY_dT_00, &
70 0 : dY_dT_out, &
71 0 : PII,dPII_dr_in,dPII_dr_00,dPII_dr_out, &
72 0 : dPII_dT_00,dPII_dT_out,DPIIZ0, &
73 0 : dsrc_dr_in,dsrc_dr_00,dsrc_dr_out,dsrc_dT_00, &
74 0 : dsrc_dT_out,dsrc_dw_00,source, &
75 0 : d_damp_dr_in,d_damp_dr_00,d_damp_dr_out,d_damp_dT_00, &
76 0 : d_damp_dT_out,d_damp_dw_00,DAMP, &
77 0 : dsrc_dr_in2,dsrc_dT_in,d_damp_dr_in2,d_damp_dT_in, &
78 0 : d_dampR_dr_in,d_dampR_dr_00,d_dampR_dr_out,d_dampR_dT_00, &
79 0 : d_dampR_dT_out,d_dampR_dw_00,DAMPR, &
80 0 : d_dampR_dr_in2,d_dampR_dT_in, &
81 0 : DLCXM,DLCX0,DLCXP,DLCY0, &
82 0 : DLCYP,DLCZ0,DLCZP, &
83 0 : DLTXM,DLTX0,DLTXP,DLTY0, &
84 0 : DLTYP,DLTZ0,DLTZP,Lt, &
85 0 : PTURB,dPtrb_dw_00,dPtrb_dr_00,dPtrb_dr_in,XXL
86 : real(dp), allocatable, dimension(:,:) :: &
87 0 : QWK,QWKEV,PHR,PHT,PHL,PHU,PHC,ZZZ
88 : complex(8), allocatable, dimension(:,:) :: &
89 0 : VRR,VRT,VRL,VRC,DVRR,DVRT,DVRL,DVRC,VRU
90 :
91 0 : real(dp) :: FFXM,FFX0,FFXP,FFY0,FFYP,FF,POM1
92 0 : real(dp) :: IGR1,IGR1XM,IGR1X0,IGR1XP,IGR1Y0,IGR1YP
93 0 : real(dp) :: IGR2,IGR2XM,IGR2X0,IGR2XP,IGR2Y0,IGR2YP
94 0 : real(dp) :: GGXM,GGX0,GGXP,GGY0,GGYP,GG,GPF
95 0 : real(dp) :: DFCX0,DFCXM,DFCXP,DFCY0,DFCYP,DFCZ0,DFCZP
96 0 : real(dp) :: DFCMX0,DFCMXM,DFCMXP,DFCMY0,DFCMYP,DFCMZ0,DFCMZP
97 0 : real(dp) :: T1,DLR,DLRP,DLRM,DLT,DLTP,DLMR,DLMRP,DLMRM,DLMT,DLMTP
98 0 : real(dp) :: W_00,W_out,BW,BK,T2,T3,DLK,DLKP,MINI
99 0 : real(dp) :: T4,POM3,POM2,POM,POM4
100 : integer :: I,J,NZN3,IG,IE,IR,IC,INFO,IMI,LD_VL,LD_VR,n,op_err
101 : real(dp) :: VRRS(15),Q(15)
102 : character (len=250) :: FILENAME
103 : character (len=1) :: NUMER1
104 : character (len=2) :: NUMER2
105 0 : complex(8):: DP_0,DV_0,VTTS(15),SCALE(15),DPEV,dP_dT_00URB
106 0 : real(dp) :: SGRP,SGRM
107 0 : real(dp) :: PSIG,TEMI,TEMM,TEM1
108 0 : real(dp) :: NORMC
109 0 : real(dp) :: QCHECK(15),ETOIEV(15),QWKPT(1000,15),ETOIPT(15),SGR(15)
110 0 : real(dp) :: EFL02
111 0 : real(dp) :: ETOI(15)
112 :
113 : !write(*,'(a55,i12,99(1pd26.16))') 'start LINA s% w(2)**2', &
114 : ! 2, s% w(2)**2, Et(NZN-1)
115 :
116 : if (.false.) then
117 : write(*,*) 'L0', L0
118 : do I=1,NZN
119 : write(*,*) 'M', I, M(I)
120 : write(*,*) 'DM', I, DM(I)
121 : write(*,*) 'DM_BAR', I, DM_BAR(I)
122 : write(*,*) 'R', I, R(I)
123 : write(*,*) 'R**3', I, R(I)**3
124 : write(*,*) 'Vol', I, Vol(I)
125 : write(*,*) 'T', I, T(I)
126 : write(*,*) 'Lr', I, Lr(I)
127 : write(*,*) 'Et', I, Et(I)
128 : end do
129 : call mesa_error(__FILE__,__LINE__,'do_LINA')
130 : end if
131 :
132 0 : ierr = 0
133 0 : EFL02 = EFL0*EFL0
134 0 : n = NZN+1
135 : allocate( &
136 : E(n), P(n), Lc(n), Hp_face(n), Y_face(n), K(n), CPS(n), QQS(n), &
137 : MX10(n),MX00(n),MX01(n),MY00(n),MY01(n), &
138 : MU10(n),MU00(n),MU01(n),MZ00(n),MZ01(n), &
139 : EVUU0(n),EVUUM(n), &
140 : EX20(n),EX10(n),EX00(n),EX01(n), &
141 : EY10(n),EY00(n),EY01(n), &
142 : EU10(n),EU00(n),EZ10(n),EZ00(n),EZ01(n), &
143 : CX10(n),CX00(n),CX01(n),CY00(n),CY01(n), &
144 : CZ00(n),CX20(n),CY10(n),CU00(n),CU10(n), &
145 : CZ01(n),CZ10(n), &
146 : dC_dr_in(n), dC_dr_00(n), dC_dr_out(n), dC_dT_out(n), &
147 : dC_dT_00(n), dC_dT_in(n),dC_dw_00(n),dC_dr_in2(n), &
148 0 : DPV(n),dP_dT_00(n),DEV(n),dE_dT_00(n),dK_dV_00(n), &
149 : dK_dT_00(n),DVR(n),DVRM(n), &
150 : ELRP(n),ELRM(n),ELR(n),ELTP(n),ELT(n), &
151 : ELZ0(n),ELZP(n), &
152 0 : CPV(n),dCp_dT_00(n),QQV(n),dQQ_dT_00(n) , &
153 0 : dP_dr_00(n), dCp_dr_00(n), dQQ_dr_00(n), &
154 0 : dP_dr_in(n),dCp_dr_in(n),dQQ_dr_in(n), &
155 : dHp_dr_in(n),dHp_dr_out(n),dHp_dr_00(n),dHp_dT_00(n),dHp_dT_out(n), &
156 : dY_dr_in(n),dY_dr_00(n),dY_dr_out(n),dY_dT_00(n), &
157 : dY_dT_out(n), &
158 : PII(n),dPII_dr_in(n),dPII_dr_00(n),dPII_dr_out(n), &
159 : dPII_dT_00(n),dPII_dT_out(n),DPIIZ0(n), &
160 : dsrc_dr_in(n),dsrc_dr_00(n),dsrc_dr_out(n),dsrc_dT_00(n), &
161 : dsrc_dT_out(n),dsrc_dw_00(n),source(n), &
162 : d_damp_dr_in(n),d_damp_dr_00(n),d_damp_dr_out(n),d_damp_dT_00(n), &
163 : d_damp_dT_out(n),d_damp_dw_00(n),DAMP(n), &
164 : dsrc_dr_in2(n),dsrc_dT_in(n),d_damp_dr_in2(n),d_damp_dT_in(n), &
165 : d_dampR_dr_in(n),d_dampR_dr_00(n),d_dampR_dr_out(n),d_dampR_dT_00(n), &
166 : d_dampR_dT_out(n),d_dampR_dw_00(n),DAMPR(n), &
167 : d_dampR_dr_in2(n),d_dampR_dT_in(n), &
168 : DLCXM(n),DLCX0(n),DLCXP(n),DLCY0(n), &
169 : DLCYP(n),DLCZ0(n),DLCZP(n), &
170 : DLTXM(n),DLTX0(n),DLTXP(n),DLTY0(n), &
171 : DLTYP(n),DLTZ0(n),DLTZP(n),Lt(n), &
172 : PTURB(n),dPtrb_dw_00(n),dPtrb_dr_00(n),dPtrb_dr_in(n), &
173 : QWK(n,15),QWKEV(n,15),VRU(n,15), &
174 : VRR(n,15), VRT(n,15), VRL(n,15), VRC(n,15), &
175 : DVRR(n,15),DVRT(n,15),DVRL(n,15),DVRC(n,15), &
176 : PHR(n,15), PHT(n,15), PHL(n,15), PHU(n,15), PHC(n,15), &
177 0 : ZZZ(4*n,4*n),XXL(4*n))
178 :
179 0 : NZN3=4*NZN
180 0 : LD_VL = LD_LLL
181 0 : LD_VR = LD_LLL
182 :
183 : ! SET ALL NECESSARY DERIVATIVES AND VALUES EQUAL 0
184 : ! FOR FROZEN-IN APPROXIMATION, OR IN CASE ALFA=0 (RADIATIVE)
185 0 : if(ALFA==0.d0)then
186 0 : do I=1,NZN
187 0 : Et(I)=0.d0
188 0 : EFL02=0.d0
189 0 : DLCX0(I)=0.d0
190 0 : DLCXM(I)=0.d0
191 0 : DLCXP(I)=0.d0
192 0 : DLCY0(I)=0.d0
193 0 : DLCYP(I)=0.d0
194 0 : DLCZ0(I)=0.d0
195 0 : DLCZP(I)=0.d0
196 0 : DLTX0(I)=0.d0
197 0 : DLTXM(I)=0.d0
198 0 : DLTXP(I)=0.d0
199 0 : DLTY0(I)=0.d0
200 0 : DLTYP(I)=0.d0
201 0 : DLTZ0(I)=0.d0
202 0 : DLTZP(I)=0.d0
203 0 : EVUU0(I)= 0.d0
204 0 : EVUUM(I)= 0.d0
205 0 : PTURB(I)= 0.d0
206 0 : dPtrb_dr_00(I)= 0.d0
207 0 : dPtrb_dr_in(I)= 0.d0
208 0 : dPtrb_dw_00(I)= 0.d0
209 0 : dC_dr_00(I) = 0.d0
210 0 : dC_dr_out(I) = 0.d0
211 0 : dC_dr_in(I) = 0.d0
212 0 : dC_dr_in2(I)= 0.d0
213 0 : dC_dT_in(I) = 0.d0
214 0 : dC_dT_00(I) = 0.d0
215 0 : dC_dT_out(I) = 0.d0
216 0 : dC_dw_00(I) = 0.d0
217 : end do
218 : end if
219 :
220 : ! LOOP 1 .. EOS
221 :
222 0 : !$OMP PARALLEL DO PRIVATE(I,T1,op_err) SCHEDULE(dynamic,2)
223 : do I=1,NZN
224 :
225 : if(Et(I)<=EFL02) Et(I)=EFL02
226 :
227 : call mesa_eos_kap(s,0, &
228 : T(I),Vol(I),P(I),DPV(I),dP_dT_00(I), &
229 : E(I),DEV(I),dE_dT_00(I),CPS(I),CPV(I),dCp_dT_00(I), &
230 : QQS(I),QQV(I),dQQ_dT_00(I),K(I),dK_dV_00(I),dK_dT_00(I),op_err)
231 : if (op_err /= 0) ierr = op_err
232 : if (ierr /= 0) cycle
233 :
234 : T1=P43/dm(I)
235 : DVR(I)=3.d0*T1*R(I)**2
236 : if(I==1) GOTO 2
237 : DVRM(I)=-3.d0*T1*R(max(1,I-1))**2
238 : ! bp: max(1,i-1) to prevent bogus warning from gfortran
239 : GOTO 3
240 : 2 DVRM(I)=-3.d0*T1*s% R_center**2
241 : 3 continue
242 : dP_dr_00(I) =DPV(I)*DVR(I)
243 : dP_dr_in(I)=DPV(I)*DVRM(I)
244 : dCp_dr_00(I) =CPV(I)*DVR(I)
245 : dCp_dr_in(I)=CPV(I)*DVRM(I)
246 : dQQ_dr_00(I) =QQV(I)*DVR(I)
247 : dQQ_dr_in(I)=QQV(I)*DVRM(I)
248 :
249 : end do
250 : !$OMP END PARALLEL DO
251 0 : if (ierr /= 0) return
252 :
253 : ! SKIP ALL DERIVATIVE CALCULATIONS IN CASE OF FROZEN-IN
254 : ! APPROXIMATION OR ALFA=0 (RADIATIVE CASE)
255 0 : if(ALFA==0.d0) GOTO 999
256 :
257 : ! SET E_T (NOW =w) BELOW AND ABOVE BOUNDARIES
258 0 : Et(NZN) = 0.d0
259 0 : do I=1,IBOTOM
260 0 : Et(I) = 0.d0
261 : end do
262 :
263 0 : do I=1,NZN-1
264 0 : POM= (R(I)**2)/(2.d0*G*M(I))
265 0 : POM2=POM*(P(I)*Vol(I)+P(I+1)*Vol(I+1))
266 0 : Hp_face(I)=POM2
267 0 : dHp_dr_in(I)=POM*(P(I)*DVRM(I)+Vol(I)*dP_dr_in(I))
268 : dHp_dr_00(I)=2.d0*Hp_face(I)/R(I)+POM*(P(I)*DVR(I)+Vol(I)*dP_dr_00(I) &
269 0 : +P(I+1)*DVRM(I+1)+Vol(I+1)*dP_dr_in(I+1))
270 0 : dHp_dr_out(I)=POM*(P(I+1)*DVR(I+1)+Vol(I+1)*dP_dr_00(I+1))
271 0 : dHp_dT_00(I)=POM*Vol(I)*dP_dT_00(I)
272 0 : dHp_dT_out(I)=POM*Vol(I+1)*dP_dT_00(I+1)
273 : end do
274 0 : POM=(R(NZN)**2)/(2.d0*G*M(NZN))
275 0 : Hp_face(NZN)=POM*P(NZN)*Vol(NZN)
276 0 : dHp_dr_in(NZN)=POM*(P(NZN)*DVRM(NZN)+Vol(NZN)*dP_dr_in(NZN))
277 : dHp_dr_00(NZN)=2.d0*Hp_face(NZN)/R(NZN)+POM* &
278 0 : (P(NZN)*DVR(NZN)+Vol(NZN)*dP_dr_00(NZN))
279 0 : dHp_dr_out(NZN)=0.d0
280 0 : dHp_dT_00(NZN)=POM*Vol(NZN)*dP_dT_00(NZN)
281 0 : dHp_dT_out(NZN)=0.d0
282 :
283 : !call mesa_error(__FILE__,__LINE__,'Hp_face')
284 :
285 0 : do I=1,NZN-1
286 0 : POM=0.5d0*(QQS(I)/CPS(I)+QQS(I+1)/CPS(I+1))
287 0 : POM2=0.5d0*(P(I+1)-P(I))
288 : IGR1=0.5d0*(QQS(I)/CPS(I)+QQS(I+1)/CPS(I+1))*(P(I+1)-P(I)) &
289 0 : -(dlog(T(I+1))-dlog(T(I)))
290 : IGR1X0=POM2*((dQQ_dr_00(I)-QQS(I)/CPS(I)*dCp_dr_00(I))/CPS(I)+ &
291 : (dQQ_dr_in(I+1)-QQS(I+1)/CPS(I+1)*dCp_dr_in(I+1))/CPS(I+1)) &
292 0 : +POM*(dP_dr_in(I+1)-dP_dr_00(I))
293 : IGR1XM=POM2*(dQQ_dr_in(I)-QQS(I)/CPS(I)*dCp_dr_in(I))/CPS(I) &
294 0 : +POM*(-dP_dr_in(I))
295 : IGR1XP=POM2*(dQQ_dr_00(I+1)-QQS(I+1)/CPS(I+1)*dCp_dr_00(I+1))/CPS(I+1) &
296 0 : +POM*(dP_dr_00(I+1))
297 : IGR1Y0=POM2*(dQQ_dT_00(I)-QQS(I)/CPS(I)*dCp_dT_00(I))/CPS(I) &
298 0 : +POM*(-dP_dT_00(I))+1.d0/T(I)
299 : IGR1YP=POM2*(dQQ_dT_00(I+1)-QQS(I+1)/CPS(I+1)*dCp_dT_00(I+1))/CPS(I+1) &
300 0 : +POM*(dP_dT_00(I+1))-1.d0/T(I+1)
301 :
302 0 : POM=2.d0/(Vol(I)+Vol(I+1))
303 0 : POM2=8.d0*PI*(R(I)**2)/dm_bar(I)*Hp_face(I)
304 0 : IGR2=4.d0*PI*(R(I)**2)*Hp_face(I)*POM/dm_bar(I)
305 : IGR2X0=2.d0*IGR2/R(I)+IGR2/Hp_face(I)*dHp_dr_00(I) &
306 0 : -POM2/(Vol(I)+Vol(I+1))**2*(DVR(I)+DVRM(I+1))
307 : IGR2XM=-POM2/(Vol(I)+Vol(I+1))**2*DVRM(I) &
308 0 : +IGR2/Hp_face(I)*dHp_dr_in(I)
309 : IGR2XP=-POM2/(Vol(I)+Vol(I+1))**2*DVR(I+1) &
310 0 : +IGR2/Hp_face(I)*dHp_dr_out(I)
311 0 : IGR2Y0=IGR2/Hp_face(I)*dHp_dT_00(I)
312 0 : IGR2YP=IGR2/Hp_face(I)*dHp_dT_out(I)
313 :
314 0 : Y_face(I)=IGR1*IGR2
315 0 : dY_dr_00(I)=IGR1*IGR2X0+IGR2*IGR1X0
316 0 : dY_dr_in(I)=IGR1*IGR2XM+IGR2*IGR1XM
317 0 : dY_dr_out(I)=IGR1*IGR2XP+IGR2*IGR1XP
318 0 : dY_dT_00(I)=IGR1*IGR2Y0+IGR2*IGR1Y0
319 0 : dY_dT_out(I)=IGR1*IGR2YP+IGR2*IGR1YP
320 : end do
321 :
322 : !call mesa_error(__FILE__,__LINE__,'IGRS')
323 :
324 0 : do I=1,NZN-1
325 0 : POM=sqrt(2.d0/3.d0)*0.5d0
326 : FF=POM*((E(I) +P(I) *Vol(I) )/T(I) &
327 0 : +(E(I+1)+P(I+1)*Vol(I+1))/T(I+1))
328 0 : FFY0=POM*(dE_dT_00(I)+dP_dT_00(I)*Vol(I)-(E(I)+P(I)*Vol(I))/T(I))/T(I)
329 : FFYP=POM*(dE_dT_00(I+1)+dP_dT_00(I+1)*Vol(I+1)-(E(I+1)+P(I+1)*Vol(I+1)) &
330 0 : /T(I+1))/T(I+1)
331 0 : FFXP=POM* ((DEV(I+1)+P(I+1))*DVR(I+1)+Vol(I+1)*dP_dr_00(I+1))/T(I+1)
332 : FFX0=POM*(((DEV(I) +P(I) )*DVR(I) +Vol(I) *dP_dr_00(I) )/T(I)+ &
333 : ((DEV(I+1)+P(I+1))*DVRM(I+1)+Vol(I+1)*dP_dr_in(I+1)) &
334 0 : /T(I+1))
335 0 : FFXM=POM* ((DEV(I) +P(I) )*DVRM(I) +Vol(I)*dP_dr_in(I))/T(I)
336 :
337 0 : POM=ALFAS*ALFA
338 0 : POM2=0.5d0*(CPS(I)+CPS(I+1))
339 0 : GG=POM*POM2*Y_face(I)
340 0 : GGXM=POM*(POM2*dY_dr_in(I)+Y_face(I)*0.5d0*dCp_dr_in(I))
341 0 : GGX0=POM*(POM2*dY_dr_00(I)+Y_face(I)*0.5d0*(dCp_dr_00(I)+dCp_dr_in(I+1)))
342 0 : GGXP=POM*(POM2*dY_dr_out(I)+Y_face(I)*0.5d0*dCp_dr_00(I+1))
343 0 : GGY0=POM*(POM2*dY_dT_00(I)+Y_face(I)*0.5d0*dCp_dT_00(I))
344 0 : GGYP=POM*(POM2*dY_dT_out(I)+Y_face(I)*0.5d0*dCp_dT_00(I+1))
345 0 : GPF=GG/FF
346 :
347 : ! corelation PI defined without e_t
348 0 : POM=1.d0
349 :
350 0 : PII(I)=POM*GG
351 0 : DPIIZ0(I)=0.d0
352 0 : dPII_dr_00(I)=POM*GGX0
353 0 : dPII_dr_in(I)=POM*GGXM
354 0 : dPII_dr_out(I)=POM*GGXP
355 0 : dPII_dT_00(I)=POM*GGY0
356 0 : dPII_dT_out(I)=POM*GGYP
357 : end do
358 :
359 : !call mesa_error(__FILE__,__LINE__,'LINA')
360 :
361 0 : do I=IBOTOM+1,NZN-1
362 : ! SOURCE TERM
363 0 : POM=0.5d0*(PII(I)/Hp_face(I)+PII(I-1)/Hp_face(I-1))
364 0 : POM2=T(I)*P(I)*QQS(I)/CPS(I)
365 0 : POM3=sqrt(Et(I))
366 0 : SOURCE(I)=POM*POM2*POM3
367 0 : TEM1=POM2*POM3*0.5d0
368 0 : TEMI=-PII(I)/Hp_face(I)**2
369 0 : TEMM=-PII(I-1)/Hp_face(I-1)**2
370 :
371 : dsrc_dr_out(I)= TEM1*(dPII_dr_out(I)/Hp_face(I) &
372 0 : +TEMI*dHp_dr_out(I))
373 : dsrc_dr_00(I)= TEM1*(dPII_dr_00(I)/Hp_face(I)+dPII_dr_out(I-1)/Hp_face(I-1) &
374 0 : +TEMI*dHp_dr_00(I)+TEMM*dHp_dr_out(I-1))
375 : dsrc_dr_in(I)= TEM1*(dPII_dr_in(I)/Hp_face(I)+dPII_dr_00(I-1)/Hp_face(I-1) &
376 0 : +TEMI*dHp_dr_in(I)+TEMM*dHp_dr_00(I-1))
377 : dsrc_dr_in2(I)=TEM1*(dPII_dr_in(I-1)/Hp_face(I-1) &
378 0 : +TEMM*dHp_dr_in(I-1))
379 :
380 : dsrc_dT_out(I)= TEM1*(dPII_dT_out(I)/Hp_face(I) &
381 0 : +TEMI*dHp_dT_out(I))
382 : dsrc_dT_00(I)= TEM1*(dPII_dT_00(I)/Hp_face(I)+dPII_dT_out(I-1)/Hp_face(I-1) &
383 0 : +TEMI*dHp_dT_00(I)+TEMM*dHp_dT_out(I-1))
384 : dsrc_dT_in(I)= TEM1*(dPII_dT_00(I-1)/Hp_face(I-1) &
385 0 : +TEMM*dHp_dT_00(I-1))
386 :
387 0 : dsrc_dw_00(I)= POM*POM2*0.5d0/sqrt(Et(I))
388 :
389 0 : POM=POM*POM3
390 :
391 : dsrc_dT_00(I)=dsrc_dT_00(I)+ &
392 : POM/CPS(I)*(P(I)*QQS(I)+T(I)*QQS(I)*dP_dT_00(I)+T(I)*P(I)*dQQ_dT_00(I) &
393 0 : -T(I)*P(I)*QQS(I)/CPS(I)*dCp_dT_00(I))
394 : dsrc_dr_00(I)=dsrc_dr_00(I)+ &
395 : POM*T(I)/CPS(I)*(QQS(I)*dP_dr_00(I)+P(I)*dQQ_dr_00(I) &
396 0 : -P(I)*QQS(I)/CPS(I)*dCp_dr_00(I))
397 : dsrc_dr_in(I)=dsrc_dr_in(I)+ &
398 : POM*T(I)/CPS(I)*(QQS(I)*dP_dr_in(I)+P(I)*dQQ_dr_in(I) &
399 0 : -P(I)*QQS(I)/CPS(I)*dCp_dr_in(I))
400 :
401 : ! SOURCE TERM ALWAYS POSITIVE
402 : ! THIS IS FORMULATION USED IN BUDAPEST-FLORIDA CODE
403 : ! IT REQUIRES E_T AS MAIN VARIABLE - OTHERWISE CONVERGENCE
404 : ! IS EXTREMELY SLOW, AND POSSIBLE ONLY IF LOW ACCURACY
405 : ! FOR CONVERGENCE CONDITION, DXXC< 1d-2 - 1.d-4, IS SET
406 : ! if(SOURCE(I)<=0.d0)then
407 : ! SOURCE(I) = 0.d0
408 : ! dsrc_dr_out(I) = 0.d0
409 : ! dsrc_dr_00(I) = 0.d0
410 : ! dsrc_dr_in(I) = 0.d0
411 : ! dsrc_dr_in2(I)= 0.d0
412 : ! dsrc_dT_out(I) = 0.d0
413 : ! dsrc_dT_in(I) = 0.d0
414 : ! dsrc_dT_00(I) = 0.d0
415 : ! dsrc_dw_00(I) = 0.d0
416 : ! end if
417 :
418 : ! DAMP TERM
419 0 : POM=(CEDE/ALFA)*(Et(I)**1.5d0-EFL02**1.5d0)
420 0 : POM2=0.5d0*(Hp_face(I)+Hp_face(I-1))
421 0 : DAMP(I)=POM/POM2
422 0 : TEM1=-0.5d0*POM/POM2**2
423 0 : d_damp_dr_out(I) =TEM1*dHp_dr_out(I)
424 0 : d_damp_dr_00(I) =TEM1*(dHp_dr_00(I)+dHp_dr_out(I-1))
425 0 : d_damp_dr_in(I) =TEM1*(dHp_dr_in(I)+dHp_dr_00(I-1))
426 0 : d_damp_dr_in2(I)=TEM1*(dHp_dr_in(I-1))
427 0 : d_damp_dT_00(I) =TEM1*(dHp_dT_00(I)+dHp_dT_out(I-1))
428 0 : d_damp_dT_out(I) =TEM1*dHp_dT_out(I)
429 0 : d_damp_dT_in(I) =TEM1*dHp_dT_00(I-1)
430 0 : d_damp_dw_00(I) =1.5d0*(CEDE/ALFA)/POM2*sqrt(Et(I))
431 :
432 : ! RADIATIVE DAMP TERM
433 0 : if(GAMMAR==0.d0)then
434 0 : DAMPR(I) = 0.d0
435 0 : d_dampR_dr_out(I) = 0.d0
436 0 : d_dampR_dr_00(I) = 0.d0
437 0 : d_dampR_dr_in(I) = 0.d0
438 0 : d_dampR_dr_in2(I)= 0.d0
439 0 : d_dampR_dT_out(I) = 0.d0
440 0 : d_dampR_dT_in(I) = 0.d0
441 0 : d_dampR_dT_00(I) = 0.d0
442 0 : d_dampR_dw_00(I) = 0.d0
443 : else
444 0 : POM=(GAMMAR**2)/(ALFA**2)*4.d0*SIG
445 0 : POM2=T(I)**3*Vol(I)**2/(CPS(I)*K(I))
446 0 : POM3=Et(I)
447 0 : POM4=0.5d0*(Hp_face(I)**2+Hp_face(I-1)**2)
448 0 : DAMPR(I)=POM*POM2*POM3/POM4
449 0 : TEM1=-DAMPR(I)/POM4
450 0 : d_dampR_dr_out(I) =TEM1*(Hp_face(I)*dHp_dr_out(I))
451 : d_dampR_dr_00(I) =TEM1*(Hp_face(I)*dHp_dr_00(I) &
452 0 : +Hp_face(I-1)*dHp_dr_out(I-1))
453 : d_dampR_dr_in(I) =TEM1*(Hp_face(I)*dHp_dr_in(I) &
454 0 : +Hp_face(I-1)*dHp_dr_00(I-1))
455 0 : d_dampR_dr_in2(I)=TEM1*(Hp_face(I-1)*dHp_dr_in(I-1))
456 0 : d_dampR_dT_out(I) =TEM1*(Hp_face(I)*dHp_dT_out(I))
457 : d_dampR_dT_00(I) =TEM1*(Hp_face(I)*dHp_dT_00(I) &
458 0 : +Hp_face(I-1)*dHp_dT_out(I-1))
459 0 : d_dampR_dT_in(I) =TEM1*(Hp_face(I-1)*dHp_dT_00(I-1))
460 :
461 0 : d_dampR_dw_00(I)=POM*POM2/POM4
462 :
463 0 : TEM1=POM*POM3/POM4
464 : d_dampR_dr_00(I)=d_dampR_dr_00(I) &
465 : +TEM1*T(I)**3*(2.d0*Vol(I)*DVR(I) &
466 : -Vol(I)**2*( 1.d0/CPS(I)*dCp_dr_00(I) &
467 : +1.d0/K(I)*dK_dV_00(I)*DVR(I))) &
468 0 : /(CPS(I)*K(I))
469 :
470 : d_dampR_dr_in(I)=d_dampR_dr_in(I) &
471 : +TEM1*T(I)**3*(2.d0*Vol(I)*DVRM(I) &
472 : -Vol(I)**2*( 1.d0/CPS(I)*dCp_dr_in(I) &
473 : +1.d0/K(I)*dK_dV_00(I)*DVRM(I))) &
474 0 : /(CPS(I)*K(I))
475 :
476 : d_dampR_dT_00(I)=d_dampR_dT_00(I) &
477 : +TEM1*Vol(I)**2*(3.d0*T(I)**2 &
478 : -T(I)**3*( 1.d0/CPS(I)*dCp_dT_00(I) &
479 : +1.d0/K(I)*dK_dT_00(I))) &
480 0 : /(CPS(I)*K(I))
481 :
482 : end if
483 :
484 0 : dC_dr_00(I) =dsrc_dr_00(I) -d_damp_dr_00(I) -d_dampR_dr_00(I)
485 0 : dC_dr_out(I) =dsrc_dr_out(I) -d_damp_dr_out(I) -d_dampR_dr_out(I)
486 0 : dC_dr_in(I) =dsrc_dr_in(I) -d_damp_dr_in(I) -d_dampR_dr_in(I)
487 0 : dC_dr_in2(I)=dsrc_dr_in2(I)-d_damp_dr_in2(I)-d_dampR_dr_in2(I)
488 0 : dC_dT_in(I) =dsrc_dT_in(I) -d_damp_dT_in(I) -d_dampR_dT_in(I)
489 0 : dC_dT_00(I) =dsrc_dT_00(I) -d_damp_dT_00(I) -d_dampR_dT_00(I)
490 0 : dC_dT_out(I) =dsrc_dT_out(I) -d_damp_dT_out(I) -d_dampR_dT_out(I)
491 0 : dC_dw_00(I) =dsrc_dw_00(I) -d_damp_dw_00(I) -d_dampR_dw_00(I)
492 :
493 : end do
494 :
495 : !call mesa_error(__FILE__,__LINE__,'LINA')
496 :
497 0 : do I=IBOTOM,NZN-1
498 : ! CONVECTIVE LUMINOSITY
499 : POM=4.d0*PI*(R(I)**2)*(T(I)/Vol(I)+T(I+1)/Vol(I+1))*0.5d0* &
500 0 : (ALFAC/ALFAS)
501 0 : POM3=(sqrt(Et(I))+sqrt(Et(I+1)))*0.5d0
502 0 : Lc(I)=POM*PII(I)*POM3
503 0 : DLCZ0(I)=0d0
504 0 : DLCZP(I)=0d0
505 0 : if (Et(I) /= 0d0) DLCZ0(I)=POM*PII(I)*0.25d0/sqrt(Et(I))
506 0 : if (Et(I+1) /= 0d0) DLCZP(I)=POM*PII(I)*0.25d0/sqrt(Et(I+1))
507 0 : POM2=4.d0*PI*(R(I)**2)*PII(I)*0.5d0*(ALFAC/ALFAS)*POM3
508 0 : POM=POM*POM3
509 : DLCX0(I)=POM*dPII_dr_00(I)+2.d0*Lc(I)/R(I) &
510 : -POM2*(T(I) /(Vol(I)**2 )*DVR(I)+ &
511 0 : T(I+1)/(Vol(I+1)**2)*DVRM(I+1))
512 0 : DLCXM(I)=POM*dPII_dr_in(I)-POM2*T(I) /(Vol(I)**2 )*DVRM(I)
513 0 : DLCXP(I)=POM*dPII_dr_out(I)-POM2*T(I+1)/(Vol(I+1)**2)*DVR(I+1)
514 0 : DLCY0(I)=POM*dPII_dT_00(I)+POM2/Vol(I)
515 0 : DLCYP(I)=POM*dPII_dT_out(I)+POM2/Vol(I+1)
516 :
517 0 : if(I==IBOTOM) DLCZ0(I)=0.d0
518 0 : if(I==NZN-1) DLCZP(I)=0.d0
519 0 : if(Et(I)<EFL02*1d-10)then
520 0 : Lc(I)=0.d0
521 0 : DLCX0(I)=0.d0
522 0 : DLCXM(I)=0.d0
523 0 : DLCXP(I)=0.d0
524 0 : DLCY0(I)=0.d0
525 0 : DLCYP(I)=0.d0
526 0 : DLCZ0(I)=0.d0
527 0 : DLCZP(I)=0.d0
528 : end if
529 :
530 : ! if(PII(I)<0.d0.or.ALFA==0.d0)then
531 : ! Lc(I)=0.d0
532 : ! DLCX0(I)=0.d0
533 : ! DLCXM(I)=0.d0
534 : ! DLCXP(I)=0.d0
535 : ! DLCY0(I)=0.d0
536 : ! DLCYP(I)=0.d0
537 : ! DLCZ0(I)=0.d0
538 : ! DLCZP(I)=0.d0
539 : ! end if
540 :
541 : ! TURBULENT LUMINOSITY
542 0 : if(ALFAT==0.d0.or.ALFA==0.d0)then
543 0 : Lt(I)=0.d0
544 0 : DLTX0(I)=0.d0
545 0 : DLTXM(I)=0.d0
546 0 : DLTXP(I)=0.d0
547 0 : DLTY0(I)=0.d0
548 0 : DLTYP(I)=0.d0
549 0 : DLTZ0(I)=0.d0
550 0 : DLTZP(I)=0.d0
551 : else
552 0 : POM=-2.d0/3.d0*ALFA*ALFAT*(4.d0*PI*(R(I)**2))**2
553 0 : POM2=Hp_face(I)*(1.d0/Vol(I)**2+1.d0/Vol(I+1)**2)*0.5d0
554 0 : POM3=(Et(I+1)**1.5d0-Et(I)**1.5d0)/dm_bar(I)
555 0 : Lt(I)=POM*POM2*POM3
556 : DLTX0(I)=4.d0*Lt(I)/R(I) &
557 : +Lt(I)/POM2*Hp_face(I)*(-1.d0/Vol(I)**3*DVR(I) &
558 : -1.d0/Vol(I+1)**3*DVRM(I+1)) &
559 0 : +Lt(I)/Hp_face(I)*dHp_dr_00(I)
560 :
561 : DLTXM(I)=Lt(I)/POM2*Hp_face(I)*(-1.d0/Vol(I)**3*DVRM(I)) &
562 0 : +Lt(I)/Hp_face(I)*dHp_dr_in(I)
563 : DLTXP(I)=Lt(I)/POM2*Hp_face(I)*(-1.d0/Vol(I+1)**3*DVR(I+1)) &
564 0 : +Lt(I)/Hp_face(I)*dHp_dr_out(I)
565 0 : DLTY0(I)=Lt(I)/Hp_face(I)*dHp_dT_00(I)
566 0 : DLTYP(I)=Lt(I)/Hp_face(I)*dHp_dT_out(I)
567 0 : DLTZ0(I)=-POM*POM2*1.5d0*sqrt(Et(I) )/dm_bar(I)
568 0 : DLTZP(I)= POM*POM2*1.5d0*sqrt(Et(I+1))/dm_bar(I)
569 : end if
570 : end do
571 :
572 : ! TURBULENT PRESSURE (ZONE)
573 0 : do I=IBOTOM+1,NZN-1
574 0 : if(ALFAP==0.d0)then
575 0 : PTURB(I) = 0.d0
576 0 : dPtrb_dw_00(I) = 0.d0
577 0 : dPtrb_dr_00(I) = 0.d0
578 0 : dPtrb_dr_in(I) = 0.d0
579 : else
580 0 : PTURB(I) = ALFAP*Et(I)/Vol(I)
581 0 : dPtrb_dw_00(I) = ALFAP/Vol(I)
582 0 : TEM1=-ALFAP*Et(I)/Vol(I)**2
583 0 : dPtrb_dr_00(I) = TEM1*DVR(I)
584 0 : dPtrb_dr_in(I) = TEM1*DVRM(I)
585 : end if
586 : end do
587 :
588 : ! EDDY VISCOSITY
589 0 : do I=IBOTOM+1,NZN-1
590 0 : if(ALFAM>=0d0) then
591 : ! Kuhfuss (1986) tensor EDDY VISCOSITY
592 0 : POM=16.d0/3.d0*PI*ALFA*ALFAM*sqrt(Et(I))
593 0 : POM1=1.d0/Vol(I)**2/dm(I)
594 0 : POM2=(R(I)**6+R(I-1)**6)*(Hp_face(I)+Hp_face(I-1))*0.25d0
595 0 : EVUU0(I)= POM*POM1*POM2/R(I)
596 0 : EVUUM(I)=-POM*POM1*POM2/R(I-1)
597 : else
598 : ! Kollath et al. 2002 EDDY VISCOSITY pressure
599 0 : POM=-(16.d0/3.d0)*PI*ALFA*abs(ALFAM)*sqrt(Et(I))
600 0 : POM1=1.d0/Vol(I)**2/dm(I)
601 0 : POM2=(R(I)**3+R(I-1)**3)*(Hp_face(I)+Hp_face(I-1))*0.25d0
602 0 : EVUU0(I)= POM*POM1*POM2/R(I)
603 0 : EVUUM(I)=-POM*POM1*POM2/R(I-1)
604 : end if
605 : end do
606 :
607 0 : do I=1,IBOTOM
608 0 : Lc(I)= 0.d0
609 0 : Lt(I)= 0.d0
610 0 : PTURB(I)= 0.d0
611 0 : dPtrb_dw_00(I)= 0.d0
612 0 : dPtrb_dr_00(I)= 0.d0
613 0 : dPtrb_dr_in(I)= 0.d0
614 0 : dC_dr_00(I) = 0.d0
615 0 : dC_dr_out(I) = 0.d0
616 0 : dC_dr_in(I) = 0.d0
617 0 : dC_dr_in2(I)= 0.d0
618 0 : dC_dT_in(I) = 0.d0
619 0 : dC_dT_00(I) = 0.d0
620 0 : dC_dT_out(I) = 0.d0
621 0 : dC_dw_00(I) = 0.d0
622 0 : EVUU0(I) = 0.d0
623 0 : EVUUM(I) = 0.d0
624 : end do
625 0 : do I=1,IBOTOM-1
626 0 : DLCX0(I)=0.d0
627 0 : DLCXM(I)=0.d0
628 0 : DLCXP(I)=0.d0
629 0 : DLCY0(I)=0.d0
630 0 : DLCYP(I)=0.d0
631 0 : DLCZ0(I)=0.d0
632 0 : DLCZP(I)=0.d0
633 0 : DLTX0(I)=0.d0
634 0 : DLTXM(I)=0.d0
635 0 : DLTXP(I)=0.d0
636 0 : DLTY0(I)=0.d0
637 0 : DLTYP(I)=0.d0
638 0 : DLTZ0(I)=0.d0
639 0 : DLTZP(I)=0.d0
640 : end do
641 0 : Lc(NZN)= 0.d0
642 0 : Lt(NZN)= 0.d0
643 0 : PTURB(NZN)= 0.d0
644 0 : dPtrb_dw_00(NZN)= 0.d0
645 0 : dPtrb_dr_00(NZN)= 0.d0
646 0 : dPtrb_dr_in(NZN)= 0.d0
647 0 : DLCX0(NZN)= 0.d0
648 0 : DLCXM(NZN)= 0.d0
649 0 : DLCXP(NZN)= 0.d0
650 0 : DLCY0(NZN)= 0.d0
651 0 : DLCYP(NZN)= 0.d0
652 0 : DLCZ0(NZN)= 0.d0
653 0 : DLCZP(NZN)= 0.d0
654 0 : DLTX0(NZN)=0.d0
655 0 : DLTXM(NZN)=0.d0
656 0 : DLTXP(NZN)=0.d0
657 0 : DLTY0(NZN)=0.d0
658 0 : DLTYP(NZN)=0.d0
659 0 : DLTZ0(NZN)=0.d0
660 0 : DLTZP(NZN)=0.d0
661 0 : dC_dr_00(NZN) = 0.d0
662 0 : dC_dr_out(NZN) = 0.d0
663 0 : dC_dr_in(NZN) = 0.d0
664 0 : dC_dr_in2(NZN)= 0.d0
665 0 : dC_dT_in(NZN) = 0.d0
666 0 : dC_dT_00(NZN) = 0.d0
667 0 : dC_dT_out(NZN) = 0.d0
668 0 : dC_dw_00(NZN) = 0.d0
669 0 : EVUU0(NZN)= 0.d0
670 0 : EVUUM(NZN)= 0.d0
671 :
672 0 : DLCZP(NZN-1)=0.d0
673 0 : DLTZP(NZN-1)=0.d0
674 :
675 : 999 continue
676 :
677 : ! LOOP 2 .. LUM PLUGS
678 : DLR = 0.d0
679 : DLRP = 0.d0
680 : DLRM = 0.d0
681 : DLT = 0.d0
682 : DLTP = 0.d0
683 : DLR = 0.d0 !-1.d0
684 :
685 : DFCX0 = 0.d0
686 : DFCXM = 0.d0
687 : DFCXP = 0.d0
688 : DFCY0 = 0.d0
689 : DFCYP = 0.d0
690 : DFCZ0 = 0.d0
691 : DFCZP = 0.d0
692 0 : do I=1,NZN
693 : ! SET LUM(I-1)
694 0 : DLMR = DLR
695 0 : DLMRP = DLRP
696 0 : DLMRM = DLRM
697 0 : DLMT = DLT
698 0 : DLMTP = DLTP
699 :
700 0 : DFCMX0 = DFCX0
701 0 : DFCMXM = DFCXM
702 0 : DFCMXP = DFCXP
703 0 : DFCMY0 = DFCY0
704 0 : DFCMYP = DFCYP
705 0 : DFCMZ0 = DFCZ0
706 0 : DFCMZP = DFCZP
707 0 : if(I==NZN) GOTO 6
708 : ! Lr(I)=Eq. A.4, Stellingwerf 1975, Appendix A
709 : ! CALC LUM(I)
710 0 : W_00=T(I)**4
711 0 : W_out=T(I+1)**4
712 0 : BW=dlog(W_out/W_00)
713 0 : BK=dlog(K(I+1)/K(I))
714 0 : T1=-CL*R(I)**4/dm_bar(I)
715 0 : T2=(W_out/K(I+1)-W_00/K(I))/(1.d0-BK/BW)
716 0 : T3=T1/(BW-BK)
717 0 : DLK= (T3/K(I)) *(W_00*BW/K(I) -T2) !dL(i)/dK(i)
718 0 : DLKP=-(T3/K(I+1))*(W_out*BW/K(I+1)-T2) !dL(i)/dK(i+1)
719 0 : DLRP= DLKP*dK_dV_00(I+1)*DVR(I+1)
720 0 : DLRM= DLK *dK_dV_00(I) *DVRM(I)
721 0 : DLR= 4.d0*T1*T2/R(I)+DLK*dK_dV_00(I)*DVR(I)+DLKP*dK_dV_00(I+1)*DVRM(I+1)
722 0 : DLTP=4.d0*(T3/T(I+1))*(W_out*BW/K(I+1)-T2*BK/BW)+DLKP*dK_dT_00(I+1)
723 0 : DLT=-4.d0*(T3/T(I))*(W_00*BW/K(I)-T2*BK/BW)+DLK*dK_dT_00(I)
724 : !-----------------------------
725 0 : DFCX0=DLCX0(I)
726 0 : DFCXM=DLCXM(I)
727 0 : DFCXP=DLCXP(I)
728 0 : DFCY0=DLCY0(I)
729 0 : DFCYP=DLCYP(I)
730 0 : DFCZ0=DLCZ0(I)
731 0 : DFCZP=DLCZP(I)
732 0 : GOTO 7
733 : ! OUTER LUM BOUNDARY CONDITION
734 : 6 continue
735 0 : DLT = 4.d0*L0/T(I) !L=4piR^2sigT^4
736 0 : DLR = 2.d0*L0/R(I) !L=4piR^2sigT^4
737 0 : DLRM = 0.d0
738 0 : DLRP = 0.d0
739 0 : DLTP = 0.d0
740 0 : DFCX0 = 0.d0
741 0 : DFCXM = 0.d0
742 0 : DFCXP = 0.d0
743 0 : DFCY0 = 0.d0
744 0 : DFCYP = 0.d0
745 0 : DFCZ0 = 0.d0
746 0 : DFCZP = 0.d0
747 : 7 continue
748 0 : ELRP(I) = DLRP +DFCXP +DLTXP(I)
749 0 : ELRM(I) = DLRM +DFCXM +DLTXM(I)
750 0 : ELR(I) = DLR +DFCX0 +DLTX0(I)
751 0 : ELTP(I) = DLTP +DFCYP +DLTYP(I)
752 0 : ELT(I) = DLT +DFCY0 +DLTY0(I)
753 0 : ELZ0(I) = DFCZ0 +DLTZ0(I)
754 0 : ELZP(I) = DFCZP +DLTZP(I)
755 :
756 : ! CALC ENERGY EQUATION(I)
757 0 : T2=P4/dm(I)
758 0 : T3=T2*(P(I)+DEV(I))
759 0 : T4=1.d0/dm(I)
760 0 : EX20(I) = -T4*( -DLMRM) -dC_dr_in2(I) -T4*( -DFCMXM)
761 0 : EX10(I) = -T4*(DLRM-DLMR ) -dC_dr_in(I) -T4*(DFCXM-DFCMX0)
762 0 : EX00(I) = -T4*(DLR -DLMRP) -dC_dr_00(I) -T4*(DFCX0-DFCMXP)
763 0 : EX01(I) = -T4*(DLRP ) -dC_dr_out(I) -T4*(DFCXP )
764 0 : EY10(I) = -T4*( -DLMT ) -dC_dT_in(I) -T4*( -DFCMY0)
765 0 : EY00(I) = -T4*(DLT -DLMTP) -dC_dT_00(I) -T4*(DFCY0-DFCMYP)
766 0 : EY01(I) = -T4*(DLTP ) -dC_dT_out(I) -T4*(DFCYP )
767 0 : EZ10(I) = -T4*( -DFCMZ0)
768 0 : EZ00(I) = -dC_dw_00(I) -T4*(DFCZ0-DFCMZP)
769 0 : EZ01(I) = -T4*(DFCZP )
770 0 : if(I==1) then
771 0 : EU10(I) = T3*s% R_center**2
772 : else
773 0 : EU10(I) = T3*R(max(1,I-1))**2
774 : ! bp: max(1,I-1) to prevent bogus warning from gfortran
775 : end if
776 0 : EU00(I) = -T3*R(I)**2
777 :
778 : ! CALC CONVECTIVE ENERGY EQUATION(I)
779 0 : T3=P4/dm(I)*PTURB(I)
780 0 : T4=1.d0/dm(I)
781 0 : if (I == 1) then
782 0 : CX20(I) = dC_dr_in2(I) -T4*( -0)
783 0 : CX10(I) = dC_dr_in(I) -T4*(DLTXM(I)-0)
784 0 : CX00(I) = dC_dr_00(I) -T4*(DLTX0(I)-0)
785 0 : CY10(I) = dC_dT_in(I) -T4*( -0)
786 0 : CY00(I) = dC_dT_00(I) -T4*(DLTY0(I)-0)
787 0 : CZ10(I) = -T4*( -0)
788 0 : CZ00(I) = dC_dw_00(I) -T4*(DLTZ0(I)-0)
789 : else
790 0 : CX20(I) = dC_dr_in2(I) -T4*( -DLTXM(I-1))
791 0 : CX10(I) = dC_dr_in(I) -T4*(DLTXM(I)-DLTX0(I-1))
792 0 : CX00(I) = dC_dr_00(I) -T4*(DLTX0(I)-DLTXP(I-1))
793 0 : CY10(I) = dC_dT_in(I) -T4*( -DLTY0(I-1))
794 0 : CY00(I) = dC_dT_00(I) -T4*(DLTY0(I)-DLTYP(I-1))
795 0 : CZ10(I) = -T4*( -DLTZ0(I-1))
796 0 : CZ00(I) = dC_dw_00(I) -T4*(DLTZ0(I)-DLTZP(I-1))
797 : end if
798 0 : CY01(I) = dC_dT_out(I) -T4*(DLTYP(I) )
799 0 : CX01(I) = dC_dr_out(I) -T4*(DLTXP(I) )
800 0 : CZ01(I) = -T4*(DLTZP(I) )
801 0 : if(I==1) then
802 0 : CU10(I) = T3*s% R_center**2
803 : else
804 0 : CU10(I) = T3*R(max(1,I-1))**2
805 : ! bp: max(1,I-1) to prevent bogus warning from gfortran
806 : end if
807 0 : CU00(I) = -T3*R(I)**2
808 :
809 : ! CALC MOMENTUM EQUATION(I)
810 0 : T1=P4*R(I)**2/dm_bar(I)
811 0 : if(ALFAM>=0.d0) T4=P4/(dm_bar(I)*R(I))
812 0 : if(ALFAM<0.d0) T4=-T1
813 0 : MU10(I) = T4*(-EVUUM(I))
814 0 : MZ00(I) = -T1*(-dPtrb_dw_00(I))
815 0 : MX10(I) = -T1*(-dP_dr_in(I)-dPtrb_dr_in(I))
816 0 : MY00(I) = -T1*(-dP_dT_00(I))
817 0 : if(I/=NZN)then
818 : MX00(I) = 4.d0*G*M(I)/R(I)**3 &
819 0 : -T1*(dP_dr_in(I+1)-dP_dr_00(I)+dPtrb_dr_in(I+1)-dPtrb_dr_00(I))
820 0 : MX01(I) = -T1*(dP_dr_00(I+1) +dPtrb_dr_00(I+1))
821 0 : MY01(I) = -T1*(dP_dT_00(I+1))
822 0 : MU00(I) = T4*(EVUUM(I+1)-EVUU0(I))
823 0 : MU01(I) = T4*(EVUU0(I+1))
824 0 : MZ01(I) = -T1*(dPtrb_dw_00(I+1))
825 : else
826 : MX00(I) = 4.d0*G*M(I)/R(I)**3 &
827 0 : -T1*(-dP_dr_00(I))
828 0 : MX01(I) = 0.d0
829 0 : MY01(I) = 0.d0
830 0 : MU00(I) = T4*(-EVUU0(I))
831 0 : MU01(I) = 0.d0
832 0 : MZ01(I) = 0.d0
833 : end if
834 :
835 : end do
836 :
837 0 : do I=1,NZN3
838 0 : do J=1,NZN3
839 0 : LLL(I,J)=0.d0
840 : end do
841 : end do
842 :
843 : !
844 : ! VELOCITY DEFINITION
845 : ! (dR/dT) = u
846 : ! MOMENTUM EQUATION
847 : ! (dU/dt) = - (4 PI R^2)(dp/dm) - (GM)/R^2
848 : ! ENERGY EQUATION
849 : ! (c_v)(dT/dt) = - (p+(de/dV)_T)(dV/dR)U - (dLr/dm)
850 : ! TURBULENT ENERGY EQUATION
851 : ! (de_t/dt) =
852 : ! SEE CODE DOCUMENTATION FOR LINEARIZATION
853 : !
854 : ! EIGENVALUE PROBLEM: LLL*X=SIGMA*X
855 : ! X={dR1,dU1,dT1,dw1,...,dRN,dUN,dTN,dwN}
856 : ! VELOCITY DEF. EQ. IN ROWS IG
857 : ! MOMENTUM EQ. IN ROWS IR
858 : ! ENERGY EQ. IN ROWS IE
859 : ! TURBULENT ENERGY EQ. IN ROWS IC
860 : ! NAMING SCHEME FOR ELEMENTS OF LLL:
861 : ! EX00(I) - d(energy==)/dR_i
862 : ! EX01(I) - d(energy==)/dR_i+1
863 : ! EX10(I) - d(energy==)/dR_i-1
864 : ! EY00(I) - d(energy==)/dT_i
865 : ! EU00(I) - d(energy==)/dU_i
866 : ! MY00(I) - d(moment==)/dT_i
867 : ! e.t.c.
868 : ! VELOCITY DEF. AND MOMENTUM EQ. AT THE INTERFACE 1....NZN
869 : ! ENERGY EQUATIONS IN THE ZONE 1....NZN
870 : !
871 :
872 0 : do I=1,NZN
873 0 : IG=4*I-3
874 0 : IR=4*I-2
875 0 : IE=4*I-1
876 0 : IC=4*I
877 :
878 0 : if(IG+1<=NZN3) LLL(IG,IG+1) = 1.d0
879 : !---------------------------------------
880 0 : if(IR-4>=1) LLL(IR,IR-4) = MU10(I)
881 0 : LLL(IR,IR) = MU00(I)
882 0 : if(IR+4<=NZN3) LLL(IR,IR+4) = MU01(I)
883 :
884 0 : if(IR-5>=1) LLL(IR,IR-5) = MX10(I)
885 0 : if(IR-1>=1) LLL(IR,IR-1) = MX00(I)
886 0 : if(IR+3<=NZN3) LLL(IR,IR+3) = MX01(I)
887 :
888 0 : if(IR+1<=NZN3) LLL(IR,IR+1) = MY00(I)
889 0 : if(IR+5<=NZN3) LLL(IR,IR+5) = MY01(I)
890 :
891 0 : if(IR+2<=NZN3) LLL(IR,IR+2) = MZ00(I)
892 0 : if(IR+6<=NZN3) LLL(IR,IR+6) = MZ01(I)
893 : !---------------------------------------
894 0 : if(IE-5>=1) LLL(IE,IE-5) = EU10(I)/dE_dT_00(I)
895 0 : if(IE-1>=1) LLL(IE,IE-1) = EU00(I)/dE_dT_00(I)
896 :
897 0 : if(IE-10>=1) LLL(IE,IE-10) = EX20(I)/dE_dT_00(I)
898 0 : if(IE-6>=1) LLL(IE,IE-6) = EX10(I)/dE_dT_00(I)
899 0 : if(IE-2>=1) LLL(IE,IE-2) = EX00(I)/dE_dT_00(I)
900 0 : if(IE+2<=NZN3) LLL(IE,IE+2) = EX01(I)/dE_dT_00(I)
901 :
902 0 : if(IE-4>=1) LLL(IE,IE-4) = EY10(I)/dE_dT_00(I)
903 0 : LLL(IE,IE) = EY00(I)/dE_dT_00(I)
904 0 : if(IE+4<=NZN3) LLL(IE,IE+4) = EY01(I)/dE_dT_00(I)
905 :
906 0 : if(IE-3>=1) LLL(IE,IE-3) = EZ10(I)/dE_dT_00(I)
907 0 : if(IE+1<=NZN3) LLL(IE,IE+1) = EZ00(I)/dE_dT_00(I)
908 0 : if(IE+5<=NZN3) LLL(IE,IE+5) = EZ01(I)/dE_dT_00(I)
909 : !---------------------------------------
910 0 : if(IC-11>=1) LLL(IC,IC-11) = CX20(I)
911 0 : if(IC-7>=1) LLL(IC,IC-7) = CX10(I)
912 0 : if(IC-3>=1) LLL(IC,IC-3) = CX00(I)
913 0 : if(IC+1<=NZN3) LLL(IC,IC+1) = CX01(I)
914 :
915 0 : if(IC-6>=1) LLL(IC,IC-6) = CU10(I)
916 0 : if(IC-2>=1) LLL(IC,IC-2) = CU00(I)
917 :
918 0 : if(IC-5>=1) LLL(IC,IC-5) = CY10(I)
919 0 : if(IC-1>=1) LLL(IC,IC-1) = CY00(I)
920 0 : if(IC+3<=NZN3) LLL(IC,IC+3) = CY01(I)
921 :
922 0 : if(IC-4>=1) LLL(IC,IC-4) = CZ10(I)
923 0 : LLL(IC,IC) = CZ00(I)
924 0 : if(IC+4<=NZN3) LLL(IC,IC+4) = CZ01(I)
925 :
926 0 : IF (IE+4 <= NZN3) then
927 : if(LLL(IE,IE+4)<0.d0)then
928 : !write(*,*) 'rerrrrrrrrrrrrrrrrrrrrrrrr',i
929 : end if
930 : end if
931 : end do
932 :
933 0 : if (s% RSP_trace_RSP_build_model) &
934 0 : write(*,*) 'waiting for DGEEV to solve eigenvalue problem....'
935 : call DGEEV('n','v',NZN3,LLL,LD_LLL,WRx,WIx,VLx,LD_VL,VRx,LD_VR, &
936 0 : WORKx,4*NZN3,INFO)
937 0 : if(INFO/=0)then
938 0 : write(*,*) 'FAILED!'
939 0 : write(*,*) 'LAPACK/DGEEV error, ierr= ',INFO
940 0 : ierr = -1
941 0 : return
942 : stop
943 : end if
944 :
945 : ! THIS IS A MODIFIED "NUMERICAL RECIPES" SORTING SUBROUTINE
946 : ! IT SORTS THE VECTOR WI (DIMENSION NZN3) IN ASCENDING ORDER
947 : ! AND IN THE SAME WAY REARRANGES THE ELEMENTS OF WR (WR IS NOT SORTED)
948 : ! THE INTEGER VECTOR ISORT CONTAINS INFORMATION ABOUT CHANGES DONE
949 : ! NAMELY ELEMENT I OF SORTED WI (AND REARRANGED WR) WAS ISORT(I) ELEMENT
950 : ! OF UNSORTED WI (AND NOT REARANGED WR)
951 0 : call SORT(NZN3,WIx,WRx,ISORTx)
952 :
953 : ! THIS LOOP FINDS THE SMALLEST BUT POSITIVE VALUE OF VECTOR WI (MINI)
954 : ! AND RETURNS THE CORRESPONDING INDEX (IMI)
955 0 : FILENAME=trim(s% log_directory) // '/' // 'LINA_period_growth_info.data'
956 : !open(15,file=trim(FILENAME),status='unknown')
957 0 : MINI=5.d0
958 0 : IMI=0
959 : !write(15,'(a)') ' I R PERIOD'
960 0 : do I=1,NZN3
961 : ! if(WIx(I)<MINI.and.WIx(I)>1d-9)then
962 0 : if((WIx(I)<MINI) .and.(WIx(I)>1.d-9))then
963 0 : if(P4*WRx(I)/WIx(I)>-.3d+1)then
964 0 : MINI=WIx(I)
965 0 : IMI=I
966 : end if
967 : end if
968 0 : if (abs(WIx(I)) > 1d-50) then
969 : !write(15,'(2(d14.8,tr2),f11.6)') &
970 : ! WIx(I),WRx(I),2.d0*PI/WIx(I)/86400d0
971 : else
972 : !write(15,'(2(d14.8,tr2),f11.6)') &
973 : ! 0d0,WRx(I),0d0
974 : end if
975 : end do
976 : !close(15)
977 :
978 0 : do J=1,NMODES
979 : 444 continue
980 0 : if(P4*WRx(IMI+J-1)/WIx(IMI+J-1)<-.5d+1)then
981 0 : IMI=IMI+1
982 0 : GOTO 444
983 : end if
984 : ! VRRS(J) IS THE MODULI OF THE SUTFACE R-EIGENVECTOR OF THE MODE J
985 : VRRS(J)=sqrt(VRx(4*NZN-3,ISORTx(IMI+J-1))**2+ &
986 0 : VRx(4*NZN-3,ISORTx(IMI+J-1)+1)**2) !surface value
987 : VTTS(J)=VRx(4*NZN-3,ISORTx(IMI+J-1))+(0.d0,1.d0)* &
988 0 : VRx(4*NZN-3,ISORTx(IMI+J-1)+1)
989 :
990 0 : SCALE(J)=R(NZN)/VTTS(J)
991 : ! PERS(J) IS THE PERIOD OF THE MODE J (IN SECONDS)
992 : OMEG(J)=WIx(IMI+J-1)
993 0 : PERS(J)=2.d0*PI/OMEG(J)
994 0 : Q(J)=PERS(J)*SQRT((M(NZN)/SUNM)*(SUNR/R(NZN))**3)/86400.d0
995 : ! ETO(J) IS THE GROWTH RATE OF MODE J
996 0 : ETO(J)= P4*WRx(IMI+J-1)/OMEG(J)
997 0 : EK(J)=0.d0
998 0 : SGRP=0.d0
999 0 : SGRM=0.d0
1000 0 : NORMC=-1.d50
1001 0 : do I=1,NZN
1002 : ! SEE LAPACK USERS GUIDE FOR CONSTRUCTION OF EIGENVECTORS
1003 :
1004 : ! VRR(I,J) IS THE dR EIGENVECTOR OF THE MODE J
1005 : ! dR_i
1006 : VRR(I,J)=VRx(4*I-3,ISORTx(IMI+J-1))+(0.d0,1.d0)* &
1007 0 : VRx(4*I-3,ISORTx(IMI+J-1)+1)
1008 : ! DVRR(I,J) IS SCALED dR/R EIGENVECTOR OF THE MODE (J)
1009 : ! (dR_i/R_i)/(dR_NZN/R_NZN)
1010 0 : DVRR(I,J)=VRR(I,J)/R(I)*SCALE(J)
1011 0 : PHR(I,J)=datan2(aimag(DVRR(I,J)),dble(DVRR(I,J)))
1012 :
1013 : ! VRT(I,J) IS THE dT EIGENVECTOR OF THE MODE J
1014 : ! dT_i
1015 : VRT(I,J)=VRx(4*I-1,ISORTx(IMI+J-1))+(0.d0,1.d0)* &
1016 0 : VRx(4*I-1,ISORTx(IMI+J-1)+1)
1017 : ! DVRT(I,J) IS SCALED dT/T EIGENVECTOR OF THE MODE (J)
1018 : ! (dT_i/T_i)/(dR_NZN/R_NZN)
1019 0 : DVRT(I,J)=VRT(I,J)/T(I)*SCALE(J)
1020 0 : PHT(I,J)=datan2(aimag(DVRT(I,J)),dble(DVRT(I,J)))
1021 :
1022 : ! VRC(I,J) IS THE dT EIGENVECTOR OF THE MODE J
1023 : ! dOMEGA_i
1024 : VRC(I,J)=VRx(4*I,ISORTx(IMI+J-1))+(0.d0,1.d0)* &
1025 0 : VRx(4*I,ISORTx(IMI+J-1)+1)
1026 : ! DVRC(I,J) IS SCALED dOMEGA/OMEGA EIGENVECTOR OF THE MODE (J)
1027 : ! (dOMEGA_i/OMEGA_i)/(dR_NZN/R_NZN)
1028 0 : DVRC(I,J)=VRC(I,J)
1029 0 : if(NORMC<abs(DVRC(I,J))) NORMC=abs(DVRC(I,J))
1030 0 : PHC(I,J)=datan2(aimag(DVRC(I,J)),dble(DVRC(I,J)))
1031 :
1032 : ! VRU(I,J) IS THE dU EIGENVECTOR OF THE MODE J
1033 : ! dU_i
1034 : VRU(I,J)=VRx(4*I-2,ISORTx(IMI+J-1))+(0.d0,1.d0)* &
1035 0 : VRx(4*I-2,ISORTx(IMI+J-1)+1)
1036 0 : PHU(I,J)=datan2(aimag(VRU(I,J)),dble(VRU(I,J)))
1037 :
1038 : ! KINETIC ENERGY OF THE MODE
1039 0 : EK(J)=EK(J)+0.5d0*dm_bar(I)*(OMEG(J)*abs(VRR(I,J)*SCALE(J)))**2
1040 :
1041 : ! WORK DONE IN ZONE I (FOR THE MODE J)
1042 0 : if(I==1) then
1043 0 : DV_0=DVR(I)*VRR(I,J)*SCALE(J)
1044 : else
1045 0 : DV_0=(DVR(I)*VRR(I,J)+DVRM(I)*VRR(I-1,J))*SCALE(J)
1046 : end if
1047 0 : DP_0=(DPV(I)*DV_0+dP_dT_00(I)*VRT(I,J)*SCALE(J))
1048 :
1049 0 : DPEV=EVUU0(I)*VRU(I,J)
1050 0 : if(I>=2)DPEV=DPEV+EVUUM(I)*VRU(I-1,J)
1051 0 : DPEV=DPEV*SCALE(J)
1052 :
1053 0 : dP_dT_00URB=dPtrb_dr_00(I)*VRR(I,J)
1054 0 : if(I>=2) dP_dT_00URB=dP_dT_00URB+dPtrb_dr_in(I)*VRR(I-1,J)
1055 0 : dP_dT_00URB=dP_dT_00URB+dPtrb_dw_00(I)*VRC(I,J)
1056 0 : dP_dT_00URB=dP_dT_00URB*SCALE(J)
1057 :
1058 0 : QWK(I,J)=-PI*dm(I)*aimag(conjg(DP_0)*DV_0)
1059 0 : if(ALFAM<0.d0)QWKEV(I,J)=-PI*dm(I)*aimag(conjg(DPEV)*DV_0)
1060 0 : QWKPT(I,J)=-PI*dm(I)*aimag(conjg(dP_dT_00URB)*DV_0)
1061 0 : if(ALFAM>0.d0)then
1062 : QWKEV(I,J)=PI*dm(I)*aimag(conjg(DPEV)*(DV_0/R(I)**3- &
1063 0 : 3.d0*Vol(I)/R(I)**4*VRR(I,J)*SCALE(J)))
1064 : end if
1065 :
1066 0 : if(QWK(I,J)+QWKEV(I,J)+QWKPT(I,J)>=0.d0) &
1067 : SGRP=SGRP+QWK(I,J)+QWKEV(I,J)+QWKPT(I,J)
1068 : if(QWK(I,J)+QWKEV(I,J)+QWKPT(I,J)<0.d0) &
1069 : SGRM=SGRM+abs(QWK(I,J)+QWKEV(I,J)+QWKPT(I,J))
1070 :
1071 0 : VEL(I,J)=abs(VRR(I,J))/VRRS(J)
1072 0 : if(abs(PHR(I,J))>1.57d0)VEL(I,J)=-VEL(I,J)
1073 : end do
1074 :
1075 : ! WRITE WORK-INTEGRALS INTO FILE
1076 0 : if(J<=9)then
1077 0 : write(NUMER1,'(I1)') J
1078 0 : FILENAME=trim(s% log_directory) // '/' // 'LINA_work'//NUMER1//'.data'
1079 : end if
1080 0 : if(J>=10)then
1081 0 : write(NUMER2,'(I2)') J
1082 0 : FILENAME=trim(s% log_directory) // '/' // 'LINA_work'//NUMER2//'.data'
1083 : end if
1084 :
1085 0 : open(57,file=trim(FILENAME),status='unknown')
1086 : write(57,'(a)') &
1087 : '#ZONE log(T) X WORK(P)' // &
1088 : ' WORK(P_NU) WORK(P_T) CWORK(P)' // &
1089 0 : ' CWORK(P_NU) CWORK(P_T)'
1090 0 : ETOI(J)=0.d0
1091 0 : ETOIEV(J)=0.d0
1092 0 : ETOIPT(J)=0.d0
1093 0 : do I=1,NZN
1094 0 : ETOI(J) =ETOI(J) +QWK(I,J)
1095 0 : ETOIEV(J)=ETOIEV(J)+QWKEV(I,J)
1096 0 : ETOIPT(J)=ETOIPT(J)+QWKPT(I,J)
1097 : write(57,'(I3,2(tr1,e15.8),6(tr1,e15.8))') &
1098 0 : NZN+1-I,dlog10(T(I)),R(I)/R(NZN), &
1099 0 : QWK(I,J)/EK(J), QWKEV(I,J)/EK(J), QWKPT(I,J)/EK(J), &
1100 0 : ETOI(J)/EK(J),ETOIEV(J)/EK(J),ETOIPT(J)/EK(J)
1101 : end do
1102 0 : write(57,'("#KINETIC ENERGY:",e15.8)') EK(J)
1103 0 : close(57)
1104 :
1105 : ! CALCULATE LUMINOSITY EIGENFUNCTIONS
1106 0 : do I=1,NZN
1107 : ! VRL(I,J) IS THE dL EIGENVECTOR OF THE MODE J
1108 : ! dL_i
1109 : VRL(I,J)= ELT(I)*VRT(I,J) &
1110 : +ELR(I)*VRR(I,J) &
1111 0 : +ELZ0(I)*VRC(I,J)
1112 0 : if(I/=1) VRL(I,J)=VRL(I,J)+ELRM(I)*VRR(I-1,J)
1113 0 : if(I/=NZN) VRL(I,J)=VRL(I,J)+ELTP(I)*VRT(I+1,J) &
1114 : +ELRP(I)*VRR(I+1,J) &
1115 0 : +ELZP(I)*VRC(I+1,J)
1116 : ! DVRL(I,J) IS SCALED dL/L EIGENVECTOR OF THE MODE (J)
1117 : ! (dL_i/L0)/(dR_NZN/R_NZN)
1118 0 : DVRL(I,J)=VRL(I,J)/L0*SCALE(J)
1119 0 : PHL(I,J)=datan2(aimag(DVRL(I,J)),dble(DVRL(I,J)))
1120 : end do
1121 :
1122 : ! WRITE EIGENVECTORS TO FILE
1123 0 : if(J<=9)then
1124 0 : write(NUMER1,'(I1)') J
1125 0 : FILENAME=trim(s% log_directory) // '/' // 'LINA_eigen'//NUMER1//'.data'
1126 : end if
1127 0 : if(J>=10)then
1128 0 : write(NUMER2,'(I2)') J
1129 0 : FILENAME=trim(s% log_directory) // '/' // 'LINA_eigen'//NUMER2//'.data'
1130 : end if
1131 0 : open(56,file=trim(FILENAME),status='unknown')
1132 : write(56,'(a)') &
1133 : '#ZONE TEMP. FRAC. RADIUS ABS(dR/R) ' // &
1134 : 'PH(dR/R) ABS(dT/T) PH(dT/T) ' // &
1135 0 : 'ABS(dL/L) PH(dL/L) ABS(dE_T) PH(dE_T)'
1136 0 : do I=1,NZN
1137 : write(56,'(I3,tr1,e15.6,tr1,e15.6,8(tr1,e15.8))') &
1138 0 : NZN+1-i,T(I),R(I)/R(NZN),abs(DVRR(I,J)),PHR(I,J), &
1139 0 : abs(DVRT(I,J)),PHT(I,J), &
1140 0 : abs(DVRL(I,J)),PHL(I,J), &
1141 0 : abs(DVRC(I,J))/NORMC,PHC(I,J)
1142 : end do
1143 0 : close(56)
1144 :
1145 0 : SGR(J)=(SGRP-SGRM)/(SGRP+SGRM)
1146 :
1147 0 : PSIG=M(NZN)/(P43*R(NZN)**3)
1148 : PSIG=sqrt(P4*G*PSIG)
1149 :
1150 : QCHECK(J)=100.d0*(ETO(J)-(ETOI(J)+ETOIEV(J)+ETOIPT(J))/EK(J)) &
1151 0 : /ETO(J)
1152 : 19 format('nonadiabatic solution, MODE: ',I2)
1153 : 20 format('eigenvalue: (',1P,D15.8,' , ',D15.8,')')
1154 : 21 format(14X,'freq.[Hz]: ',D11.5,' , sigma: ',D12.6)
1155 : 22 format(14X,'period[d]:',F11.6,1X,' , Q:',F9.6)
1156 : 23 format(9X,'KE growte rate:',F13.9,2X,', KE: ',D12.6)
1157 : 24 format('control: KE growte rate:',F13.9)
1158 : 25 format(' f: (',1P,D15.8,' , ',D15.8,'), abs(f):',D15.8)
1159 : 26 format('control: f: (',1P,D15.8,' , ',D15.8,'), abs(f):',D15.8)
1160 :
1161 : end do
1162 :
1163 : !write(*,'(a55,i12,99(1pd26.16))') 'end LINA s% w(2)**2 Et', &
1164 : ! 2, s% w(2)**2, Et(NZN-1)
1165 :
1166 : return
1167 0 : end subroutine do_LINA
1168 :
1169 :
1170 0 : SUBROUTINE mesa_eos_kap (s,k,G,H, & !input: temp,volume &
1171 : P,PV,PT,E,EV,ET,CP,CPV,dCp_dT_00, &
1172 : Q,QV,QT,OP,OPV,OPT,ierr)
1173 : use rsp_eval_eos_and_kap, only : eval_mesa_eos_and_kap
1174 : type (star_info), pointer :: s
1175 : integer, intent(out) :: ierr
1176 : integer :: k, j
1177 : real(8) :: G,H,P,PV,PT, &
1178 : E,EV,ET,CP,CPV,dCp_dT_00, &
1179 0 : Q,QV,QT,OP,OPV,OPT,cs, &
1180 0 : Pgas,d_Pg_dV,d_Pg_dT,Prad,d_Pr_dT, &
1181 0 : egas,d_egas_dV,d_egas_dT,erad,d_erad_dV,d_erad_dT
1182 : include 'formats'
1183 0 : if (k <= 0 .or. k > NZN) then
1184 0 : j = 0
1185 : else
1186 0 : j = NZN+1-k
1187 : end if
1188 0 : if (is_bad(G+H)) then
1189 0 : write(*,2) 'LINA mesa_eos_kap G H', k, G, H
1190 0 : call mesa_error(__FILE__,__LINE__,'mesa_eos_kap')
1191 : end if
1192 : call eval_mesa_eos_and_kap(s,j,G,H, &
1193 : Pgas,d_Pg_dV,d_Pg_dT,Prad,d_Pr_dT,&
1194 : egas,d_egas_dV,d_egas_dT,erad,d_erad_dV,d_erad_dT, &
1195 : cs,CP,CPV,dCp_dT_00, &
1196 0 : Q,QV,QT,OP,OPV,OPT,ierr)
1197 0 : if (ierr /= 0) return
1198 0 : E = egas + erad
1199 0 : EV = d_egas_dV + d_erad_dV
1200 0 : ET = d_egas_dT + d_erad_dT
1201 0 : P = Pgas + Prad
1202 0 : PV = d_Pg_dV
1203 0 : PT = d_Pg_dT + d_Pr_dT
1204 0 : end SUBROUTINE mesa_eos_kap
1205 :
1206 :
1207 : ! THIS IS A MODIFIED "NUMERICAL RECIPES" SORTING SUBROUTINE
1208 : ! IT SORTS THE VECTOR RA (DIMENSION N) IN ASCENDING ORDER
1209 : ! AND IN THE SAME WAY REARRANGES THE ELEMENTS OF RB (RB IS NOT SORTED)
1210 : ! THE INTEGER VECTOR ISORT CONTAINS INFORMATION ABOUT CHANGES DONE
1211 : ! NAMELY ELEMENT I OF SORTED RA (AND REARRANGED RB) WAS ISORT(I) ELEMENT
1212 : ! OF UNSORTED RA (AND NOT REARANGED RB)
1213 0 : subroutine SORT(N,RA,RB,ISORT)
1214 : integer :: N,ISORT(N)
1215 : real(dp) :: RA(N),RB(N)
1216 : integer :: L,IR,I,J,RRI
1217 0 : real(dp) :: RRA,RRB
1218 :
1219 0 : do I=1,N
1220 0 : ISORT(I)=I
1221 : end do
1222 :
1223 0 : L=N/2+1
1224 0 : IR=N
1225 : 10 continue
1226 0 : if(L>1)then
1227 0 : L=L-1
1228 0 : RRA=RA(L)
1229 0 : RRB=RB(L)
1230 0 : RRI=ISORT(L)
1231 : else
1232 0 : RRA=RA(IR)
1233 0 : RRB=RB(IR)
1234 0 : RRI=ISORT(IR)
1235 0 : RA(IR)=RA(1)
1236 0 : RB(IR)=RB(1)
1237 0 : ISORT(IR)=ISORT(1)
1238 0 : IR=IR-1
1239 0 : if(IR==1)then
1240 0 : RA(1)=RRA
1241 0 : RB(1)=RRB
1242 0 : ISORT(1)=RRI
1243 0 : return
1244 : end if
1245 : end if
1246 0 : I=L
1247 0 : J=L+L
1248 0 : 20 if(J<=IR)then
1249 0 : if(J<IR)then
1250 0 : if(RA(J)<RA(J+1))J=J+1
1251 : end if
1252 0 : if(RRA<RA(J))then
1253 0 : RA(I)=RA(J)
1254 0 : RB(I)=RB(J)
1255 0 : ISORT(I)=ISORT(J)
1256 0 : I=J
1257 0 : J=J+J
1258 : else
1259 0 : J=IR+1
1260 : end if
1261 : GOTO 20
1262 : end if
1263 0 : RA(I)=RRA
1264 0 : RB(I)=RRB
1265 0 : ISORT(I)=RRI
1266 0 : GOTO 10
1267 0 : end subroutine SORT
1268 :
1269 : end module rsp_lina
|
