Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2020 Adam Jermyn & 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 auto_diff_support
21 :
22 : use star_private_def
23 : use const_def, only: dp, sqrt_2_div_3
24 : use auto_diff
25 :
26 : implicit none
27 :
28 : ! current use of xtra's
29 : ! xtra1 is w_div_wc
30 :
31 : public
32 :
33 : contains
34 :
35 52348 : type(auto_diff_real_star_order1) function shift_p1(val_00) result(val_p1)
36 : type(auto_diff_real_star_order1), intent(in) :: val_00
37 : integer :: j
38 52348 : val_p1%val = val_00%val
39 52348 : do j=auto_diff_star_num_vars-2,1,-3 ! p1 gets 00, 00 gets m1, m1 gets 0d0
40 575828 : val_p1%d1Array(j+2) = val_00%d1Array(j+1)
41 575828 : val_p1%d1Array(j+1) = val_00%d1Array(j)
42 575828 : val_p1%d1Array(j) = 0d0
43 : end do
44 52348 : end function shift_p1
45 :
46 52304 : type(auto_diff_real_star_order1) function shift_m1(val_00) result(val_m1)
47 : type(auto_diff_real_star_order1), intent(in) :: val_00
48 : integer :: j
49 52304 : val_m1%val = val_00%val
50 52304 : do j=1,auto_diff_star_num_vars,3 ! m1 gets 00, 00 gets p1, p1 gets 0d0
51 575344 : val_m1%d1Array(j) = val_00%d1Array(j+1)
52 575344 : val_m1%d1Array(j+1) = val_00%d1Array(j+2)
53 575344 : val_m1%d1Array(j+2) = 0d0
54 : end do
55 52304 : end function shift_m1
56 :
57 209392 : subroutine unwrap(var, val, &
58 : dlnd_m1, dlnd_00, dlnd_p1, &
59 : dlnT_m1, dlnT_00, dlnT_p1, &
60 : dw_m1, dw_00, dw_p1, &
61 : dlnR_m1, dlnR_00, dlnR_p1, &
62 : dv_m1, dv_00, dv_p1, &
63 : dL_m1, dL_00, dL_p1, &
64 : dHp_m1, dHp_00, dHp_p1, &
65 : dw_div_wc_m1, dw_div_wc_00, dw_div_wc_p1, &
66 : djrot_m1, djrot_00, djrot_p1, &
67 : dxtra1_m1, dxtra1_00, dxtra1_p1, &
68 : dxtra2_m1, dxtra2_00, dxtra2_p1)
69 : type(auto_diff_real_star_order1), intent(in) :: var
70 : real(dp), intent(out) :: &
71 : val, dlnd_m1, dlnd_00, dlnd_p1, dlnT_m1, dlnT_00, dlnT_p1, &
72 : dw_m1, dw_00, dw_p1, dlnR_m1, dlnR_00, dlnR_p1, &
73 : dv_m1, dv_00, dv_p1, dL_m1, dL_00, dL_p1, &
74 : dHp_m1, dHp_00, dHp_p1, &
75 : dw_div_wc_m1, dw_div_wc_00, dw_div_wc_p1, &
76 : djrot_m1, djrot_00, djrot_p1, &
77 : dxtra1_m1, dxtra1_00, dxtra1_p1, &
78 : dxtra2_m1, dxtra2_00, dxtra2_p1
79 209392 : val = var%val
80 209392 : dlnd_m1 = var%d1Array(i_lnd_m1)
81 209392 : dlnd_00 = var%d1Array(i_lnd_00)
82 209392 : dlnd_p1 = var%d1Array(i_lnd_p1)
83 209392 : dlnT_m1 = var%d1Array(i_lnT_m1)
84 209392 : dlnT_00 = var%d1Array(i_lnT_00)
85 209392 : dlnT_p1 = var%d1Array(i_lnT_p1)
86 209392 : dw_m1 = var%d1Array(i_w_m1)
87 209392 : dw_00 = var%d1Array(i_w_00)
88 209392 : dw_p1 = var%d1Array(i_w_p1)
89 209392 : dlnR_m1 = var%d1Array(i_lnR_m1)
90 209392 : dlnR_00 = var%d1Array(i_lnR_00)
91 209392 : dlnR_p1 = var%d1Array(i_lnR_p1)
92 209392 : dv_m1 = var%d1Array(i_v_m1)
93 209392 : dv_00 = var%d1Array(i_v_00)
94 209392 : dv_p1 = var%d1Array(i_v_p1)
95 209392 : dL_m1 = var%d1Array(i_L_m1)
96 209392 : dL_00 = var%d1Array(i_L_00)
97 209392 : dL_p1 = var%d1Array(i_L_p1)
98 209392 : dHp_m1 = var%d1Array(i_Hp_m1)
99 209392 : dHp_00 = var%d1Array(i_Hp_00)
100 209392 : dHp_p1 = var%d1Array(i_Hp_p1)
101 209392 : dw_div_wc_m1 = var%d1Array(i_w_div_wc_m1)
102 209392 : dw_div_wc_00 = var%d1Array(i_w_div_wc_00)
103 209392 : dw_div_wc_p1 = var%d1Array(i_w_div_wc_p1)
104 209392 : djrot_m1 = var%d1Array(i_jrot_m1)
105 209392 : djrot_00 = var%d1Array(i_jrot_00)
106 209392 : djrot_p1 = var%d1Array(i_jrot_p1)
107 209392 : dxtra1_m1 = var%d1Array(i_xtra1_m1)
108 209392 : dxtra1_00 = var%d1Array(i_xtra1_00)
109 209392 : dxtra1_p1 = var%d1Array(i_xtra1_p1)
110 209392 : dxtra2_m1 = var%d1Array(i_xtra2_m1)
111 209392 : dxtra2_00 = var%d1Array(i_xtra2_00)
112 209392 : dxtra2_p1 = var%d1Array(i_xtra2_p1)
113 209392 : end subroutine unwrap
114 :
115 176 : subroutine wrap(var, val, &
116 : dlnd_m1, dlnd_00, dlnd_p1, &
117 : dlnT_m1, dlnT_00, dlnT_p1, &
118 : dw_m1, dw_00, dw_p1, &
119 : dlnR_m1, dlnR_00, dlnR_p1, &
120 : dv_m1, dv_00, dv_p1, &
121 : dL_m1, dL_00, dL_p1, &
122 : dHp_m1, dHp_00, dHp_p1, &
123 : dw_div_wc_m1, dw_div_wc_00, dw_div_wc_p1, &
124 : djrot_m1, djrot_00, djrot_p1, &
125 : dxtra1_m1, dxtra1_00, dxtra1_p1, &
126 : dxtra2_m1, dxtra2_00, dxtra2_p1)
127 : type(auto_diff_real_star_order1), intent(out) :: var
128 : real(dp), intent(in) :: &
129 : val, dlnd_m1, dlnd_00, dlnd_p1, dlnT_m1, dlnT_00, dlnT_p1, &
130 : dw_m1, dw_00, dw_p1, dlnR_m1, dlnR_00, dlnR_p1, &
131 : dv_m1, dv_00, dv_p1, dL_m1, dL_00, dL_p1, &
132 : dHp_m1, dHp_00, dHp_p1, &
133 : dw_div_wc_m1, dw_div_wc_00, dw_div_wc_p1, &
134 : djrot_m1, djrot_00, djrot_p1, &
135 : dxtra1_m1, dxtra1_00, dxtra1_p1, &
136 : dxtra2_m1, dxtra2_00, dxtra2_p1
137 176 : var%val = val
138 176 : var%d1Array(i_lnd_m1) = dlnd_m1
139 176 : var%d1Array(i_lnd_00) = dlnd_00
140 176 : var%d1Array(i_lnd_p1) = dlnd_p1
141 176 : var%d1Array(i_lnT_m1) = dlnT_m1
142 176 : var%d1Array(i_lnT_00) = dlnT_00
143 176 : var%d1Array(i_lnT_p1) = dlnT_p1
144 176 : var%d1Array(i_w_m1) = dw_m1
145 176 : var%d1Array(i_w_00) = dw_00
146 176 : var%d1Array(i_w_p1) = dw_p1
147 176 : var%d1Array(i_lnR_m1) = dlnR_m1
148 176 : var%d1Array(i_lnR_00) = dlnR_00
149 176 : var%d1Array(i_lnR_p1) = dlnR_p1
150 176 : var%d1Array(i_v_m1) = dv_m1
151 176 : var%d1Array(i_v_00) = dv_00
152 176 : var%d1Array(i_v_p1) = dv_p1
153 176 : var%d1Array(i_L_m1) = dL_m1
154 176 : var%d1Array(i_L_00) = dL_00
155 176 : var%d1Array(i_L_p1) = dL_p1
156 176 : var%d1Array(i_Hp_m1) = dHp_m1
157 176 : var%d1Array(i_Hp_00) = dHp_00
158 176 : var%d1Array(i_Hp_p1) = dHp_p1
159 176 : var%d1Array(i_w_div_wc_m1) = dw_div_wc_m1
160 176 : var%d1Array(i_w_div_wc_00) = dw_div_wc_00
161 176 : var%d1Array(i_w_div_wc_p1) = dw_div_wc_p1
162 176 : var%d1Array(i_jrot_m1) = djrot_m1
163 176 : var%d1Array(i_jrot_00) = djrot_00
164 176 : var%d1Array(i_jrot_p1) = djrot_p1
165 176 : var%d1Array(i_xtra1_m1) = dxtra1_m1
166 176 : var%d1Array(i_xtra1_00) = dxtra1_00
167 176 : var%d1Array(i_xtra1_p1) = dxtra1_p1
168 176 : var%d1Array(i_xtra2_m1) = dxtra2_m1
169 176 : var%d1Array(i_xtra2_00) = dxtra2_00
170 176 : var%d1Array(i_xtra2_p1) = dxtra2_p1
171 176 : end subroutine wrap
172 :
173 :
174 : ! The following routines turn regular star variables into auto_diff_real_star_order1 variables.
175 : ! For independent variables this is a straightforward wrapping. For dependent variables like eos and kap
176 : ! outputs we also pull in information about their partials from the relevant module.
177 :
178 : !----------------------------------------------------------------------------------------------------
179 : !
180 : ! We begin with quantities defined on cells (T, rho, e_turb, kap, pressure, energy, entropy).
181 : ! We need these on cells k-1, k and k+1.
182 : !
183 : !----------------------------------------------------------------------------------------------------
184 :
185 :
186 212160 : function wrap_T_m1(s, k) result(T_m1)
187 : type (star_info), pointer :: s
188 : type(auto_diff_real_star_order1) :: T_m1
189 : integer, intent(in) :: k
190 212160 : T_m1 = 0d0
191 212160 : if (k > 1) then
192 212160 : T_m1 % val = s%T(k-1)
193 212160 : T_m1 % d1Array(i_lnT_m1) = s%T(k-1)
194 : end if
195 212160 : end function wrap_T_m1
196 :
197 212292 : function wrap_T_00(s, k) result(T_00)
198 : type (star_info), pointer :: s
199 : type(auto_diff_real_star_order1) :: T_00
200 : integer, intent(in) :: k
201 212292 : T_00 = 0d0
202 212292 : T_00 % val = s%T(k)
203 212292 : T_00 % d1Array(i_lnT_00) = s%T(k)
204 212292 : end function wrap_T_00
205 :
206 0 : function wrap_T_p1(s, k) result(T_p1)
207 : type (star_info), pointer :: s
208 : type(auto_diff_real_star_order1) :: T_p1
209 : integer, intent(in) :: k
210 0 : T_p1 = 0d0
211 0 : if (k < s%nz) then
212 0 : T_p1 % val = s%T(k+1)
213 0 : T_p1 % d1Array(i_lnT_p1) = s%T(k+1)
214 : end if
215 0 : end function wrap_T_p1
216 :
217 0 : function wrap_lnT_m1(s, k) result(lnT_m1)
218 : type (star_info), pointer :: s
219 : type(auto_diff_real_star_order1) :: lnT_m1
220 : integer, intent(in) :: k
221 0 : lnT_m1 = 0d0
222 0 : if (k > 1) then
223 0 : lnT_m1 % val = s%lnT(k-1)
224 0 : lnT_m1 % d1Array(i_lnT_m1) = 1d0
225 : end if
226 0 : end function wrap_lnT_m1
227 :
228 44 : function wrap_lnT_00(s, k) result(lnT_00)
229 : type (star_info), pointer :: s
230 : type(auto_diff_real_star_order1) :: lnT_00
231 : integer, intent(in) :: k
232 44 : lnT_00 = 0d0
233 44 : lnT_00 % val = s%lnT(k)
234 44 : lnT_00 % d1Array(i_lnT_00) = 1d0
235 44 : end function wrap_lnT_00
236 :
237 0 : function wrap_lnT_p1(s, k) result(lnT_p1)
238 : type (star_info), pointer :: s
239 : type(auto_diff_real_star_order1) :: lnT_p1
240 : integer, intent(in) :: k
241 0 : lnT_p1 = 0d0
242 0 : if (k < s%nz) then
243 0 : lnT_p1 % val = s%lnT(k+1)
244 0 : lnT_p1 % d1Array(i_lnT_p1) = 1d0
245 : end if
246 0 : end function wrap_lnT_p1
247 :
248 0 : function wrap_dxh_lnT(s, k) result(dxh_lnT)
249 : type (star_info), pointer :: s
250 : type(auto_diff_real_star_order1) :: dxh_lnT
251 : integer, intent(in) :: k
252 0 : dxh_lnT = 0d0
253 0 : dxh_lnT % val = s%dxh_lnT(k)
254 0 : dxh_lnT % d1Array(i_lnT_00) = 1d0
255 0 : end function wrap_dxh_lnT
256 :
257 239784 : function wrap_d_m1(s, k) result(d_m1)
258 : type (star_info), pointer :: s
259 : type(auto_diff_real_star_order1) :: d_m1
260 : integer, intent(in) :: k
261 239784 : d_m1 = 0d0
262 239784 : if (k > 1) then
263 239784 : d_m1 % val = s%rho(k-1)
264 239784 : d_m1 % d1Array(i_lnd_m1) = s%rho(k-1)
265 : end if
266 239784 : end function wrap_d_m1
267 :
268 292330 : function wrap_d_00(s, k) result(d_00)
269 : type (star_info), pointer :: s
270 : type(auto_diff_real_star_order1) :: d_00
271 : integer, intent(in) :: k
272 292330 : d_00 = 0d0
273 292330 : d_00 % val = s%rho(k)
274 292330 : d_00 % d1Array(i_lnd_00) = s%rho(k)
275 292330 : end function wrap_d_00
276 :
277 0 : function wrap_d_p1(s, k) result(d_p1)
278 : type (star_info), pointer :: s
279 : type(auto_diff_real_star_order1) :: d_p1
280 : integer, intent(in) :: k
281 0 : d_p1 = 0d0
282 0 : if (k < s%nz) then
283 0 : d_p1 % val = s%rho(k+1)
284 0 : d_p1 % d1Array(i_lnd_p1) = s%rho(k+1)
285 : end if
286 0 : end function wrap_d_p1
287 :
288 0 : function wrap_lnd_m1(s, k) result(lnd_m1)
289 : type (star_info), pointer :: s
290 : type(auto_diff_real_star_order1) :: lnd_m1
291 : integer, intent(in) :: k
292 0 : lnd_m1 = 0d0
293 0 : if (k > 1) then
294 0 : lnd_m1 % val = s%lnd(k-1)
295 0 : lnd_m1 % d1Array(i_lnd_m1) = 1d0
296 : end if
297 0 : end function wrap_lnd_m1
298 :
299 0 : function wrap_lnd_00(s, k) result(lnd_00)
300 : type (star_info), pointer :: s
301 : type(auto_diff_real_star_order1) :: lnd_00
302 : integer, intent(in) :: k
303 0 : lnd_00 = 0d0
304 0 : lnd_00 % val = s%lnd(k)
305 0 : lnd_00 % d1Array(i_lnd_00) = 1d0
306 0 : end function wrap_lnd_00
307 :
308 0 : function wrap_lnd_p1(s, k) result(lnd_p1)
309 : type (star_info), pointer :: s
310 : type(auto_diff_real_star_order1) :: lnd_p1
311 : integer, intent(in) :: k
312 0 : lnd_p1 = 0d0
313 0 : if (k < s%nz) then
314 0 : lnd_p1 % val = s%lnd(k+1)
315 0 : lnd_p1 % d1Array(i_lnd_p1) = 1d0
316 : end if
317 0 : end function wrap_lnd_p1
318 :
319 0 : function wrap_dxh_lnd(s, k) result(dxh_lnd)
320 : type (star_info), pointer :: s
321 : type(auto_diff_real_star_order1) :: dxh_lnd
322 : integer, intent(in) :: k
323 0 : dxh_lnd = 0d0
324 0 : dxh_lnd % val = s%dxh_lnd(k)
325 0 : dxh_lnd % d1Array(i_lnd_00) = 1d0
326 0 : end function wrap_dxh_lnd
327 :
328 0 : function wrap_w_m1(s, k) result(w_m1)
329 : type (star_info), pointer :: s
330 : type(auto_diff_real_star_order1) :: w_m1
331 : integer, intent(in) :: k
332 0 : w_m1 = 0d0
333 0 : if (k > 1) then
334 0 : w_m1 % val = s%w(k-1)
335 0 : w_m1 % d1Array(i_w_m1) = 1d0
336 : end if
337 0 : end function wrap_w_m1
338 :
339 0 : function wrap_w_00(s, k) result(w_00)
340 : type (star_info), pointer :: s
341 : type(auto_diff_real_star_order1) :: w_00
342 : integer, intent(in) :: k
343 0 : w_00 = 0d0
344 0 : w_00 % val = s%w(k)
345 0 : w_00 % d1Array(i_w_00) = 1d0
346 0 : end function wrap_w_00
347 :
348 0 : function wrap_w_p1(s, k) result(w_p1)
349 : type (star_info), pointer :: s
350 : type(auto_diff_real_star_order1) :: w_p1
351 : integer, intent(in) :: k
352 0 : w_p1 = 0d0
353 0 : if (k < s%nz) then
354 0 : w_p1 % val = s%w(k+1)
355 0 : w_p1 % d1Array(i_w_p1) = 1d0
356 : end if
357 0 : end function wrap_w_p1
358 :
359 0 : real(dp) function get_etrb(s,k)
360 : type (star_info), pointer :: s
361 : integer, intent(in) :: k
362 0 : get_etrb = pow2(s% w(k))
363 0 : end function get_etrb
364 :
365 0 : real(dp) function get_etrb_start(s,k)
366 : type (star_info), pointer :: s
367 : integer, intent(in) :: k
368 0 : get_etrb_start = pow2(s% w_start(k))
369 0 : end function get_etrb_start
370 :
371 0 : real(dp) function get_RSP2_conv_velocity(s,k) result (cv) ! at face k
372 : type (star_info), pointer :: s
373 : integer, intent(in) :: k
374 0 : real(dp) :: alfa, beta
375 0 : if (k == 1) then
376 : cv = 0d0
377 : else
378 0 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
379 0 : beta = 1d0 - alfa
380 0 : cv = sqrt_2_div_3*(alfa*s% w(k) + beta*s% w(k-1))
381 : end if
382 0 : end function get_RSP2_conv_velocity
383 :
384 0 : real(dp) function get_w(s,k)
385 : type (star_info), pointer :: s
386 : integer, intent(in) :: k
387 0 : get_w = s% w(k)
388 0 : end function get_w
389 :
390 0 : real(dp) function get_w_start(s,k)
391 : type (star_info), pointer :: s
392 : integer, intent(in) :: k
393 0 : get_w_start = s% w_start(k)
394 0 : end function get_w_start
395 :
396 0 : function wrap_etrb_m1(s, k) result(etrb_m1)
397 : type (star_info), pointer :: s
398 : type(auto_diff_real_star_order1) :: etrb_m1
399 : integer, intent(in) :: k
400 0 : etrb_m1 = pow2(wrap_w_m1(s,k))
401 0 : end function wrap_etrb_m1
402 :
403 0 : function wrap_etrb_00(s, k) result(etrb_00)
404 : type (star_info), pointer :: s
405 : type(auto_diff_real_star_order1) :: etrb_00
406 : integer, intent(in) :: k
407 0 : etrb_00 = pow2(wrap_w_00(s,k))
408 0 : end function wrap_etrb_00
409 :
410 0 : function wrap_etrb_p1(s, k) result(etrb_p1)
411 : type (star_info), pointer :: s
412 : type(auto_diff_real_star_order1) :: etrb_p1
413 : integer, intent(in) :: k
414 0 : etrb_p1 = pow2(wrap_w_p1(s,k))
415 0 : end function wrap_etrb_p1
416 :
417 159856 : function wrap_kap_m1(s, k) result(kap_m1)
418 : type (star_info), pointer :: s
419 : type(auto_diff_real_star_order1) :: kap_m1
420 : integer, intent(in) :: k
421 159856 : kap_m1 = 0d0
422 159856 : if (k > 1) then
423 159856 : kap_m1 % val = s%opacity(k-1)
424 159856 : kap_m1 % d1Array(i_lnd_m1) = s%d_opacity_dlnd(k-1)
425 159856 : kap_m1 % d1Array(i_lnT_m1) = s%d_opacity_dlnT(k-1)
426 : end if
427 159856 : end function wrap_kap_m1
428 :
429 159988 : function wrap_kap_00(s, k) result(kap_00)
430 : type (star_info), pointer :: s
431 : type(auto_diff_real_star_order1) :: kap_00
432 : integer, intent(in) :: k
433 159988 : kap_00 = 0d0
434 159988 : kap_00 % val = s%opacity(k)
435 159988 : kap_00 % d1Array(i_lnd_00) = s%d_opacity_dlnd(k)
436 159988 : kap_00 % d1Array(i_lnT_00) = s%d_opacity_dlnT(k)
437 159988 : end function wrap_kap_00
438 :
439 0 : function wrap_kap_p1(s, k) result(kap_p1)
440 : type (star_info), pointer :: s
441 : type(auto_diff_real_star_order1) :: kap_p1
442 : integer, intent(in) :: k
443 0 : kap_p1 = 0d0
444 0 : if (k < s%nz) then
445 0 : kap_p1 % val = s%opacity(k+1)
446 0 : kap_p1 % d1Array(i_lnd_p1) = s%d_opacity_dlnd(k+1)/s% rho(k+1)
447 0 : kap_p1 % d1Array(i_lnT_p1) = s%d_opacity_dlnT(k+1)
448 : end if
449 0 : end function wrap_kap_p1
450 :
451 0 : function wrap_s_m1(s, k) result(s_m1)
452 : type (star_info), pointer :: s
453 : type(auto_diff_real_star_order1) :: s_m1
454 : integer, intent(in) :: k
455 0 : s_m1 = 0d0
456 0 : if (k > 1) then
457 0 : s_m1%val = s% entropy(k-1)
458 0 : s_m1%d1Array(i_lnd_m1) = s% dS_dRho_for_partials(k-1)*s% rho(k-1)
459 0 : s_m1%d1Array(i_lnT_m1) = s% dS_dT_for_partials(k-1)*s% T(k-1)
460 : end if
461 0 : end function wrap_s_m1
462 :
463 0 : function wrap_s_00(s, k) result(s_00)
464 : type (star_info), pointer :: s
465 : type(auto_diff_real_star_order1) :: s_00
466 : integer, intent(in) :: k
467 0 : s_00 = 0d0
468 0 : s_00%val = s% entropy(k)
469 0 : s_00%d1Array(i_lnd_00) = s% dS_dRho_for_partials(k)*s% rho(k)
470 0 : s_00%d1Array(i_lnT_00) = s% dS_dT_for_partials(k)*s% T(k)
471 0 : end function wrap_s_00
472 :
473 0 : function wrap_s_p1(s, k) result(s_p1)
474 : type (star_info), pointer :: s
475 : type(auto_diff_real_star_order1) :: s_p1
476 : integer, intent(in) :: k
477 0 : s_p1 = 0d0
478 0 : if (k < s%nz) then
479 0 : s_p1%val = s% entropy(k+1)
480 0 : s_p1%d1Array(i_lnd_p1) = s% dS_dRho_for_partials(k+1)*s% rho(k+1)
481 0 : s_p1%d1Array(i_lnT_p1) = s% dS_dT_for_partials(k+1)*s% T(k+1)
482 : end if
483 0 : end function wrap_s_p1
484 :
485 0 : function wrap_e_m1(s, k) result(e_m1)
486 : type (star_info), pointer :: s
487 : type(auto_diff_real_star_order1) :: e_m1
488 : integer, intent(in) :: k
489 0 : e_m1 = 0d0
490 0 : if (k > 1) then
491 0 : e_m1%val = s% energy(k-1)
492 0 : e_m1%d1Array(i_lnd_m1) = s% dE_dRho_for_partials(k-1)*s% rho(k-1)
493 0 : e_m1%d1Array(i_lnT_m1) = s% Cv_for_partials(k-1)*s% T(k-1)
494 : end if
495 0 : end function wrap_e_m1
496 :
497 0 : function wrap_e_00(s, k) result(e_00)
498 : type (star_info), pointer :: s
499 : type(auto_diff_real_star_order1) :: e_00
500 : integer, intent(in) :: k
501 0 : e_00 = 0d0
502 0 : e_00%val = s% energy(k)
503 0 : e_00%d1Array(i_lnd_00) = s% dE_dRho_for_partials(k)*s% rho(k)
504 0 : e_00%d1Array(i_lnT_00) = s% Cv_for_partials(k)*s% T(k)
505 0 : end function wrap_e_00
506 :
507 0 : function wrap_e_p1(s, k) result(e_p1)
508 : type (star_info), pointer :: s
509 : type(auto_diff_real_star_order1) :: e_p1
510 : integer, intent(in) :: k
511 0 : e_p1 = 0d0
512 0 : if (k < s%nz) then
513 0 : e_p1%val = s% energy(k+1)
514 0 : e_p1%d1Array(i_lnd_p1) = s% dE_dRho_for_partials(k+1)*s% rho(k+1)
515 0 : e_p1%d1Array(i_lnT_p1) = s% Cv_for_partials(k+1)*s% T(k+1)
516 : end if
517 0 : end function wrap_e_p1
518 :
519 396696 : function wrap_Peos_m1(s, k) result(Peos_m1)
520 : type (star_info), pointer :: s
521 : type(auto_diff_real_star_order1) :: Peos_m1
522 : integer, intent(in) :: k
523 396696 : Peos_m1 = 0d0
524 396696 : if (k > 1) then
525 396696 : Peos_m1%val = s% Peos(k-1)
526 396696 : Peos_m1%d1Array(i_lnd_m1) = s%Peos(k-1) * s% chiRho_for_partials(k-1)
527 396696 : Peos_m1%d1Array(i_lnT_m1) = s%Peos(k-1) * s% chiT_for_partials(k-1)
528 : end if
529 396696 : end function wrap_Peos_m1
530 :
531 501546 : function wrap_Peos_00(s, k) result(Peos_00)
532 : type (star_info), pointer :: s
533 : type(auto_diff_real_star_order1) :: Peos_00
534 : integer, intent(in) :: k
535 501546 : Peos_00 = 0d0
536 501546 : Peos_00%val = s% Peos(k)
537 501546 : Peos_00%d1Array(i_lnd_00) = s%Peos(k) * s% chiRho_for_partials(k)
538 501546 : Peos_00%d1Array(i_lnT_00) = s%Peos(k) * s% chiT_for_partials(k)
539 501546 : end function wrap_Peos_00
540 :
541 0 : function wrap_Peos_p1(s, k) result(Peos_p1)
542 : type (star_info), pointer :: s
543 : type(auto_diff_real_star_order1) :: Peos_p1
544 : integer, intent(in) :: k
545 0 : Peos_p1 = 0d0
546 0 : if (k < s%nz) then
547 0 : Peos_p1%val = s% Peos(k+1)
548 0 : Peos_p1%d1Array(i_lnd_p1) = s%Peos(k+1) * s% chiRho_for_partials(k+1)
549 0 : Peos_p1%d1Array(i_lnT_p1) = s%Peos(k+1) * s% chiT_for_partials(k+1)
550 : end if
551 0 : end function wrap_Peos_p1
552 :
553 0 : function wrap_lnPeos_m1(s, k) result(lnPeos_m1)
554 : type (star_info), pointer :: s
555 : type(auto_diff_real_star_order1) :: lnPeos_m1
556 : integer, intent(in) :: k
557 0 : lnPeos_m1 = 0d0
558 0 : if (k > 1) then
559 0 : lnPeos_m1%val = s% lnPeos(k-1)
560 0 : lnPeos_m1%d1Array(i_lnd_m1) = s% chiRho_for_partials(k-1)
561 0 : lnPeos_m1%d1Array(i_lnT_m1) = s% chiT_for_partials(k-1)
562 : end if
563 0 : end function wrap_lnPeos_m1
564 :
565 44 : function wrap_lnPeos_00(s, k) result(lnPeos_00)
566 : type (star_info), pointer :: s
567 : type(auto_diff_real_star_order1) :: lnPeos_00
568 : integer, intent(in) :: k
569 44 : lnPeos_00 = 0d0
570 44 : lnPeos_00%val = s% lnPeos(k)
571 44 : lnPeos_00%d1Array(i_lnd_00) = s% chiRho_for_partials(k)
572 44 : lnPeos_00%d1Array(i_lnT_00) = s% chiT_for_partials(k)
573 44 : end function wrap_lnPeos_00
574 :
575 0 : function wrap_lnPeos_p1(s, k) result(lnPeos_p1)
576 : type (star_info), pointer :: s
577 : type(auto_diff_real_star_order1) :: lnPeos_p1
578 : integer, intent(in) :: k
579 0 : lnPeos_p1 = 0d0
580 0 : if (k < s%nz) then
581 0 : lnPeos_p1%val = s% lnPeos(k+1)
582 0 : lnPeos_p1%d1Array(i_lnd_p1) = s% chiRho_for_partials(k+1)
583 0 : lnPeos_p1%d1Array(i_lnT_p1) = s% chiT_for_partials(k+1)
584 : end if
585 0 : end function wrap_lnPeos_p1
586 :
587 79928 : function wrap_ChiRho_m1(s, k) result(ChiRho_m1)
588 : use eos_def, only: i_ChiRho
589 : type (star_info), pointer :: s
590 : type(auto_diff_real_star_order1) :: ChiRho_m1
591 : integer, intent(in) :: k
592 79928 : ChiRho_m1 = 0d0
593 79928 : if (k > 1) then
594 79928 : ChiRho_m1%val = s% ChiRho(k-1)
595 79928 : ChiRho_m1%d1Array(i_lnd_m1) = s% d_eos_dlnd(i_ChiRho,k-1)
596 79928 : ChiRho_m1%d1Array(i_lnT_m1) = s% d_eos_dlnT(i_ChiRho,k-1)
597 : end if
598 79928 : end function wrap_ChiRho_m1
599 :
600 79994 : function wrap_ChiRho_00(s, k) result(ChiRho_00)
601 79928 : use eos_def, only: i_ChiRho
602 : type (star_info), pointer :: s
603 : type(auto_diff_real_star_order1) :: ChiRho_00
604 : integer, intent(in) :: k
605 79994 : ChiRho_00 = 0d0
606 79994 : ChiRho_00%val = s% ChiRho(k)
607 79994 : ChiRho_00%d1Array(i_lnd_00) = s% d_eos_dlnd(i_ChiRho,k)
608 79994 : ChiRho_00%d1Array(i_lnT_00) = s% d_eos_dlnT(i_ChiRho,k)
609 79994 : end function wrap_ChiRho_00
610 :
611 0 : function wrap_ChiRho_p1(s, k) result(ChiRho_p1)
612 79994 : use eos_def, only: i_ChiRho
613 : type (star_info), pointer :: s
614 : type(auto_diff_real_star_order1) :: ChiRho_p1
615 : integer, intent(in) :: k
616 0 : ChiRho_p1 = 0d0
617 0 : if (k < s% nz) then
618 0 : ChiRho_p1%val = s% ChiRho(k+1)
619 0 : ChiRho_p1%d1Array(i_lnd_p1) = s% d_eos_dlnd(i_ChiRho,k+1)
620 0 : ChiRho_p1%d1Array(i_lnT_p1) = s% d_eos_dlnT(i_ChiRho,k+1)
621 : end if
622 0 : end function wrap_ChiRho_p1
623 :
624 79928 : function wrap_ChiT_m1(s, k) result(ChiT_m1)
625 0 : use eos_def, only: i_ChiT
626 : type (star_info), pointer :: s
627 : type(auto_diff_real_star_order1) :: ChiT_m1
628 : integer, intent(in) :: k
629 79928 : ChiT_m1 = 0d0
630 79928 : if (k > 1) then
631 79928 : ChiT_m1%val = s% ChiT(k-1)
632 79928 : ChiT_m1%d1Array(i_lnd_m1) = s% d_eos_dlnd(i_ChiT,k-1)
633 79928 : ChiT_m1%d1Array(i_lnT_m1) = s% d_eos_dlnT(i_ChiT,k-1)
634 : end if
635 79928 : end function wrap_ChiT_m1
636 :
637 79994 : function wrap_ChiT_00(s, k) result(ChiT_00)
638 79928 : use eos_def, only: i_ChiT
639 : type (star_info), pointer :: s
640 : type(auto_diff_real_star_order1) :: ChiT_00
641 : integer, intent(in) :: k
642 79994 : ChiT_00 = 0d0
643 79994 : ChiT_00%val = s% ChiT(k)
644 79994 : ChiT_00%d1Array(i_lnd_00) = s% d_eos_dlnd(i_ChiT,k)
645 79994 : ChiT_00%d1Array(i_lnT_00) = s% d_eos_dlnT(i_ChiT,k)
646 79994 : end function wrap_ChiT_00
647 :
648 0 : function wrap_ChiT_p1(s, k) result(ChiT_p1)
649 79994 : use eos_def, only: i_ChiT
650 : type (star_info), pointer :: s
651 : type(auto_diff_real_star_order1) :: ChiT_p1
652 : integer, intent(in) :: k
653 0 : ChiT_p1 = 0d0
654 0 : if (k < s% nz) then
655 0 : ChiT_p1%val = s% ChiT(k+1)
656 0 : ChiT_p1%d1Array(i_lnd_p1) = s% d_eos_dlnd(i_ChiT,k+1)
657 0 : ChiT_p1%d1Array(i_lnT_p1) = s% d_eos_dlnT(i_ChiT,k+1)
658 : end if
659 0 : end function wrap_ChiT_p1
660 :
661 79928 : function wrap_Cp_m1(s, k) result(Cp_m1)
662 0 : use eos_def, only: i_Cp
663 : type (star_info), pointer :: s
664 : type(auto_diff_real_star_order1) :: Cp_m1
665 : integer, intent(in) :: k
666 79928 : Cp_m1 = 0d0
667 79928 : if (k > 1) then
668 79928 : Cp_m1%val = s% Cp(k-1)
669 79928 : Cp_m1%d1Array(i_lnd_m1) = s% d_eos_dlnd(i_Cp,k-1)
670 79928 : Cp_m1%d1Array(i_lnT_m1) = s% d_eos_dlnT(i_Cp,k-1)
671 : end if
672 79928 : end function wrap_Cp_m1
673 :
674 79994 : function wrap_Cp_00(s, k) result(Cp_00)
675 79928 : use eos_def, only: i_Cp
676 : type (star_info), pointer :: s
677 : type(auto_diff_real_star_order1) :: Cp_00
678 : integer, intent(in) :: k
679 79994 : Cp_00 = 0d0
680 79994 : Cp_00%val = s% Cp(k)
681 79994 : Cp_00%d1Array(i_lnd_00) = s% d_eos_dlnd(i_Cp,k)
682 79994 : Cp_00%d1Array(i_lnT_00) = s% d_eos_dlnT(i_Cp,k)
683 79994 : end function wrap_Cp_00
684 :
685 0 : function wrap_Cp_p1(s, k) result(Cp_p1)
686 79994 : use eos_def, only: i_Cp
687 : type (star_info), pointer :: s
688 : type(auto_diff_real_star_order1) :: Cp_p1
689 : integer, intent(in) :: k
690 0 : Cp_p1 = 0d0
691 0 : if (k < s% nz) then
692 0 : Cp_p1%val = s% Cp(k+1)
693 0 : Cp_p1%d1Array(i_lnd_p1) = s% d_eos_dlnd(i_Cp,k+1)
694 0 : Cp_p1%d1Array(i_lnT_p1) = s% d_eos_dlnT(i_Cp,k+1)
695 : end if
696 0 : end function wrap_Cp_p1
697 :
698 79928 : function wrap_grad_ad_m1(s, k) result(grad_ad_m1)
699 0 : use eos_def, only: i_grad_ad
700 : type (star_info), pointer :: s
701 : type(auto_diff_real_star_order1) :: grad_ad_m1
702 : integer, intent(in) :: k
703 79928 : grad_ad_m1 = 0d0
704 79928 : if (k > 1) then
705 79928 : grad_ad_m1%val = s% grada(k-1)
706 79928 : grad_ad_m1%d1Array(i_lnd_m1) = s% d_eos_dlnd(i_grad_ad,k-1)
707 79928 : grad_ad_m1%d1Array(i_lnT_m1) = s% d_eos_dlnT(i_grad_ad,k-1)
708 : end if
709 79928 : end function wrap_grad_ad_m1
710 :
711 79994 : function wrap_grad_ad_00(s, k) result(grad_ad_00)
712 79928 : use eos_def, only: i_grad_ad
713 : type (star_info), pointer :: s
714 : type(auto_diff_real_star_order1) :: grad_ad_00
715 : integer, intent(in) :: k
716 79994 : grad_ad_00 = 0d0
717 79994 : grad_ad_00%val = s% grada(k)
718 79994 : grad_ad_00%d1Array(i_lnd_00) = s% d_eos_dlnd(i_grad_ad,k)
719 79994 : grad_ad_00%d1Array(i_lnT_00) = s% d_eos_dlnT(i_grad_ad,k)
720 79994 : end function wrap_grad_ad_00
721 :
722 0 : function wrap_grad_ad_p1(s, k) result(grad_ad_p1)
723 79994 : use eos_def, only: i_grad_ad
724 : type (star_info), pointer :: s
725 : type(auto_diff_real_star_order1) :: grad_ad_p1
726 : integer, intent(in) :: k
727 0 : grad_ad_p1 = 0d0
728 0 : if (k < s% nz) then
729 0 : grad_ad_p1%val = s% grada(k+1)
730 0 : grad_ad_p1%d1Array(i_lnd_p1) = s% d_eos_dlnd(i_grad_ad,k+1)
731 0 : grad_ad_p1%d1Array(i_lnT_p1) = s% d_eos_dlnT(i_grad_ad,k+1)
732 : end if
733 0 : end function wrap_grad_ad_p1
734 :
735 0 : function wrap_gamma1_m1(s, k) result(gamma1_m1)
736 0 : use eos_def, only: i_gamma1
737 : type (star_info), pointer :: s
738 : type(auto_diff_real_star_order1) :: gamma1_m1
739 : integer, intent(in) :: k
740 0 : gamma1_m1 = 0d0
741 0 : if (k > 1) then
742 0 : gamma1_m1%val = s% gamma1(k-1)
743 0 : gamma1_m1%d1Array(i_lnd_m1) = s% d_eos_dlnd(i_gamma1,k-1)
744 0 : gamma1_m1%d1Array(i_lnT_m1) = s% d_eos_dlnT(i_gamma1,k-1)
745 : end if
746 0 : end function wrap_gamma1_m1
747 :
748 0 : function wrap_gamma1_00(s, k) result(gamma1_00)
749 0 : use eos_def, only: i_gamma1
750 : type (star_info), pointer :: s
751 : type(auto_diff_real_star_order1) :: gamma1_00
752 : integer, intent(in) :: k
753 0 : gamma1_00 = 0d0
754 0 : gamma1_00%val = s% gamma1(k)
755 0 : gamma1_00%d1Array(i_lnd_00) = s% d_eos_dlnd(i_gamma1,k)
756 0 : gamma1_00%d1Array(i_lnT_00) = s% d_eos_dlnT(i_gamma1,k)
757 0 : end function wrap_gamma1_00
758 :
759 0 : function wrap_gamma1_p1(s, k) result(gamma1_p1)
760 0 : use eos_def, only: i_gamma1
761 : type (star_info), pointer :: s
762 : type(auto_diff_real_star_order1) :: gamma1_p1
763 : integer, intent(in) :: k
764 0 : gamma1_p1 = 0d0
765 0 : if (k < s% nz) then
766 0 : gamma1_p1%val = s% gamma1(k+1)
767 0 : gamma1_p1%d1Array(i_lnd_p1) = s% d_eos_dlnd(i_gamma1,k+1)
768 0 : gamma1_p1%d1Array(i_lnT_p1) = s% d_eos_dlnT(i_gamma1,k+1)
769 : end if
770 0 : end function wrap_gamma1_p1
771 :
772 0 : function wrap_latent_ddlnT_m1(s, k) result(latent_ddlnT_m1)
773 0 : use eos_def, only: i_latent_ddlnT
774 : type (star_info), pointer :: s
775 : type(auto_diff_real_star_order1) :: latent_ddlnT_m1
776 : integer, intent(in) :: k
777 0 : latent_ddlnT_m1 = 0d0
778 0 : if (k > 1) then
779 0 : latent_ddlnT_m1%val = s% latent_ddlnT(k-1)
780 0 : latent_ddlnT_m1%d1Array(i_lnd_m1) = s% d_eos_dlnd(i_latent_ddlnT,k-1)
781 0 : latent_ddlnT_m1%d1Array(i_lnT_m1) = s% d_eos_dlnT(i_latent_ddlnT,k-1)
782 : end if
783 0 : end function wrap_latent_ddlnT_m1
784 :
785 0 : function wrap_latent_ddlnT_00(s, k) result(latent_ddlnT_00)
786 0 : use eos_def, only: i_latent_ddlnT
787 : type (star_info), pointer :: s
788 : type(auto_diff_real_star_order1) :: latent_ddlnT_00
789 : integer, intent(in) :: k
790 0 : latent_ddlnT_00 = 0d0
791 0 : latent_ddlnT_00%val = s% latent_ddlnT(k)
792 0 : latent_ddlnT_00%d1Array(i_lnd_00) = s% d_eos_dlnd(i_latent_ddlnT,k)
793 0 : latent_ddlnT_00%d1Array(i_lnT_00) = s% d_eos_dlnT(i_latent_ddlnT,k)
794 0 : end function wrap_latent_ddlnT_00
795 :
796 0 : function wrap_latent_ddlnT_p1(s, k) result(latent_ddlnT_p1)
797 0 : use eos_def, only: i_latent_ddlnT
798 : type (star_info), pointer :: s
799 : type(auto_diff_real_star_order1) :: latent_ddlnT_p1
800 : integer, intent(in) :: k
801 0 : latent_ddlnT_p1 = 0d0
802 0 : if (k < s% nz) then
803 0 : latent_ddlnT_p1%val = s% latent_ddlnT(k+1)
804 0 : latent_ddlnT_p1%d1Array(i_lnd_p1) = s% d_eos_dlnd(i_latent_ddlnT,k+1)
805 0 : latent_ddlnT_p1%d1Array(i_lnT_p1) = s% d_eos_dlnT(i_latent_ddlnT,k+1)
806 : end if
807 0 : end function wrap_latent_ddlnT_p1
808 :
809 0 : function wrap_latent_ddlnRho_m1(s, k) result(latent_ddlnRho_m1)
810 0 : use eos_def, only: i_latent_ddlnRho
811 : type (star_info), pointer :: s
812 : type(auto_diff_real_star_order1) :: latent_ddlnRho_m1
813 : integer, intent(in) :: k
814 0 : latent_ddlnRho_m1 = 0d0
815 0 : if (k > 1) then
816 0 : latent_ddlnRho_m1%val = s% latent_ddlnRho(k-1)
817 0 : latent_ddlnRho_m1%d1Array(i_lnd_m1) = s% d_eos_dlnd(i_latent_ddlnRho,k-1)
818 0 : latent_ddlnRho_m1%d1Array(i_lnT_m1) = s% d_eos_dlnT(i_latent_ddlnRho,k-1)
819 : end if
820 0 : end function wrap_latent_ddlnRho_m1
821 :
822 0 : function wrap_latent_ddlnRho_00(s, k) result(latent_ddlnRho_00)
823 0 : use eos_def, only: i_latent_ddlnRho
824 : type (star_info), pointer :: s
825 : type(auto_diff_real_star_order1) :: latent_ddlnRho_00
826 : integer, intent(in) :: k
827 0 : latent_ddlnRho_00 = 0d0
828 0 : latent_ddlnRho_00%val = s% latent_ddlnRho(k)
829 0 : latent_ddlnRho_00%d1Array(i_lnd_00) = s% d_eos_dlnd(i_latent_ddlnRho,k)
830 0 : latent_ddlnRho_00%d1Array(i_lnT_00) = s% d_eos_dlnT(i_latent_ddlnRho,k)
831 0 : end function wrap_latent_ddlnRho_00
832 :
833 0 : function wrap_latent_ddlnRho_p1(s, k) result(latent_ddlnRho_p1)
834 0 : use eos_def, only: i_latent_ddlnRho
835 : type (star_info), pointer :: s
836 : type(auto_diff_real_star_order1) :: latent_ddlnRho_p1
837 : integer, intent(in) :: k
838 0 : latent_ddlnRho_p1 = 0d0
839 0 : if (k < s% nz) then
840 0 : latent_ddlnRho_p1%val = s% latent_ddlnRho(k+1)
841 0 : latent_ddlnRho_p1%d1Array(i_lnd_p1) = s% d_eos_dlnd(i_latent_ddlnRho,k+1)
842 0 : latent_ddlnRho_p1%d1Array(i_lnT_p1) = s% d_eos_dlnT(i_latent_ddlnRho,k+1)
843 : end if
844 0 : end function wrap_latent_ddlnRho_p1
845 :
846 0 : function wrap_L_m1(s, k) result(L_m1)
847 : type (star_info), pointer :: s
848 : type(auto_diff_real_star_order1) :: L_m1
849 : integer, intent(in) :: k
850 0 : L_m1 = 0d0
851 0 : if (k > 1) then
852 0 : L_m1 % val = s%L(k-1)
853 0 : L_m1 % d1Array(i_L_m1) = 1d0
854 : end if
855 0 : end function wrap_L_m1
856 :
857 212336 : function wrap_L_00(s, k) result(L_00)
858 : type (star_info), pointer :: s
859 : type(auto_diff_real_star_order1) :: L_00
860 : integer, intent(in) :: k
861 212336 : L_00 = 0d0
862 212336 : L_00 % val = s%L(k)
863 212336 : L_00 % d1Array(i_L_00) = 1d0
864 212336 : end function wrap_L_00
865 :
866 52348 : function wrap_L_p1(s, k) result(L_p1)
867 : type (star_info), pointer :: s
868 : type(auto_diff_real_star_order1) :: L_p1
869 : integer, intent(in) :: k
870 52348 : L_p1 = 0d0
871 52348 : if (k < s%nz) then
872 52304 : L_p1 % val = s%L(k+1)
873 52304 : L_p1 % d1Array(i_L_p1) = 1d0
874 : else
875 44 : L_p1 %val = s%L_center
876 : ! L_center is a constant
877 : end if
878 52348 : end function wrap_L_p1
879 :
880 0 : function wrap_r_m1(s, k) result(r_m1)
881 : type (star_info), pointer :: s
882 : type(auto_diff_real_star_order1) :: r_m1
883 : integer, intent(in) :: k
884 0 : r_m1 = 0d0
885 0 : if (k > 1) then
886 0 : r_m1 % val = s%r(k-1)
887 0 : r_m1 % d1Array(i_lnR_m1) = s%r(k-1)
888 : end if
889 0 : end function wrap_r_m1
890 :
891 369248 : function wrap_r_00(s, k) result(r_00)
892 : type (star_info), pointer :: s
893 : type(auto_diff_real_star_order1) :: r_00
894 : integer, intent(in) :: k
895 369248 : r_00 = 0d0
896 369248 : r_00 % val = s%r(k)
897 369248 : r_00 % d1Array(i_lnR_00) = s%r(k)
898 369248 : end function wrap_r_00
899 :
900 52348 : function wrap_r_p1(s, k) result(r_p1)
901 : type (star_info), pointer :: s
902 : type(auto_diff_real_star_order1) :: r_p1
903 : integer, intent(in) :: k
904 52348 : r_p1 = 0d0
905 52348 : if (k < s%nz) then
906 52304 : r_p1 % val = s%r(k+1)
907 52304 : r_p1 % d1Array(i_lnR_p1) = s%r(k+1)
908 : else
909 44 : r_p1 % val = s%r_center
910 : end if
911 52348 : end function wrap_r_p1
912 :
913 0 : function wrap_lnR_m1(s, k) result(lnR_m1)
914 : type (star_info), pointer :: s
915 : type(auto_diff_real_star_order1) :: lnR_m1
916 : integer, intent(in) :: k
917 0 : lnR_m1 = 0d0
918 0 : if (k > 1) then
919 0 : lnR_m1 % val = s%lnR(k-1)
920 0 : lnR_m1 % d1Array(i_lnR_m1) = 1d0
921 : end if
922 0 : end function wrap_lnR_m1
923 :
924 52348 : function wrap_lnR_00(s, k) result(lnR_00)
925 : type (star_info), pointer :: s
926 : type(auto_diff_real_star_order1) :: lnR_00
927 : integer, intent(in) :: k
928 52348 : lnR_00 = 0d0
929 52348 : lnR_00 % val = s%lnR(k)
930 52348 : lnR_00 % d1Array(i_lnR_00) = 1d0
931 52348 : end function wrap_lnR_00
932 :
933 0 : function wrap_lnR_p1(s, k) result(lnR_p1)
934 : type (star_info), pointer :: s
935 : type(auto_diff_real_star_order1) :: lnR_p1
936 : integer, intent(in) :: k
937 0 : lnR_p1 = 0d0
938 0 : if (k < s%nz) then
939 0 : lnR_p1 % val = s%lnR(k+1)
940 0 : lnR_p1 % d1Array(i_lnR_p1) = 1d0
941 : else
942 0 : lnR_p1 % val = log(max(1d0,s%r_center))
943 : end if
944 0 : end function wrap_lnR_p1
945 :
946 0 : function wrap_dxh_lnR(s, k) result(dxh_lnR)
947 : type (star_info), pointer :: s
948 : type(auto_diff_real_star_order1) :: dxh_lnR
949 : integer, intent(in) :: k
950 0 : dxh_lnR = 0d0
951 0 : dxh_lnR % val = s%dxh_lnR(k)
952 0 : dxh_lnR % d1Array(i_lnR_00) = 1d0
953 0 : end function wrap_dxh_lnR
954 :
955 0 : function wrap_u_face_m1(s, k) result(v_m1)
956 : type (star_info), pointer :: s
957 : type(auto_diff_real_star_order1) :: v_m1
958 : integer, intent(in) :: k
959 0 : v_m1 = 0
960 0 : if (k > 1) v_m1 = shift_m1(s% u_face_ad(k-1))
961 0 : end function wrap_u_face_m1
962 :
963 0 : function wrap_u_face_00(s, k) result(v_00)
964 : type (star_info), pointer :: s
965 : type(auto_diff_real_star_order1) :: v_00
966 : integer, intent(in) :: k
967 0 : v_00 = s% u_face_ad(k)
968 0 : end function wrap_u_face_00
969 :
970 0 : function wrap_u_face_p1(s, k) result(v_p1)
971 : type (star_info), pointer :: s
972 : type(auto_diff_real_star_order1) :: v_p1
973 : integer, intent(in) :: k
974 0 : v_p1 = 0
975 0 : if (k < s% nz) v_p1 = shift_p1(s% u_face_ad(k+1))
976 0 : end function wrap_u_face_p1
977 :
978 0 : function wrap_v_m1(s, k) result(v_m1)
979 : type (star_info), pointer :: s
980 : type(auto_diff_real_star_order1) :: v_m1
981 : integer, intent(in) :: k
982 0 : if (s% u_flag) then
983 0 : v_m1 = wrap_u_face_m1(s,k)
984 0 : return
985 : end if
986 0 : v_m1 = 0d0
987 0 : if (k > 1) then
988 0 : v_m1 % val = s%v(k-1)
989 0 : v_m1 % d1Array(i_v_m1) = 1d0
990 : end if
991 0 : end function wrap_v_m1
992 :
993 0 : function wrap_v_00(s, k) result(v_00)
994 : type (star_info), pointer :: s
995 : type(auto_diff_real_star_order1) :: v_00
996 : integer, intent(in) :: k
997 0 : if (s% u_flag) then
998 0 : v_00 = wrap_u_face_00(s,k)
999 0 : return
1000 : end if
1001 0 : v_00 = 0d0
1002 0 : v_00 % val = s%v(k)
1003 0 : v_00 % d1Array(i_v_00) = 1d0
1004 0 : end function wrap_v_00
1005 :
1006 0 : function wrap_v_p1(s, k) result(v_p1)
1007 : type (star_info), pointer :: s
1008 : type(auto_diff_real_star_order1) :: v_p1
1009 : integer, intent(in) :: k
1010 0 : if (s% u_flag) then
1011 0 : v_p1 = wrap_u_face_p1(s,k)
1012 0 : return
1013 : end if
1014 0 : v_p1 = 0d0
1015 0 : if (k < s%nz) then
1016 0 : v_p1 % val = s%v(k+1)
1017 0 : v_p1 % d1Array(i_v_p1) = 1d0
1018 : else
1019 0 : v_p1 % val = s%v_center
1020 : ! v_center is a constant
1021 : end if
1022 0 : end function wrap_v_p1
1023 :
1024 0 : function wrap_opt_time_center_r_m1(s, k) result(r_tc)
1025 : type (star_info), pointer :: s
1026 : integer, intent(in) :: k
1027 : type(auto_diff_real_star_order1) :: r_tc
1028 0 : r_tc = wrap_r_m1(s,k)
1029 0 : if (s% using_velocity_time_centering) then
1030 0 : if (k > 1) r_tc = 0.5d0*(r_tc + s% r_start(k-1))
1031 : end if
1032 0 : end function wrap_opt_time_center_r_m1
1033 :
1034 0 : function wrap_opt_time_center_r_00(s, k) result(r_tc)
1035 : type (star_info), pointer :: s
1036 : integer, intent(in) :: k
1037 : type(auto_diff_real_star_order1) :: r_tc
1038 0 : r_tc = wrap_r_00(s,k)
1039 0 : if (s% using_velocity_time_centering) &
1040 0 : r_tc = 0.5d0*(r_tc + s% r_start(k))
1041 0 : end function wrap_opt_time_center_r_00
1042 :
1043 0 : function wrap_opt_time_center_r_p1(s, k) result(r_tc)
1044 : type (star_info), pointer :: s
1045 : integer, intent(in) :: k
1046 : type(auto_diff_real_star_order1) :: r_tc
1047 0 : r_tc = wrap_r_p1(s,k)
1048 0 : if (s% using_velocity_time_centering) then
1049 0 : if (k < s% nz) then
1050 0 : r_tc = 0.5d0*(r_tc + s% r_start(k+1))
1051 : else
1052 0 : r_tc = 0.5d0*(r_tc + s% r_center)
1053 : end if
1054 : end if
1055 0 : end function wrap_opt_time_center_r_p1
1056 :
1057 0 : function wrap_opt_time_center_v_m1(s, k) result(v_tc)
1058 : type (star_info), pointer :: s
1059 : integer, intent(in) :: k
1060 : type(auto_diff_real_star_order1) :: v_tc
1061 0 : v_tc = 0
1062 0 : if (k == 1) return
1063 0 : if (s% v_flag) then
1064 0 : v_tc = wrap_v_m1(s,k)
1065 0 : if (s% using_velocity_time_centering) &
1066 0 : v_tc = 0.5d0*(v_tc + s% v_start(k-1))
1067 0 : else if (s% u_flag) then
1068 0 : v_tc = wrap_u_face_m1(s,k)
1069 0 : if (s% using_velocity_time_centering) &
1070 0 : v_tc = 0.5d0*(v_tc + s% u_face_start(k-1))
1071 : end if
1072 0 : end function wrap_opt_time_center_v_m1
1073 :
1074 0 : function wrap_opt_time_center_v_00(s, k) result(v_tc)
1075 : type (star_info), pointer :: s
1076 : integer, intent(in) :: k
1077 : type(auto_diff_real_star_order1) :: v_tc
1078 0 : v_tc = 0
1079 0 : if (s% v_flag) then
1080 0 : v_tc = wrap_v_00(s,k)
1081 0 : if (s% using_velocity_time_centering) &
1082 0 : v_tc = 0.5d0*(v_tc + s% v_start(k))
1083 0 : else if (s% u_flag) then
1084 0 : v_tc = wrap_u_face_00(s,k)
1085 0 : if (s% using_velocity_time_centering) &
1086 0 : v_tc = 0.5d0*(v_tc + s% u_face_start(k))
1087 : end if
1088 0 : end function wrap_opt_time_center_v_00
1089 :
1090 0 : function wrap_opt_time_center_v_p1(s, k) result(v_tc)
1091 : type (star_info), pointer :: s
1092 : integer, intent(in) :: k
1093 : type(auto_diff_real_star_order1) :: v_tc
1094 0 : v_tc = 0
1095 0 : if (k == s% nz) return
1096 0 : if (s% v_flag) then
1097 0 : v_tc = wrap_v_p1(s,k)
1098 0 : if (s% using_velocity_time_centering) &
1099 0 : v_tc = 0.5d0*(v_tc + s% v_start(k+1))
1100 0 : else if (s% u_flag) then
1101 0 : v_tc = wrap_u_face_p1(s,k)
1102 0 : if (s% using_velocity_time_centering) &
1103 0 : v_tc = 0.5d0*(v_tc + s% u_face_start(k+1))
1104 : end if
1105 0 : end function wrap_opt_time_center_v_p1
1106 :
1107 0 : function wrap_u_m1(s, k) result(v_m1)
1108 : type (star_info), pointer :: s
1109 : type(auto_diff_real_star_order1) :: v_m1
1110 : integer, intent(in) :: k
1111 0 : v_m1 = 0d0
1112 0 : if (k > 1) then
1113 0 : v_m1 % val = s%u(k-1)
1114 0 : v_m1 % d1Array(i_v_m1) = 1d0
1115 : end if
1116 0 : end function wrap_u_m1
1117 :
1118 0 : function wrap_u_00(s, k) result(v_00)
1119 : type (star_info), pointer :: s
1120 : type(auto_diff_real_star_order1) :: v_00
1121 : integer, intent(in) :: k
1122 0 : v_00 = 0d0
1123 0 : v_00 % val = s%u(k)
1124 0 : v_00 % d1Array(i_v_00) = 1d0
1125 0 : end function wrap_u_00
1126 :
1127 0 : function wrap_u_p1(s, k) result(v_p1)
1128 : type (star_info), pointer :: s
1129 : type(auto_diff_real_star_order1) :: v_p1
1130 : integer, intent(in) :: k
1131 0 : v_p1 = 0d0
1132 0 : if (k < s%nz) then
1133 0 : v_p1 % val = s%u(k+1)
1134 0 : v_p1 % d1Array(i_v_p1) = 1d0
1135 : else
1136 0 : v_p1 % val = 0d0
1137 : ! v_center is a constant
1138 : end if
1139 0 : end function wrap_u_p1
1140 :
1141 0 : function wrap_Hp_m1(s, k) result(Hp_m1)
1142 : type (star_info), pointer :: s
1143 : type(auto_diff_real_star_order1) :: Hp_m1
1144 : integer, intent(in) :: k
1145 0 : Hp_m1 = 0d0
1146 0 : if (k > 1) then
1147 0 : Hp_m1 % val = s%Hp_face(k-1)
1148 0 : Hp_m1 % d1Array(i_Hp_m1) = 1d0
1149 : end if
1150 0 : end function wrap_Hp_m1
1151 :
1152 0 : function wrap_Hp_00(s, k) result(Hp_00)
1153 : type (star_info), pointer :: s
1154 : type(auto_diff_real_star_order1) :: Hp_00
1155 : integer, intent(in) :: k
1156 0 : Hp_00 = 0d0
1157 0 : Hp_00 % val = s%Hp_face(k)
1158 0 : Hp_00 % d1Array(i_Hp_00) = 1d0
1159 0 : end function wrap_Hp_00
1160 :
1161 0 : function wrap_Hp_p1(s, k) result(Hp_p1)
1162 : type (star_info), pointer :: s
1163 : type(auto_diff_real_star_order1) :: Hp_p1
1164 : integer, intent(in) :: k
1165 0 : Hp_p1 = 0d0
1166 0 : if (k < s%nz) then
1167 0 : Hp_p1 % val = s%Hp_face(k+1)
1168 0 : Hp_p1 % d1Array(i_Hp_p1) = 1d0
1169 : end if
1170 0 : end function wrap_Hp_p1
1171 :
1172 :
1173 0 : function wrap_w_div_wc_m1(s, k) result(w_div_wc_m1)
1174 : type (star_info), pointer :: s
1175 : type(auto_diff_real_star_order1) :: w_div_wc_m1
1176 : integer, intent(in) :: k
1177 0 : w_div_wc_m1 = 0d0
1178 0 : if (k > 1) then
1179 0 : w_div_wc_m1 % val = s% w_div_w_crit_roche(k-1)
1180 0 : w_div_wc_m1 % d1Array(i_w_div_wc_m1) = 1d0
1181 : end if
1182 0 : end function wrap_w_div_wc_m1
1183 :
1184 0 : function wrap_w_div_wc_00(s, k) result(w_div_wc_00)
1185 : type (star_info), pointer :: s
1186 : type(auto_diff_real_star_order1) :: w_div_wc_00
1187 : integer, intent(in) :: k
1188 0 : w_div_wc_00 = 0d0
1189 0 : w_div_wc_00 % val = s% w_div_w_crit_roche(k)
1190 0 : w_div_wc_00 % d1Array(i_w_div_wc_00) = 1d0
1191 0 : end function wrap_w_div_wc_00
1192 :
1193 0 : function wrap_w_div_wc_p1(s, k) result(w_div_wc_p1)
1194 : type (star_info), pointer :: s
1195 : type(auto_diff_real_star_order1) :: w_div_wc_p1
1196 : integer, intent(in) :: k
1197 0 : w_div_wc_p1 = 0d0
1198 0 : if (k < s%nz) then
1199 0 : w_div_wc_p1 % val = s% w_div_w_crit_roche(k+1)
1200 0 : w_div_wc_p1 % d1Array(i_w_div_wc_p1) = 1d0
1201 : end if
1202 0 : end function wrap_w_div_wc_p1
1203 :
1204 :
1205 0 : function wrap_jrot_m1(s, k) result(jrot_m1)
1206 : type (star_info), pointer :: s
1207 : type(auto_diff_real_star_order1) :: jrot_m1
1208 : integer, intent(in) :: k
1209 0 : jrot_m1 = 0d0
1210 0 : if (k > 1) then
1211 0 : jrot_m1 % val = s% j_rot(k-1)
1212 0 : jrot_m1 % d1Array(i_jrot_m1) = 1d0
1213 : end if
1214 0 : end function wrap_jrot_m1
1215 :
1216 0 : function wrap_jrot_00(s, k) result(jrot_00)
1217 : type (star_info), pointer :: s
1218 : type(auto_diff_real_star_order1) :: jrot_00
1219 : integer, intent(in) :: k
1220 0 : jrot_00 = 0d0
1221 0 : jrot_00 % val = s% j_rot(k)
1222 0 : jrot_00 % d1Array(i_jrot_00) = 1d0
1223 0 : end function wrap_jrot_00
1224 :
1225 0 : function wrap_jrot_p1(s, k) result(jrot_p1)
1226 : type (star_info), pointer :: s
1227 : type(auto_diff_real_star_order1) :: jrot_p1
1228 : integer, intent(in) :: k
1229 0 : jrot_p1 = 0d0
1230 0 : if (k < s%nz) then
1231 0 : jrot_p1 % val = s% j_rot(k+1)
1232 0 : jrot_p1 % d1Array(i_jrot_p1) = 1d0
1233 : end if
1234 0 : end function wrap_jrot_p1
1235 :
1236 :
1237 0 : function wrap_omega_m1(s, k) result(omega_m1)
1238 : type (star_info), pointer :: s
1239 : type(auto_diff_real_star_order1) :: omega_m1, jrot_m1
1240 : integer, intent(in) :: k
1241 0 : omega_m1 = 0d0
1242 0 : jrot_m1 = wrap_jrot_m1(s,k)
1243 0 : if (k > 1) then
1244 0 : omega_m1 = jrot_m1/shift_m1(s% i_rot(k-1))
1245 : end if
1246 0 : end function wrap_omega_m1
1247 :
1248 0 : function wrap_omega_00(s, k) result(omega_00)
1249 : type (star_info), pointer :: s
1250 : type(auto_diff_real_star_order1) :: omega_00, jrot_00
1251 : integer, intent(in) :: k
1252 0 : jrot_00 = wrap_jrot_00(s,k)
1253 0 : omega_00 = jrot_00/s% i_rot(k)
1254 0 : end function wrap_omega_00
1255 :
1256 0 : function wrap_omega_p1(s, k) result(omega_p1)
1257 : type (star_info), pointer :: s
1258 : type(auto_diff_real_star_order1) :: omega_p1, jrot_p1
1259 : integer, intent(in) :: k
1260 0 : omega_p1 = 0d0
1261 0 : jrot_p1 = wrap_jrot_p1(s,k)
1262 0 : if (k < s%nz) then
1263 0 : omega_p1 = jrot_p1/shift_p1(s% i_rot(k+1))
1264 : end if
1265 0 : end function wrap_omega_p1
1266 :
1267 :
1268 0 : function wrap_xtra1_m1(s, k) result(xtra1_m1)
1269 : type (star_info), pointer :: s
1270 : type(auto_diff_real_star_order1) :: xtra1_m1
1271 : integer, intent(in) :: k
1272 0 : xtra1_m1 = 0d0
1273 0 : if (k > 1) then
1274 0 : xtra1_m1 % val = 0d0
1275 0 : xtra1_m1 % d1Array(i_xtra1_m1) = 1d0
1276 : end if
1277 0 : end function wrap_xtra1_m1
1278 :
1279 0 : function wrap_xtra1_00(s, k) result(xtra1_00)
1280 : type (star_info), pointer :: s
1281 : type(auto_diff_real_star_order1) :: xtra1_00
1282 : integer, intent(in) :: k
1283 0 : xtra1_00 = 0d0
1284 0 : xtra1_00 % val = 0d0
1285 0 : xtra1_00 % d1Array(i_xtra1_00) = 1d0
1286 0 : end function wrap_xtra1_00
1287 :
1288 0 : function wrap_xtra1_p1(s, k) result(xtra1_p1)
1289 : type (star_info), pointer :: s
1290 : type(auto_diff_real_star_order1) :: xtra1_p1
1291 : integer, intent(in) :: k
1292 0 : xtra1_p1 = 0d0
1293 0 : if (k < s%nz) then
1294 0 : xtra1_p1 % val = 0d0
1295 0 : xtra1_p1 % d1Array(i_xtra1_p1) = 1d0
1296 : end if
1297 0 : end function wrap_xtra1_p1
1298 :
1299 :
1300 0 : function wrap_xtra2_m1(s, k) result(xtra2_m1)
1301 : type (star_info), pointer :: s
1302 : type(auto_diff_real_star_order1) :: xtra2_m1
1303 : integer, intent(in) :: k
1304 0 : xtra2_m1 = 0d0
1305 0 : if (k > 1) then
1306 0 : xtra2_m1 % val = 0d0
1307 0 : xtra2_m1 % d1Array(i_xtra2_m1) = 1d0
1308 : end if
1309 0 : end function wrap_xtra2_m1
1310 :
1311 0 : function wrap_xtra2_00(s, k) result(xtra2_00)
1312 : type (star_info), pointer :: s
1313 : type(auto_diff_real_star_order1) :: xtra2_00
1314 : integer, intent(in) :: k
1315 0 : xtra2_00 = 0d0
1316 0 : xtra2_00 % val = 0d0
1317 0 : xtra2_00 % d1Array(i_xtra2_00) = 1d0
1318 0 : end function wrap_xtra2_00
1319 :
1320 0 : function wrap_xtra2_p1(s, k) result(xtra2_p1)
1321 : type (star_info), pointer :: s
1322 : type(auto_diff_real_star_order1) :: xtra2_p1
1323 : integer, intent(in) :: k
1324 0 : xtra2_p1 = 0d0
1325 0 : if (k < s%nz) then
1326 0 : xtra2_p1 % val = 0d0
1327 0 : xtra2_p1 % d1Array(i_xtra2_p1) = 1d0
1328 : end if
1329 0 : end function wrap_xtra2_p1
1330 :
1331 : end module auto_diff_support
|
